Actual source code: ex8.c
petsc-3.5.4 2015-05-23
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,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: return(0);
76: }
80: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
81: {
83: PetscScalar *x;
86: if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
87: VecGetArray(X,&x);
88: x[0] = 1;
89: x[1] = 0;
90: x[2] = 0;
91: VecRestoreArray(X,&x);
92: return(0);
93: }
97: static PetscErrorCode RoberCreate(Problem p)
98: {
101: p->destroy = 0;
102: p->function = &RoberFunction;
103: p->jacobian = &RoberJacobian;
104: p->solution = &RoberSolution;
105: p->final_time = 1e11;
106: p->n = 3;
107: return(0);
108: }
110: /*
111: Stiff scalar valued problem
112: */
114: typedef struct {
115: PetscReal lambda;
116: } CECtx;
120: static PetscErrorCode CEDestroy(Problem p)
121: {
125: PetscFree(p->data);
126: return(0);
127: }
131: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
132: {
134: PetscReal l = ((CECtx*)ctx)->lambda;
135: PetscScalar *x,*xdot,*f;
138: VecGetArray(X,&x);
139: VecGetArray(Xdot,&xdot);
140: VecGetArray(F,&f);
141: f[0] = xdot[0] + l*(x[0] - PetscCosReal(t));
142: #if 0
143: PetscPrintf(PETSC_COMM_WORLD," f(t=%g,x=%g,xdot=%g) = %g\n",(double)t,(double)x[0],(double)xdot[0],(double)f[0]);
144: #endif
145: VecRestoreArray(X,&x);
146: VecRestoreArray(Xdot,&xdot);
147: VecRestoreArray(F,&f);
148: return(0);
149: }
153: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
154: {
155: PetscReal l = ((CECtx*)ctx)->lambda;
157: PetscInt rowcol[] = {0};
158: PetscScalar *x,*xdot,J[1][1];
161: VecGetArray(X,&x);
162: VecGetArray(Xdot,&xdot);
163: J[0][0] = a + l;
164: MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
165: VecRestoreArray(X,&x);
166: VecRestoreArray(Xdot,&xdot);
168: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
169: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
170: if (A != B) {
171: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
172: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
173: }
174: return(0);
175: }
179: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
180: {
181: PetscReal l = ((CECtx*)ctx)->lambda;
183: PetscScalar *x;
186: VecGetArray(X,&x);
187: x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t);
188: VecRestoreArray(X,&x);
189: return(0);
190: }
194: static PetscErrorCode CECreate(Problem p)
195: {
197: CECtx *ce;
200: PetscMalloc(sizeof(CECtx),&ce);
201: p->data = (void*)ce;
203: p->destroy = &CEDestroy;
204: p->function = &CEFunction;
205: p->jacobian = &CEJacobian;
206: p->solution = &CESolution;
207: p->final_time = 10;
208: p->n = 1;
209: p->hasexact = PETSC_TRUE;
211: ce->lambda = 10;
212: PetscOptionsBegin(p->comm,NULL,"CE options","");
213: {
214: PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);
215: }
216: PetscOptionsEnd();
217: return(0);
218: }
220: /*
221: * Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
222: */
225: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
226: {
228: PetscScalar *x,*xdot,*f;
231: VecGetArray(X,&x);
232: VecGetArray(Xdot,&xdot);
233: VecGetArray(F,&f);
234: f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
235: f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
236: f[2] = xdot[2] - 0.161*(x[0] - x[2]);
237: VecRestoreArray(X,&x);
238: VecRestoreArray(Xdot,&xdot);
239: VecRestoreArray(F,&f);
240: return(0);
241: }
245: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
246: {
248: PetscInt rowcol[] = {0,1,2};
249: PetscScalar *x,*xdot,J[3][3];
252: VecGetArray(X,&x);
253: VecGetArray(Xdot,&xdot);
254: J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
255: J[0][1] = -77.27*(1. - x[0]);
256: J[0][2] = 0;
257: J[1][0] = 1./77.27*x[1];
258: J[1][1] = a + 1./77.27*(1. + x[0]);
259: J[1][2] = -1./77.27;
260: J[2][0] = -0.161;
261: J[2][1] = 0;
262: J[2][2] = a + 0.161;
263: MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
264: VecRestoreArray(X,&x);
265: VecRestoreArray(Xdot,&xdot);
267: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
268: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
269: if (A != B) {
270: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
271: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
272: }
273: return(0);
274: }
278: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
279: {
281: PetscScalar *x;
284: if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
285: VecGetArray(X,&x);
286: x[0] = 1;
287: x[1] = 2;
288: x[2] = 3;
289: VecRestoreArray(X,&x);
290: return(0);
291: }
295: static PetscErrorCode OregoCreate(Problem p)
296: {
299: p->destroy = 0;
300: p->function = &OregoFunction;
301: p->jacobian = &OregoJacobian;
302: p->solution = &OregoSolution;
303: p->final_time = 360;
304: p->n = 3;
305: return(0);
306: }
309: /*
310: * User-defined monitor for comparing to exact solutions when possible
311: */
312: typedef struct {
313: MPI_Comm comm;
314: Problem problem;
315: Vec x;
316: } MonitorCtx;
320: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
321: {
323: MonitorCtx *mon = (MonitorCtx*)ctx;
324: PetscReal h,nrm_x,nrm_exact,nrm_diff;
327: if (!mon->problem->solution) return(0);
328: (*mon->problem->solution)(t,mon->x,mon->problem->data);
329: VecNorm(x,NORM_2,&nrm_x);
330: VecNorm(mon->x,NORM_2,&nrm_exact);
331: VecAYPX(mon->x,-1,x);
332: VecNorm(mon->x,NORM_2,&nrm_diff);
333: TSGetTimeStep(ts,&h);
334: 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);
335: return(0);
336: }
341: int main(int argc,char **argv)
342: {
343: PetscFunctionList plist = NULL;
344: char pname[256];
345: TS ts; /* nonlinear solver */
346: Vec x,r; /* solution, residual vectors */
347: Mat A; /* Jacobian matrix */
348: Problem problem;
349: PetscBool use_monitor;
350: PetscInt steps,maxsteps = 1000,nonlinits,linits,snesfails,rejects;
351: PetscReal ftime;
352: MonitorCtx mon;
353: PetscErrorCode ierr;
354: PetscMPIInt size;
356: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
357: Initialize program
358: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
359: PetscInitialize(&argc,&argv,(char*)0,help);
360: MPI_Comm_size(PETSC_COMM_WORLD,&size);
361: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
363: /* Register the available problems */
364: PetscFunctionListAdd(&plist,"rober",&RoberCreate);
365: PetscFunctionListAdd(&plist,"ce",&CECreate);
366: PetscFunctionListAdd(&plist,"orego",&OregoCreate);
367: PetscStrcpy(pname,"ce");
369: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
370: Set runtime options
371: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
372: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
373: {
374: PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
375: use_monitor = PETSC_FALSE;
376: PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
377: }
378: PetscOptionsEnd();
380: /* Create the new problem */
381: PetscNew(&problem);
382: problem->comm = MPI_COMM_WORLD;
383: {
384: PetscErrorCode (*pcreate)(Problem);
386: PetscFunctionListFind(plist,pname,&pcreate);
387: if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
388: (*pcreate)(problem);
389: }
391: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
392: Create necessary matrix and vectors
393: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
394: MatCreate(PETSC_COMM_WORLD,&A);
395: MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
396: MatSetFromOptions(A);
397: MatSetUp(A);
399: MatGetVecs(A,&x,NULL);
400: VecDuplicate(x,&r);
402: mon.comm = PETSC_COMM_WORLD;
403: mon.problem = problem;
404: VecDuplicate(x,&mon.x);
406: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
407: Create timestepping solver context
408: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
409: TSCreate(PETSC_COMM_WORLD,&ts);
410: TSSetProblemType(ts,TS_NONLINEAR);
411: TSSetType(ts,TSROSW); /* Rosenbrock-W */
412: TSSetIFunction(ts,NULL,problem->function,problem->data);
413: TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
414: TSSetDuration(ts,maxsteps,problem->final_time);
415: TSSetMaxStepRejections(ts,10);
416: TSSetMaxSNESFailures(ts,-1); /* unlimited */
417: if (use_monitor) {
418: TSMonitorSet(ts,&MonitorError,&mon,NULL);
419: }
421: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
422: Set initial conditions
423: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
424: (*problem->solution)(0,x,problem->data);
425: TSSetInitialTimeStep(ts,0.0,.001);
426: TSSetSolution(ts,x);
428: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
429: Set runtime options
430: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
431: TSSetFromOptions(ts);
433: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
434: Solve nonlinear system
435: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
436: TSSolve(ts,x);
437: TSGetSolveTime(ts,&ftime);
438: TSGetTimeStepNumber(ts,&steps);
439: TSGetSNESFailures(ts,&snesfails);
440: TSGetStepRejections(ts,&rejects);
441: TSGetSNESIterations(ts,&nonlinits);
442: TSGetKSPIterations(ts,&linits);
443: 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);
444: if (problem->hasexact) {
445: MonitorError(ts,steps,ftime,x,&mon);
446: }
448: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
449: Free work space. All PETSc objects should be destroyed when they
450: are no longer needed.
451: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
452: MatDestroy(&A);
453: VecDestroy(&x);
454: VecDestroy(&r);
455: VecDestroy(&mon.x);
456: TSDestroy(&ts);
457: if (problem->destroy) {
458: (*problem->destroy)(problem);
459: }
460: PetscFree(problem);
461: PetscFunctionListDestroy(&plist);
463: PetscFinalize();
464: return(0);
465: }