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