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