1 static char help[] = "Nonlinear DAE benchmark problems.\n";
2
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>
13
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 };
26
27 /*
28 Stiff 3-variable system from chemical reactions, due to Robertson (1966), problem ROBER in Hairer&Wanner, ODE 2, 1996
29 https://archimede.uniba.it/~testset/report/rober.pdf
30 */
RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)31 static PetscErrorCode RoberFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
32 {
33 PetscScalar *f;
34 const PetscScalar *x, *xdot;
35
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 }
48
RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,PetscCtx ctx)49 static PetscErrorCode RoberJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
50 {
51 PetscInt rowcol[] = {0, 1, 2};
52 PetscScalar J[3][3];
53 const PetscScalar *x, *xdot;
54
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));
70
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 }
79
RoberSolution(PetscReal t,Vec X,PetscCtx ctx)80 static PetscErrorCode RoberSolution(PetscReal t, Vec X, PetscCtx ctx)
81 {
82 PetscScalar *x;
83
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 }
93
RoberCreate(Problem p)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 }
105
106 /*
107 Stiff scalar valued problem
108 */
109
110 typedef struct {
111 PetscReal lambda;
112 } CECtx;
113
CEDestroy(Problem p)114 static PetscErrorCode CEDestroy(Problem p)
115 {
116 PetscFunctionBeginUser;
117 PetscCall(PetscFree(p->data));
118 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120
CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)121 static PetscErrorCode CEFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
122 {
123 PetscReal l = ((CECtx *)ctx)->lambda;
124 PetscScalar *f;
125 const PetscScalar *x, *xdot;
126
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 }
140
CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,PetscCtx ctx)141 static PetscErrorCode CEJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
142 {
143 PetscReal l = ((CECtx *)ctx)->lambda;
144 PetscInt rowcol[] = {0};
145 PetscScalar J[1][1];
146 const PetscScalar *x, *xdot;
147
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));
155
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 }
164
CESolution(PetscReal t,Vec X,PetscCtx ctx)165 static PetscErrorCode CESolution(PetscReal t, Vec X, PetscCtx ctx)
166 {
167 PetscReal l = ((CECtx *)ctx)->lambda;
168 PetscScalar *x;
169
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 }
176
CECreate(Problem p)177 static PetscErrorCode CECreate(Problem p)
178 {
179 CECtx *ce;
180
181 PetscFunctionBeginUser;
182 PetscCall(PetscMalloc(sizeof(CECtx), &ce));
183 p->data = (void *)ce;
184
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;
192
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 }
201
202 /*
203 Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
204 */
OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)205 static PetscErrorCode OregoFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
206 {
207 PetscScalar *f;
208 const PetscScalar *x, *xdot;
209
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 }
222
OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,PetscCtx ctx)223 static PetscErrorCode OregoJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
224 {
225 PetscInt rowcol[] = {0, 1, 2};
226 PetscScalar J[3][3];
227 const PetscScalar *x, *xdot;
228
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));
244
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 }
253
OregoSolution(PetscReal t,Vec X,PetscCtx ctx)254 static PetscErrorCode OregoSolution(PetscReal t, Vec X, PetscCtx ctx)
255 {
256 PetscScalar *x;
257
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 }
267
OregoCreate(Problem p)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 }
279
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;
288
MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,PetscCtx ctx)289 static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal t, Vec x, PetscCtx ctx)
290 {
291 MonitorCtx *mon = (MonitorCtx *)ctx;
292 PetscReal h, nrm_x, nrm_exact, nrm_diff;
293
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 }
306
main(int argc,char ** argv)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;
321
322 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
323 Initialize program
324 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
325 PetscFunctionBeginUser;
326 PetscCall(PetscInitialize(&argc, &argv, NULL, 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");
329
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)));
335
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();
347
348 /* Create the new problem */
349 PetscCall(PetscNew(&problem));
350 problem->comm = MPI_COMM_WORLD;
351 {
352 PetscErrorCode (*pcreate)(Problem);
353
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 }
358
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));
366
367 PetscCall(MatCreateVecs(A, &x, NULL));
368 PetscCall(VecDuplicate(x, &r));
369
370 mon.comm = PETSC_COMM_WORLD;
371 mon.problem = problem;
372 PetscCall(VecDuplicate(x, &mon.x));
373
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));
387
388 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
389 Set initial conditions
390 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
391 PetscCall((*problem->solution)(0, x, problem->data));
392 PetscCall(TSSetTimeStep(ts, .001));
393 PetscCall(TSSetSolution(ts, x));
394
395 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
396 Set runtime options
397 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
398 PetscCall(TSSetFromOptions(ts));
399
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 }
413
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));
426
427 PetscCall(PetscFinalize());
428 return 0;
429 }
430
431 /*TEST
432
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
436
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
441
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
446
447 test:
448 suffix: 4
449 output_file: output/empty.out
450
451 test:
452 suffix: 5
453 args: -snes_lag_jacobian 20 -snes_lag_jacobian_persists
454 output_file: output/empty.out
455
456 TEST*/
457