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:   PetscInitialize(&argc,&argv,NULL,help);
 55:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

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

 62:   /*
 63:      Allow user to set the grid dimensions and nonlinearity parameter at run-time
 64:   */
 65:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
 66:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
 67:   N  = user.mx*user.my;
 68:   dt = .5/PetscMax(user.mx,user.my);
 69:   PetscOptionsGetReal(NULL,NULL,"-param",&user.param,NULL);

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

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

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

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

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

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

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

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

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

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

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

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

148:   TSSetUp(ts);

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

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

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

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

177: /* --------------------  Form initial approximation ----------------- */

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

186:   mx     = user->mx;
187:   my     = user->my;
188:   lambda = user->param;

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

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

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

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

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

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

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

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

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

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

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

307: /*TEST

309:     test:
310:       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

312:     test:
313:       suffix: 2
314:       args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5

316: TEST*/