Actual source code: ex8.c
2: /* Program usage: ./ex8 [-help] [all PETSc options] */
4: static char help[] = "Nonlinear DAE.\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: void *data;
29: };
31: /*
32: * User-defined routines
33: */
35: /*
36: * Stiff 3-variable system from chemical reactions, due to Robertson (1966), problem ROBER in Hairer&Wanner, ODE 2, 1996
37: */
40: static PetscErrorCode RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
41: {
43: PetscScalar *x,*xdot,*f;
46: VecGetArray(X,&x);
47: VecGetArray(Xdot,&xdot);
48: VecGetArray(F,&f);
49: f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2];
50: f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]);
51: f[2] = xdot[2] - 3e7*PetscSqr(x[1]);
52: VecRestoreArray(X,&x);
53: VecRestoreArray(Xdot,&xdot);
54: VecRestoreArray(F,&f);
55: return(0);
56: }
60: static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
61: {
63: PetscInt rowcol[] = {0,1,2};
64: PetscScalar *x,*xdot,J[3][3];
67: VecGetArray(X,&x);
68: VecGetArray(Xdot,&xdot);
69: J[0][0] = a + 0.04; J[0][1] = -1e4*x[2]; J[0][2] = -1e4*x[1];
70: J[1][0] = -0.04; J[1][1] = a + 1e4*x[2] + 3e7*2*x[1]; J[1][2] = 1e4*x[1];
71: J[2][0] = 0; J[2][1] = -3e7*2*x[1]; J[2][2] = a;
72: MatSetValues(*B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
73: VecRestoreArray(X,&x);
74: VecRestoreArray(Xdot,&xdot);
76: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
78: if (*A != *B) {
79: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
81: }
82: *flag = SAME_NONZERO_PATTERN;
83: return(0);
84: }
88: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
89: {
91: PetscScalar *x;
94: if (t != 0) SETERRQ(1,"not implemented");
95: VecGetArray(X,&x);
96: x[0] = 1;
97: x[1] = 0;
98: x[2] = 0;
99: VecRestoreArray(X,&x);
100: return(0);
101: }
105: static PetscErrorCode RoberCreate(Problem p)
106: {
109: p->destroy = 0;
110: p->function = &RoberFunction;
111: p->jacobian = &RoberJacobian;
112: p->solution = &RoberSolution;
113: p->final_time = 1e11;
114: p->n = 3;
115: return(0);
116: }
119: typedef struct {
120: PetscReal lambda;
121: } CECtx;
123: /*
124: * Stiff scalar valued problem with an exact solution
125: */
128: static PetscErrorCode CEDestroy(Problem p)
129: {
133: PetscFree(p->data);
134: return(0);
135: }
139: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
140: {
142: PetscReal l = ((CECtx*)ctx)->lambda;
143: PetscScalar *x,*xdot,*f;
146: VecGetArray(X,&x);
147: VecGetArray(Xdot,&xdot);
148: VecGetArray(F,&f);
149: f[0] = xdot[0] + l*(x[0] - cos(t));
150: #if 0
151: PetscPrintf(PETSC_COMM_WORLD," f(t=%G,x=%G,xdot=%G) = %G\n",t,x[0],xdot[0],f[0]);
152: #endif
153: VecRestoreArray(X,&x);
154: VecRestoreArray(Xdot,&xdot);
155: VecRestoreArray(F,&f);
156: return(0);
157: }
161: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
162: {
163: PetscReal l = ((CECtx*)ctx)->lambda;
165: PetscInt rowcol[] = {0};
166: PetscScalar *x,*xdot,J[1][1];
169: VecGetArray(X,&x);
170: VecGetArray(Xdot,&xdot);
171: J[0][0] = a + l;
172: MatSetValues(*B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
173: VecRestoreArray(X,&x);
174: VecRestoreArray(Xdot,&xdot);
176: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
177: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
178: if (*A != *B) {
179: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
180: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
181: }
182: *flag = SAME_NONZERO_PATTERN;
183: return(0);
184: }
188: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
189: {
190: PetscReal l = ((CECtx*)ctx)->lambda;
192: PetscScalar *x;
195: VecGetArray(X,&x);
196: x[0] = l/(l*l+1)*(l*cos(t)+sin(t)) - l*l/(l*l+1)*exp(-l*t);
197: VecRestoreArray(X,&x);
198: return(0);
199: }
203: static PetscErrorCode CECreate(Problem p)
204: {
206: CECtx *ce;
209: PetscMalloc(sizeof(CECtx),&ce);
210: p->data = (void*)ce;
212: p->destroy = &CEDestroy;
213: p->function = &CEFunction;
214: p->jacobian = &CEJacobian;
215: p->solution = &CESolution;
216: p->final_time = 10;
217: p->n = 1;
219: ce->lambda = 10;
220: PetscOptionsBegin(p->comm,PETSC_NULL,"CE options","");
221: {
222: PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,PETSC_NULL);
223: }
224: PetscOptionsEnd();
225: return(0);
226: }
229: /*
230: * User-defined monitor for comparing to exact solutions when possible
231: */
232: typedef struct {
233: MPI_Comm comm;
234: Problem problem;
235: Vec x;
236: } MonitorCtx;
240: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
241: {
243: MonitorCtx *mon = (MonitorCtx*)ctx;
244: PetscReal h,nrm_x,nrm_exact,nrm_diff;
247: if (!mon->problem->solution) return(0);
248: (*mon->problem->solution)(t,mon->x,mon->problem->data);
249: VecNorm(x,NORM_2,&nrm_x);
250: VecNorm(mon->x,NORM_2,&nrm_exact);
251: VecAYPX(mon->x,-1,x);
252: VecNorm(mon->x,NORM_2,&nrm_diff);
253: TSGetTimeStep(ts,&h);
254: 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);
255: return(0);
256: }
262: int main(int argc,char **argv)
263: {
264: PetscFList plist = PETSC_NULL;
265: char pname[256];
266: TS ts; /* nonlinear solver */
267: Vec x,r; /* solution, residual vectors */
268: Mat A; /* Jacobian matrix */
269: Problem problem;
270: PetscTruth use_monitor;
271: PetscInt steps,maxsteps = 100;
272: PetscReal ftime;
273: MonitorCtx mon;
274: PetscErrorCode ierr;
276: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277: Initialize program
278: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
279: PetscInitialize(&argc,&argv,(char *)0,help);
281: /* Register the available problems */
282: PetscFListAdd(&plist,"rober","",(void(*)(void))&RoberCreate);
283: PetscFListAdd(&plist,"ce", "",(void(*)(void))&CECreate );
284: PetscStrcpy(pname,"ce");
286: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
287: Set runtime options
288: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
289: PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Timestepping benchmark options","");
290: {
291: PetscOptionsList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),PETSC_NULL);
292: use_monitor = PETSC_FALSE;
293: PetscOptionsTruth("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,PETSC_NULL);
294: }
295: PetscOptionsEnd();
297: /* Create the new problem */
298: PetscNew(struct _Problem,&problem);
299: problem->comm = MPI_COMM_WORLD;
300: {
301: PetscErrorCode (*pcreate)(Problem);
303: PetscFListFind(plist,MPI_COMM_WORLD,pname,(void(**)(void))&pcreate);
304: if (!pcreate) SETERRQ1(1,"No problem '%s'",pname);
305: (*pcreate)(problem);
306: }
309: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
310: Create necessary matrix and vectors, solve same ODE on every process
311: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
312: MatCreate(PETSC_COMM_WORLD,&A);
313: MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
314: MatSetFromOptions(A);
316: MatGetVecs(A,&x,PETSC_NULL);
317: VecDuplicate(x,&r);
319: mon.comm = PETSC_COMM_WORLD;
320: mon.problem = problem;
321: VecDuplicate(x,&mon.x);
323: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
324: Create timestepping solver context
325: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
326: TSCreate(PETSC_COMM_WORLD,&ts);
327: TSSetProblemType(ts,TS_NONLINEAR);
328: TSSetType(ts,TSGL); /* General Linear method, TSTHETA can also solve DAE */
329: TSSetIFunction(ts,problem->function,problem->data);
330: TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
331: TSSetDuration(ts,maxsteps,problem->final_time);
332: if (use_monitor) {
333: TSMonitorSet(ts,&MonitorError,&mon,PETSC_NULL);
334: }
336: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
337: Set initial conditions
338: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
339: (*problem->solution)(0,x,problem->data);
340: TSSetInitialTimeStep(ts,0.0,.001);
341: TSSetSolution(ts,x);
343: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
344: Set runtime options
345: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
346: TSSetFromOptions(ts);
348: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349: Solve nonlinear system
350: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
351: TSStep(ts,&steps,&ftime);
352: PetscPrintf(PETSC_COMM_WORLD,"steps %D, ftime %G\n",steps,ftime);
354: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
355: Free work space. All PETSc objects should be destroyed when they
356: are no longer needed.
357: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
358: MatDestroy(A);
359: VecDestroy(x);
360: VecDestroy(r);
361: VecDestroy(mon.x);
362: TSDestroy(ts);
363: if (problem->destroy) {
364: (*problem->destroy)(problem);
365: }
366: PetscFree(problem);
367: PetscFListDestroy(&plist);
369: PetscFinalize();
370: return(0);
371: }