Actual source code: ex8.c


  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:   PetscScalar       *f;
 34:   const PetscScalar *x, *xdot;

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

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

 56:   VecGetArrayRead(X, &x);
 57:   VecGetArrayRead(Xdot, &xdot);
 58:   J[0][0] = a + 0.04;
 59:   J[0][1] = -1e4 * x[2];
 60:   J[0][2] = -1e4 * x[1];
 61:   J[1][0] = -0.04;
 62:   J[1][1] = a + 1e4 * x[2] + 3e7 * 2 * x[1];
 63:   J[1][2] = 1e4 * x[1];
 64:   J[2][0] = 0;
 65:   J[2][1] = -3e7 * 2 * x[1];
 66:   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: }

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

 86:   VecGetArray(X, &x);
 87:   x[0] = 1;
 88:   x[1] = 0;
 89:   x[2] = 0;
 90:   VecRestoreArray(X, &x);
 91:   return 0;
 92: }

 94: static PetscErrorCode RoberCreate(Problem p)
 95: {
 97:   p->destroy    = 0;
 98:   p->function   = &RoberFunction;
 99:   p->jacobian   = &RoberJacobian;
100:   p->solution   = &RoberSolution;
101:   p->final_time = 1e11;
102:   p->n          = 3;
103:   return 0;
104: }

106: /*
107:      Stiff scalar valued problem
108: */

110: typedef struct {
111:   PetscReal lambda;
112: } CECtx;

