Actual source code: ex1.c

  1: /*$Id: ex1.c,v 1.31 2001/08/07 21:31:30 bsmith Exp $*/

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

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

 11: */

 13: /* ------------------------------------------------------------------------

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

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

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

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

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

 47: /* 
 48:    User-defined routines
 49: */
 50: extern int  FormJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*),
 51:      FormFunction(TS,PetscReal,Vec,Vec,void*),
 52:      FormInitialGuess(Vec,AppCtx*);

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

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

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

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

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

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

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

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

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

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

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

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

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

159:   TSSetUp(ts);

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

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

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

183: /* --------------------  Form initial approximation ----------------- */

185: #undef __FUNCT__
187: int FormInitialGuess(Vec X,AppCtx *user)
188: {
189:   int         i,j,row,mx,my,ierr;
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) --------------------- */

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

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

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

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

264: #undef __FUNCT__
266: int FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
267: {
268:   AppCtx      *user = (AppCtx*)ptr;
269:   Mat         jac = *B;
270:   int         i,j,row,mx,my,col[5],ierr;
271:   PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc,*x;
272:   PetscReal   hx,hy,hxdhy,hydhx;


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:   VecGetArray(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:         MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
291:         continue;
292:       }
293:       v[0] = hxdhy; col[0] = row - mx;
294:       v[1] = hydhx; col[1] = row - 1;
295:       v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row;
296:       v[3] = hydhx; col[3] = row + 1;
297:       v[4] = hxdhy; col[4] = row + mx;
298:       MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
299:     }
300:   }
301:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
302:   VecRestoreArray(X,&x);
303:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
304:   *flag = SAME_NONZERO_PATTERN;
305:   return 0;
306: }