Actual source code: ex22.c

      1: static const char help[] = "Time-dependent advection-reaction PDE in 1d, demonstrates IMEX methods.\n";
      2: /*
      3:    u_t + a1*u_x = -k1*u + k2*v + s1
      4:    v_t + a2*v_x = k1*u - k2*v + s2
      5:    0 < x < 1;
      6:    a1 = 1, k1 = 10^6, s1 = 0,
      7:    a2 = 0, k2 = 2*k1, s2 = 1

      9:    Initial conditions:
     10:    u(x,0) = 1 + s2*x
     11:    v(x,0) = k0/k1*u(x,0) + s1/k1

     13:    Upstream boundary conditions:
     14:    u(0,t) = 1-sin(12*t)^4

     16:    Useful command-line parameters:
     17:    -ts_arkimex_fully_implicit - treats advection implicitly, only relevant with -ts_type arkimex (default)
     18:    -ts_type rosw - use Rosenbrock-W method (linearly implicit IMEX, amortizes assembly and preconditioner setup)
     19: */

     21: #include <petscdm.h>
     22: #include <petscdmda.h>
     23: #include <petscts.h>

     25: typedef PetscScalar Field[2];

     27: typedef struct _User *User;
     28: struct _User {
     29:   PetscReal a[2]; /* Advection speeds */
     30:   PetscReal k[2]; /* Reaction coefficients */
     31:   PetscReal s[2]; /* Source terms */
     32: };

     34: static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
     35: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
     36: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
     37: static PetscErrorCode FormInitialSolution(TS, Vec, void *);

     39: int main(int argc, char **argv)
     40: {
     41:   TS                ts;         /* time integrator */
     42:   SNES              snes;       /* nonlinear solver */
     43:   SNESLineSearch    linesearch; /* line search */
     44:   Vec               X;          /* solution, residual vectors */
     45:   Mat               J;          /* Jacobian matrix */
     46:   PetscInt          steps, mx;
     47:   DM                da;
     48:   PetscReal         ftime, dt;
     49:   struct _User      user; /* user-defined work context */
     50:   TSConvergedReason reason;

     52:   PetscFunctionBeginUser;
     53:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

     55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     56:      Create distributed array (DMDA) to manage parallel grid and vectors
     57:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     58:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da));
     59:   PetscCall(DMSetFromOptions(da));
     60:   PetscCall(DMSetUp(da));

     62:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     63:      Extract global vectors from DMDA;
     64:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     65:   PetscCall(DMCreateGlobalVector(da, &X));

     67:   /* Initialize user application context */
     68:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "");
     69:   {
     70:     user.a[0] = 1;
     71:     PetscCall(PetscOptionsReal("-a0", "Advection rate 0", "", user.a[0], &user.a[0], NULL));
     72:     user.a[1] = 0;
     73:     PetscCall(PetscOptionsReal("-a1", "Advection rate 1", "", user.a[1], &user.a[1], NULL));
     74:     user.k[0] = 1e6;
     75:     PetscCall(PetscOptionsReal("-k0", "Reaction rate 0", "", user.k[0], &user.k[0], NULL));
     76:     user.k[1] = 2 * user.k[0];
     77:     PetscCall(PetscOptionsReal("-k1", "Reaction rate 1", "", user.k[1], &user.k[1], NULL));
     78:     user.s[0] = 0;
     79:     PetscCall(PetscOptionsReal("-s0", "Source 0", "", user.s[0], &user.s[0], NULL));
     80:     user.s[1] = 1;
     81:     PetscCall(PetscOptionsReal("-s1", "Source 1", "", user.s[1], &user.s[1], NULL));
     82:   }
     83:   PetscOptionsEnd();

     85:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     86:      Create timestepping solver context
     87:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     88:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     89:   PetscCall(TSSetDM(ts, da));
     90:   PetscCall(TSSetType(ts, TSARKIMEX));
     91:   PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
     92:   PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user));
     93:   PetscCall(DMSetMatType(da, MATAIJ));
     94:   PetscCall(DMCreateMatrix(da, &J));
     95:   PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user));

     97:   /* A line search in the nonlinear solve can fail due to ill-conditioning unless an absolute tolerance is set. Since
     98:    * this problem is linear, we deactivate the line search. For a linear problem, it is usually recommended to also use
     99:    * SNESSetType(snes,SNESKSPONLY). */
    100:   PetscCall(TSGetSNES(ts, &snes));
    101:   PetscCall(SNESGetLineSearch(snes, &linesearch));
    102:   PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC));

    104:   ftime = .1;
    105:   PetscCall(TSSetMaxTime(ts, ftime));
    106:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    108:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    109:      Set initial conditions
    110:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    111:   PetscCall(FormInitialSolution(ts, X, &user));
    112:   PetscCall(TSSetSolution(ts, X));
    113:   PetscCall(VecGetSize(X, &mx));
    114:   dt = .1 * PetscMax(user.a[0], user.a[1]) / mx; /* Advective CFL, I don't know why it needs so much safety factor. */
    115:   PetscCall(TSSetTimeStep(ts, dt));

    117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    118:      Set runtime options
    119:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    120:   PetscCall(TSSetFromOptions(ts));

    122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    123:      Solve nonlinear system
    124:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    125:   PetscCall(TSSolve(ts, X));
    126:   PetscCall(TSGetSolveTime(ts, &ftime));
    127:   PetscCall(TSGetStepNumber(ts, &steps));
    128:   PetscCall(TSGetConvergedReason(ts, &reason));
    129:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));

    131:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    132:      Free work space.
    133:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    134:   PetscCall(MatDestroy(&J));
    135:   PetscCall(VecDestroy(&X));
    136:   PetscCall(TSDestroy(&ts));
    137:   PetscCall(DMDestroy(&da));
    138:   PetscCall(PetscFinalize());
    139:   return 0;
    140: }

    142: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
    143: {
    144:   User          user = (User)ptr;
    145:   DM            da;
    146:   DMDALocalInfo info;
    147:   PetscInt      i;
    148:   Field        *f;
    149:   const Field  *x, *xdot;

    151:   PetscFunctionBeginUser;
    152:   PetscCall(TSGetDM(ts, &da));
    153:   PetscCall(DMDAGetLocalInfo(da, &info));

    155:   /* Get pointers to vector data */
    156:   PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x));
    157:   PetscCall(DMDAVecGetArrayRead(da, Xdot, (void *)&xdot));
    158:   PetscCall(DMDAVecGetArray(da, F, &f));

    160:   /* Compute function over the locally owned part of the grid */
    161:   for (i = info.xs; i < info.xs + info.xm; i++) {
    162:     f[i][0] = xdot[i][0] + user->k[0] * x[i][0] - user->k[1] * x[i][1] - user->s[0];
    163:     f[i][1] = xdot[i][1] - user->k[0] * x[i][0] + user->k[1] * x[i][1] - user->s[1];
    164:   }

    166:   /* Restore vectors */
    167:   PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x));
    168:   PetscCall(DMDAVecRestoreArrayRead(da, Xdot, (void *)&xdot));
    169:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
    174: {
    175:   User          user = (User)ptr;
    176:   DM            da;
    177:   Vec           Xloc;
    178:   DMDALocalInfo info;
    179:   PetscInt      i, j;
    180:   PetscReal     hx;
    181:   Field        *f;
    182:   const Field  *x;

    184:   PetscFunctionBeginUser;
    185:   PetscCall(TSGetDM(ts, &da));
    186:   PetscCall(DMDAGetLocalInfo(da, &info));
    187:   hx = 1.0 / (PetscReal)info.mx;

    189:   /*
    190:      Scatter ghost points to local vector,using the 2-step process
    191:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    192:      By placing code between these two statements, computations can be
    193:      done while messages are in transition.
    194:   */
    195:   PetscCall(DMGetLocalVector(da, &Xloc));
    196:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    197:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));

    199:   /* Get pointers to vector data */
    200:   PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
    201:   PetscCall(DMDAVecGetArray(da, F, &f));

    203:   /* Compute function over the locally owned part of the grid */
    204:   for (i = info.xs; i < info.xs + info.xm; i++) {
    205:     const PetscReal *a = user->a;
    206:     PetscReal        u0t[2];
    207:     u0t[0] = 1.0 - PetscPowRealInt(PetscSinReal(12 * t), 4);
    208:     u0t[1] = 0.0;
    209:     for (j = 0; j < 2; j++) {
    210:       if (i == 0) f[i][j] = a[j] / hx * (1. / 3 * u0t[j] + 0.5 * x[i][j] - x[i + 1][j] + 1. / 6 * x[i + 2][j]);
    211:       else if (i == 1) f[i][j] = a[j] / hx * (-1. / 12 * u0t[j] + 2. / 3 * x[i - 1][j] - 2. / 3 * x[i + 1][j] + 1. / 12 * x[i + 2][j]);
    212:       else if (i == info.mx - 2) f[i][j] = a[j] / hx * (-1. / 6 * x[i - 2][j] + x[i - 1][j] - 0.5 * x[i][j] - 1. / 3 * x[i + 1][j]);
    213:       else if (i == info.mx - 1) f[i][j] = a[j] / hx * (-x[i][j] + x[i - 1][j]);
    214:       else f[i][j] = a[j] / hx * (-1. / 12 * x[i - 2][j] + 2. / 3 * x[i - 1][j] - 2. / 3 * x[i + 1][j] + 1. / 12 * x[i + 2][j]);
    215:     }
    216:   }

    218:   /* Restore vectors */
    219:   PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
    220:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    221:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    222:   PetscFunctionReturn(PETSC_SUCCESS);
    223: }

    225: /* --------------------------------------------------------------------- */
    226: /*
    227:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
    228: */
    229: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
    230: {
    231:   User          user = (User)ptr;
    232:   DMDALocalInfo info;
    233:   PetscInt      i;
    234:   DM            da;
    235:   const Field  *x, *xdot;

    237:   PetscFunctionBeginUser;
    238:   PetscCall(TSGetDM(ts, &da));
    239:   PetscCall(DMDAGetLocalInfo(da, &info));

    241:   /* Get pointers to vector data */
    242:   PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x));
    243:   PetscCall(DMDAVecGetArrayRead(da, Xdot, (void *)&xdot));

    245:   /* Compute function over the locally owned part of the grid */
    246:   for (i = info.xs; i < info.xs + info.xm; i++) {
    247:     const PetscReal *k = user->k;
    248:     PetscScalar      v[2][2];

    250:     v[0][0] = a + k[0];
    251:     v[0][1] = -k[1];
    252:     v[1][0] = -k[0];
    253:     v[1][1] = a + k[1];
    254:     PetscCall(MatSetValuesBlocked(Jpre, 1, &i, 1, &i, &v[0][0], INSERT_VALUES));
    255:   }

    257:   /* Restore vectors */
    258:   PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x));
    259:   PetscCall(DMDAVecRestoreArrayRead(da, Xdot, (void *)&xdot));

    261:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
    262:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
    263:   if (J != Jpre) {
    264:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    265:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    266:   }
    267:   PetscFunctionReturn(PETSC_SUCCESS);
    268: }

    270: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
    271: {
    272:   User          user = (User)ctx;
    273:   DM            da;
    274:   PetscInt      i;
    275:   DMDALocalInfo info;
    276:   Field        *x;
    277:   PetscReal     hx;

    279:   PetscFunctionBeginUser;
    280:   PetscCall(TSGetDM(ts, &da));
    281:   PetscCall(DMDAGetLocalInfo(da, &info));
    282:   hx = 1.0 / (PetscReal)info.mx;

    284:   /* Get pointers to vector data */
    285:   PetscCall(DMDAVecGetArray(da, X, &x));

    287:   /* Compute function over the locally owned part of the grid */
    288:   for (i = info.xs; i < info.xs + info.xm; i++) {
    289:     PetscReal r = (i + 1) * hx, ik = user->k[1] != 0.0 ? 1.0 / user->k[1] : 1.0;
    290:     x[i][0] = 1 + user->s[1] * r;
    291:     x[i][1] = user->k[0] * ik * x[i][0] + user->s[1] * ik;
    292:   }
    293:   PetscCall(DMDAVecRestoreArray(da, X, &x));
    294:   PetscFunctionReturn(PETSC_SUCCESS);
    295: }

    297: /*TEST

    299:     test:
    300:       args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_arkimex_type 4 -ts_adapt_type none -ts_dt .005 -ts_max_time .1
    301:       requires: !single

    303:     test:
    304:       suffix: 2
    305:       args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_type rosw -ts_dt 1e-3 -ts_adapt_type none -ts_dt .005 -ts_max_time .1
    306:       nsize: 2

    308:     test:
    309:       suffix: 3
    310:       args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type ra34pw2 -ts_dt 5e-3 -ts_max_time .1 -ts_adapt_type none
    311:       nsize: 2

    313:     test:
    314:       suffix: 4
    315:       args: -ts_type eimex -da_grid_x 200 -ts_eimex_order_adapt true -ts_dt 0.001 -ts_monitor -ts_max_steps 100
    316:       filter: sed "s/ITS/TIME/g"
    317:       nsize: 2

    319: TEST*/