Actual source code: ex8.c
1: static char help[] = "Nonlinear DAE benchmark problems.\n";
3: /*
4: Include "petscts.h" so that we can use TS solvers. Note that this
5: file automatically includes:
6: petscsys.h - base PETSc routines petscvec.h - vectors
7: petscmat.h - matrices
8: petscis.h - index sets petscksp.h - Krylov subspace methods
9: petscviewer.h - viewers petscpc.h - preconditioners
10: petscksp.h - linear solvers
11: */
12: #include <petscts.h>
14: typedef struct _Problem *Problem;
15: struct _Problem {
16: PetscErrorCode (*destroy)(Problem);
17: TSIFunctionFn *function;
18: TSIJacobianFn *jacobian;
19: PetscErrorCode (*solution)(PetscReal, Vec, void *);
20: MPI_Comm comm;
21: PetscReal final_time;
22: PetscInt n;
23: PetscBool hasexact;
24: void *data;
25: };
27: /*
28: Stiff 3-variable system from chemical reactions, due to Robertson (1966), problem ROBER in Hairer&Wanner, ODE 2, 1996
29: */
30: static PetscErrorCode RoberFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
31: {
32: PetscScalar *f;
33: const PetscScalar *x, *xdot;
35: PetscFunctionBeginUser;
36: PetscCall(VecGetArrayRead(X, &x));
37: PetscCall(VecGetArrayRead(Xdot, &xdot));
38: PetscCall(VecGetArray(F, &f));
39: f[0] = xdot[0] + 0.04 * x[0] - 1e4 * x[1] * x[2];
40: f[1] = xdot[1] - 0.04 * x[0] + 1e4 * x[1] * x[2] + 3e7 * PetscSqr(x[1]);
41: f[2] = xdot[2] - 3e7 * PetscSqr(x[1]);
42: PetscCall(VecRestoreArrayRead(X, &x));
43: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
44: PetscCall(VecRestoreArray(F, &f));
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
48: static PetscErrorCode RoberJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
49: {
50: PetscInt rowcol[] = {0, 1, 2};
51: PetscScalar J[3][3];
52: const PetscScalar *x, *xdot;
54: PetscFunctionBeginUser;
55: PetscCall(VecGetArrayRead(X, &x));
56: PetscCall(VecGetArrayRead(Xdot, &xdot));
57: J[0][0] = a + 0.04;
58: J[0][1] = -1e4 * x[2];
59: J[0][2] = -1e4 * x[1];
60: J[1][0] = -0.04;
61: J[1][1] = a + 1e4 * x[2] + 3e7 * 2 * x[1];
62: J[1][2] = 1e4 * x[1];
63: J[2][0] = 0;
64: J[2][1] = -3e7 * 2 * x[1];
65: J[2][2] = a;
66: PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
67: PetscCall(VecRestoreArrayRead(X, &x));
68: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
70: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
71: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
72: if (A != B) {
73: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
74: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
75: }
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode RoberSolution(PetscReal t, Vec X, void *ctx)
80: {
81: PetscScalar *x;
83: PetscFunctionBeginUser;
84: PetscCheck(t == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP, "not implemented");
85: PetscCall(VecGetArray(X, &x));
86: x[0] = 1;
87: x[1] = 0;
88: x[2] = 0;
89: PetscCall(VecRestoreArray(X, &x));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode RoberCreate(Problem p)
94: {
95: PetscFunctionBeginUser;
96: p->destroy = 0;
97: p->function = &RoberFunction;
98: p->jacobian = &RoberJacobian;
99: p->solution = &RoberSolution;
100: p->final_time = 1e11;
101: p->n = 3;
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: /*
106: Stiff scalar valued problem
107: */
109: typedef struct {
110: PetscReal lambda;
111: } CECtx;
113: static PetscErrorCode CEDestroy(Problem p)
114: {
115: PetscFunctionBeginUser;
116: PetscCall(PetscFree(p->data));
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: static PetscErrorCode CEFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
121: {
122: PetscReal l = ((CECtx *)ctx)->lambda;
123: PetscScalar *f;
124: const PetscScalar *x, *xdot;
126: PetscFunctionBeginUser;
127: PetscCall(VecGetArrayRead(X, &x));
128: PetscCall(VecGetArrayRead(Xdot, &xdot));
129: PetscCall(VecGetArray(F, &f));
130: f[0] = xdot[0] + l * (x[0] - PetscCosReal(t));
131: #if 0
132: 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]));
133: #endif
134: PetscCall(VecRestoreArrayRead(X, &x));
135: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
136: PetscCall(VecRestoreArray(F, &f));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: static PetscErrorCode CEJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
141: {
142: PetscReal l = ((CECtx *)ctx)->lambda;
143: PetscInt rowcol[] = {0};
144: PetscScalar J[1][1];
145: const PetscScalar *x, *xdot;
147: PetscFunctionBeginUser;
148: PetscCall(VecGetArrayRead(X, &x));
149: PetscCall(VecGetArrayRead(Xdot, &xdot));
150: J[0][0] = a + l;
151: PetscCall(MatSetValues(B, 1, rowcol, 1, rowcol, &J[0][0], INSERT_VALUES));
152: PetscCall(VecRestoreArrayRead(X, &x));
153: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
155: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
156: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
157: if (A != B) {
158: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
159: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
160: }
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: static PetscErrorCode CESolution(PetscReal t, Vec X, void *ctx)
165: {
166: PetscReal l = ((CECtx *)ctx)->lambda;
167: PetscScalar *x;
169: PetscFunctionBeginUser;
170: PetscCall(VecGetArray(X, &x));
171: x[0] = l / (l * l + 1) * (l * PetscCosReal(t) + PetscSinReal(t)) - l * l / (l * l + 1) * PetscExpReal(-l * t);
172: PetscCall(VecRestoreArray(X, &x));
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: static PetscErrorCode CECreate(Problem p)
177: {
178: CECtx *ce;
180: PetscFunctionBeginUser;
181: PetscCall(PetscMalloc(sizeof(CECtx), &ce));
182: p->data = (void *)ce;
184: p->destroy = &CEDestroy;
185: p->function = &CEFunction;
186: p->jacobian = &CEJacobian;
187: p->solution = &CESolution;
188: p->final_time = 10;
189: p->n = 1;
190: p->hasexact = PETSC_TRUE;
192: ce->lambda = 10;
193: PetscOptionsBegin(p->comm, NULL, "CE options", "");
194: {
195: PetscCall(PetscOptionsReal("-problem_ce_lambda", "Parameter controlling stiffness: xdot + lambda*(x - cos(t))", "", ce->lambda, &ce->lambda, NULL));
196: }
197: PetscOptionsEnd();
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: /*
202: Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
203: */
204: static PetscErrorCode OregoFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
205: {
206: PetscScalar *f;
207: const PetscScalar *x, *xdot;
209: PetscFunctionBeginUser;
210: PetscCall(VecGetArrayRead(X, &x));
211: PetscCall(VecGetArrayRead(Xdot, &xdot));
212: PetscCall(VecGetArray(F, &f));
213: f[0] = xdot[0] - 77.27 * (x[1] + x[0] * (1. - 8.375e-6 * x[0] - x[1]));
214: f[1] = xdot[1] - 1 / 77.27 * (x[2] - (1. + x[0]) * x[1]);
215: f[2] = xdot[2] - 0.161 * (x[0] - x[2]);
216: PetscCall(VecRestoreArrayRead(X, &x));
217: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
218: PetscCall(VecRestoreArray(F, &f));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode OregoJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
223: {
224: PetscInt rowcol[] = {0, 1, 2};
225: PetscScalar J[3][3];
226: const PetscScalar *x, *xdot;
228: PetscFunctionBeginUser;
229: PetscCall(VecGetArrayRead(X, &x));
230: PetscCall(VecGetArrayRead(Xdot, &xdot));
231: J[0][0] = a - 77.27 * ((1. - 8.375e-6 * x[0] - x[1]) - 8.375e-6 * x[0]);
232: J[0][1] = -77.27 * (1. - x[0]);
233: J[0][2] = 0;
234: J[1][0] = 1. / 77.27 * x[1];
235: J[1][1] = a + 1. / 77.27 * (1. + x[0]);
236: J[1][2] = -1. / 77.27;
237: J[2][0] = -0.161;
238: J[2][1] = 0;
239: J[2][2] = a + 0.161;
240: PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
241: PetscCall(VecRestoreArrayRead(X, &x));
242: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
244: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
245: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
246: if (A != B) {
247: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
248: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
249: }
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: static PetscErrorCode OregoSolution(PetscReal t, Vec X, void *ctx)
254: {
255: PetscScalar *x;
257: PetscFunctionBeginUser;
258: PetscCheck(t == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP, "not implemented");
259: PetscCall(VecGetArray(X, &x));
260: x[0] = 1;
261: x[1] = 2;
262: x[2] = 3;
263: PetscCall(VecRestoreArray(X, &x));
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static PetscErrorCode OregoCreate(Problem p)
268: {
269: PetscFunctionBeginUser;
270: p->destroy = 0;
271: p->function = &OregoFunction;
272: p->jacobian = &OregoJacobian;
273: p->solution = &OregoSolution;
274: p->final_time = 360;
275: p->n = 3;
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*
280: User-defined monitor for comparing to exact solutions when possible
281: */
282: typedef struct {
283: MPI_Comm comm;
284: Problem problem;
285: Vec x;
286: } MonitorCtx;
288: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal t, Vec x, void *ctx)
289: {
290: MonitorCtx *mon = (MonitorCtx *)ctx;
291: PetscReal h, nrm_x, nrm_exact, nrm_diff;
293: PetscFunctionBeginUser;
294: if (!mon->problem->solution) PetscFunctionReturn(PETSC_SUCCESS);
295: PetscCall((*mon->problem->solution)(t, mon->x, mon->problem->data));
296: PetscCall(VecNorm(x, NORM_2, &nrm_x));
297: PetscCall(VecNorm(mon->x, NORM_2, &nrm_exact));
298: PetscCall(VecAYPX(mon->x, -1, x));
299: PetscCall(VecNorm(mon->x, NORM_2, &nrm_diff));
300: PetscCall(TSGetTimeStep(ts, &h));
301: if (step < 0) PetscCall(PetscPrintf(mon->comm, "Interpolated final solution "));
302: 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));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: int main(int argc, char **argv)
307: {
308: PetscFunctionList plist = NULL;
309: char pname[256];
310: TS ts; /* nonlinear solver */
311: Vec x, r; /* solution, residual vectors */
312: Mat A; /* Jacobian matrix */
313: Problem problem;
314: PetscBool use_monitor = PETSC_FALSE;
315: PetscBool use_result = PETSC_FALSE;
316: PetscInt steps, nonlinits, linits, snesfails, rejects;
317: PetscReal ftime;
318: MonitorCtx mon;
319: PetscMPIInt size;
321: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
322: Initialize program
323: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
324: PetscFunctionBeginUser;
325: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
326: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
327: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
329: /* Register the available problems */
330: PetscCall(PetscFunctionListAdd(&plist, "rober", &RoberCreate));
331: PetscCall(PetscFunctionListAdd(&plist, "ce", &CECreate));
332: PetscCall(PetscFunctionListAdd(&plist, "orego", &OregoCreate));
333: PetscCall(PetscStrncpy(pname, "ce", sizeof(pname)));
335: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
336: Set runtime options
337: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
338: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Timestepping benchmark options", "");
339: {
340: PetscCall(PetscOptionsFList("-problem_type", "Name of problem to run", "", plist, pname, pname, sizeof(pname), NULL));
341: use_monitor = PETSC_FALSE;
342: PetscCall(PetscOptionsBool("-monitor_error", "Display errors relative to exact solutions", "", use_monitor, &use_monitor, NULL));
343: PetscCall(PetscOptionsBool("-monitor_result", "Display result", "", use_result, &use_result, NULL));
344: }
345: PetscOptionsEnd();
347: /* Create the new problem */
348: PetscCall(PetscNew(&problem));
349: problem->comm = MPI_COMM_WORLD;
350: {
351: PetscErrorCode (*pcreate)(Problem);
353: PetscCall(PetscFunctionListFind(plist, pname, &pcreate));
354: PetscCheck(pcreate, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "No problem '%s'", pname);
355: PetscCall((*pcreate)(problem));
356: }
358: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
359: Create necessary matrix and vectors
360: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
361: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
362: PetscCall(MatSetSizes(A, problem->n, problem->n, PETSC_DETERMINE, PETSC_DETERMINE));
363: PetscCall(MatSetFromOptions(A));
364: PetscCall(MatSetUp(A));
366: PetscCall(MatCreateVecs(A, &x, NULL));
367: PetscCall(VecDuplicate(x, &r));
369: mon.comm = PETSC_COMM_WORLD;
370: mon.problem = problem;
371: PetscCall(VecDuplicate(x, &mon.x));
373: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
374: Create timestepping solver context
375: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
376: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
377: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
378: PetscCall(TSSetType(ts, TSROSW)); /* Rosenbrock-W */
379: PetscCall(TSSetIFunction(ts, NULL, problem->function, problem->data));
380: PetscCall(TSSetIJacobian(ts, A, A, problem->jacobian, problem->data));
381: PetscCall(TSSetMaxTime(ts, problem->final_time));
382: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
383: PetscCall(TSSetMaxStepRejections(ts, 10));
384: PetscCall(TSSetMaxSNESFailures(ts, -1)); /* unlimited */
385: if (use_monitor) PetscCall(TSMonitorSet(ts, &MonitorError, &mon, NULL));
387: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: Set initial conditions
389: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
390: PetscCall((*problem->solution)(0, x, problem->data));
391: PetscCall(TSSetTimeStep(ts, .001));
392: PetscCall(TSSetSolution(ts, x));
394: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
395: Set runtime options
396: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
397: PetscCall(TSSetFromOptions(ts));
399: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
400: Solve nonlinear system
401: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
402: PetscCall(TSSolve(ts, x));
403: PetscCall(TSGetSolveTime(ts, &ftime));
404: PetscCall(TSGetStepNumber(ts, &steps));
405: PetscCall(TSGetSNESFailures(ts, &snesfails));
406: PetscCall(TSGetStepRejections(ts, &rejects));
407: PetscCall(TSGetSNESIterations(ts, &nonlinits));
408: PetscCall(TSGetKSPIterations(ts, &linits));
409: if (use_result) {
410: 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));
411: }
413: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414: Free work space. All PETSc objects should be destroyed when they
415: are no longer needed.
416: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
417: PetscCall(MatDestroy(&A));
418: PetscCall(VecDestroy(&x));
419: PetscCall(VecDestroy(&r));
420: PetscCall(VecDestroy(&mon.x));
421: PetscCall(TSDestroy(&ts));
422: if (problem->destroy) PetscCall((*problem->destroy)(problem));
423: PetscCall(PetscFree(problem));
424: PetscCall(PetscFunctionListDestroy(&plist));
426: PetscCall(PetscFinalize());
427: return 0;
428: }
430: /*TEST
432: test:
433: requires: !complex
434: args: -monitor_result -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex
436: test:
437: suffix: 2
438: requires: !single !complex
439: 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
441: test:
442: suffix: 3
443: requires: !single !complex
444: 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
446: test:
447: suffix: 4
449: test:
450: suffix: 5
451: args: -snes_lag_jacobian 20 -snes_lag_jacobian_persists
453: TEST*/