Actual source code: ex8.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: static char help[] = "Nonlinear DAE benchmark problems.\n";

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

 15: typedef struct _Problem* Problem;
 16: struct _Problem {
 17:   PetscErrorCode (*destroy)(Problem);
 18:   TSIFunction    function;
 19:   TSIJacobian    jacobian;
 20:   PetscErrorCode (*solution)(PetscReal,Vec,void*);
 21:   MPI_Comm       comm;
 22:   PetscReal      final_time;
 23:   PetscInt       n;
 24:   PetscBool      hasexact;
 25:   void           *data;
 26: };

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

 40:   VecGetArrayRead(X,&x);
 41:   VecGetArrayRead(Xdot,&xdot);
 42:   VecGetArray(F,&f);
 43:   f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2];
 44:   f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]);
 45:   f[2] = xdot[2] - 3e7*PetscSqr(x[1]);
 46:   VecRestoreArrayRead(X,&x);
 47:   VecRestoreArrayRead(Xdot,&xdot);
 48:   VecRestoreArray(F,&f);
 49:   return(0);
 50: }

 54: static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
 55: {
 56:   PetscErrorCode    ierr;
 57:   PetscInt          rowcol[] = {0,1,2};
 58:   PetscScalar       J[3][3];
 59:   const PetscScalar *x,*xdot;

 62:   VecGetArrayRead(X,&x);
 63:   VecGetArrayRead(Xdot,&xdot);
 64:   J[0][0] = a + 0.04;     J[0][1] = -1e4*x[2];                   J[0][2] = -1e4*x[1];
 65:   J[1][0] = -0.04;        J[1][1] = a + 1e4*x[2] + 3e7*2*x[1];   J[1][2] = 1e4*x[1];
 66:   J[2][0] = 0;            J[2][1] = -3e7*2*x[1];                 J[2][2] = a;
 67:   MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
 68:   VecRestoreArrayRead(X,&x);
 69:   VecRestoreArrayRead(Xdot,&xdot);

 71:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 72:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 73:   if (A != B) {
 74:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 75:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 76:   }
 77:   return(0);
 78: }

 82: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
 83: {
 85:   PetscScalar    *x;

 88:   if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
 89:   VecGetArray(X,&x);
 90:   x[0] = 1;
 91:   x[1] = 0;
 92:   x[2] = 0;
 93:   VecRestoreArray(X,&x);
 94:   return(0);
 95: }

 99: static PetscErrorCode RoberCreate(Problem p)
100: {

103:   p->destroy    = 0;
104:   p->function   = &RoberFunction;
105:   p->jacobian   = &RoberJacobian;
106:   p->solution   = &RoberSolution;
107:   p->final_time = 1e11;
108:   p->n          = 3;
109:   return(0);
110: }

112: /*
113:      Stiff scalar valued problem
114: */

116: typedef struct {
117:   PetscReal lambda;
118: } CECtx;

122: static PetscErrorCode CEDestroy(Problem p)
123: {

127:   PetscFree(p->data);
128:   return(0);
129: }

133: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
134: {
135:   PetscErrorCode    ierr;
136:   PetscReal         l = ((CECtx*)ctx)->lambda;
137:   PetscScalar       *f;
138:   const PetscScalar *x,*xdot;

141:   VecGetArrayRead(X,&x);
142:   VecGetArrayRead(Xdot,&xdot);
143:   VecGetArray(F,&f);
144:   f[0] = xdot[0] + l*(x[0] - PetscCosReal(t));
145: #if 0
146:   PetscPrintf(PETSC_COMM_WORLD," f(t=%g,x=%g,xdot=%g) = %g\n",(double)t,(double)x[0],(double)xdot[0],(double)f[0]);
147: #endif
148:   VecRestoreArrayRead(X,&x);
149:   VecRestoreArrayRead(Xdot,&xdot);
150:   VecRestoreArray(F,&f);
151:   return(0);
152: }

156: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
157: {
158:   PetscReal         l = ((CECtx*)ctx)->lambda;
159:   PetscErrorCode    ierr;
160:   PetscInt          rowcol[] = {0};
161:   PetscScalar       J[1][1];
162:   const PetscScalar *x,*xdot;

165:   VecGetArrayRead(X,&x);
166:   VecGetArrayRead(Xdot,&xdot);
167:   J[0][0] = a + l;
168:   MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
169:   VecRestoreArrayRead(X,&x);
170:   VecRestoreArrayRead(Xdot,&xdot);

172:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
173:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
174:   if (A != B) {
175:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
176:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
177:   }
178:   return(0);
179: }

