Actual source code: ex8.c

petsc-3.12.5 2020-03-29
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: */
 31: static PetscErrorCode RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
 32: {
 33:   PetscErrorCode    ierr;
 34:   PetscScalar       *f;
 35:   const PetscScalar *x,*xdot;

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

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

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

 67:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 68:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 69:   if (A != B) {
 70:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 71:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 72:   }
 73:   return(0);
 74: }

 76: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
 77: {
 79:   PetscScalar    *x;

 82:   if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
 83:   VecGetArray(X,&x);
 84:   x[0] = 1;
 85:   x[1] = 0;
 86:   x[2] = 0;
 87:   VecRestoreArray(X,&x);
 88:   return(0);
 89: }

 91: static PetscErrorCode RoberCreate(Problem p)
 92: {

 95:   p->destroy    = 0;
 96:   p->function   = &RoberFunction;
 97:   p->jacobian   = &RoberJacobian;
 98:   p->solution   = &RoberSolution;
 99:   p->final_time = 1e11;
100:   p->n          = 3;
101:   return(0);
102: }

104: /*
105:      Stiff scalar valued problem
106: */

108: typedef struct {
109:   PetscReal lambda;
110: } CECtx;

112: static PetscErrorCode CEDestroy(Problem p)
113: {

117:   PetscFree(p->data);
118:   return(0);
119: }

121: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
122: {
123:   PetscErrorCode    ierr;
124:   PetscReal         l = ((CECtx*)ctx)->lambda;
125:   PetscScalar       *f;
126:   const PetscScalar *x,*xdot;

129:   VecGetArrayRead(X,&x);
130:   VecGetArrayRead(Xdot,&xdot);
131:   VecGetArray(F,&f);
132:   f[0] = xdot[0] + l*(x[0] - PetscCosReal(t));
133: #if 0
134:   PetscPrintf(PETSC_COMM_WORLD," f(t=%g,x=%g,xdot=%g) = %g\n",(double)t,(double)x[0],(double)xdot[0],(double)f[0]);
135: #endif
136:   VecRestoreArrayRead(X,&x);
137:   VecRestoreArrayRead(Xdot,&xdot);
138:   VecRestoreArray(F,&f);
139:   return(0);
140: }

142: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
143: {
144:   PetscReal         l = ((CECtx*)ctx)->lambda;
145:   PetscErrorCode    ierr;
146:   PetscInt          rowcol[] = {0};
147:   PetscScalar       J[1][1];
148:   const PetscScalar *x,*xdot;

151:   VecGetArrayRead(X,&x);
152:   VecGetArrayRead(Xdot,&xdot);
153:   J[0][0] = a + l;
154:   MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
155:   VecRestoreArrayRead(X,&x);
156:   VecRestoreArrayRead(Xdot,&xdot);

158:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
159:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
160:   if (A != B) {
161:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
162:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
163:   }
164:   return(0);
165: }

167: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
168: {
169:   PetscReal      l = ((CECtx*)ctx)->lambda;
171:   PetscScalar    *x;

174:   VecGetArray(X,&x);
175:   x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t);
176:   VecRestoreArray(X,&x);
177:   return(0);
178: }

180: static PetscErrorCode CECreate(Problem p)
181: {
183:   CECtx          *ce;

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

189:   p->destroy    = &CEDestroy;
190:   p->function   = &CEFunction;
191:   p->jacobian   = &CEJacobian;
192:   p->solution   = &CESolution;
193:   p->final_time = 10;
194:   p->n          = 1;
195:   p->hasexact   = PETSC_TRUE;

197:   ce->lambda = 10;
198:   PetscOptionsBegin(p->comm,NULL,"CE options","");
199:   {
200:     PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);
201:   }
202:   PetscOptionsEnd();
203:   return(0);
204: }

