Actual source code: ex8.c

  2: /* Program usage:  ./ex8 [-help] [all PETSc options] */

  4: static char help[] = "Nonlinear DAE.\n";


  7: /*
  8:    Include "petscts.h" so that we can use TS solvers.  Note that this
  9:    file automatically includes:
 10:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 11:      petscmat.h - matrices
 12:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 13:      petscviewer.h - viewers               petscpc.h  - preconditioners
 14:      petscksp.h   - linear solvers
 15: */
 16:  #include petscts.h

 18: typedef struct _Problem *Problem;
 19: struct _Problem {
 20:   PetscErrorCode (*destroy)(Problem);
 21:   TSIFunction function;
 22:   TSIJacobian jacobian;
 23:   PetscErrorCode (*solution)(PetscReal,Vec,void*);

 25:   MPI_Comm comm;
 26:   PetscReal final_time;
 27:   PetscInt n;
 28:   void *data;
 29: };

 31: /*
 32: *  User-defined routines
 33: */

 35: /*
 36: *  Stiff 3-variable system from chemical reactions, due to Robertson (1966), problem ROBER in Hairer&Wanner, ODE 2, 1996
 37: */
 40: static PetscErrorCode RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
 41: {
 43:   PetscScalar *x,*xdot,*f;

 46:   VecGetArray(X,&x);
 47:   VecGetArray(Xdot,&xdot);
 48:   VecGetArray(F,&f);
 49:   f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2];
 50:   f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]);
 51:   f[2] = xdot[2] - 3e7*PetscSqr(x[1]);
 52:   VecRestoreArray(X,&x);
 53:   VecRestoreArray(Xdot,&xdot);
 54:   VecRestoreArray(F,&f);
 55:   return(0);
 56: }

 60: static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
 61: {
 63:   PetscInt rowcol[] = {0,1,2};
 64:   PetscScalar *x,*xdot,J[3][3];

 67:   VecGetArray(X,&x);
 68:   VecGetArray(Xdot,&xdot);
 69:   J[0][0] = a + 0.04;     J[0][1] = -1e4*x[2];                   J[0][2] = -1e4*x[1];
 70:   J[1][0] = -0.04;        J[1][1] = a + 1e4*x[2] + 3e7*2*x[1];   J[1][2] = 1e4*x[1];
 71:   J[2][0] = 0;            J[2][1] = -3e7*2*x[1];                 J[2][2] = a;
 72:   MatSetValues(*B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
 73:   VecRestoreArray(X,&x);
 74:   VecRestoreArray(Xdot,&xdot);

 76:   MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
 77:   MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
 78:   if (*A != *B) {
 79:     MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
 80:     MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
 81:   }
 82:   *flag = SAME_NONZERO_PATTERN;
 83:   return(0);
 84: }

 88: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
 89: {
 91:   PetscScalar *x;

 94:   if (t != 0) SETERRQ(1,"not implemented");
 95:   VecGetArray(X,&x);
 96:   x[0] = 1;
 97:   x[1] = 0;
 98:   x[2] = 0;
 99:   VecRestoreArray(X,&x);
100:   return(0);
101: }

105: static PetscErrorCode RoberCreate(Problem p)
106: {

109:   p->destroy    = 0;
110:   p->function   = &RoberFunction;
111:   p->jacobian   = &RoberJacobian;
112:   p->solution   = &RoberSolution;
113:   p->final_time = 1e11;
114:   p->n          = 3;
115:   return(0);
116: }


119: typedef struct {
120:   PetscReal lambda;
121: } CECtx;

123: /*
124: * Stiff scalar valued problem with an exact solution
125: */
128: static PetscErrorCode CEDestroy(Problem p)
129: {

133:   PetscFree(p->data);
134:   return(0);
135: }

139: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
140: {
142:   PetscReal l = ((CECtx*)ctx)->lambda;
143:   PetscScalar *x,*xdot,*f;

146:   VecGetArray(X,&x);
147:   VecGetArray(Xdot,&xdot);
148:   VecGetArray(F,&f);
149:   f[0] = xdot[0] + l*(x[0] - cos(t));
150: #if 0
151:   PetscPrintf(PETSC_COMM_WORLD," f(t=%G,x=%G,xdot=%G) = %G\n",t,x[0],xdot[0],f[0]);
152: #endif
153:   VecRestoreArray(X,&x);
154:   VecRestoreArray(Xdot,&xdot);
155:   VecRestoreArray(F,&f);
156:   return(0);
157: }

161: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
162: {
163:   PetscReal l = ((CECtx*)ctx)->lambda;
165:   PetscInt rowcol[] = {0};
166:   PetscScalar *x,*xdot,J[1][1];

169:   VecGetArray(X,&x);
170:   VecGetArray(Xdot,&xdot);
171:   J[0][0] = a + l;
172:   MatSetValues(*B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
173:   VecRestoreArray(X,&x);
174:   VecRestoreArray(Xdot,&xdot);

176:   MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
177:   MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
178:   if (*A != *B) {
179:     MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
180:     MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
181:   }
182:   *flag = SAME_NONZERO_PATTERN;
183:   return(0);
184: }

188: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
189: {
190:   PetscReal l = ((CECtx*)ctx)->lambda;
192:   PetscScalar *x;

195:   VecGetArray(X,&x);
196:   x[0] = l/(l*l+1)*(l*cos(t)+sin(t)) - l*l/(l*l+1)*exp(-l*t);
197:   VecRestoreArray(X,&x);
198:   return(0);
199: }

203: static PetscErrorCode CECreate(Problem p)
204: {
206:   CECtx         *ce;

209:   PetscMalloc(sizeof(CECtx),&ce);
210:   p->data = (void*)ce;

212:   p->destroy    = &CEDestroy;
213:   p->function   = &CEFunction;
214:   p->jacobian   = &CEJacobian;
215:   p->solution   = &CESolution;
216:   p->final_time = 10;
217:   p->n          = 1;

219:   ce->lambda = 10;
220:   PetscOptionsBegin(p->comm,PETSC_NULL,"CE options","");
221:   {
222:     PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,PETSC_NULL);
223:   }
224:   PetscOptionsEnd();
225:   return(0);
226: }


229: /*
230: *  User-defined monitor for comparing to exact solutions when possible
231: */
232: typedef struct {
233:   MPI_Comm  comm;
234:   Problem   problem;
235:   Vec       x;
236: } MonitorCtx;

240: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
241: {
243:   MonitorCtx *mon = (MonitorCtx*)ctx;
244:   PetscReal h,nrm_x,nrm_exact,nrm_diff;

247:   if (!mon->problem->solution) return(0);
248:   (*mon->problem->solution)(t,mon->x,mon->problem->data);
249:   VecNorm(x,NORM_2,&nrm_x);
250:   VecNorm(mon->x,NORM_2,&nrm_exact);
251:   VecAYPX(mon->x,-1,x);
252:   VecNorm(mon->x,NORM_2,&nrm_diff);
253:   TSGetTimeStep(ts,&h);
254:   PetscPrintf(mon->comm,"step %4D t=%12.8e h=% 8.2e  |x|=%9.2e  |x_e|=%9.2e  |x-x_e|=%9.2e\n",step,t,h,nrm_x,nrm_exact,nrm_diff);
255:   return(0);
256: }



262: int main(int argc,char **argv)
263: {
264:   PetscFList      plist = PETSC_NULL;
265:   char            pname[256];
266:   TS              ts;           /* nonlinear solver */
267:   Vec             x,r;          /* solution, residual vectors */
268:   Mat             A;            /* Jacobian matrix */
269:   Problem         problem;
270:   PetscTruth      use_monitor;
271:   PetscInt        steps,maxsteps = 100;
272:   PetscReal       ftime;
273:   MonitorCtx      mon;
274:   PetscErrorCode  ierr;

276:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277:      Initialize program
278:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
279:   PetscInitialize(&argc,&argv,(char *)0,help);

281:   /* Register the available problems */
282:   PetscFListAdd(&plist,"rober","",(void(*)(void))&RoberCreate);
283:   PetscFListAdd(&plist,"ce",   "",(void(*)(void))&CECreate   );
284:   PetscStrcpy(pname,"ce");

286:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
287:     Set runtime options
288:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
289:   PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Timestepping benchmark options","");
290:   {
291:     PetscOptionsList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),PETSC_NULL);
292:     use_monitor = PETSC_FALSE;
293:     PetscOptionsTruth("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,PETSC_NULL);
294:   }
295:   PetscOptionsEnd();

297:   /* Create the new problem */
298:   PetscNew(struct _Problem,&problem);
299:   problem->comm = MPI_COMM_WORLD;
300:   {
301:     PetscErrorCode (*pcreate)(Problem);

303:     PetscFListFind(plist,MPI_COMM_WORLD,pname,(void(**)(void))&pcreate);
304:     if (!pcreate) SETERRQ1(1,"No problem '%s'",pname);
305:     (*pcreate)(problem);
306:   }


309:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
310:     Create necessary matrix and vectors, solve same ODE on every process
311:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
312:   MatCreate(PETSC_COMM_WORLD,&A);
313:   MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
314:   MatSetFromOptions(A);

316:   MatGetVecs(A,&x,PETSC_NULL);
317:   VecDuplicate(x,&r);

319:   mon.comm = PETSC_COMM_WORLD;
320:   mon.problem = problem;
321:   VecDuplicate(x,&mon.x);

323:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
324:      Create timestepping solver context
325:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
326:   TSCreate(PETSC_COMM_WORLD,&ts);
327:   TSSetProblemType(ts,TS_NONLINEAR);
328:   TSSetType(ts,TSGL); /* General Linear method, TSTHETA can also solve DAE */
329:   TSSetIFunction(ts,problem->function,problem->data);
330:   TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
331:   TSSetDuration(ts,maxsteps,problem->final_time);
332:   if (use_monitor) {
333:     TSMonitorSet(ts,&MonitorError,&mon,PETSC_NULL);
334:   }

336:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
337:      Set initial conditions
338:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
339:   (*problem->solution)(0,x,problem->data);
340:   TSSetInitialTimeStep(ts,0.0,.001);
341:   TSSetSolution(ts,x);

343:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
344:      Set runtime options
345:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
346:   TSSetFromOptions(ts);

348:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349:      Solve nonlinear system
350:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
351:   TSStep(ts,&steps,&ftime);
352:   PetscPrintf(PETSC_COMM_WORLD,"steps %D, ftime %G\n",steps,ftime);

354:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
355:      Free work space.  All PETSc objects should be destroyed when they
356:      are no longer needed.
357:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
358:   MatDestroy(A);
359:   VecDestroy(x);
360:   VecDestroy(r);
361:   VecDestroy(mon.x);
362:   TSDestroy(ts);
363:   if (problem->destroy) {
364:     (*problem->destroy)(problem);
365:   }
366:   PetscFree(problem);
367:   PetscFListDestroy(&plist);

369:   PetscFinalize();
370:   return(0);
371: }