Actual source code: ex8.c
petsc-3.8.4 2018-03-24
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_SELF,1,"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_SELF,1,"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;
329: PetscInt steps,nonlinits,linits,snesfails,rejects;
330: PetscReal ftime;
331: MonitorCtx mon;
332: PetscErrorCode ierr;
333: PetscMPIInt size;
335: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
336: Initialize program
337: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
338: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
339: MPI_Comm_size(PETSC_COMM_WORLD,&size);
340: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
342: /* Register the available problems */
343: PetscFunctionListAdd(&plist,"rober",&RoberCreate);
344: PetscFunctionListAdd(&plist,"ce",&CECreate);
345: PetscFunctionListAdd(&plist,"orego",&OregoCreate);
346: PetscStrcpy(pname,"ce");
348: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349: Set runtime options
350: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
351: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
352: {
353: PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
354: use_monitor = PETSC_FALSE;
355: PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
356: }
357: PetscOptionsEnd();
359: /* Create the new problem */
360: PetscNew(&problem);
361: problem->comm = MPI_COMM_WORLD;
362: {
363: PetscErrorCode (*pcreate)(Problem);
365: PetscFunctionListFind(plist,pname,&pcreate);
366: if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
367: (*pcreate)(problem);
368: }
370: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371: Create necessary matrix and vectors
372: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
373: MatCreate(PETSC_COMM_WORLD,&A);
374: MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
375: MatSetFromOptions(A);
376: MatSetUp(A);
378: MatCreateVecs(A,&x,NULL);
379: VecDuplicate(x,&r);
381: mon.comm = PETSC_COMM_WORLD;
382: mon.problem = problem;
383: VecDuplicate(x,&mon.x);
385: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
386: Create timestepping solver context
387: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
388: TSCreate(PETSC_COMM_WORLD,&ts);
389: TSSetProblemType(ts,TS_NONLINEAR);
390: TSSetType(ts,TSROSW); /* Rosenbrock-W */
391: TSSetIFunction(ts,NULL,problem->function,problem->data);
392: TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
393: TSSetMaxTime(ts,problem->final_time);
394: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
395: TSSetMaxStepRejections(ts,10);
396: TSSetMaxSNESFailures(ts,-1); /* unlimited */
397: if (use_monitor) {
398: TSMonitorSet(ts,&MonitorError,&mon,NULL);
399: }
401: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
402: Set initial conditions
403: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
404: (*problem->solution)(0,x,problem->data);
405: TSSetTimeStep(ts,.001);
406: TSSetSolution(ts,x);
408: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
409: Set runtime options
410: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
411: TSSetFromOptions(ts);
413: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414: Solve nonlinear system
415: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
416: TSSolve(ts,x);
417: TSGetSolveTime(ts,&ftime);
418: TSGetStepNumber(ts,&steps);
419: TSGetSNESFailures(ts,&snesfails);
420: TSGetStepRejections(ts,&rejects);
421: TSGetSNESIterations(ts,&nonlinits);
422: TSGetKSPIterations(ts,&linits);
423: 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);
425: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
426: Free work space. All PETSc objects should be destroyed when they
427: are no longer needed.
428: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
429: MatDestroy(&A);
430: VecDestroy(&x);
431: VecDestroy(&r);
432: VecDestroy(&mon.x);
433: TSDestroy(&ts);
434: if (problem->destroy) {
435: (*problem->destroy)(problem);
436: }
437: PetscFree(problem);
438: PetscFunctionListDestroy(&plist);
440: PetscFinalize();
441: return ierr;
442: }