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;

 55:   PetscInitialize(&argc, &argv, NULL, help);
 56:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 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:   PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL);
 67:   PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL);
 68:   N  = user.mx * user.my;
 69:   dt = .5 / PetscMax(user.mx, user.my);
 70:   PetscOptionsGetReal(NULL, NULL, "-param", &user.param, NULL);

 73:   /*
 74:       Create vectors to hold the solution and function value
 75:   */
 76:   VecCreateSeq(PETSC_COMM_SELF, N, &x);
 77:   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:   MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 5, 0, &J);

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

 93:   /*
 94:      Tell the timestepper context where to compute solutions
 95:   */
 96:   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:   TSSetRHSFunction(ts, NULL, FormFunction, &user);

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

112:   /*
113:        Form the initial guess for the problem
114:   */
115:   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:   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:   TSSetTimeStep(ts, dt);

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

137:   /*
138:       Use the default strategy for increasing the timestep
139:   */
140:   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:   TSSetFromOptions(ts);

149:   TSSetUp(ts);

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

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

162:   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:   VecDestroy(&x);
168:   VecDestroy(&r);
169:   MatDestroy(&J);
170:   TSDestroy(&ts);
171:   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:   mx     = user->mx;
188:   my     = user->my;
189:   lambda = user->param;

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

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

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

221:   mx     = user->mx;
222:   my     = user->my;
223:   lambda = user->param;

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

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

256: /*
257:    Calculate the Jacobian matrix J(X,t).

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

272:   mx     = user->mx;
273:   my     = user->my;
274:   lambda = user->param;

276:   hx    = 1.0 / (PetscReal)(mx - 1);
277:   hy    = 1.0 / (PetscReal)(my - 1);
278:   sc    = hx * hy;
279:   hxdhy = hx / hy;
280:   hydhx = hy / hx;

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

313: /*TEST

315:     test:
316:       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

318:     test:
319:       suffix: 2
320:       args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5

322: TEST*/