183: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
184: {
185:   PetscReal      l = ((CECtx*)ctx)->lambda;
187:   PetscScalar    *x;

190:   VecGetArray(X,&x);
191:   x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t);
192:   VecRestoreArray(X,&x);
193:   return(0);
194: }

198: static PetscErrorCode CECreate(Problem p)
199: {
201:   CECtx          *ce;

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

207:   p->destroy    = &CEDestroy;
208:   p->function   = &CEFunction;
209:   p->jacobian   = &CEJacobian;
210:   p->solution   = &CESolution;
211:   p->final_time = 10;
212:   p->n          = 1;
213:   p->hasexact   = PETSC_TRUE;

215:   ce->lambda = 10;
216:   PetscOptionsBegin(p->comm,NULL,"CE options","");
217:   {
218:     PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);
219:   }
220:   PetscOptionsEnd();
221:   return(0);
222: }

224: /*
225: *  Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
226: */
229: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
230: {
231:   PetscErrorCode    ierr;
232:   PetscScalar       *f;
233:   const PetscScalar *x,*xdot;

236:   VecGetArrayRead(X,&x);
237:   VecGetArrayRead(Xdot,&xdot);
238:   VecGetArray(F,&f);
239:   f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
240:   f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
241:   f[2] = xdot[2] - 0.161*(x[0] - x[2]);
242:   VecRestoreArrayRead(X,&x);
243:   VecRestoreArrayRead(Xdot,&xdot);
244:   VecRestoreArray(F,&f);
245:   return(0);
246: }

250: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
251: {
252:   PetscErrorCode    ierr;
253:   PetscInt          rowcol[] = {0,1,2};
254:   PetscScalar       J[3][3];
255:   const PetscScalar *x,*xdot;

258:   VecGetArrayRead(X,&x);
259:   VecGetArrayRead(Xdot,&xdot);
260:   J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
261:   J[0][1] = -77.27*(1. - x[0]);
262:   J[0][2] = 0;
263:   J[1][0] = 1./77.27*x[1];
264:   J[1][1] = a + 1./77.27*(1. + x[0]);
265:   J[1][2] = -1./77.27;
266:   J[2][0] = -0.161;
267:   J[2][1] = 0;
268:   J[2][2] = a + 0.161;
269:   MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
270:   VecRestoreArrayRead(X,&x);
271:   VecRestoreArrayRead(Xdot,&xdot);

273:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
274:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
275:   if (A != B) {
276:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
277:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
278:   }
279:   return(0);
280: }

284: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
285: {
287:   PetscScalar    *x;

290:   if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
291:   VecGetArray(X,&x);
292:   x[0] = 1;
293:   x[1] = 2;
294:   x[2] = 3;
295:   VecRestoreArray(X,&x);
296:   return(0);
297: }

301: static PetscErrorCode OregoCreate(Problem p)
302: {

305:   p->destroy    = 0;
306:   p->function   = &OregoFunction;
307:   p->jacobian   = &OregoJacobian;
308:   p->solution   = &OregoSolution;
309:   p->final_time = 360;
310:   p->n          = 3;
311:   return(0);
312: }


315: /*
316: *  User-defined monitor for comparing to exact solutions when possible
317: */
318: typedef struct {
319:   MPI_Comm comm;
320:   Problem  problem;
321:   Vec      x;
322: } MonitorCtx;

326: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
327: {
329:   MonitorCtx     *mon = (MonitorCtx*)ctx;
330:   PetscReal      h,nrm_x,nrm_exact,nrm_diff;

333:   if (!mon->problem->solution) return(0);
334:   (*mon->problem->solution)(t,mon->x,mon->problem->data);
335:   VecNorm(x,NORM_2,&nrm_x);
336:   VecNorm(mon->x,NORM_2,&nrm_exact);
337:   VecAYPX(mon->x,-1,x);
338:   VecNorm(mon->x,NORM_2,&nrm_diff);
339:   TSGetTimeStep(ts,&h);
340:   if (step < 0) {
341:     PetscPrintf(mon->comm,"Interpolated final solution ");
342:   }
343:   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,(double)t,(double)h,(double)nrm_x,(double)nrm_exact,(double)nrm_diff);
344:   return(0);
345: }


