Actual source code: ex1.c

petsc-master 2017-06-28
Report Typos and Errors

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

  4: /*
  5:    Concepts: TS^pseudo-timestepping
  6:    Concepts: pseudo-timestepping
  7:    Concepts: TS^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 "petscsys.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: */
 49: extern PetscErrorCode  FormJacobian(TS,PetscReal,Vec,Mat,Mat,void*), FormFunction(TS,PetscReal,Vec,Vec,void*), FormInitialGuess(Vec,AppCtx*);

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

 63:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
 64:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 65:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only");

 67:   user.mx    = 4;
 68:   user.my    = 4;
 69:   user.param = 6.0;

 71:   /*
 72:      Allow user to set the grid dimensions and nonlinearity parameter at run-time
 73:   */
 74:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
 75:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
 76:   N  = user.mx*user.my;
 77:   dt = .5/PetscMax(user.mx,user.my);
 78:   PetscOptionsGetReal(NULL,NULL,"-param",&user.param,NULL);
 79:   if (user.param >= param_max || user.param <= param_min) SETERRQ(PETSC_COMM_SELF,1,"Parameter is out of range");

 81:   /*
 82:       Create vectors to hold the solution and function value
 83:   */
 84:   VecCreateSeq(PETSC_COMM_SELF,N,&x);
 85:   VecDuplicate(x,&r);

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

 95:   /*
 96:      Create timestepper context
 97:   */
 98:   TSCreate(PETSC_COMM_WORLD,&ts);
 99:   TSSetProblemType(ts,TS_NONLINEAR);

101:   /*
102:      Tell the timestepper context where to compute solutions
103:   */
104:   TSSetSolution(ts,x);

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

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

120:   /*
121:        Form the initial guess for the problem
122:   */
123:   FormInitialGuess(x,&user);

125:   /*
126:        This indicates that we are using pseudo timestepping to
127:      find a steady state solution to the nonlinear problem.
128:   */
129:   TSSetType(ts,TSPSEUDO);

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

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

144:   /*
145:       Use the default strategy for increasing the timestep
146:   */
147:   TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,0);

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

156:   TSSetUp(ts);

158:   /*
159:       Perform the solve. This is where the timestepping takes place.
160:   */
161:   TSSolve(ts,x);
162:   TSGetSolveTime(ts,&ftime);

164:   /*
165:       Get the number of steps
166:   */
167:   TSGetTimeStepNumber(ts,&its);

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

171:   /*
172:      Free the data structures constructed above
173:   */
174:   VecDestroy(&x);
175:   VecDestroy(&r);
176:   MatDestroy(&J);
177:   TSDestroy(&ts);
178:   PetscFinalize();
179:   return ierr;
180: }
181: /* ------------------------------------------------------------------ */
182: /*           Bratu (Solid Fuel Ignition) Test Problem                 */
183: /* ------------------------------------------------------------------ */

185: /* --------------------  Form initial approximation ----------------- */

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

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

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

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

220: PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
221: {
222:   AppCtx            *user = (AppCtx*)ptr;
223:   PetscErrorCode    ierr;
224:   PetscInt          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,*f;
228:   const PetscScalar *x;

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:   VecGetArrayRead(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:   VecRestoreArrayRead(X,&x);
260:   VecRestoreArray(F,&f);
261:   return 0;
262: }
263: /* --------------------  Evaluate Jacobian F'(x) -------------------- */

265: /*
266:    Calculate the Jacobian matrix J(X,t).

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


283:   mx     = user->mx;
284:   my     = user->my;
285:   lambda = user->param;

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

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