Actual source code: ex8.c
petsc-3.3-p7 2013-05-11
2: /* Program usage: ./ex8 [-help] [all PETSc options] */
4: static char help[] = "Nonlinear DAE benchmark problems.\n";
7: /*
8: Include "petscts.h" so that we can use TS solvers. Note that this
9: file automatically includes:
10: petscsys.h - base PETSc routines petscvec.h - vectors
11: petscmat.h - matrices
12: petscis.h - index sets petscksp.h - Krylov subspace methods
13: petscviewer.h - viewers petscpc.h - preconditioners
14: petscksp.h - linear solvers
15: */
16: #include <petscts.h>
18: typedef struct _Problem *Problem;
19: struct _Problem {
20: PetscErrorCode (*destroy)(Problem);
21: TSIFunction function;
22: TSIJacobian jacobian;
23: PetscErrorCode (*solution)(PetscReal,Vec,void*);
25: MPI_Comm comm;
26: PetscReal final_time;
27: PetscInt n;
28: PetscBool hasexact;
29: void *data;
30: };
32: /*
33: * User-defined routines
34: */
36: /*
37: * Stiff 3-variable system from chemical reactions, due to Robertson (1966), problem ROBER in Hairer&Wanner, ODE 2, 1996
38: */
41: static PetscErrorCode RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
42: {
44: PetscScalar *x,*xdot,*f;
47: VecGetArray(X,&x);
48: VecGetArray(Xdot,&xdot);
49: VecGetArray(F,&f);
50: f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2];
51: f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]);
52: f[2] = xdot[2] - 3e7*PetscSqr(x[1]);
53: VecRestoreArray(X,&x);
54: VecRestoreArray(Xdot,&xdot);
55: VecRestoreArray(F,&f);
56: return(0);
57: }
61: static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
62: {
64: PetscInt rowcol[] = {0,1,2};
65: PetscScalar *x,*xdot,J[3][3];
68: VecGetArray(X,&x);
69: VecGetArray(Xdot,&xdot);
70: J[0][0] = a + 0.04; J[0][1] = -1e4*x[2]; J[0][2] = -1e4*x[1];
71: J[1][0] = -0.04; J[1][1] = a + 1e4*x[2] + 3e7*2*x[1]; J[1][2] = 1e4*x[1];
72: J[2][0] = 0; J[2][1] = -3e7*2*x[1]; J[2][2] = a;
73: MatSetValues(*B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
74: VecRestoreArray(X,&x);
75: VecRestoreArray(Xdot,&xdot);
77: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
79: if (*A != *B) {
80: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
82: }
83: *flag = SAME_NONZERO_PATTERN;
84: return(0);
85: }
89: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
90: {
92: PetscScalar *x;
95: if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
96: VecGetArray(X,&x);
97: x[0] = 1;
98: x[1] = 0;
99: x[2] = 0;
100: VecRestoreArray(X,&x);
101: return(0);
102: }
106: static PetscErrorCode RoberCreate(Problem p)
107: {
110: p->destroy = 0;
111: p->function = &RoberFunction;
112: p->jacobian = &RoberJacobian;
113: p->solution = &RoberSolution;
114: p->final_time = 1e11;
115: p->n = 3;
116: return(0);
117: }
120: typedef struct {
121: PetscReal lambda;
122: } CECtx;
124: /*
125: * Stiff scalar valued problem with an exact solution
126: */
129: static PetscErrorCode CEDestroy(Problem p)
130: {
134: PetscFree(p->data);
135: return(0);
136: }
140: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
141: {
143: PetscReal l = ((CECtx*)ctx)->lambda;
144: PetscScalar *x,*xdot,*f;
147: VecGetArray(X,&x);
148: VecGetArray(Xdot,&xdot);
149: VecGetArray(F,&f);
150: f[0] = xdot[0] + l*(x[0] - cos(t));
151: #if 0
152: PetscPrintf(PETSC_COMM_WORLD," f(t=%G,x=%G,xdot=%G) = %G\n",t,x[0],xdot[0],f[0]);
153: #endif
154: VecRestoreArray(X,&x);
155: VecRestoreArray(Xdot,&xdot);
156: VecRestoreArray(F,&f);
157: return(0);
158: }
162: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
163: {
164: PetscReal l = ((CECtx*)ctx)->lambda;
166: PetscInt rowcol[] = {0};
167: PetscScalar *x,*xdot,J[1][1];
170: VecGetArray(X,&x);
171: VecGetArray(Xdot,&xdot);
172: J[0][0] = a + l;
173: MatSetValues(*B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
174: VecRestoreArray(X,&x);
175: VecRestoreArray(Xdot,&xdot);
177: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
178: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
179: if (*A != *B) {
180: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
181: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
182: }
183: *flag = SAME_NONZERO_PATTERN;
184: return(0);
185: }
189: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
190: {
191: PetscReal l = ((CECtx*)ctx)->lambda;
193: PetscScalar *x;
196: VecGetArray(X,&x);
197: x[0] = l/(l*l+1)*(l*cos(t)+sin(t)) - l*l/(l*l+1)*exp(-l*t);
198: VecRestoreArray(X,&x);
199: return(0);
200: }
204: static PetscErrorCode CECreate(Problem p)
205: {
207: CECtx *ce;
210: PetscMalloc(sizeof(CECtx),&ce);
211: p->data = (void*)ce;
213: p->destroy = &CEDestroy;
214: p->function = &CEFunction;
215: p->jacobian = &CEJacobian;
216: p->solution = &CESolution;
217: p->final_time = 10;
218: p->n = 1;
219: p->hasexact = PETSC_TRUE;
221: ce->lambda = 10;
222: PetscOptionsBegin(p->comm,PETSC_NULL,"CE options","");
223: {
224: PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,PETSC_NULL);
225: }
226: PetscOptionsEnd();
227: return(0);
228: }
230: /*
231: * Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
232: */
235: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
236: {
238: PetscScalar *x,*xdot,*f;
241: VecGetArray(X,&x);
242: VecGetArray(Xdot,&xdot);
243: VecGetArray(F,&f);
244: f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
245: f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
246: f[2] = xdot[2] - 0.161*(x[0] - x[2]);
247: VecRestoreArray(X,&x);
248: VecRestoreArray(Xdot,&xdot);
249: VecRestoreArray(F,&f);
250: return(0);
251: }
255: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
256: {
258: PetscInt rowcol[] = {0,1,2};
259: PetscScalar *x,*xdot,J[3][3];
262: VecGetArray(X,&x);
263: VecGetArray(Xdot,&xdot);
264: J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
265: J[0][1] = -77.27*(1. - x[0]);
266: J[0][2] = 0;
267: J[1][0] = 1./77.27*x[1];
268: J[1][1] = a + 1./77.27*(1. + x[0]);
269: J[1][2] = -1./77.27;
270: J[2][0] = -0.161;
271: J[2][1] = 0;
272: J[2][2] = a + 0.161;
273: MatSetValues(*B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
274: VecRestoreArray(X,&x);
275: VecRestoreArray(Xdot,&xdot);
277: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
278: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
279: if (*A != *B) {
280: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
281: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
282: }
283: *flag = SAME_NONZERO_PATTERN;
284: return(0);
285: }
289: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
290: {
292: PetscScalar *x;
295: if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
296: VecGetArray(X,&x);
297: x[0] = 1;
298: x[1] = 2;
299: x[2] = 3;
300: VecRestoreArray(X,&x);
301: return(0);
302: }
306: static PetscErrorCode OregoCreate(Problem p)
307: {
310: p->destroy = 0;
311: p->function = &OregoFunction;
312: p->jacobian = &OregoJacobian;
313: p->solution = &OregoSolution;
314: p->final_time = 360;
315: p->n = 3;
316: return(0);
317: }
320: /*
321: * User-defined monitor for comparing to exact solutions when possible
322: */
323: typedef struct {
324: MPI_Comm comm;
325: Problem problem;
326: Vec x;
327: } MonitorCtx;
331: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
332: {
334: MonitorCtx *mon = (MonitorCtx*)ctx;
335: PetscReal h,nrm_x,nrm_exact,nrm_diff;
338: if (!mon->problem->solution) return(0);
339: (*mon->problem->solution)(t,mon->x,mon->problem->data);
340: VecNorm(x,NORM_2,&nrm_x);
341: VecNorm(mon->x,NORM_2,&nrm_exact);
342: VecAYPX(mon->x,-1,x);
343: VecNorm(mon->x,NORM_2,&nrm_diff);
344: TSGetTimeStep(ts,&h);
345: 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);
346: return(0);
347: }
353: int main(int argc,char **argv)
354: {
355: PetscFList plist = PETSC_NULL;
356: char pname[256];
357: TS ts; /* nonlinear solver */
358: Vec x,r; /* solution, residual vectors */
359: Mat A; /* Jacobian matrix */
360: Problem problem;
361: PetscBool use_monitor;
362: PetscInt steps,maxsteps = 1000,nonlinits,linits,snesfails,rejects;
363: PetscReal ftime;
364: MonitorCtx mon;
365: PetscErrorCode ierr;
367: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
368: Initialize program
369: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
370: PetscInitialize(&argc,&argv,(char *)0,help);
372: /* Register the available problems */
373: PetscFListAdd(&plist,"rober","",(void(*)(void))&RoberCreate);
374: PetscFListAdd(&plist,"ce", "",(void(*)(void))&CECreate );
375: PetscFListAdd(&plist,"orego","",(void(*)(void))&OregoCreate);
376: PetscStrcpy(pname,"ce");
378: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
379: Set runtime options
380: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
381: PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Timestepping benchmark options","");
382: {
383: PetscOptionsList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),PETSC_NULL);
384: use_monitor = PETSC_FALSE;
385: PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,PETSC_NULL);
386: }
387: PetscOptionsEnd();
389: /* Create the new problem */
390: PetscNew(struct _Problem,&problem);
391: problem->comm = MPI_COMM_WORLD;
392: {
393: PetscErrorCode (*pcreate)(Problem);
395: PetscFListFind(plist,MPI_COMM_WORLD,pname,PETSC_FALSE,(void(**)(void))&pcreate);
396: if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
397: (*pcreate)(problem);
398: }
401: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
402: Create necessary matrix and vectors, solve same ODE on every process
403: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
404: MatCreate(PETSC_COMM_WORLD,&A);
405: MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
406: MatSetFromOptions(A);
407: MatSetUp(A);
409: MatGetVecs(A,&x,PETSC_NULL);
410: VecDuplicate(x,&r);
412: mon.comm = PETSC_COMM_WORLD;
413: mon.problem = problem;
414: VecDuplicate(x,&mon.x);
416: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
417: Create timestepping solver context
418: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
419: TSCreate(PETSC_COMM_WORLD,&ts);
420: TSSetProblemType(ts,TS_NONLINEAR);
421: TSSetType(ts,TSROSW); /* Rosenbrock-W */
422: TSSetIFunction(ts,PETSC_NULL,problem->function,problem->data);
423: TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
424: TSSetDuration(ts,maxsteps,problem->final_time);
425: TSSetMaxStepRejections(ts,10);
426: TSSetMaxSNESFailures(ts,-1); /* unlimited */
427: if (use_monitor) {
428: TSMonitorSet(ts,&MonitorError,&mon,PETSC_NULL);
429: }
431: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
432: Set initial conditions
433: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
434: (*problem->solution)(0,x,problem->data);
435: TSSetInitialTimeStep(ts,0.0,.001);
436: TSSetSolution(ts,x);
438: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
439: Set runtime options
440: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
441: TSSetFromOptions(ts);
443: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
444: Solve nonlinear system
445: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
446: TSSolve(ts,x,&ftime);
447: TSGetTimeStepNumber(ts,&steps);
448: TSGetSNESFailures(ts,&snesfails);
449: TSGetStepRejections(ts,&rejects);
450: TSGetSNESIterations(ts,&nonlinits);
451: TSGetKSPIterations(ts,&linits);
452: PetscPrintf(PETSC_COMM_WORLD,"steps %D (%D rejected, %D SNES fails), ftime %G, nonlinits %D, linits %D\n",steps,rejects,snesfails,ftime,nonlinits,linits);
453: if (problem->hasexact) {MonitorError(ts,steps,ftime,x,&mon);}
455: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
456: Free work space. All PETSc objects should be destroyed when they
457: are no longer needed.
458: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
459: MatDestroy(&A);
460: VecDestroy(&x);
461: VecDestroy(&r);
462: VecDestroy(&mon.x);
463: TSDestroy(&ts);
464: if (problem->destroy) {
465: (*problem->destroy)(problem);
466: }
467: PetscFree(problem);
468: PetscFListDestroy(&plist);
470: PetscFinalize();
471: return(0);
472: }