Actual source code: ex1.c


  2: static char help[] = "Solves the time independent Bratu problem using pseudo-timestepping.";

  4: /* ------------------------------------------------------------------------

  6:     This code demonstrates how one may solve a nonlinear problem
  7:     with pseudo-timestepping. In this simple example, the pseudo-timestep
  8:     is the same for all grid points, i.e., this is equivalent to using
  9:     the backward Euler method with a variable timestep.

 11:     Note: This example does not require pseudo-timestepping since it
 12:     is an easy nonlinear problem, but it is included to demonstrate how
 13:     the pseudo-timestepping may be done.

 15:     See snes/tutorials/ex4.c[ex4f.F] and
 16:     snes/tutorials/ex5.c[ex5f.F] where the problem is described
 17:     and solved using Newton's method alone.

 19:   ----------------------------------------------------------------------------- */
 20: /*
 21:     Include "petscts.h" to use the PETSc timestepping routines. Note that
 22:     this file automatically includes "petscsys.h" and other lower-level
 23:     PETSc include files.
 24: */
 25: #include <petscts.h>

 27: /*
 28:   Create an application context to contain data needed by the
 29:   application-provided call-back routines, FormJacobian() and
 30:   FormFunction().
 31: */
 32: typedef struct {
 33:   PetscReal param; /* test problem parameter */
 34:   PetscInt  mx;    /* Discretization in x-direction */
 35:   PetscInt  my;    /* Discretization in y-direction */
 36: } AppCtx;

 38: /*
 39:    User-defined routines
 40: */
 41: extern PetscErrorCode FormJacobian(TS, PetscReal, Vec, Mat, Mat, void *), FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialGuess(Vec, AppCtx *);

 43: int main(int argc, char **argv)
 44: {
 45:   TS          ts;     /* timestepping context */
 46:   Vec         x, r;   /* solution, residual vectors */
 47:   Mat         J;      /* Jacobian matrix */
 48:   AppCtx      user;   /* user-defined work context */
 49:   PetscInt    its, N; /* iterations for convergence */
 50:   PetscReal   param_max = 6.81, param_min = 0., dt;
 51:   PetscReal   ftime;
 52:   PetscMPIInt size;

 54:   PetscFunctionBeginUser;
 55:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 56:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 57:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");

 59:   user.mx    = 4;
 60:   user.my    = 4;
 61:   user.param = 6.0;

 63:   /*
 64:      Allow user to set the grid dimensions and nonlinearity parameter at run-time
 65:   */
 66:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL));
 67:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL));
 68:   N  = user.mx * user.my;
 69:   dt = .5 / PetscMax(user.mx, user.my);
 70:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-param", &user.param, NULL));
 71:   PetscCheck(user.param < param_max && user.param >= param_min, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Parameter is out of range");

 73:   /*
 74:       Create vectors to hold the solution and function value
 75:   */
 76:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
 77:   PetscCall(VecDuplicate(x, &r));

 79:   /*
 80:     Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
 81:     in the sparse matrix. Note that this is not the optimal strategy; see
 82:     the Performance chapter of the users manual for information on
 83:     preallocating memory in sparse matrices.
 84:   */
 85:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 5, 0, &J));

 87:   /*
 88:      Create timestepper context
 89:   */
 90:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
 91:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));

 93:   /*
 94:      Tell the timestepper context where to compute solutions
 95:   */
 96:   PetscCall(TSSetSolution(ts, x));

 98:   /*
 99:      Provide the call-back for the nonlinear function we are
100:      evaluating. Thus whenever the timestepping routines need the
101:      function they will call this routine. Note the final argument
102:      is the application context used by the call-back functions.
103:   */
104:   PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &user));

106:   /*
107:      Set the Jacobian matrix and the function used to compute
108:      Jacobians.
109:   */
110:   PetscCall(TSSetRHSJacobian(ts, J, J, FormJacobian, &user));

112:   /*
113:        Form the initial guess for the problem
114:   */
115:   PetscCall(FormInitialGuess(x, &user));

117:   /*
118:        This indicates that we are using pseudo timestepping to
119:      find a steady state solution to the nonlinear problem.
120:   */
121:   PetscCall(TSSetType(ts, TSPSEUDO));

123:   /*
124:        Set the initial time to start at (this is arbitrary for
125:      steady state problems); and the initial timestep given above
126:   */
127:   PetscCall(TSSetTimeStep(ts, dt));

129:   /*
130:       Set a large number of timesteps and final duration time
131:      to insure convergence to steady state.
132:   */
133:   PetscCall(TSSetMaxSteps(ts, 10000));
134:   PetscCall(TSSetMaxTime(ts, 1e12));
135:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

137:   /*
138:       Use the default strategy for increasing the timestep
139:   */
140:   PetscCall(TSPseudoSetTimeStep(ts, TSPseudoTimeStepDefault, 0));

142:   /*
143:       Set any additional options from the options database. This
144:      includes all options for the nonlinear and linear solvers used
145:      internally the timestepping routines.
146:   */
147:   PetscCall(TSSetFromOptions(ts));

149:   PetscCall(TSSetUp(ts));

151:   /*
152:       Perform the solve. This is where the timestepping takes place.
153:   */
154:   PetscCall(TSSolve(ts, x));
155:   PetscCall(TSGetSolveTime(ts, &ftime));

157:   /*
158:       Get the number of steps
159:   */
160:   PetscCall(TSGetStepNumber(ts, &its));

162:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of pseudo timesteps = %" PetscInt_FMT " final time %4.2e\n", its, (double)ftime));

