Actual source code: ex8.c

petsc-3.5.0 2014-06-30
Report Typos and Errors
  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: }