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;
37: VecGetArrayRead(X,&x);
38: VecGetArrayRead(Xdot,&xdot);
39: 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: VecRestoreArrayRead(X,&x);
44: VecRestoreArrayRead(Xdot,&xdot);
45: VecRestoreArray(F,&f);
46: return 0;
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;
56: VecGetArrayRead(X,&x);
57: VecGetArrayRead(Xdot,&xdot);
58: J[0][0] = a + 0.04; J[0][1] = -1e4*x[2]; J[0][2] = -1e4*x[1];
59: J[1][0] = -0.04; J[1][1] = a + 1e4*x[2] + 3e7*2*x[1]; J[1][2] = 1e4*x[1];
60: J[2][0] = 0; J[2][1] = -3e7*2*x[1]; J[2][2] = a;
61: MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
62: VecRestoreArrayRead(X,&x);
63: VecRestoreArrayRead(Xdot,&xdot);
65: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
67: if (A != B) {
68: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
70: }
71: return 0;
72: }
74: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
75: {
76: PetscScalar *x;
80: VecGetArray(X,&x);
81: x[0] = 1;
82: x[1] = 0;
83: x[2] = 0;
84: VecRestoreArray(X,&x);
85: return 0;
86: }
88: static PetscErrorCode RoberCreate(Problem p)
89: {
91: p->destroy = 0;
92: p->function = &RoberFunction;
93: p->jacobian = &RoberJacobian;
94: p->solution = &RoberSolution;
95: p->final_time = 1e11;
96: p->n = 3;
97: return 0;
98: }
100: /*
101: Stiff scalar valued problem
102: */
104: typedef struct {
105: PetscReal lambda;
106: } CECtx;
108: static PetscErrorCode CEDestroy(Problem p)
109: {
111: PetscFree(p->data);
112: return 0;
113: }
115: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
116: {
117: PetscReal l = ((CECtx*)ctx)->lambda;
118: PetscScalar *f;
119: const PetscScalar *x,*xdot;
122: VecGetArrayRead(X,&x);
123: VecGetArrayRead(Xdot,&xdot);
124: VecGetArray(F,&f);
125: f[0] = xdot[0] + l*(x[0] - PetscCosReal(t));
126: #if 0
127: PetscPrintf(PETSC_COMM_WORLD," f(t=%g,x=%g,xdot=%g) = %g\n",(double)t,(double)x[0],(double)xdot[0],(double)f[0]);
128: #endif
129: VecRestoreArrayRead(X,&x);
130: VecRestoreArrayRead(Xdot,&xdot);
131: VecRestoreArray(F,&f);
132: return 0;
133: }
135: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
136: {
137: PetscReal l = ((CECtx*)ctx)->lambda;
138: PetscInt rowcol[] = {0};
139: PetscScalar J[1][1];
140: const PetscScalar *x,*xdot;
143: VecGetArrayRead(X,&x);
144: VecGetArrayRead(Xdot,&xdot);
145: J[0][0] = a + l;
146: MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
147: VecRestoreArrayRead(X,&x);
148: VecRestoreArrayRead(Xdot,&xdot);
150: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
151: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
152: if (A != B) {
153: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
154: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
155: }
156: return 0;
157: }
159: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
160: {
161: PetscReal l = ((CECtx*)ctx)->lambda;
162: PetscScalar *x;
165: VecGetArray(X,&x);
166: x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t);
167: VecRestoreArray(X,&x);
168: return 0;
169: }
171: static PetscErrorCode CECreate(Problem p)
172: {
174: CECtx *ce;
177: PetscMalloc(sizeof(CECtx),&ce);
178: p->data = (void*)ce;
180: p->destroy = &CEDestroy;
181: p->function = &CEFunction;
182: p->jacobian = &CEJacobian;
183: p->solution = &CESolution;
184: p->final_time = 10;
185: p->n = 1;
186: p->hasexact = PETSC_TRUE;
188: ce->lambda = 10;
189: PetscOptionsBegin(p->comm,NULL,"CE options","");
190: {
191: PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);
192: }
193: PetscOptionsEnd();
194: return 0;
195: }
197: /*
198: Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
199: */
200: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
201: {
202: PetscScalar *f;
203: const PetscScalar *x,*xdot;
206: VecGetArrayRead(X,&x);
207: VecGetArrayRead(Xdot,&xdot);
208: VecGetArray(F,&f);
209: f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
210: f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
211: f[2] = xdot[2] - 0.161*(x[0] - x[2]);
212: VecRestoreArrayRead(X,&x);
213: VecRestoreArrayRead(Xdot,&xdot);
214: VecRestoreArray(F,&f);
215: return 0;
216: }
218: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
219: {
220: PetscInt rowcol[] = {0,1,2};
221: PetscScalar J[3][3];
222: const PetscScalar *x,*xdot;
225: VecGetArrayRead(X,&x);
226: VecGetArrayRead(Xdot,&xdot);
227: J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
228: J[0][1] = -77.27*(1. - x[0]);
229: J[0][2] = 0;
230: J[1][0] = 1./77.27*x[1];
231: J[1][1] = a + 1./77.27*(1. + x[0]);
232: J[1][2] = -1./77.27;
233: J[2][0] = -0.161;
234: J[2][1] = 0;
235: J[2][2] = a + 0.161;
236: MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
237: VecRestoreArrayRead(X,&x);
238: VecRestoreArrayRead(Xdot,&xdot);
240: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
241: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
242: if (A != B) {
243: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
244: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
245: }
246: return 0;
247: }
249: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
250: {
251: PetscScalar *x;
255: VecGetArray(X,&x);
256: x[0] = 1;
257: x[1] = 2;
258: x[2] = 3;
259: VecRestoreArray(X,&x);
260: return 0;
261: }
263: static PetscErrorCode OregoCreate(Problem p)
264: {
266: p->destroy = 0;
267: p->function = &OregoFunction;
268: p->jacobian = &OregoJacobian;
269: p->solution = &OregoSolution;
270: p->final_time = 360;
271: p->n = 3;
272: return 0;
273: }
275: /*
276: User-defined monitor for comparing to exact solutions when possible
277: */
278: typedef struct {
279: MPI_Comm comm;
280: Problem problem;
281: Vec x;
282: } MonitorCtx;
284: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
285: {
286: MonitorCtx *mon = (MonitorCtx*)ctx;
287: PetscReal h,nrm_x,nrm_exact,nrm_diff;
290: if (!mon->problem->solution) return 0;
291: (*mon->problem->solution)(t,mon->x,mon->problem->data);
292: VecNorm(x,NORM_2,&nrm_x);
293: VecNorm(mon->x,NORM_2,&nrm_exact);
294: VecAYPX(mon->x,-1,x);
295: VecNorm(mon->x,NORM_2,&nrm_diff);
296: TSGetTimeStep(ts,&h);
297: if (step < 0) {
298: PetscPrintf(mon->comm,"Interpolated final solution ");
299: }
300: 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);
301: return 0;
302: }
304: int main(int argc,char **argv)
305: {
306: PetscFunctionList plist = NULL;
307: char pname[256];
308: TS ts; /* nonlinear solver */
309: Vec x,r; /* solution, residual vectors */
310: Mat A; /* Jacobian matrix */
311: Problem problem;
312: PetscBool use_monitor = PETSC_FALSE;
313: PetscBool use_result = PETSC_FALSE;
314: PetscInt steps,nonlinits,linits,snesfails,rejects;
315: PetscReal ftime;
316: MonitorCtx mon;
317: PetscErrorCode ierr;
318: PetscMPIInt size;
320: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
321: Initialize program
322: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
323: PetscInitialize(&argc,&argv,(char*)0,help);
324: MPI_Comm_size(PETSC_COMM_WORLD,&size);
327: /* Register the available problems */
328: PetscFunctionListAdd(&plist,"rober",&RoberCreate);
329: PetscFunctionListAdd(&plist,"ce",&CECreate);
330: PetscFunctionListAdd(&plist,"orego",&OregoCreate);
331: PetscStrcpy(pname,"ce");
333: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
334: Set runtime options
335: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
336: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
337: {
338: PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
339: use_monitor = PETSC_FALSE;
340: PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
341: PetscOptionsBool("-monitor_result","Display result","",use_result,&use_result,NULL);
342: }
343: PetscOptionsEnd();
345: /* Create the new problem */
346: PetscNew(&problem);
347: problem->comm = MPI_COMM_WORLD;
348: {
349: PetscErrorCode (*pcreate)(Problem);
351: PetscFunctionListFind(plist,pname,&pcreate);
353: (*pcreate)(problem);
354: }
356: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
357: Create necessary matrix and vectors
358: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
359: MatCreate(PETSC_COMM_WORLD,&A);
360: MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
361: MatSetFromOptions(A);
362: MatSetUp(A);
364: MatCreateVecs(A,&x,NULL);
365: VecDuplicate(x,&r);
367: mon.comm = PETSC_COMM_WORLD;
368: mon.problem = problem;
369: VecDuplicate(x,&mon.x);
371: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
372: Create timestepping solver context
373: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
374: TSCreate(PETSC_COMM_WORLD,&ts);
375: TSSetProblemType(ts,TS_NONLINEAR);
376: TSSetType(ts,TSROSW); /* Rosenbrock-W */
377: TSSetIFunction(ts,NULL,problem->function,problem->data);
378: TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
379: TSSetMaxTime(ts,problem->final_time);
380: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
381: TSSetMaxStepRejections(ts,10);
382: TSSetMaxSNESFailures(ts,-1); /* unlimited */
383: if (use_monitor) {
384: TSMonitorSet(ts,&MonitorError,&mon,NULL);
385: }
387: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: Set initial conditions
389: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
390: (*problem->solution)(0,x,problem->data);
391: TSSetTimeStep(ts,.001);
392: TSSetSolution(ts,x);
394: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
395: Set runtime options
396: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
397: TSSetFromOptions(ts);
399: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
400: Solve nonlinear system
401: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
402: TSSolve(ts,x);
403: TSGetSolveTime(ts,&ftime);
404: TSGetStepNumber(ts,&steps);
405: TSGetSNESFailures(ts,&snesfails);
406: TSGetStepRejections(ts,&rejects);
407: TSGetSNESIterations(ts,&nonlinits);
408: TSGetKSPIterations(ts,&linits);
409: if (use_result) {
410: 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);
411: }
413: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414: Free work space. All PETSc objects should be destroyed when they
415: are no longer needed.
416: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
417: MatDestroy(&A);
418: VecDestroy(&x);
419: VecDestroy(&r);
420: VecDestroy(&mon.x);
421: TSDestroy(&ts);
422: if (problem->destroy) {
423: (*problem->destroy)(problem);
424: }
425: PetscFree(problem);
426: PetscFunctionListDestroy(&plist);
428: 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
451: test:
452: suffix: 5
453: args: -snes_lag_jacobian 20 -snes_lag_jacobian_persists
455: TEST*/