Actual source code: ex8.c
petsc-3.12.4 2020-02-04
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: */
     31: static PetscErrorCode RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
     32: {
     33:   PetscErrorCode    ierr;
     34:   PetscScalar       *f;
     35:   const PetscScalar *x,*xdot;

     38:   VecGetArrayRead(X,&x);
     39:   VecGetArrayRead(Xdot,&xdot);
     40:   VecGetArray(F,&f);
     41:   f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2];
     42:   f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]);
     43:   f[2] = xdot[2] - 3e7*PetscSqr(x[1]);
     44:   VecRestoreArrayRead(X,&x);
     45:   VecRestoreArrayRead(Xdot,&xdot);
     46:   VecRestoreArray(F,&f);
     47:   return(0);
     48: }

     50: static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
     51: {
     52:   PetscErrorCode    ierr;
     53:   PetscInt          rowcol[] = {0,1,2};
     54:   PetscScalar       J[3][3];
     55:   const PetscScalar *x,*xdot;

     58:   VecGetArrayRead(X,&x);
     59:   VecGetArrayRead(Xdot,&xdot);
     60:   J[0][0] = a + 0.04;     J[0][1] = -1e4*x[2];                   J[0][2] = -1e4*x[1];
     61:   J[1][0] = -0.04;        J[1][1] = a + 1e4*x[2] + 3e7*2*x[1];   J[1][2] = 1e4*x[1];
     62:   J[2][0] = 0;            J[2][1] = -3e7*2*x[1];                 J[2][2] = a;
     63:   MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
     64:   VecRestoreArrayRead(X,&x);
     65:   VecRestoreArrayRead(Xdot,&xdot);

     67:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     68:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
     69:   if (A != B) {
     70:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
     71:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
     72:   }
     73:   return(0);
     74: }

     76: static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
     77: {
     79:   PetscScalar    *x;

     82:   if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
     83:   VecGetArray(X,&x);
     84:   x[0] = 1;
     85:   x[1] = 0;
     86:   x[2] = 0;
     87:   VecRestoreArray(X,&x);
     88:   return(0);
     89: }

     91: static PetscErrorCode RoberCreate(Problem p)
     92: {

     95:   p->destroy    = 0;
     96:   p->function   = &RoberFunction;
     97:   p->jacobian   = &RoberJacobian;
     98:   p->solution   = &RoberSolution;
     99:   p->final_time = 1e11;
    100:   p->n          = 3;
    101:   return(0);
    102: }

    104: /*
    105:      Stiff scalar valued problem
    106: */

    108: typedef struct {
    109:   PetscReal lambda;
    110: } CECtx;

    112: static PetscErrorCode CEDestroy(Problem p)
    113: {

    117:   PetscFree(p->data);
    118:   return(0);
    119: }

    121: static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
    122: {
    123:   PetscErrorCode    ierr;
    124:   PetscReal         l = ((CECtx*)ctx)->lambda;
    125:   PetscScalar       *f;
    126:   const PetscScalar *x,*xdot;

    129:   VecGetArrayRead(X,&x);
    130:   VecGetArrayRead(Xdot,&xdot);
    131:   VecGetArray(F,&f);
    132:   f[0] = xdot[0] + l*(x[0] - PetscCosReal(t));
    133: #if 0
    134:   PetscPrintf(PETSC_COMM_WORLD," f(t=%g,x=%g,xdot=%g) = %g\n",(double)t,(double)x[0],(double)xdot[0],(double)f[0]);
    135: #endif
    136:   VecRestoreArrayRead(X,&x);
    137:   VecRestoreArrayRead(Xdot,&xdot);
    138:   VecRestoreArray(F,&f);
    139:   return(0);
    140: }

    142: static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
    143: {
    144:   PetscReal         l = ((CECtx*)ctx)->lambda;
    145:   PetscErrorCode    ierr;
    146:   PetscInt          rowcol[] = {0};
    147:   PetscScalar       J[1][1];
    148:   const PetscScalar *x,*xdot;

    151:   VecGetArrayRead(X,&x);
    152:   VecGetArrayRead(Xdot,&xdot);
    153:   J[0][0] = a + l;
    154:   MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);
    155:   VecRestoreArrayRead(X,&x);
    156:   VecRestoreArrayRead(Xdot,&xdot);

    158:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    159:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    160:   if (A != B) {
    161:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    162:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    163:   }
    164:   return(0);
    165: }

    167: static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
    168: {
    169:   PetscReal      l = ((CECtx*)ctx)->lambda;
    171:   PetscScalar    *x;

    174:   VecGetArray(X,&x);
    175:   x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t);
    176:   VecRestoreArray(X,&x);
    177:   return(0);
    178: }

    180: static PetscErrorCode CECreate(Problem p)
    181: {
    183:   CECtx          *ce;

    186:   PetscMalloc(sizeof(CECtx),&ce);
    187:   p->data = (void*)ce;

    189:   p->destroy    = &CEDestroy;
    190:   p->function   = &CEFunction;
    191:   p->jacobian   = &CEJacobian;
    192:   p->solution   = &CESolution;
    193:   p->final_time = 10;
    194:   p->n          = 1;
    195:   p->hasexact   = PETSC_TRUE;

    197:   ce->lambda = 10;
    198:   PetscOptionsBegin(p->comm,NULL,"CE options","");
    199:   {
    200:     PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);
    201:   }
    202:   PetscOptionsEnd();
    203:   return(0);
    204: }

    206: /*
    207: *  Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
    208: */
    209: static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
    210: {
    211:   PetscErrorCode    ierr;
    212:   PetscScalar       *f;
    213:   const PetscScalar *x,*xdot;

    216:   VecGetArrayRead(X,&x);
    217:   VecGetArrayRead(Xdot,&xdot);
    218:   VecGetArray(F,&f);
    219:   f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
    220:   f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
    221:   f[2] = xdot[2] - 0.161*(x[0] - x[2]);
    222:   VecRestoreArrayRead(X,&x);
    223:   VecRestoreArrayRead(Xdot,&xdot);
    224:   VecRestoreArray(F,&f);
    225:   return(0);
    226: }

    228: static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
    229: {
    230:   PetscErrorCode    ierr;
    231:   PetscInt          rowcol[] = {0,1,2};
    232:   PetscScalar       J[3][3];
    233:   const PetscScalar *x,*xdot;

    236:   VecGetArrayRead(X,&x);
    237:   VecGetArrayRead(Xdot,&xdot);
    238:   J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
    239:   J[0][1] = -77.27*(1. - x[0]);
    240:   J[0][2] = 0;
    241:   J[1][0] = 1./77.27*x[1];
    242:   J[1][1] = a + 1./77.27*(1. + x[0]);
    243:   J[1][2] = -1./77.27;
    244:   J[2][0] = -0.161;
    245:   J[2][1] = 0;
    246:   J[2][2] = a + 0.161;
    247:   MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);
    248:   VecRestoreArrayRead(X,&x);
    249:   VecRestoreArrayRead(Xdot,&xdot);

    251:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    252:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    253:   if (A != B) {
    254:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    255:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    256:   }
    257:   return(0);
    258: }

    260: static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
    261: {
    263:   PetscScalar    *x;

    266:   if (t != 0) SETERRQ(PETSC_COMM_SELF,1,"not implemented");
    267:   VecGetArray(X,&x);
    268:   x[0] = 1;
    269:   x[1] = 2;
    270:   x[2] = 3;
    271:   VecRestoreArray(X,&x);
    272:   return(0);
    273: }

    275: static PetscErrorCode OregoCreate(Problem p)
    276: {

    279:   p->destroy    = 0;
    280:   p->function   = &OregoFunction;
    281:   p->jacobian   = &OregoJacobian;
    282:   p->solution   = &OregoSolution;
    283:   p->final_time = 360;
    284:   p->n          = 3;
    285:   return(0);
    286: }

    289: /*
    290: *  User-defined monitor for comparing to exact solutions when possible
    291: */
    292: typedef struct {
    293:   MPI_Comm comm;
    294:   Problem  problem;
    295:   Vec      x;
    296: } MonitorCtx;

    298: static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
    299: {
    301:   MonitorCtx     *mon = (MonitorCtx*)ctx;
    302:   PetscReal      h,nrm_x,nrm_exact,nrm_diff;

    305:   if (!mon->problem->solution) return(0);
    306:   (*mon->problem->solution)(t,mon->x,mon->problem->data);
    307:   VecNorm(x,NORM_2,&nrm_x);
    308:   VecNorm(mon->x,NORM_2,&nrm_exact);
    309:   VecAYPX(mon->x,-1,x);
    310:   VecNorm(mon->x,NORM_2,&nrm_diff);
    311:   TSGetTimeStep(ts,&h);
    312:   if (step < 0) {
    313:     PetscPrintf(mon->comm,"Interpolated final solution ");
    314:   }
    315:   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);
    316:   return(0);
    317: }

    320: int main(int argc,char **argv)
    321: {
    322:   PetscFunctionList plist = NULL;
    323:   char              pname[256];
    324:   TS                ts;            /* nonlinear solver */
    325:   Vec               x,r;           /* solution, residual vectors */
    326:   Mat               A;             /* Jacobian matrix */
    327:   Problem           problem;
    328:   PetscBool         use_monitor;
    329:   PetscInt          steps,nonlinits,linits,snesfails,rejects;
    330:   PetscReal         ftime;
    331:   MonitorCtx        mon;
    332:   PetscErrorCode    ierr;
    333:   PetscMPIInt       size;

    335:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    336:      Initialize program
    337:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    338:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
    339:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
    340:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

    342:   /* Register the available problems */
    343:   PetscFunctionListAdd(&plist,"rober",&RoberCreate);
    344:   PetscFunctionListAdd(&plist,"ce",&CECreate);
    345:   PetscFunctionListAdd(&plist,"orego",&OregoCreate);
    346:   PetscStrcpy(pname,"ce");

    348:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    349:     Set runtime options
    350:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    351:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");
    352:   {
    353:     PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);
    354:     use_monitor = PETSC_FALSE;
    355:     PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);
    356:   }
    357:   PetscOptionsEnd();

    359:   /* Create the new problem */
    360:   PetscNew(&problem);
    361:   problem->comm = MPI_COMM_WORLD;
    362:   {
    363:     PetscErrorCode (*pcreate)(Problem);

    365:     PetscFunctionListFind(plist,pname,&pcreate);
    366:     if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname);
    367:     (*pcreate)(problem);
    368:   }

    370:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    371:     Create necessary matrix and vectors
    372:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    373:   MatCreate(PETSC_COMM_WORLD,&A);
    374:   MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);
    375:   MatSetFromOptions(A);
    376:   MatSetUp(A);

    378:   MatCreateVecs(A,&x,NULL);
    379:   VecDuplicate(x,&r);

    381:   mon.comm    = PETSC_COMM_WORLD;
    382:   mon.problem = problem;
    383:   VecDuplicate(x,&mon.x);

    385:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    386:      Create timestepping solver context
    387:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    388:   TSCreate(PETSC_COMM_WORLD,&ts);
    389:   TSSetProblemType(ts,TS_NONLINEAR);
    390:   TSSetType(ts,TSROSW); /* Rosenbrock-W */
    391:   TSSetIFunction(ts,NULL,problem->function,problem->data);
    392:   TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);
    393:   TSSetMaxTime(ts,problem->final_time);
    394:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
    395:   TSSetMaxStepRejections(ts,10);
    396:   TSSetMaxSNESFailures(ts,-1); /* unlimited */
    397:   if (use_monitor) {
    398:     TSMonitorSet(ts,&MonitorError,&mon,NULL);
    399:   }

    401:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    402:      Set initial conditions
    403:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    404:   (*problem->solution)(0,x,problem->data);
    405:   TSSetTimeStep(ts,.001);
    406:   TSSetSolution(ts,x);

    408:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    409:      Set runtime options
    410:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    411:   TSSetFromOptions(ts);

    413:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    414:      Solve nonlinear system
    415:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    416:   TSSolve(ts,x);
    417:   TSGetSolveTime(ts,&ftime);
    418:   TSGetStepNumber(ts,&steps);
    419:   TSGetSNESFailures(ts,&snesfails);
    420:   TSGetStepRejections(ts,&rejects);
    421:   TSGetSNESIterations(ts,&nonlinits);
    422:   TSGetKSPIterations(ts,&linits);
    423:   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);

    425:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    426:      Free work space.  All PETSc objects should be destroyed when they
    427:      are no longer needed.
    428:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    429:   MatDestroy(&A);
    430:   VecDestroy(&x);
    431:   VecDestroy(&r);
    432:   VecDestroy(&mon.x);
    433:   TSDestroy(&ts);
    434:   if (problem->destroy) {
    435:     (*problem->destroy)(problem);
    436:   }
    437:   PetscFree(problem);
    438:   PetscFunctionListDestroy(&plist);

    440:   PetscFinalize();
    441:   return ierr;
    442: }

    444: /*TEST

    446:     test:
    447:       requires: !complex
    448:       args: -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex

    450:     test:
    451:       suffix: 2
    452:       requires: !single !complex
    453:       args: -ts_atol 1e-2 -ts_rtol 1e-2 -ts_max_time 15 -ts_type arkimex -ts_arkimex_type 2e -problem_type orego -ts_arkimex_initial_guess_extrapolate 0 -ts_adapt_time_step_increase_delay 4

    455:     test:
    456:       suffix: 3
    457:       requires: !single !complex
    458:       args: -ts_atol 1e-2 -ts_rtol 1e-2 -ts_max_time 15 -ts_type arkimex -ts_arkimex_type 2e -problem_type orego -ts_arkimex_initial_guess_extrapolate 1

    460: TEST*/