Actual source code: ex8.c
petsc-3.14.6 2021-03-30
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: PetscErrorCode ierr;
34: PetscScalar *f;
35: const PetscScalar *x,*xdot;
38: VecGetArrayRead(X,&x);
39: VecGetArrayRead(Xdot,&xdot);
40: 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: VecRestoreArrayRead(X,&x);
45: VecRestoreArrayRead(Xdot,&xdot);
46: VecRestoreArray(F,&f);
47: return(0);
48: }
50: static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
51: {
52: PetscErrorCode ierr;
53: PetscInt rowcol[] = {0,1,2};
54: PetscScalar J[3][3];
55: const PetscScalar *x,*xdot;
58: VecGetArrayRead(X,&x);
59: VecGetArrayRead(Xdot,&xdot);
60: J[0][0] = a + 0.04; J[0][1] = -1e4*x[2]; J[0][2] = -1e4*x[1];
61: J[1][0] = -0.04; J[1][1] = a + 1e4*x[2] + 3e7*2*x[1]; J[1][2] = 1e4*x[1];
62: J[2][0] = 0; J[2][1] = -3e7*2*x[1]; J[2][2] = a;
63: MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
64: VecRestoreArrayRead(X,&x);
65: VecRestoreArrayRead(Xdot,&xdot);
67: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
69: if (A != B) {
70: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
71: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
72: }
73: return(0);
74: }
76: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
77: {
79: PetscScalar *x;
82: if (t != 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not implemented");
83: VecGetArray(X,&x);
84: x[0] = 1;
85: x[1] = 0;
86: x[2] = 0;
87: VecRestoreArray(X,&x);
88: return(0);
89: }
91: static PetscErrorCode RoberCreate(Problem p)
92: {
95: p->destroy = 0;
96: p->function = &RoberFunction;
97: p->jacobian = &RoberJacobian;
98: p->solution = &RoberSolution;
99: p->final_time = 1e11;
100: p->n = 3;
101: return(0);
102: }
104: /*
105: Stiff scalar valued problem
106: */
108: typedef struct {
109: PetscReal lambda;
110: } CECtx;
112: static PetscErrorCode CEDestroy(Problem p)
113: {
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: PetscErrorCode ierr;
124: PetscReal l = ((CECtx*)ctx)->lambda;
125: PetscScalar *f;
126: const PetscScalar *x,*xdot;
129: VecGetArrayRead(X,&x);
130: VecGetArrayRead(Xdot,&xdot);
131: VecGetArray(F,&f);
132: f[0] = xdot[0] + l*(x[0] - PetscCosReal(t));
133: #if 0
134: 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: VecRestoreArrayRead(X,&x);
137: VecRestoreArrayRead(Xdot,&xdot);
138: VecRestoreArray(F,&f);
139: return(0);
140: }
142: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
143: {
144: PetscReal l = ((CECtx*)ctx)->lambda;
145: PetscErrorCode ierr;
146: PetscInt rowcol[] = {0};
147: PetscScalar J[1][1];
148: const PetscScalar *x,*xdot;
151: VecGetArrayRead(X,&x);
152: VecGetArrayRead(Xdot,&xdot);
153: J[0][0] = a + l;
154: MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
155: VecRestoreArrayRead(X,&x);
156: VecRestoreArrayRead(Xdot,&xdot);
158: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
159: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
160: if (A != B) {
161: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
162: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
163: }
164: return(0);
165: }
167: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
168: {
169: PetscReal l = ((CECtx*)ctx)->lambda;
171: PetscScalar *x;
174: VecGetArray(X,&x);
175: x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t);
176: VecRestoreArray(X,&x);
177: return(0);
178: }
180: static PetscErrorCode CECreate(Problem p)
181: {
183: CECtx *ce;
186: PetscMalloc(sizeof(CECtx),&ce);
187: p->data = (void*)ce;
189: p->destroy = &CEDestroy;
190: p->function = &CEFunction;
191: p->jacobian = &CEJacobian;
192: p->solution = &CESolution;
193: p->final_time = 10;
194: p->n = 1;
195: p->hasexact = PETSC_TRUE;
197: ce->lambda = 10;
198: PetscOptionsBegin(p->comm,NULL,"CE options","");
199: {
200: PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);
201: }
202: PetscOptionsEnd();
203: return(0);
204: }
206: /*
207: * Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
208: */
209: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
210: {
211: PetscErrorCode ierr;
212: PetscScalar *f;
213: const PetscScalar *x,*xdot;
216: VecGetArrayRead(X,&x);
217: VecGetArrayRead(Xdot,&xdot);
218: VecGetArray(F,&f);
219: f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
220: f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
221: f[2] = xdot[2] - 0.161*(x[0] - x[2]);
222: VecRestoreArrayRead(X,&x);
223: VecRestoreArrayRead(Xdot,&xdot);
224: VecRestoreArray(F,&f);
225: return(0);
226: }
228: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
229: {
230: PetscErrorCode ierr;
231: PetscInt rowcol[] = {0,1,2};
232: PetscScalar J[3][3];
233: const PetscScalar *x,*xdot;
236: VecGetArrayRead(X,&x);
237: VecGetArrayRead(Xdot,&xdot);
238: J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
239: J[0][1] = -77.27*(1. - x[0]);
240: J[0][2] = 0;
241: J[1][0] = 1./77.27*x[1];
242: J[1][1] = a + 1./77.27*(1. + x[0]);
243: J[1][2] = -1./77.27;
244: J[2][0] = -0.161;
245: J[2][1] = 0;
246: J[2][2] = a + 0.161;
247: MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
248: VecRestoreArrayRead(X,&x);
249: VecRestoreArrayRead(Xdot,&xdot);
251: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
252: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
253: if (A != B) {
254: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
255: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
256: }
257: return(0);
258: }
260: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
261: {
263: PetscScalar *x;
266: if (t != 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not implemented");
267: VecGetArray(X,&x);
268: x[0] = 1;
269: x[1] = 2;
270: x[2] = 3;
271: VecRestoreArray(X,&x);
272: return(0);
273: }
275: static PetscErrorCode OregoCreate(Problem p)
276: {
279: p->destroy = 0;
280: p->function = &OregoFunction;
281: p->jacobian = &OregoJacobian;
282: p->solution = &OregoSolution;
283: p->final_time = 360;
284: p->n = 3;
285: return(0);
286: }
289: /*
290: * User-defined monitor for comparing to exact solutions when possible
291: */
292: typedef struct {
293: MPI_Comm comm;
294: Problem problem;
295: Vec x;
296: } MonitorCtx;
298: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
299: {
301: MonitorCtx *mon = (MonitorCtx*)ctx;
302: PetscReal h,nrm_x,nrm_exact,nrm_diff;
305: if (!mon->problem->solution) return(0);
306: (*mon->problem->solution)(t,mon->x,mon->problem->data);
307: VecNorm(x,NORM_2,&nrm_x);
308: VecNorm(mon->x,NORM_2,&nrm_exact);
309: VecAYPX(mon->x,-1,x);
310: VecNorm(mon->x,NORM_2,&nrm_diff);
311: TSGetTimeStep(ts,&h);
312: if (step < 0) {
313: PetscPrintf(mon->comm,"Interpolated final solution ");
314: }
315: PetscPrintf(mon->comm,"step %4D 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);
316: return(0);
317: }
320: int main(int argc,char **argv)
321: {
322: PetscFunctionList plist = NULL;
323: char pname[256];
324: TS ts; /* nonlinear solver */
325: Vec x,r; /* solution, residual vectors */
326: Mat A; /* Jacobian matrix */
327: Problem problem;
328: PetscBool use_monitor = PETSC_FALSE;
329: PetscBool use_result = PETSC_FALSE;
330: PetscInt steps,nonlinits,linits,snesfails,rejects;
331: PetscReal ftime;
332: MonitorCtx mon;
333: PetscErrorCode ierr;
334: PetscMPIInt size;
336: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
337: Initialize program
338: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
339: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
340: MPI_Comm_size(PETSC_COMM_WORLD,&size);
341: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs");
343: /* Register the available problems */
344: PetscFunctionListAdd(&plist,"rober",&RoberCreate);
345: PetscFunctionListAdd(&plist,"ce",&CECreate);
346: PetscFunctionListAdd(&plist,"orego",&OregoCreate);
347: PetscStrcpy(pname,"ce");
349: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
350: Set runtime options
351: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
352: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
353: {
354: PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
355: use_monitor = PETSC_FALSE;
356: PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
357: PetscOptionsBool("-monitor_result","Display result","",use_result,&use_result,NULL);
358: }
359: PetscOptionsEnd();
361: /* Create the new problem */
362: PetscNew(&problem);
363: problem->comm = MPI_COMM_WORLD;
364: {
365: PetscErrorCode (*pcreate)(Problem);
367: PetscFunctionListFind(plist,pname,&pcreate);
368: if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
369: (*pcreate)(problem);
370: }
372: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
373: Create necessary matrix and vectors
374: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
375: MatCreate(PETSC_COMM_WORLD,&A);
376: MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
377: MatSetFromOptions(A);
378: MatSetUp(A);
380: MatCreateVecs(A,&x,NULL);
381: VecDuplicate(x,&r);
383: mon.comm = PETSC_COMM_WORLD;
384: mon.problem = problem;
385: VecDuplicate(x,&mon.x);
387: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: Create timestepping solver context
389: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
390: TSCreate(PETSC_COMM_WORLD,&ts);
391: TSSetProblemType(ts,TS_NONLINEAR);
392: TSSetType(ts,TSROSW); /* Rosenbrock-W */
393: TSSetIFunction(ts,NULL,problem->function,problem->data);
394: TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
395: TSSetMaxTime(ts,problem->final_time);
396: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
397: TSSetMaxStepRejections(ts,10);
398: TSSetMaxSNESFailures(ts,-1); /* unlimited */
399: if (use_monitor) {
400: TSMonitorSet(ts,&MonitorError,&mon,NULL);
401: }
403: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
404: Set initial conditions
405: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
406: (*problem->solution)(0,x,problem->data);
407: TSSetTimeStep(ts,.001);
408: TSSetSolution(ts,x);
410: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
411: Set runtime options
412: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
413: TSSetFromOptions(ts);
415: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
416: Solve nonlinear system
417: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
418: TSSolve(ts,x);
419: TSGetSolveTime(ts,&ftime);
420: TSGetStepNumber(ts,&steps);
421: TSGetSNESFailures(ts,&snesfails);
422: TSGetStepRejections(ts,&rejects);
423: TSGetSNESIterations(ts,&nonlinits);
424: TSGetKSPIterations(ts,&linits);
425: if (use_result) {
426: PetscPrintf(PETSC_COMM_WORLD,"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",steps,rejects,snesfails,(double)ftime,nonlinits,linits);
427: }
429: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
430: Free work space. All PETSc objects should be destroyed when they
431: are no longer needed.
432: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
433: MatDestroy(&A);
434: VecDestroy(&x);
435: VecDestroy(&r);
436: VecDestroy(&mon.x);
437: TSDestroy(&ts);
438: if (problem->destroy) {
439: (*problem->destroy)(problem);
440: }
441: PetscFree(problem);
442: PetscFunctionListDestroy(&plist);
444: PetscFinalize();
445: return ierr;
446: }
448: /*TEST
450: test:
451: requires: !complex
452: args: -monitor_result -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex
454: test:
455: suffix: 2
456: requires: !single !complex
457: 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
459: test:
460: suffix: 3
461: requires: !single !complex
462: 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
464: test:
465: suffix: 4
467: test:
468: suffix: 5
469: args: -snes_lag_jacobian 20 -snes_lag_jacobian_persists
471: TEST*/