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;

 36:   PetscFunctionBeginUser;
 37:   PetscCall(VecGetArrayRead(X, &x));
 38:   PetscCall(VecGetArrayRead(Xdot, &xdot));
 39:   PetscCall(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:   PetscCall(VecRestoreArrayRead(X, &x));
 44:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
 45:   PetscCall(VecRestoreArray(F, &f));
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

 55:   PetscFunctionBeginUser;
 56:   PetscCall(VecGetArrayRead(X, &x));
 57:   PetscCall(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:   PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
 68:   PetscCall(VecRestoreArrayRead(X, &x));
 69:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));

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

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

 84:   PetscFunctionBeginUser;
 85:   PetscCheck(t == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP, "not implemented");
 86:   PetscCall(VecGetArray(X, &x));
 87:   x[0] = 1;
 88:   x[1] = 0;
 89:   x[2] = 0;
 90:   PetscCall(VecRestoreArray(X, &x));
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

 94: static PetscErrorCode RoberCreate(Problem p)
 95: {
 96:   PetscFunctionBeginUser;
 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:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

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

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

114: static PetscErrorCode CEDestroy(Problem p)
115: {
116:   PetscFunctionBeginUser;
117:   PetscCall(PetscFree(p->data));
118:   PetscFunctionReturn(PETSC_SUCCESS);
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;

127:   PetscFunctionBeginUser;
128:   PetscCall(VecGetArrayRead(X, &x));
129:   PetscCall(VecGetArrayRead(Xdot, &xdot));
130:   PetscCall(VecGetArray(F, &f));
131:   f[0] = xdot[0] + l * (x[0] - PetscCosReal(t));
132: #if 0
133:   PetscCall(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:   PetscCall(VecRestoreArrayRead(X, &x));
136:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
137:   PetscCall(VecRestoreArray(F, &f));
138:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

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

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

181:   PetscFunctionBeginUser;
182:   PetscCall(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:     PetscCall(PetscOptionsReal("-problem_ce_lambda", "Parameter controlling stiffness: xdot + lambda*(x - cos(t))", "", ce->lambda, &ce->lambda, NULL));
197:   }
198:   PetscOptionsEnd();
199:   PetscFunctionReturn(PETSC_SUCCESS);
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;

210:   PetscFunctionBeginUser;
211:   PetscCall(VecGetArrayRead(X, &x));
212:   PetscCall(VecGetArrayRead(Xdot, &xdot));
213:   PetscCall(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:   PetscCall(VecRestoreArrayRead(X, &x));
218:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
219:   PetscCall(VecRestoreArray(F, &f));
220:   PetscFunctionReturn(PETSC_SUCCESS);
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;

229:   PetscFunctionBeginUser;
230:   PetscCall(VecGetArrayRead(X, &x));
231:   PetscCall(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:   PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
242:   PetscCall(VecRestoreArrayRead(X, &x));
243:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));

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

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

258:   PetscFunctionBeginUser;
259:   PetscCheck(t == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP, "not implemented");
260:   PetscCall(VecGetArray(X, &x));
261:   x[0] = 1;
262:   x[1] = 2;
263:   x[2] = 3;
264:   PetscCall(VecRestoreArray(X, &x));
265:   PetscFunctionReturn(PETSC_SUCCESS);
266: }

268: static PetscErrorCode OregoCreate(Problem p)
269: {
270:   PetscFunctionBeginUser;
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:   PetscFunctionReturn(PETSC_SUCCESS);
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;

294:   PetscFunctionBeginUser;
295:   if (!mon->problem->solution) PetscFunctionReturn(PETSC_SUCCESS);
296:   PetscCall((*mon->problem->solution)(t, mon->x, mon->problem->data));
297:   PetscCall(VecNorm(x, NORM_2, &nrm_x));
298:   PetscCall(VecNorm(mon->x, NORM_2, &nrm_exact));
299:   PetscCall(VecAYPX(mon->x, -1, x));
300:   PetscCall(VecNorm(mon->x, NORM_2, &nrm_diff));
301:   PetscCall(TSGetTimeStep(ts, &h));
302:   if (step < 0) PetscCall(PetscPrintf(mon->comm, "Interpolated final solution "));
303:   PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
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:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
325:   PetscFunctionBeginUser;
326:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
327:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
328:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

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

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

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

354:     PetscCall(PetscFunctionListFind(plist, pname, &pcreate));
355:     PetscCheck(pcreate, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "No problem '%s'", pname);
356:     PetscCall((*pcreate)(problem));
357:   }

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

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

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

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

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

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

400:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
401:      Solve nonlinear system
402:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
403:   PetscCall(TSSolve(ts, x));
404:   PetscCall(TSGetSolveTime(ts, &ftime));
405:   PetscCall(TSGetStepNumber(ts, &steps));
406:   PetscCall(TSGetSNESFailures(ts, &snesfails));
407:   PetscCall(TSGetStepRejections(ts, &rejects));
408:   PetscCall(TSGetSNESIterations(ts, &nonlinits));
409:   PetscCall(TSGetKSPIterations(ts, &linits));
410:   if (use_result) {
411:     PetscCall(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:   PetscCall(MatDestroy(&A));
419:   PetscCall(VecDestroy(&x));
420:   PetscCall(VecDestroy(&r));
421:   PetscCall(VecDestroy(&mon.x));
422:   PetscCall(TSDestroy(&ts));
423:   if (problem->destroy) PetscCall((*problem->destroy)(problem));
424:   PetscCall(PetscFree(problem));
425:   PetscCall(PetscFunctionListDestroy(&plist));

427:   PetscCall(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*/