Actual source code: ex1.c

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

  4: /*
  5:    Concepts: TS^pseudo-timestepping
  6:    Concepts: pseudo-timestepping
  7:    Concepts: nonlinear problems
  8:    Processors: 1

 10: */

 12: /* ------------------------------------------------------------------------

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

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

 23:     See snes/examples/tutorials/ex4.c[ex4f.F] and 
 24:     snes/examples/tutorials/ex5.c[ex5f.F] where the problem is described
 25:     and solved using Newton's method alone.

 27:   ----------------------------------------------------------------------------- */
 28: /*
 29:     Include "petscts.h" to use the PETSc timestepping routines. Note that
 30:     this file automatically includes "petsc.h" and other lower-level
 31:     PETSc include files.
 32: */
 33:  #include petscts.h

 35: /*
 36:   Create an application context to contain data needed by the 
 37:   application-provided call-back routines, FormJacobian() and
 38:   FormFunction().
 39: */
 40: typedef struct {
 41:   PetscReal   param;        /* test problem parameter */
 42:   PetscInt    mx;           /* Discretization in x-direction */
 43:   PetscInt    my;           /* Discretization in y-direction */
 44: } AppCtx;

 46: /* 
 47:    User-defined routines
 48: */
 50:      FormFunction(TS,PetscReal,Vec,Vec,void*),
 51:      FormInitialGuess(Vec,AppCtx*);

 55: int main(int argc,char **argv)
 56: {
 57:   TS             ts;                 /* timestepping context */
 58:   Vec            x,r;               /* solution, residual vectors */
 59:   Mat            J;                  /* Jacobian matrix */
 60:   AppCtx         user;               /* user-defined work context */
 61:   PetscInt       its,N;                /* iterations for convergence */
 63:   PetscReal      param_max = 6.81,param_min = 0.,dt;
 64:   PetscReal      ftime;

 66:   PetscInitialize(&argc,&argv,PETSC_NULL,help);
 67:   user.mx        = 4;
 68:   user.my        = 4;
 69:   user.param     = 6.0;
 70: 
 71:   /*
 72:      Allow user to set the grid dimensions and nonlinearity parameter at run-time
 73:   */
 74:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
 75:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
 76:   PetscOptionsGetReal(PETSC_NULL,"-param",&user.param,PETSC_NULL);
 77:   if (user.param >= param_max || user.param <= param_min) {
 78:     SETERRQ(1,"Parameter is out of range");
 79:   }
 80:   dt = .5/PetscMax(user.mx,user.my);
 81:   PetscOptionsGetReal(PETSC_NULL,"-dt",&dt,PETSC_NULL);
 82:   N          = user.mx*user.my;
 83: 
 84:   /* 
 85:       Create vectors to hold the solution and function value
 86:   */
 87:   VecCreateSeq(PETSC_COMM_SELF,N,&x);
 88:   VecDuplicate(x,&r);

 90:   /*
 91:     Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
 92:     in the sparse matrix. Note that this is not the optimal strategy; see
 93:     the Performance chapter of the users manual for information on 
 94:     preallocating memory in sparse matrices.
 95:   */
 96:   MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J);

 98:   /* 
 99:      Create timestepper context 
100:   */
101:   TSCreate(PETSC_COMM_WORLD,&ts);
102:   TSSetProblemType(ts,TS_NONLINEAR);

104:   /*
105:      Tell the timestepper context where to compute solutions
106:   */
107:   TSSetSolution(ts,x);

109:   /*
110:      Provide the call-back for the nonlinear function we are 
111:      evaluating. Thus whenever the timestepping routines need the
112:      function they will call this routine. Note the final argument
113:      is the application context used by the call-back functions.
114:   */
115:   TSSetRHSFunction(ts,FormFunction,&user);

117:   /*
118:      Set the Jacobian matrix and the function used to compute 
119:      Jacobians.
120:   */
121:   TSSetRHSJacobian(ts,J,J,FormJacobian,&user);

123:   /*
124:        For the initial guess for the problem
125:   */
126:   FormInitialGuess(x,&user);

128:   /*
129:        This indicates that we are using pseudo timestepping to 
130:      find a steady state solution to the nonlinear problem.
131:   */
132:   TSSetType(ts,TS_PSEUDO);

134:   /*
135:        Set the initial time to start at (this is arbitrary for 
136:      steady state problems; and the initial timestep given above
137:   */
138:   TSSetInitialTimeStep(ts,0.0,dt);

140:   /*
141:       Set a large number of timesteps and final duration time
142:      to insure convergence to steady state.
143:   */
144:   TSSetDuration(ts,1000,1.e12);

146:   /*
147:       Use the default strategy for increasing the timestep
148:   */
149:   TSPseudoSetTimeStep(ts,TSPseudoDefaultTimeStep,0);

151:   /*
152:       Set any additional options from the options database. This
153:      includes all options for the nonlinear and linear solvers used
154:      internally the the timestepping routines.
155:   */
156:   TSSetFromOptions(ts);