206: /*
207: *  Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
208: */
209: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
210: {
211:   PetscErrorCode    ierr;
212:   PetscScalar       *f;
213:   const PetscScalar *x,*xdot;

216:   VecGetArrayRead(X,&x);
217:   VecGetArrayRead(Xdot,&xdot);
218:   VecGetArray(F,&f);
219:   f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
220:   f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
221:   f[2] = xdot[2] - 0.161*(x[0] - x[2]);
222:   VecRestoreArrayRead(X,&x);
223:   VecRestoreArrayRead(Xdot,&xdot);
224:   VecRestoreArray(F,&f);
225:   return(0);
226: }

228: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
229: {
230:   PetscErrorCode    ierr;
231:   PetscInt          rowcol[] = {0,1,2};
232:   PetscScalar       J[3][3];
233:   const PetscScalar *x,*xdot;

236:   VecGetArrayRead(X,&x);
237:   VecGetArrayRead(Xdot,&xdot);
238:   J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
239:   J[0][1] = -77.27*(1. - x[0]);
240:   J[0][2] = 0;
241:   J[1][0] = 1./77.27*x[1];
242:   J[1][1] = a + 1./77.27*(1. + x[0]);
243:   J[1][2] = -1./77.27;
244:   J[2][0] = -0.161;
245:   J[2][1] = 0;
246:   J[2][2] = a + 0.161;
247:   MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
248:   VecRestoreArrayRead(X,&x);
249:   VecRestoreArrayRead(Xdot,&xdot);

251:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
252:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
253:   if (A != B) {
254:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
255:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
256:   }
257:   return(0);
258: }

260: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
261: {
263:   PetscScalar    *x;

266:   if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
267:   VecGetArray(X,&x);
268:   x[0] = 1;
269:   x[1] = 2;
270:   x[2] = 3;
271:   VecRestoreArray(X,&x);
272:   return(0);
273: }

275: static PetscErrorCode OregoCreate(Problem p)
276: {

279:   p->destroy    = 0;
280:   p->function   = &OregoFunction;
281:   p->jacobian   = &OregoJacobian;
282:   p->solution   = &OregoSolution;
283:   p->final_time = 360;
284:   p->n          = 3;
285:   return(0);
286: }


289: /*
290: *  User-defined monitor for comparing to exact solutions when possible
291: */
292: typedef struct {
293:   MPI_Comm comm;
294:   Problem  problem;
295:   Vec      x;
296: } MonitorCtx;

298: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
299: {
301:   MonitorCtx     *mon = (MonitorCtx*)ctx;
302:   PetscReal      h,nrm_x,nrm_exact,nrm_diff;

305:   if (!mon->problem->solution) return(0);
306:   (*mon->problem->solution)(t,mon->x,mon->problem->data);
307:   VecNorm(x,NORM_2,&nrm_x);
308:   VecNorm(mon->x,NORM_2,&nrm_exact);
309:   VecAYPX(mon->x,-1,x);
310:   VecNorm(mon->x,NORM_2,&nrm_diff);
311:   TSGetTimeStep(ts,&h);
312:   if (step < 0) {
313:     PetscPrintf(mon->comm,"Interpolated final solution ");
314:   }
315:   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);
316:   return(0);
317: }


