Actual source code: ex8.c

petsc-3.6.1 2015-07-22
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:   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);
341:   return(0);
342: }


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

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

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

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

386:   /* Create the new problem */
387:   PetscNew(&problem);
388:   problem->comm = MPI_COMM_WORLD;
389:   {
390:     PetscErrorCode (*pcreate)(Problem);

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

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

405:   MatCreateVecs(A,&x,NULL);
406:   VecDuplicate(x,&r);

408:   mon.comm    = PETSC_COMM_WORLD;
409:   mon.problem = problem;
410:   VecDuplicate(x,&mon.x);

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

427:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
428:      Set initial conditions
429:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
430:   (*problem->solution)(0,x,problem->data);
431:   TSSetInitialTimeStep(ts,0.0,.001);
432:   TSSetSolution(ts,x);

434:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
435:      Set runtime options
436:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
437:   TSSetFromOptions(ts);

439:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
440:      Solve nonlinear system
441:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
442:   TSSolve(ts,x);
443:   TSGetSolveTime(ts,&ftime);
444:   TSGetTimeStepNumber(ts,&steps);
445:   TSGetSNESFailures(ts,&snesfails);
446:   TSGetStepRejections(ts,&rejects);
447:   TSGetSNESIterations(ts,&nonlinits);
448:   TSGetKSPIterations(ts,&linits);
449:   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);
450:   if (problem->hasexact) {
451:     MonitorError(ts,steps,ftime,x,&mon);
452:   }

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

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