164:   /*
165:      Free the data structures constructed above
166:   */
167:   PetscCall(VecDestroy(&x));
168:   PetscCall(VecDestroy(&r));
169:   PetscCall(MatDestroy(&J));
170:   PetscCall(TSDestroy(&ts));
171:   PetscCall(PetscFinalize());
172:   return 0;
173: }
174: /* ------------------------------------------------------------------ */
175: /*           Bratu (Solid Fuel Ignition) Test Problem                 */
176: /* ------------------------------------------------------------------ */

178: /* --------------------  Form initial approximation ----------------- */

180: PetscErrorCode FormInitialGuess(Vec X, AppCtx *user)
181: {
182:   PetscInt     i, j, row, mx, my;
183:   PetscReal    one = 1.0, lambda;
184:   PetscReal    temp1, temp, hx, hy;
185:   PetscScalar *x;

187:   PetscFunctionBeginUser;
188:   mx     = user->mx;
189:   my     = user->my;
190:   lambda = user->param;

192:   hx = one / (PetscReal)(mx - 1);
193:   hy = one / (PetscReal)(my - 1);

195:   PetscCall(VecGetArray(X, &x));
196:   temp1 = lambda / (lambda + one);
197:   for (j = 0; j < my; j++) {
198:     temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy;
199:     for (i = 0; i < mx; i++) {
200:       row = i + j * mx;
201:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
202:         x[row] = 0.0;
203:         continue;
204:       }
205:       x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp));
206:     }
207:   }
208:   PetscCall(VecRestoreArray(X, &x));
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }
211: /* --------------------  Evaluate Function F(x) --------------------- */

213: PetscErrorCode FormFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
214: {
215:   AppCtx            *user = (AppCtx *)ptr;
216:   PetscInt           i, j, row, mx, my;
217:   PetscReal          two = 2.0, one = 1.0, lambda;
218:   PetscReal          hx, hy, hxdhy, hydhx;
219:   PetscScalar        ut, ub, ul, ur, u, uxx, uyy, sc, *f;
220:   const PetscScalar *x;

222:   PetscFunctionBeginUser;
223:   mx     = user->mx;
224:   my     = user->my;
225:   lambda = user->param;

227:   hx    = one / (PetscReal)(mx - 1);
228:   hy    = one / (PetscReal)(my - 1);
229:   sc    = hx * hy;
230:   hxdhy = hx / hy;
231:   hydhx = hy / hx;

233:   PetscCall(VecGetArrayRead(X, &x));
234:   PetscCall(VecGetArray(F, &f));
235:   for (j = 0; j < my; j++) {
236:     for (i = 0; i < mx; i++) {
237:       row = i + j * mx;
238:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
239:         f[row] = x[row];
240:         continue;
241:       }
242:       u      = x[row];
243:       ub     = x[row - mx];
244:       ul     = x[row - 1];
245:       ut     = x[row + mx];
246:       ur     = x[row + 1];
247:       uxx    = (-ur + two * u - ul) * hydhx;
248:       uyy    = (-ut + two * u - ub) * hxdhy;
249:       f[row] = -uxx + -uyy + sc * lambda * PetscExpScalar(u);
250:     }
251:   }
252:   PetscCall(VecRestoreArrayRead(X, &x));
253:   PetscCall(VecRestoreArray(F, &f));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }
256: /* --------------------  Evaluate Jacobian F'(x) -------------------- */

258: /*
259:    Calculate the Jacobian matrix J(X,t).

261:    Note: We put the Jacobian in the preconditioner storage B instead of J. This
262:    way we can give the -snes_mf_operator flag to check our work. This replaces
263:    J with a finite difference approximation, using our analytic Jacobian B for
264:    the preconditioner.
265: */
266: PetscErrorCode FormJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat B, void *ptr)
267: {
268:   AppCtx            *user = (AppCtx *)ptr;
269:   PetscInt           i, j, row, mx, my, col[5];
270:   PetscScalar        two = 2.0, one = 1.0, lambda, v[5], sc;
271:   const PetscScalar *x;
272:   PetscReal          hx, hy, hxdhy, hydhx;

274:   PetscFunctionBeginUser;
275:   mx     = user->mx;
276:   my     = user->my;
277:   lambda = user->param;

279:   hx    = 1.0 / (PetscReal)(mx - 1);
280:   hy    = 1.0 / (PetscReal)(my - 1);
281:   sc    = hx * hy;
282:   hxdhy = hx / hy;
283:   hydhx = hy / hx;

285:   PetscCall(VecGetArrayRead(X, &x));
286:   for (j = 0; j < my; j++) {
287:     for (i = 0; i < mx; i++) {
288:       row = i + j * mx;
289:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
290:         PetscCall(MatSetValues(B, 1, &row, 1, &row, &one, INSERT_VALUES));
291:         continue;
292:       }
293:       v[0]   = hxdhy;
294:       col[0] = row - mx;
295:       v[1]   = hydhx;
296:       col[1] = row - 1;
297:       v[2]   = -two * (hydhx + hxdhy) + sc * lambda * PetscExpScalar(x[row]);
298:       col[2] = row;
299:       v[3]   = hydhx;
300:       col[3] = row + 1;
301:       v[4]   = hxdhy;
302:       col[4] = row + mx;
303:       PetscCall(MatSetValues(B, 1, &row, 5, col, v, INSERT_VALUES));
304:     }
305:   }
306:   PetscCall(VecRestoreArrayRead(X, &x));
307:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
308:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
309:   if (J != B) {
310:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
311:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
312:   }
313:   PetscFunctionReturn(PETSC_SUCCESS);
314: }

316: /*TEST

318:     test:
319:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_monitor_pseudo -snes_atol 1.e-7 -ts_pseudo_frtol 1.e-5 -ts_view draw:tikz:fig.tex

321:     test:
322:       suffix: 2
323:       args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5

325: TEST*/