Actual source code: ex8.c
petsc-3.4.5 2014-06-29
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: */
33: static PetscErrorCode RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
34: {
36: PetscScalar *x,*xdot,*f;
39: VecGetArray(X,&x);
40: VecGetArray(Xdot,&xdot);
41: VecGetArray(F,&f);
42: f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2];
43: f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]);
44: f[2] = xdot[2] - 3e7*PetscSqr(x[1]);
45: VecRestoreArray(X,&x);
46: VecRestoreArray(Xdot,&xdot);
47: VecRestoreArray(F,&f);
48: return(0);
49: }
53: static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
54: {
56: PetscInt rowcol[] = {0,1,2};
57: PetscScalar *x,*xdot,J[3][3];
60: VecGetArray(X,&x);
61: VecGetArray(Xdot,&xdot);
62: J[0][0] = a + 0.04; J[0][1] = -1e4*x[2]; J[0][2] = -1e4*x[1];
63: J[1][0] = -0.04; J[1][1] = a + 1e4*x[2] + 3e7*2*x[1]; J[1][2] = 1e4*x[1];
64: J[2][0] = 0; J[2][1] = -3e7*2*x[1]; J[2][2] = a;
65: MatSetValues(*B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
66: VecRestoreArray(X,&x);
67: VecRestoreArray(Xdot,&xdot);
69: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
71: if (*A != *B) {
72: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
74: }
75: *flag = SAME_NONZERO_PATTERN;
76: return(0);
77: }
81: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
82: {
84: PetscScalar *x;
87: if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
88: VecGetArray(X,&x);
89: x[0] = 1;
90: x[1] = 0;
91: x[2] = 0;
92: VecRestoreArray(X,&x);
93: return(0);
94: }
98: static PetscErrorCode RoberCreate(Problem p)
99: {
102: p->destroy = 0;
103: p->function = &RoberFunction;
104: p->jacobian = &RoberJacobian;
105: p->solution = &RoberSolution;
106: p->final_time = 1e11;
107: p->n = 3;
108: return(0);
109: }
111: /*
112: Stiff scalar valued problem
113: */
115: typedef struct {
116: PetscReal lambda;
117: } CECtx;
121: static PetscErrorCode CEDestroy(Problem p)
122: {
126: PetscFree(p->data);
127: return(0);
128: }
132: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
133: {
135: PetscReal l = ((CECtx*)ctx)->lambda;
136: PetscScalar *x,*xdot,*f;
139: VecGetArray(X,&x);
140: VecGetArray(Xdot,&xdot);
141: VecGetArray(F,&f);
142: f[0] = xdot[0] + l*(x[0] - cos(t));
143: #if 0
144: PetscPrintf(PETSC_COMM_WORLD," f(t=%G,x=%G,xdot=%G) = %G\n",t,x[0],xdot[0],f[0]);
145: #endif
146: VecRestoreArray(X,&x);
147: VecRestoreArray(Xdot,&xdot);
148: VecRestoreArray(F,&f);
149: return(0);
150: }
154: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
155: {
156: PetscReal l = ((CECtx*)ctx)->lambda;
158: PetscInt rowcol[] = {0};
159: PetscScalar *x,*xdot,J[1][1];
162: VecGetArray(X,&x);
163: VecGetArray(Xdot,&xdot);
164: J[0][0] = a + l;
165: MatSetValues(*B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
166: VecRestoreArray(X,&x);
167: VecRestoreArray(Xdot,&xdot);
169: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
170: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
171: if (*A != *B) {
172: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
173: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
174: }
175: *flag = SAME_NONZERO_PATTERN;
176: return(0);
177: }
181: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
182: {
183: PetscReal l = ((CECtx*)ctx)->lambda;
185: PetscScalar *x;
188: VecGetArray(X,&x);
189: x[0] = l/(l*l+1)*(l*cos(t)+sin(t)) - l*l/(l*l+1)*exp(-l*t);
190: VecRestoreArray(X,&x);
191: return(0);
192: }
196: static PetscErrorCode CECreate(Problem p)
197: {
199: CECtx *ce;
202: PetscMalloc(sizeof(CECtx),&ce);
203: p->data = (void*)ce;
205: p->destroy = &CEDestroy;
206: p->function = &CEFunction;
207: p->jacobian = &CEJacobian;
208: p->solution = &CESolution;
209: p->final_time = 10;
210: p->n = 1;
211: p->hasexact = PETSC_TRUE;
213: ce->lambda = 10;
214: PetscOptionsBegin(p->comm,NULL,"CE options","");
215: {
216: PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);
217: }
218: PetscOptionsEnd();
219: return(0);
220: }
222: /*
223: * Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
224: */
227: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
228: {
230: PetscScalar *x,*xdot,*f;
233: VecGetArray(X,&x);
234: VecGetArray(Xdot,&xdot);
235: VecGetArray(F,&f);
236: f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
237: f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
238: f[2] = xdot[2] - 0.161*(x[0] - x[2]);
239: VecRestoreArray(X,&x);
240: VecRestoreArray(Xdot,&xdot);
241: VecRestoreArray(F,&f);
242: return(0);
243: }
247: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
248: {
250: PetscInt rowcol[] = {0,1,2};
251: PetscScalar *x,*xdot,J[3][3];
254: VecGetArray(X,&x);
255: VecGetArray(Xdot,&xdot);
256: J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
257: J[0][1] = -77.27*(1. - x[0]);
258: J[0][2] = 0;
259: J[1][0] = 1./77.27*x[1];
260: J[1][1] = a + 1./77.27*(1. + x[0]);
261: J[1][2] = -1./77.27;
262: J[2][0] = -0.161;
263: J[2][1] = 0;
264: J[2][2] = a + 0.161;
265: MatSetValues(*B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
266: VecRestoreArray(X,&x);
267: VecRestoreArray(Xdot,&xdot);
269: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
270: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
271: if (*A != *B) {
272: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
273: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
274: }
275: *flag = SAME_NONZERO_PATTERN;
276: return(0);
277: }
281: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
282: {
284: PetscScalar *x;
287: if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
288: VecGetArray(X,&x);
289: x[0] = 1;
290: x[1] = 2;
291: x[2] = 3;
292: VecRestoreArray(X,&x);
293: return(0);
294: }
298: static PetscErrorCode OregoCreate(Problem p)
299: {
302: p->destroy = 0;
303: p->function = &OregoFunction;
304: p->jacobian = &OregoJacobian;
305: p->solution = &OregoSolution;
306: p->final_time = 360;
307: p->n = 3;
308: return(0);
309: }
312: /*
313: * User-defined monitor for comparing to exact solutions when possible
314: */
315: typedef struct {
316: MPI_Comm comm;
317: Problem problem;
318: Vec x;
319: } MonitorCtx;
323: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
324: {
326: MonitorCtx *mon = (MonitorCtx*)ctx;
327: PetscReal h,nrm_x,nrm_exact,nrm_diff;
330: if (!mon->problem->solution) return(0);
331: (*mon->problem->solution)(t,mon->x,mon->problem->data);
332: VecNorm(x,NORM_2,&nrm_x);
333: VecNorm(mon->x,NORM_2,&nrm_exact);
334: VecAYPX(mon->x,-1,x);
335: VecNorm(mon->x,NORM_2,&nrm_diff);
336: TSGetTimeStep(ts,&h);
337: 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,t,h,nrm_x,nrm_exact,nrm_diff);
338: return(0);
339: }
344: int main(int argc,char **argv)
345: {
346: PetscFunctionList plist = NULL;
347: char pname[256];
348: TS ts; /* nonlinear solver */
349: Vec x,r; /* solution, residual vectors */
350: Mat A; /* Jacobian matrix */
351: Problem problem;
352: PetscBool use_monitor;
353: PetscInt steps,maxsteps = 1000,nonlinits,linits,snesfails,rejects;
354: PetscReal ftime;
355: MonitorCtx mon;
356: PetscErrorCode ierr;
357: PetscMPIInt size;
359: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
360: Initialize program
361: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
362: PetscInitialize(&argc,&argv,(char*)0,help);
363: MPI_Comm_size(PETSC_COMM_WORLD,&size);
364: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
366: /* Register the available problems */
367: PetscFunctionListAdd(&plist,"rober",&RoberCreate);
368: PetscFunctionListAdd(&plist,"ce",&CECreate);
369: PetscFunctionListAdd(&plist,"orego",&OregoCreate);
370: PetscStrcpy(pname,"ce");
372: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
373: Set runtime options
374: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
375: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
376: {
377: PetscOptionsList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
378: use_monitor = PETSC_FALSE;
379: PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
380: }
381: PetscOptionsEnd();
383: /* Create the new problem */
384: PetscNew(struct _Problem,&problem);
385: problem->comm = MPI_COMM_WORLD;
386: {
387: PetscErrorCode (*pcreate)(Problem);
389: PetscFunctionListFind(plist,pname,&pcreate);
390: if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
391: (*pcreate)(problem);
392: }
394: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
395: Create necessary matrix and vectors
396: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
397: MatCreate(PETSC_COMM_WORLD,&A);
398: MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
399: MatSetFromOptions(A);
400: MatSetUp(A);
402: MatGetVecs(A,&x,NULL);
403: VecDuplicate(x,&r);
405: mon.comm = PETSC_COMM_WORLD;
406: mon.problem = problem;
407: VecDuplicate(x,&mon.x);
409: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
410: Create timestepping solver context
411: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
412: TSCreate(PETSC_COMM_WORLD,&ts);
413: TSSetProblemType(ts,TS_NONLINEAR);
414: TSSetType(ts,TSROSW); /* Rosenbrock-W */
415: TSSetIFunction(ts,NULL,problem->function,problem->data);
416: TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
417: TSSetDuration(ts,maxsteps,problem->final_time);
418: TSSetMaxStepRejections(ts,10);
419: TSSetMaxSNESFailures(ts,-1); /* unlimited */
420: if (use_monitor) {
421: TSMonitorSet(ts,&MonitorError,&mon,NULL);
422: }
424: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
425: Set initial conditions
426: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
427: (*problem->solution)(0,x,problem->data);
428: TSSetInitialTimeStep(ts,0.0,.001);
429: TSSetSolution(ts,x);
431: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
432: Set runtime options
433: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
434: TSSetFromOptions(ts);
436: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
437: Solve nonlinear system
438: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
439: TSSolve(ts,x);
440: TSGetSolveTime(ts,&ftime);
441: TSGetTimeStepNumber(ts,&steps);
442: TSGetSNESFailures(ts,&snesfails);
443: TSGetStepRejections(ts,&rejects);
444: TSGetSNESIterations(ts,&nonlinits);
445: TSGetKSPIterations(ts,&linits);
446: PetscPrintf(PETSC_COMM_WORLD,"steps %D (%D rejected, %D SNES fails), ftime %G, nonlinits %D, linits %D\n",steps,rejects,snesfails,ftime,nonlinits,linits);
447: if (problem->hasexact) {
448: MonitorError(ts,steps,ftime,x,&mon);
449: }
451: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
452: Free work space. All PETSc objects should be destroyed when they
453: are no longer needed.
454: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
455: MatDestroy(&A);
456: VecDestroy(&x);
457: VecDestroy(&r);
458: VecDestroy(&mon.x);
459: TSDestroy(&ts);
460: if (problem->destroy) {
461: (*problem->destroy)(problem);
462: }
463: PetscFree(problem);
464: PetscFunctionListDestroy(&plist);
466: PetscFinalize();
467: return(0);
468: }