Actual source code: ex8.c
petsc-3.6.1 2015-08-06
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: {
35: PetscErrorCode ierr;
36: PetscScalar *f;
37: const PetscScalar *x,*xdot;
40: VecGetArrayRead(X,&x);
41: VecGetArrayRead(Xdot,&xdot);
42: VecGetArray(F,&f);
43: f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2];
44: f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]);
45: f[2] = xdot[2] - 3e7*PetscSqr(x[1]);
46: VecRestoreArrayRead(X,&x);
47: VecRestoreArrayRead(Xdot,&xdot);
48: VecRestoreArray(F,&f);
49: return(0);
50: }
54: static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
55: {
56: PetscErrorCode ierr;
57: PetscInt rowcol[] = {0,1,2};
58: PetscScalar J[3][3];
59: const PetscScalar *x,*xdot;
62: VecGetArrayRead(X,&x);
63: VecGetArrayRead(Xdot,&xdot);
64: J[0][0] = a + 0.04; J[0][1] = -1e4*x[2]; J[0][2] = -1e4*x[1];
65: J[1][0] = -0.04; J[1][1] = a + 1e4*x[2] + 3e7*2*x[1]; J[1][2] = 1e4*x[1];
66: J[2][0] = 0; J[2][1] = -3e7*2*x[1]; J[2][2] = a;
67: MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
68: VecRestoreArrayRead(X,&x);
69: VecRestoreArrayRead(Xdot,&xdot);
71: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
72: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
73: if (A != B) {
74: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
76: }
77: return(0);
78: }
82: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
83: {
85: PetscScalar *x;
88: if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
89: VecGetArray(X,&x);
90: x[0] = 1;
91: x[1] = 0;
92: x[2] = 0;
93: VecRestoreArray(X,&x);
94: return(0);
95: }
99: static PetscErrorCode RoberCreate(Problem p)
100: {
103: p->destroy = 0;
104: p->function = &RoberFunction;
105: p->jacobian = &RoberJacobian;
106: p->solution = &RoberSolution;
107: p->final_time = 1e11;
108: p->n = 3;
109: return(0);
110: }
112: /*
113: Stiff scalar valued problem
114: */
116: typedef struct {
117: PetscReal lambda;
118: } CECtx;
122: static PetscErrorCode CEDestroy(Problem p)
123: {
127: PetscFree(p->data);
128: return(0);
129: }
133: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
134: {
135: PetscErrorCode ierr;
136: PetscReal l = ((CECtx*)ctx)->lambda;
137: PetscScalar *f;
138: const PetscScalar *x,*xdot;
141: VecGetArrayRead(X,&x);
142: VecGetArrayRead(Xdot,&xdot);
143: VecGetArray(F,&f);
144: f[0] = xdot[0] + l*(x[0] - PetscCosReal(t));
145: #if 0
146: PetscPrintf(PETSC_COMM_WORLD," f(t=%g,x=%g,xdot=%g) = %g\n",(double)t,(double)x[0],(double)xdot[0],(double)f[0]);
147: #endif
148: VecRestoreArrayRead(X,&x);
149: VecRestoreArrayRead(Xdot,&xdot);
150: VecRestoreArray(F,&f);
151: return(0);
152: }
156: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
157: {
158: PetscReal l = ((CECtx*)ctx)->lambda;
159: PetscErrorCode ierr;
160: PetscInt rowcol[] = {0};
161: PetscScalar J[1][1];
162: const PetscScalar *x,*xdot;
165: VecGetArrayRead(X,&x);
166: VecGetArrayRead(Xdot,&xdot);
167: J[0][0] = a + l;
168: MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
169: VecRestoreArrayRead(X,&x);
170: VecRestoreArrayRead(Xdot,&xdot);
172: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
173: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
174: if (A != B) {
175: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
176: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
177: }
178: return(0);
179: }
183: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
184: {
185: PetscReal l = ((CECtx*)ctx)->lambda;
187: PetscScalar *x;
190: VecGetArray(X,&x);
191: x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t);
192: VecRestoreArray(X,&x);
193: return(0);
194: }
198: static PetscErrorCode CECreate(Problem p)
199: {
201: CECtx *ce;
204: PetscMalloc(sizeof(CECtx),&ce);
205: p->data = (void*)ce;
207: p->destroy = &CEDestroy;
208: p->function = &CEFunction;
209: p->jacobian = &CEJacobian;
210: p->solution = &CESolution;
211: p->final_time = 10;
212: p->n = 1;
213: p->hasexact = PETSC_TRUE;
215: ce->lambda = 10;
216: PetscOptionsBegin(p->comm,NULL,"CE options","");
217: {
218: PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);
219: }
220: PetscOptionsEnd();
221: return(0);
222: }
224: /*
225: * Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
226: */
229: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
230: {
231: PetscErrorCode ierr;
232: PetscScalar *f;
233: const PetscScalar *x,*xdot;
236: VecGetArrayRead(X,&x);
237: VecGetArrayRead(Xdot,&xdot);
238: VecGetArray(F,&f);
239: f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
240: f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
241: f[2] = xdot[2] - 0.161*(x[0] - x[2]);
242: VecRestoreArrayRead(X,&x);
243: VecRestoreArrayRead(Xdot,&xdot);
244: VecRestoreArray(F,&f);
245: return(0);
246: }
250: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
251: {
252: PetscErrorCode ierr;
253: PetscInt rowcol[] = {0,1,2};
254: PetscScalar J[3][3];
255: const PetscScalar *x,*xdot;
258: VecGetArrayRead(X,&x);
259: VecGetArrayRead(Xdot,&xdot);
260: J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
261: J[0][1] = -77.27*(1. - x[0]);
262: J[0][2] = 0;
263: J[1][0] = 1./77.27*x[1];
264: J[1][1] = a + 1./77.27*(1. + x[0]);
265: J[1][2] = -1./77.27;
266: J[2][0] = -0.161;
267: J[2][1] = 0;
268: J[2][2] = a + 0.161;
269: MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
270: VecRestoreArrayRead(X,&x);
271: VecRestoreArrayRead(Xdot,&xdot);
273: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
274: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
275: if (A != B) {
276: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
277: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
278: }
279: return(0);
280: }
284: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
285: {
287: PetscScalar *x;
290: if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
291: VecGetArray(X,&x);
292: x[0] = 1;
293: x[1] = 2;
294: x[2] = 3;
295: VecRestoreArray(X,&x);
296: return(0);
297: }
301: static PetscErrorCode OregoCreate(Problem p)
302: {
305: p->destroy = 0;
306: p->function = &OregoFunction;
307: p->jacobian = &OregoJacobian;
308: p->solution = &OregoSolution;
309: p->final_time = 360;
310: p->n = 3;
311: return(0);
312: }
315: /*
316: * User-defined monitor for comparing to exact solutions when possible
317: */
318: typedef struct {
319: MPI_Comm comm;
320: Problem problem;
321: Vec x;
322: } MonitorCtx;
326: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
327: {
329: MonitorCtx *mon = (MonitorCtx*)ctx;
330: PetscReal h,nrm_x,nrm_exact,nrm_diff;
333: if (!mon->problem->solution) return(0);
334: (*mon->problem->solution)(t,mon->x,mon->problem->data);
335: VecNorm(x,NORM_2,&nrm_x);
336: VecNorm(mon->x,NORM_2,&nrm_exact);
337: VecAYPX(mon->x,-1,x);
338: VecNorm(mon->x,NORM_2,&nrm_diff);
339: TSGetTimeStep(ts,&h);
340: 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);
341: return(0);
342: }
347: int main(int argc,char **argv)
348: {
349: PetscFunctionList plist = NULL;
350: char pname[256];
351: TS ts; /* nonlinear solver */
352: Vec x,r; /* solution, residual vectors */
353: Mat A; /* Jacobian matrix */
354: Problem problem;
355: PetscBool use_monitor;
356: PetscInt steps,maxsteps = 1000,nonlinits,linits,snesfails,rejects;
357: PetscReal ftime;
358: MonitorCtx mon;
359: PetscErrorCode ierr;
360: PetscMPIInt size;
362: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
363: Initialize program
364: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
365: PetscInitialize(&argc,&argv,(char*)0,help);
366: MPI_Comm_size(PETSC_COMM_WORLD,&size);
367: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
369: /* Register the available problems */
370: PetscFunctionListAdd(&plist,"rober",&RoberCreate);
371: PetscFunctionListAdd(&plist,"ce",&CECreate);
372: PetscFunctionListAdd(&plist,"orego",&OregoCreate);
373: PetscStrcpy(pname,"ce");
375: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
376: Set runtime options
377: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
378: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
379: {
380: PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
381: use_monitor = PETSC_FALSE;
382: PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
383: }
384: PetscOptionsEnd();
386: /* Create the new problem */
387: PetscNew(&problem);
388: problem->comm = MPI_COMM_WORLD;
389: {
390: PetscErrorCode (*pcreate)(Problem);
392: PetscFunctionListFind(plist,pname,&pcreate);
393: if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
394: (*pcreate)(problem);
395: }
397: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
398: Create necessary matrix and vectors
399: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
400: MatCreate(PETSC_COMM_WORLD,&A);
401: MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
402: MatSetFromOptions(A);
403: MatSetUp(A);
405: MatCreateVecs(A,&x,NULL);
406: VecDuplicate(x,&r);
408: mon.comm = PETSC_COMM_WORLD;
409: mon.problem = problem;
410: VecDuplicate(x,&mon.x);
412: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
413: Create timestepping solver context
414: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
415: TSCreate(PETSC_COMM_WORLD,&ts);
416: TSSetProblemType(ts,TS_NONLINEAR);
417: TSSetType(ts,TSROSW); /* Rosenbrock-W */
418: TSSetIFunction(ts,NULL,problem->function,problem->data);
419: TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
420: TSSetDuration(ts,maxsteps,problem->final_time);
421: TSSetMaxStepRejections(ts,10);
422: TSSetMaxSNESFailures(ts,-1); /* unlimited */
423: if (use_monitor) {
424: TSMonitorSet(ts,&MonitorError,&mon,NULL);
425: }
427: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
428: Set initial conditions
429: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
430: (*problem->solution)(0,x,problem->data);
431: TSSetInitialTimeStep(ts,0.0,.001);
432: TSSetSolution(ts,x);
434: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
435: Set runtime options
436: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
437: TSSetFromOptions(ts);
439: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
440: Solve nonlinear system
441: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
442: TSSolve(ts,x);
443: TSGetSolveTime(ts,&ftime);
444: TSGetTimeStepNumber(ts,&steps);
445: TSGetSNESFailures(ts,&snesfails);
446: TSGetStepRejections(ts,&rejects);
447: TSGetSNESIterations(ts,&nonlinits);
448: TSGetKSPIterations(ts,&linits);
449: 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);
450: if (problem->hasexact) {
451: MonitorError(ts,steps,ftime,x,&mon);
452: }
454: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
455: Free work space. All PETSc objects should be destroyed when they
456: are no longer needed.
457: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
458: MatDestroy(&A);
459: VecDestroy(&x);
460: VecDestroy(&r);
461: VecDestroy(&mon.x);
462: TSDestroy(&ts);
463: if (problem->destroy) {
464: (*problem->destroy)(problem);
465: }
466: PetscFree(problem);
467: PetscFunctionListDestroy(&plist);
469: PetscFinalize();
470: return(0);
471: }