114: static PetscErrorCode CEDestroy(Problem p)
115: {
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:   PetscReal          l = ((CECtx *)ctx)->lambda;
124:   PetscScalar       *f;
125:   const PetscScalar *x, *xdot;

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

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

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

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

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

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

177: static PetscErrorCode CECreate(Problem p)
178: {
179:   CECtx *ce;

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

185:   p->destroy    = &CEDestroy;
186:   p->function   = &CEFunction;
187:   p->jacobian   = &CEJacobian;
188:   p->solution   = &CESolution;
189:   p->final_time = 10;
190:   p->n          = 1;
191:   p->hasexact   = PETSC_TRUE;

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

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

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

223: static PetscErrorCode OregoJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
224: {
225:   PetscInt           rowcol[] = {0, 1, 2};
226:   PetscScalar        J[3][3];
227:   const PetscScalar *x, *xdot;

230:   VecGetArrayRead(X, &x);
231:   VecGetArrayRead(Xdot, &xdot);
232:   J[0][0] = a - 77.27 * ((1. - 8.375e-6 * x[0] - x[1]) - 8.375e-6 * x[0]);
233:   J[0][1] = -77.27 * (1. - x[0]);
234:   J[0][2] = 0;
235:   J[1][0] = 1. / 77.27 * x[1];
236:   J[1][1] = a + 1. / 77.27 * (1. + x[0]);
237:   J[1][2] = -1. / 77.27;
238:   J[2][0] = -0.161;
239:   J[2][1] = 0;
240:   J[2][2] = a + 0.161;
241:   MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES);
242:   VecRestoreArrayRead(X, &x);
243:   VecRestoreArrayRead(Xdot, &xdot);

245:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
246:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
247:   if (A != B) {
248:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
249:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
250:   }
251:   return 0;
252: }

254: static PetscErrorCode OregoSolution(PetscReal t, Vec X, void *ctx)
255: {
256:   PetscScalar *x;

260:   VecGetArray(X, &x);
261:   x[0] = 1;
262:   x[1] = 2;
263:   x[2] = 3;
264:   VecRestoreArray(X, &x);
265:   return 0;
266: }

268: static PetscErrorCode OregoCreate(Problem p)
269: {
271:   p->destroy    = 0;
272:   p->function   = &OregoFunction;
273:   p->jacobian   = &OregoJacobian;
274:   p->solution   = &OregoSolution;
275:   p->final_time = 360;
276:   p->n          = 3;
277:   return 0;
278: }

280: /*
281:    User-defined monitor for comparing to exact solutions when possible
282: */
283: typedef struct {
284:   MPI_Comm comm;
285:   Problem  problem;
286:   Vec      x;
287: } MonitorCtx;

289: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal t, Vec x, void *ctx)
290: {
291:   MonitorCtx *mon = (MonitorCtx *)ctx;
292:   PetscReal   h, nrm_x, nrm_exact, nrm_diff;

295:   if (!mon->problem->solution) return 0;
296:   (*mon->problem->solution)(t, mon->x, mon->problem->data);
297:   VecNorm(x, NORM_2, &nrm_x);
298:   VecNorm(mon->x, NORM_2, &nrm_exact);
299:   VecAYPX(mon->x, -1, x);
300:   VecNorm(mon->x, NORM_2, &nrm_diff);
301:   TSGetTimeStep(ts, &h);
302:   if (step < 0) PetscPrintf(mon->comm, "Interpolated final solution ");
303:   PetscPrintf(mon->comm, "step %4" PetscInt_FMT " 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);
304:   return 0;
305: }

307: int main(int argc, char **argv)
308: {
309:   PetscFunctionList plist = NULL;
310:   char              pname[256];
311:   TS                ts;   /* nonlinear solver */
312:   Vec               x, r; /* solution, residual vectors */
313:   Mat               A;    /* Jacobian matrix */
314:   Problem           problem;
315:   PetscBool         use_monitor = PETSC_FALSE;
316:   PetscBool         use_result  = PETSC_FALSE;
317:   PetscInt          steps, nonlinits, linits, snesfails, rejects;
318:   PetscReal         ftime;
319:   MonitorCtx        mon;
320:   PetscMPIInt       size;

322:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
323:      Initialize program
324:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
326:   PetscInitialize(&argc, &argv, (char *)0, help);
327:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

330:   /* Register the available problems */
331:   PetscFunctionListAdd(&plist, "rober", &RoberCreate);
332:   PetscFunctionListAdd(&plist, "ce", &CECreate);
333:   PetscFunctionListAdd(&plist, "orego", &OregoCreate);
334:   PetscStrcpy(pname, "ce");

336:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
337:     Set runtime options
338:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
339:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Timestepping benchmark options", "");
340:   {
341:     PetscOptionsFList("-problem_type", "Name of problem to run", "", plist, pname, pname, sizeof(pname), NULL);
342:     use_monitor = PETSC_FALSE;
343:     PetscOptionsBool("-monitor_error", "Display errors relative to exact solutions", "", use_monitor, &use_monitor, NULL);
344:     PetscOptionsBool("-monitor_result", "Display result", "", use_result, &use_result, NULL);
345:   }
346:   PetscOptionsEnd();

348:   /* Create the new problem */
349:   PetscNew(&problem);
350:   problem->comm = MPI_COMM_WORLD;
351:   {
352:     PetscErrorCode (*pcreate)(Problem);

354:     PetscFunctionListFind(plist, pname, &pcreate);
356:     (*pcreate)(problem);
357:   }

359:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
360:     Create necessary matrix and vectors
361:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
362:   MatCreate(PETSC_COMM_WORLD, &A);
363:   MatSetSizes(A, problem->n, problem->n, PETSC_DETERMINE, PETSC_DETERMINE);
364:   MatSetFromOptions(A);
365:   MatSetUp(A);

367:   MatCreateVecs(A, &x, NULL);
368:   VecDuplicate(x, &r);

370:   mon.comm    = PETSC_COMM_WORLD;
371:   mon.problem = problem;
372:   VecDuplicate(x, &mon.x);

374:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
375:      Create timestepping solver context
376:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
377:   TSCreate(PETSC_COMM_WORLD, &ts);
378:   TSSetProblemType(ts, TS_NONLINEAR);
379:   TSSetType(ts, TSROSW); /* Rosenbrock-W */
380:   TSSetIFunction(ts, NULL, problem->function, problem->data);
381:   TSSetIJacobian(ts, A, A, problem->jacobian, problem->data);
382:   TSSetMaxTime(ts, problem->final_time);
383:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
384:   TSSetMaxStepRejections(ts, 10);
385:   TSSetMaxSNESFailures(ts, -1); /* unlimited */
386:   if (use_monitor) TSMonitorSet(ts, &MonitorError, &mon, NULL);

388:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
389:      Set initial conditions
390:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
391:   (*problem->solution)(0, x, problem->data);
392:   TSSetTimeStep(ts, .001);
393:   TSSetSolution(ts, x);

395:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
396:      Set runtime options
397:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
398:   TSSetFromOptions(ts);

400:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
401:      Solve nonlinear system
402:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
403:   TSSolve(ts, x);
404:   TSGetSolveTime(ts, &ftime);
405:   TSGetStepNumber(ts, &steps);
406:   TSGetSNESFailures(ts, &snesfails);
407:   TSGetStepRejections(ts, &rejects);
408:   TSGetSNESIterations(ts, &nonlinits);
409:   TSGetKSPIterations(ts, &linits);
410:   if (use_result) {
411:     PetscPrintf(PETSC_COMM_WORLD, "steps %" PetscInt_FMT " (%" PetscInt_FMT " rejected, %" PetscInt_FMT " SNES fails), ftime %g, nonlinits %" PetscInt_FMT ", linits %" PetscInt_FMT "\n", steps, rejects, snesfails, (double)ftime, nonlinits, linits);
412:   }

414:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
415:      Free work space.  All PETSc objects should be destroyed when they
416:      are no longer needed.
417:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
418:   MatDestroy(&A);
419:   VecDestroy(&x);
420:   VecDestroy(&r);
421:   VecDestroy(&mon.x);
422:   TSDestroy(&ts);
423:   if (problem->destroy) (*problem->destroy)(problem);
424:   PetscFree(problem);
425:   PetscFunctionListDestroy(&plist);

427:   PetscFinalize();
428:   return 0;
429: }

431: /*TEST

433:     test:
434:       requires: !complex
435:       args:  -monitor_result -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex

437:     test:
438:       suffix: 2
439:       requires: !single !complex
440:       args: -monitor_result -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

442:     test:
443:       suffix: 3
444:       requires: !single !complex
445:       args: -monitor_result -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

447:     test:
448:       suffix: 4

450:     test:
451:       suffix: 5
452:       args: -snes_lag_jacobian 20 -snes_lag_jacobian_persists

454: TEST*/