158:   TSSetUp(ts);

160:   /*
161:       Perform the solve. This is where the timestepping takes place.
162:   */
163:   TSStep(ts,&its,&ftime);
164: 
165:   printf("Number of pseudo timesteps = %d final time %4.2e\n",(int)its,ftime);

167:   /* 
168:      Free the data structures constructed above
169:   */
170:   VecDestroy(x);
171:   VecDestroy(r);
172:   MatDestroy(J);
173:   TSDestroy(ts);
174:   PetscFinalize();

176:   return 0;
177: }
178: /* ------------------------------------------------------------------ */
179: /*           Bratu (Solid Fuel Ignition) Test Problem                 */
180: /* ------------------------------------------------------------------ */

182: /* --------------------  Form initial approximation ----------------- */

186: PetscErrorCode FormInitialGuess(Vec X,AppCtx *user)
187: {
188:   PetscInt       i,j,row,mx,my;
190:   PetscReal      one = 1.0,lambda;
191:   PetscReal      temp1,temp,hx,hy;
192:   PetscScalar    *x;

194:   mx         = user->mx;
195:   my         = user->my;
196:   lambda = user->param;

198:   hx    = one / (PetscReal)(mx-1);
199:   hy    = one / (PetscReal)(my-1);

201:   VecGetArray(X,&x);
202:   temp1 = lambda/(lambda + one);
203:   for (j=0; j<my; j++) {
204:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
205:     for (i=0; i<mx; i++) {
206:       row = i + j*mx;
207:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
208:         x[row] = 0.0;
209:         continue;
210:       }
211:       x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
212:     }
213:   }
214:   VecRestoreArray(X,&x);
215:   return 0;
216: }
217: /* --------------------  Evaluate Function F(x) --------------------- */

221: PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
222: {
223:   AppCtx         *user = (AppCtx*)ptr;
225:   PetscInt       i,j,row,mx,my;
226:   PetscReal      two = 2.0,one = 1.0,lambda;
227:   PetscReal      hx,hy,hxdhy,hydhx;
228:   PetscScalar    ut,ub,ul,ur,u,uxx,uyy,sc,*x,*f;

230:   mx         = user->mx;
231:   my         = user->my;
232:   lambda = user->param;

234:   hx    = one / (PetscReal)(mx-1);
235:   hy    = one / (PetscReal)(my-1);
236:   sc    = hx*hy;
237:   hxdhy = hx/hy;
238:   hydhx = hy/hx;

240:   VecGetArray(X,&x);
241:   VecGetArray(F,&f);
242:   for (j=0; j<my; j++) {
243:     for (i=0; i<mx; i++) {
244:       row = i + j*mx;
245:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
246:         f[row] = x[row];
247:         continue;
248:       }
249:       u = x[row];
250:       ub = x[row - mx];
251:       ul = x[row - 1];
252:       ut = x[row + mx];
253:       ur = x[row + 1];
254:       uxx = (-ur + two*u - ul)*hydhx;
255:       uyy = (-ut + two*u - ub)*hxdhy;
256:       f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u);
257:     }
258:   }
259:   VecRestoreArray(X,&x);
260:   VecRestoreArray(F,&f);
261:   return 0;
262: }
263: /* --------------------  Evaluate Jacobian F'(x) -------------------- */

267: /*
268:    Calculate the Jacobian matrix J(X,t).

270:    Note: We put the Jacobian in the preconditioner storage B instead of J. This
271:    way we can give the -snes_mf_operator flag to check our work. This replaces
272:    J with a finite difference approximation, using our analytic Jacobian B for
273:    the preconditioner.
274: */
275: PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
276: {
277:   AppCtx         *user = (AppCtx*)ptr;
278:   Mat            jac = *B;
279:   PetscInt       i,j,row,mx,my,col[5];
281:   PetscScalar    two = 2.0,one = 1.0,lambda,v[5],sc,*x;
282:   PetscReal      hx,hy,hxdhy,hydhx;


285:   mx         = user->mx;
286:   my         = user->my;
287:   lambda = user->param;

289:   hx    = 1.0 / (PetscReal)(mx-1);
290:   hy    = 1.0 / (PetscReal)(my-1);
291:   sc    = hx*hy;
292:   hxdhy = hx/hy;
293:   hydhx = hy/hx;

295:   VecGetArray(X,&x);
296:   for (j=0; j<my; j++) {
297:     for (i=0; i<mx; i++) {
298:       row = i + j*mx;
299:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
300:         MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
301:         continue;
302:       }
303:       v[0] = hxdhy; col[0] = row - mx;
304:       v[1] = hydhx; col[1] = row - 1;
305:       v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row;
306:       v[3] = hydhx; col[3] = row + 1;
307:       v[4] = hxdhy; col[4] = row + mx;
308:       MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
309:     }
310:   }
311:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
312:   VecRestoreArray(X,&x);
313:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
314:   *flag = SAME_NONZERO_PATTERN;
315:   return 0;
316: }