350: int main(int argc,char **argv)
351: {
352:   PetscFunctionList plist = NULL;
353:   char              pname[256];
354:   TS                ts;            /* nonlinear solver */
355:   Vec               x,r;           /* solution, residual vectors */
356:   Mat               A;             /* Jacobian matrix */
357:   Problem           problem;
358:   PetscBool         use_monitor;
359:   PetscInt          steps,maxsteps = 1000,nonlinits,linits,snesfails,rejects;
360:   PetscReal         ftime;
361:   MonitorCtx        mon;
362:   PetscErrorCode    ierr;
363:   PetscMPIInt       size;

365:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
366:      Initialize program
367:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
368:   PetscInitialize(&argc,&argv,(char*)0,help);
369:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
370:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

372:   /* Register the available problems */
373:   PetscFunctionListAdd(&plist,"rober",&RoberCreate);
374:   PetscFunctionListAdd(&plist,"ce",&CECreate);
375:   PetscFunctionListAdd(&plist,"orego",&OregoCreate);
376:   PetscStrcpy(pname,"ce");

378:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
379:     Set runtime options
380:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
381:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
382:   {
383:     PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
384:     use_monitor = PETSC_FALSE;
385:     PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
386:   }
387:   PetscOptionsEnd();

389:   /* Create the new problem */
390:   PetscNew(&problem);
391:   problem->comm = MPI_COMM_WORLD;
392:   {
393:     PetscErrorCode (*pcreate)(Problem);

395:     PetscFunctionListFind(plist,pname,&pcreate);
396:     if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
397:     (*pcreate)(problem);
398:   }

400:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
401:     Create necessary matrix and vectors
402:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
403:   MatCreate(PETSC_COMM_WORLD,&A);
404:   MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
405:   MatSetFromOptions(A);
406:   MatSetUp(A);

408:   MatCreateVecs(A,&x,NULL);
409:   VecDuplicate(x,&r);

411:   mon.comm    = PETSC_COMM_WORLD;
412:   mon.problem = problem;
413:   VecDuplicate(x,&mon.x);

415:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
416:      Create timestepping solver context
417:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
418:   TSCreate(PETSC_COMM_WORLD,&ts);
419:   TSSetProblemType(ts,TS_NONLINEAR);
420:   TSSetType(ts,TSROSW); /* Rosenbrock-W */
421:   TSSetIFunction(ts,NULL,problem->function,problem->data);
422:   TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
423:   TSSetDuration(ts,maxsteps,problem->final_time);
424:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
425:   TSSetMaxStepRejections(ts,10);
426:   TSSetMaxSNESFailures(ts,-1); /* unlimited */
427:   if (use_monitor) {
428:     TSMonitorSet(ts,&MonitorError,&mon,NULL);
429:   }

431:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
432:      Set initial conditions
433:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
434:   (*problem->solution)(0,x,problem->data);
435:   TSSetInitialTimeStep(ts,0.0,.001);
436:   TSSetSolution(ts,x);

438:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
439:      Set runtime options
440:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
441:   TSSetFromOptions(ts);

443:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
444:      Solve nonlinear system
445:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
446:   TSSolve(ts,x);
447:   TSGetSolveTime(ts,&ftime);
448:   TSGetTimeStepNumber(ts,&steps);
449:   TSGetSNESFailures(ts,&snesfails);
450:   TSGetStepRejections(ts,&rejects);
451:   TSGetSNESIterations(ts,&nonlinits);
452:   TSGetKSPIterations(ts,&linits);
453:   PetscPrintf(PETSC_COMM_WORLD,"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",steps,rejects,snesfails,(double)ftime,nonlinits,linits);

455:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
456:      Free work space.  All PETSc objects should be destroyed when they
457:      are no longer needed.
458:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
459:   MatDestroy(&A);
460:   VecDestroy(&x);
461:   VecDestroy(&r);
462:   VecDestroy(&mon.x);
463:   TSDestroy(&ts);
464:   if (problem->destroy) {
465:     (*problem->destroy)(problem);
466:   }
467:   PetscFree(problem);
468:   PetscFunctionListDestroy(&plist);

470:   PetscFinalize();
471:   return(0);
472: }