320: int main(int argc,char **argv)
321: {
322:   PetscFunctionList plist = NULL;
323:   char              pname[256];
324:   TS                ts;            /* nonlinear solver */
325:   Vec               x,r;           /* solution, residual vectors */
326:   Mat               A;             /* Jacobian matrix */
327:   Problem           problem;
328:   PetscBool         use_monitor;
329:   PetscInt          steps,nonlinits,linits,snesfails,rejects;
330:   PetscReal         ftime;
331:   MonitorCtx        mon;
332:   PetscErrorCode    ierr;
333:   PetscMPIInt       size;

335:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
336:      Initialize program
337:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
338:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
339:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
340:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

342:   /* Register the available problems */
343:   PetscFunctionListAdd(&plist,"rober",&RoberCreate);
344:   PetscFunctionListAdd(&plist,"ce",&CECreate);
345:   PetscFunctionListAdd(&plist,"orego",&OregoCreate);
346:   PetscStrcpy(pname,"ce");

348:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349:     Set runtime options
350:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
351:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
352:   {
353:     PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
354:     use_monitor = PETSC_FALSE;
355:     PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
356:   }
357:   PetscOptionsEnd();

359:   /* Create the new problem */
360:   PetscNew(&problem);
361:   problem->comm = MPI_COMM_WORLD;
362:   {
363:     PetscErrorCode (*pcreate)(Problem);

365:     PetscFunctionListFind(plist,pname,&pcreate);
366:     if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
367:     (*pcreate)(problem);
368:   }

370:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371:     Create necessary matrix and vectors
372:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
373:   MatCreate(PETSC_COMM_WORLD,&A);
374:   MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
375:   MatSetFromOptions(A);
376:   MatSetUp(A);

378:   MatCreateVecs(A,&x,NULL);
379:   VecDuplicate(x,&r);

381:   mon.comm    = PETSC_COMM_WORLD;
382:   mon.problem = problem;
383:   VecDuplicate(x,&mon.x);

385:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
386:      Create timestepping solver context
387:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
388:   TSCreate(PETSC_COMM_WORLD,&ts);
389:   TSSetProblemType(ts,TS_NONLINEAR);
390:   TSSetType(ts,TSROSW); /* Rosenbrock-W */
391:   TSSetIFunction(ts,NULL,problem->function,problem->data);
392:   TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
393:   TSSetMaxTime(ts,problem->final_time);
394:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
395:   TSSetMaxStepRejections(ts,10);
396:   TSSetMaxSNESFailures(ts,-1); /* unlimited */
397:   if (use_monitor) {
398:     TSMonitorSet(ts,&MonitorError,&mon,NULL);
399:   }

401:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
402:      Set initial conditions
403:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
404:   (*problem->solution)(0,x,problem->data);
405:   TSSetTimeStep(ts,.001);
406:   TSSetSolution(ts,x);

408:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
409:      Set runtime options
410:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
411:   TSSetFromOptions(ts);

413:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414:      Solve nonlinear system
415:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
416:   TSSolve(ts,x);
417:   TSGetSolveTime(ts,&ftime);
418:   TSGetStepNumber(ts,&steps);
419:   TSGetSNESFailures(ts,&snesfails);
420:   TSGetStepRejections(ts,&rejects);
421:   TSGetSNESIterations(ts,&nonlinits);
422:   TSGetKSPIterations(ts,&linits);
423:   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);

425:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
426:      Free work space.  All PETSc objects should be destroyed when they
427:      are no longer needed.
428:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
429:   MatDestroy(&A);
430:   VecDestroy(&x);
431:   VecDestroy(&r);
432:   VecDestroy(&mon.x);
433:   TSDestroy(&ts);
434:   if (problem->destroy) {
435:     (*problem->destroy)(problem);
436:   }
437:   PetscFree(problem);
438:   PetscFunctionListDestroy(&plist);

440:   PetscFinalize();
441:   return ierr;
442: }

444: /*TEST

446:     test:
447:       requires: !complex
448:       args: -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex

450:     test:
451:       suffix: 2
452:       requires: !single !complex
453:       args: -ts_atol 1e-2 -ts_rtol 1e-2 -ts_max_time 15 -ts_type arkimex -ts_arkimex_type 2e -problem_type orego -ts_arkimex_initial_guess_extrapolate 0 -ts_adapt_time_step_increase_delay 4

455:     test:
456:       suffix: 3
457:       requires: !single !complex
458:       args: -ts_atol 1e-2 -ts_rtol 1e-2 -ts_max_time 15 -ts_type arkimex -ts_arkimex_type 2e -problem_type orego -ts_arkimex_initial_guess_extrapolate 1

460: TEST*/