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*/