Actual source code: ex8.c

      1: static char help[] = "Nonlinear DAE benchmark problems.\n";

      3: /*
      4:    Include "petscts.h" so that we can use TS solvers.  Note that this
      5:    file automatically includes:
      6:      petscsys.h       - base PETSc routines   petscvec.h - vectors
      7:      petscmat.h - matrices
      8:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
      9:      petscviewer.h - viewers               petscpc.h  - preconditioners
     10:      petscksp.h   - linear solvers
     11: */
     12: #include <petscts.h>

     14: typedef struct _Problem *Problem;
     15: struct _Problem {
     16:   PetscErrorCode (*destroy)(Problem);
     17:   TSIFunctionFn *function;
     18:   TSIJacobianFn *jacobian;
     19:   PetscErrorCode (*solution)(PetscReal, Vec, void *);
     20:   MPI_Comm  comm;
     21:   PetscReal final_time;
     22:   PetscInt  n;
     23:   PetscBool hasexact;
     24:   void     *data;
     25: };

     27: /*
     28:       Stiff 3-variable system from chemical reactions, due to Robertson (1966), problem ROBER in Hairer&Wanner, ODE 2, 1996
     29:       https://archimede.uniba.it/~testset/report/rober.pdf
     30: */
     31: static PetscErrorCode RoberFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
     32: {
     33:   PetscScalar       *f;
     34:   const PetscScalar *x, *xdot;

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

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

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

     71:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     72:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     73:   if (A != B) {
     74:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     75:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
     76:   }
     77:   PetscFunctionReturn(PETSC_SUCCESS);
     78: }

     80: static PetscErrorCode RoberSolution(PetscReal t, Vec X, void *ctx)
     81: {
     82:   PetscScalar *x;

     84:   PetscFunctionBeginUser;
     85:   PetscCheck(t == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP, "not implemented");
     86:   PetscCall(VecGetArray(X, &x));
     87:   x[0] = 1;
     88:   x[1] = 0;
     89:   x[2] = 0;
     90:   PetscCall(VecRestoreArray(X, &x));
     91:   PetscFunctionReturn(PETSC_SUCCESS);
     92: }

     94: static PetscErrorCode RoberCreate(Problem p)
     95: {
     96:   PetscFunctionBeginUser;
     97:   p->destroy    = 0;
     98:   p->function   = &RoberFunction;
     99:   p->jacobian   = &RoberJacobian;
    100:   p->solution   = &RoberSolution;
    101:   p->final_time = 1e11;
    102:   p->n          = 3;
    103:   PetscFunctionReturn(PETSC_SUCCESS);
    104: }

    106: /*
    107:      Stiff scalar valued problem
    108: */

    110: typedef struct {
    111:   PetscReal lambda;
    112: } CECtx;

    114: static PetscErrorCode CEDestroy(Problem p)
    115: {
    116:   PetscFunctionBeginUser;
    117:   PetscCall(PetscFree(p->data));
    118:   PetscFunctionReturn(PETSC_SUCCESS);
    119: }

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

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

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

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

    156:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    157:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    158:   if (A != B) {
    159:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    160:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    161:   }
    162:   PetscFunctionReturn(PETSC_SUCCESS);
    163: }

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

    170:   PetscFunctionBeginUser;
    171:   PetscCall(VecGetArray(X, &x));
    172:   x[0] = l / (l * l + 1) * (l * PetscCosReal(t) + PetscSinReal(t)) - l * l / (l * l + 1) * PetscExpReal(-l * t);
    173:   PetscCall(VecRestoreArray(X, &x));
    174:   PetscFunctionReturn(PETSC_SUCCESS);
    175: }

    177: static PetscErrorCode CECreate(Problem p)
    178: {
    179:   CECtx *ce;

    181:   PetscFunctionBeginUser;
    182:   PetscCall(PetscMalloc(sizeof(CECtx), &ce));
    183:   p->data = (void *)ce;

    185:   p->destroy    = &CEDestroy;
    186:   p->function   = &CEFunction;
    187:   p->jacobian   = &CEJacobian;
    188:   p->solution   = &CESolution;
    189:   p->final_time = 10;
    190:   p->n          = 1;
    191:   p->hasexact   = PETSC_TRUE;

    193:   ce->lambda = 10;
    194:   PetscOptionsBegin(p->comm, NULL, "CE options", "");
    195:   {
    196:     PetscCall(PetscOptionsReal("-problem_ce_lambda", "Parameter controlling stiffness: xdot + lambda*(x - cos(t))", "", ce->lambda, &ce->lambda, NULL));
    197:   }
    198:   PetscOptionsEnd();
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

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

    210:   PetscFunctionBeginUser;
    211:   PetscCall(VecGetArrayRead(X, &x));
    212:   PetscCall(VecGetArrayRead(Xdot, &xdot));
    213:   PetscCall(VecGetArray(F, &f));
    214:   f[0] = xdot[0] - 77.27 * (x[1] + x[0] * (1. - 8.375e-6 * x[0] - x[1]));
    215:   f[1] = xdot[1] - 1 / 77.27 * (x[2] - (1. + x[0]) * x[1]);
    216:   f[2] = xdot[2] - 0.161 * (x[0] - x[2]);
    217:   PetscCall(VecRestoreArrayRead(X, &x));
    218:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
    219:   PetscCall(VecRestoreArray(F, &f));
    220:   PetscFunctionReturn(PETSC_SUCCESS);
    221: }

    223: static PetscErrorCode OregoJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
    224: {
    225:   PetscInt           rowcol[] = {0, 1, 2};
    226:   PetscScalar        J[3][3];
    227:   const PetscScalar *x, *xdot;

    229:   PetscFunctionBeginUser;
    230:   PetscCall(VecGetArrayRead(X, &x));
    231:   PetscCall(VecGetArrayRead(Xdot, &xdot));
    232:   J[0][0] = a - 77.27 * ((1. - 8.375e-6 * x[0] - x[1]) - 8.375e-6 * x[0]);
    233:   J[0][1] = -77.27 * (1. - x[0]);
    234:   J[0][2] = 0;
    235:   J[1][0] = 1. / 77.27 * x[1];
    236:   J[1][1] = a + 1. / 77.27 * (1. + x[0]);
    237:   J[1][2] = -1. / 77.27;
    238:   J[2][0] = -0.161;
    239:   J[2][1] = 0;
    240:   J[2][2] = a + 0.161;
    241:   PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
    242:   PetscCall(VecRestoreArrayRead(X, &x));
    243:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));

    245:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    246:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    247:   if (A != B) {
    248:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    249:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    250:   }
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

    254: static PetscErrorCode OregoSolution(PetscReal t, Vec X, void *ctx)
    255: {
    256:   PetscScalar *x;

    258:   PetscFunctionBeginUser;
    259:   PetscCheck(t == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP, "not implemented");
    260:   PetscCall(VecGetArray(X, &x));
    261:   x[0] = 1;
    262:   x[1] = 2;
    263:   x[2] = 3;
    264:   PetscCall(VecRestoreArray(X, &x));
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: static PetscErrorCode OregoCreate(Problem p)
    269: {
    270:   PetscFunctionBeginUser;
    271:   p->destroy    = 0;
    272:   p->function   = &OregoFunction;
    273:   p->jacobian   = &OregoJacobian;
    274:   p->solution   = &OregoSolution;
    275:   p->final_time = 360;
    276:   p->n          = 3;
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }

    280: /*
    281:    User-defined monitor for comparing to exact solutions when possible
    282: */
    283: typedef struct {
    284:   MPI_Comm comm;
    285:   Problem  problem;
    286:   Vec      x;
    287: } MonitorCtx;

    289: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal t, Vec x, void *ctx)
    290: {
    291:   MonitorCtx *mon = (MonitorCtx *)ctx;
    292:   PetscReal   h, nrm_x, nrm_exact, nrm_diff;

    294:   PetscFunctionBeginUser;
    295:   if (!mon->problem->solution) PetscFunctionReturn(PETSC_SUCCESS);
    296:   PetscCall((*mon->problem->solution)(t, mon->x, mon->problem->data));
    297:   PetscCall(VecNorm(x, NORM_2, &nrm_x));
    298:   PetscCall(VecNorm(mon->x, NORM_2, &nrm_exact));
    299:   PetscCall(VecAYPX(mon->x, -1, x));
    300:   PetscCall(VecNorm(mon->x, NORM_2, &nrm_diff));
    301:   PetscCall(TSGetTimeStep(ts, &h));
    302:   if (step < 0) PetscCall(PetscPrintf(mon->comm, "Interpolated final solution "));
    303:   PetscCall(PetscPrintf(mon->comm, "step %4" PetscInt_FMT " 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));
    304:   PetscFunctionReturn(PETSC_SUCCESS);
    305: }

    307: int main(int argc, char **argv)
    308: {
    309:   PetscFunctionList plist = NULL;
    310:   char              pname[256];
    311:   TS                ts;   /* nonlinear solver */
    312:   Vec               x, r; /* solution, residual vectors */
    313:   Mat               A;    /* Jacobian matrix */
    314:   Problem           problem;
    315:   PetscBool         use_monitor = PETSC_FALSE;
    316:   PetscBool         use_result  = PETSC_FALSE;
    317:   PetscInt          steps, nonlinits, linits, snesfails, rejects;
    318:   PetscReal         ftime;
    319:   MonitorCtx        mon;
    320:   PetscMPIInt       size;

    322:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    323:      Initialize program
    324:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    325:   PetscFunctionBeginUser;
    326:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    327:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    328:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

    330:   /* Register the available problems */
    331:   PetscCall(PetscFunctionListAdd(&plist, "rober", &RoberCreate));
    332:   PetscCall(PetscFunctionListAdd(&plist, "ce", &CECreate));
    333:   PetscCall(PetscFunctionListAdd(&plist, "orego", &OregoCreate));
    334:   PetscCall(PetscStrncpy(pname, "ce", sizeof(pname)));

    336:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    337:     Set runtime options
    338:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    339:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Timestepping benchmark options", "");
    340:   {
    341:     PetscCall(PetscOptionsFList("-problem_type", "Name of problem to run", "", plist, pname, pname, sizeof(pname), NULL));
    342:     use_monitor = PETSC_FALSE;
    343:     PetscCall(PetscOptionsBool("-monitor_error", "Display errors relative to exact solutions", "", use_monitor, &use_monitor, NULL));
    344:     PetscCall(PetscOptionsBool("-monitor_result", "Display result", "", use_result, &use_result, NULL));
    345:   }
    346:   PetscOptionsEnd();

    348:   /* Create the new problem */
    349:   PetscCall(PetscNew(&problem));
    350:   problem->comm = MPI_COMM_WORLD;
    351:   {
    352:     PetscErrorCode (*pcreate)(Problem);

    354:     PetscCall(PetscFunctionListFind(plist, pname, &pcreate));
    355:     PetscCheck(pcreate, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "No problem '%s'", pname);
    356:     PetscCall((*pcreate)(problem));
    357:   }

    359:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    360:     Create necessary matrix and vectors
    361:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    362:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    363:   PetscCall(MatSetSizes(A, problem->n, problem->n, PETSC_DETERMINE, PETSC_DETERMINE));
    364:   PetscCall(MatSetFromOptions(A));
    365:   PetscCall(MatSetUp(A));

    367:   PetscCall(MatCreateVecs(A, &x, NULL));
    368:   PetscCall(VecDuplicate(x, &r));

    370:   mon.comm    = PETSC_COMM_WORLD;
    371:   mon.problem = problem;
    372:   PetscCall(VecDuplicate(x, &mon.x));

    374:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    375:      Create timestepping solver context
    376:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    377:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    378:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    379:   PetscCall(TSSetType(ts, TSROSW)); /* Rosenbrock-W */
    380:   PetscCall(TSSetIFunction(ts, NULL, problem->function, problem->data));
    381:   PetscCall(TSSetIJacobian(ts, A, A, problem->jacobian, problem->data));
    382:   PetscCall(TSSetMaxTime(ts, problem->final_time));
    383:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    384:   PetscCall(TSSetMaxStepRejections(ts, 10));
    385:   PetscCall(TSSetMaxSNESFailures(ts, -1)); /* unlimited */
    386:   if (use_monitor) PetscCall(TSMonitorSet(ts, &MonitorError, &mon, NULL));

    388:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    389:      Set initial conditions
    390:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    391:   PetscCall((*problem->solution)(0, x, problem->data));
    392:   PetscCall(TSSetTimeStep(ts, .001));
    393:   PetscCall(TSSetSolution(ts, x));

    395:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    396:      Set runtime options
    397:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    398:   PetscCall(TSSetFromOptions(ts));

    400:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    401:      Solve nonlinear system
    402:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    403:   PetscCall(TSSolve(ts, x));
    404:   PetscCall(TSGetSolveTime(ts, &ftime));
    405:   PetscCall(TSGetStepNumber(ts, &steps));
    406:   PetscCall(TSGetSNESFailures(ts, &snesfails));
    407:   PetscCall(TSGetStepRejections(ts, &rejects));
    408:   PetscCall(TSGetSNESIterations(ts, &nonlinits));
    409:   PetscCall(TSGetKSPIterations(ts, &linits));
    410:   if (use_result) {
    411:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "steps %" PetscInt_FMT " (%" PetscInt_FMT " rejected, %" PetscInt_FMT " SNES fails), ftime %g, nonlinits %" PetscInt_FMT ", linits %" PetscInt_FMT "\n", steps, rejects, snesfails, (double)ftime, nonlinits, linits));
    412:   }

    414:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    415:      Free work space.  All PETSc objects should be destroyed when they
    416:      are no longer needed.
    417:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    418:   PetscCall(MatDestroy(&A));
    419:   PetscCall(VecDestroy(&x));
    420:   PetscCall(VecDestroy(&r));
    421:   PetscCall(VecDestroy(&mon.x));
    422:   PetscCall(TSDestroy(&ts));
    423:   if (problem->destroy) PetscCall((*problem->destroy)(problem));
    424:   PetscCall(PetscFree(problem));
    425:   PetscCall(PetscFunctionListDestroy(&plist));

    427:   PetscCall(PetscFinalize());
    428:   return 0;
    429: }

    431: /*TEST

    433:     test:
    434:       requires: !complex
    435:       args: -monitor_result -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex

    437:     test:
    438:       suffix: 2
    439:       requires: !single !complex
    440:       args: -monitor_result -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

    442:     test:
    443:       suffix: 3
    444:       requires: !single !complex
    445:       args: -monitor_result -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

    447:     test:
    448:       suffix: 4
    449:       output_file: output/empty.out

    451:     test:
    452:       suffix: 5
    453:       args: -snes_lag_jacobian 20 -snes_lag_jacobian_persists
    454:       output_file: output/empty.out

    456: TEST*/