Actual source code: ex27.c

      1: static char help[] = "Time-Dependent Reactive Flow example in 2D with Darcy Flow";

      3: /*

      5: This example solves the elementary chemical reaction:

      7: SP_A + 2SP_B = SP_C

      9: Subject to predetermined flow modeled as though it were in a porous media.
     10: This flow is modeled as advection and diffusion of the three species as

     12: v = porosity*saturation*grad(q)

     14: and the time-dependent equation solved for each species as
     15: advection + diffusion + reaction:

     17: v dot grad SP_i + dSP_i / dt + dispersivity*div*grad(SP_i) + R(SP_i) = 0

     19: The following options are available:

     21: -length_x - The length of the domain in the direction of the flow
     22: -length_y - The length of the domain in the direction orthogonal to the flow

     24: -gradq_inflow - The inflow speed as if the saturation and porosity were 1.
     25: -saturation - The saturation of the porous medium.
     26: -porosity - The porosity of the medium.
     27: -dispersivity - The dispersivity of the flow.
     28: -rate_constant - The rate constant for the chemical reaction
     29: -stoich - The stoichiometry matrix for the reaction
     30: -sp_inflow - The species concentrations at the inflow
     31: -sp_0 - The species concentrations at time 0

     33:  */

     35: #include <petscdm.h>
     36: #include <petscdmda.h>
     37: #include <petscsnes.h>
     38: #include <petscts.h>

     40: #define N_SPECIES   3
     41: #define N_REACTIONS 1
     42: #define DIM         2

     44: #define stoich(i, j) ctx->stoichiometry[N_SPECIES * i + j]

     46: typedef struct {
     47:   PetscScalar sp[N_SPECIES];
     48: } Field;

     50: typedef struct {
     51:   Field     x_inflow;
     52:   Field     x_0;
     53:   PetscReal stoichiometry[N_SPECIES * N_REACTIONS];
     54:   PetscReal porosity;
     55:   PetscReal dispersivity;
     56:   PetscReal saturation;
     57:   PetscReal rate_constant[N_REACTIONS];
     58:   PetscReal gradq_inflow;
     59:   PetscReal length[DIM];
     60: } AppCtx;

     62: extern PetscErrorCode FormInitialGuess(DM da, AppCtx *ctx, Vec X);
     63: extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo *, PetscReal, Field **, Field **, Field **, AppCtx *);
     64: extern PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
     65: extern PetscErrorCode ReactingFlowPostCheck(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);

     67: PetscErrorCode SetFromOptions(AppCtx *ctx)
     68: {
     69:   PetscInt i, j;

     71:   PetscFunctionBeginUser;
     72:   ctx->dispersivity     = 0.5;
     73:   ctx->porosity         = 0.25;
     74:   ctx->saturation       = 1.0;
     75:   ctx->gradq_inflow     = 1.0;
     76:   ctx->rate_constant[0] = 0.5;

     78:   ctx->length[0] = 100.0;
     79:   ctx->length[1] = 100.0;

     81:   ctx->x_0.sp[0] = 0.0;
     82:   ctx->x_0.sp[1] = 0.0;
     83:   ctx->x_0.sp[2] = 0.0;

     85:   ctx->x_inflow.sp[0] = 0.05;
     86:   ctx->x_inflow.sp[1] = 0.05;
     87:   ctx->x_inflow.sp[2] = 0.0;

     89:   for (i = 0; i < N_REACTIONS; i++) {
     90:     for (j = 0; j < N_SPECIES; j++) stoich(i, j) = 0.;
     91:   }

     93:   /* set up a pretty easy example */
     94:   stoich(0, 0) = -1.;
     95:   stoich(0, 1) = -2.;
     96:   stoich(0, 2) = 1.;

     98:   PetscInt as = N_SPECIES;
     99:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-length_x", &ctx->length[0], NULL));
    100:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-length_y", &ctx->length[1], NULL));
    101:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-porosity", &ctx->porosity, NULL));
    102:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-saturation", &ctx->saturation, NULL));
    103:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-dispersivity", &ctx->dispersivity, NULL));
    104:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-gradq_inflow", &ctx->gradq_inflow, NULL));
    105:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-rate_constant", &ctx->rate_constant[0], NULL));
    106:   PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-sp_inflow", ctx->x_inflow.sp, &as, NULL));
    107:   PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-sp_0", ctx->x_0.sp, &as, NULL));
    108:   as = N_SPECIES;
    109:   PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-stoich", ctx->stoichiometry, &as, NULL));
    110:   PetscFunctionReturn(PETSC_SUCCESS);
    111: }

    113: int main(int argc, char **argv)
    114: {
    115:   TS             ts;
    116:   SNES           snes;
    117:   SNESLineSearch linesearch;
    118:   Vec            x;
    119:   AppCtx         ctx;
    120:   DM             da;

    122:   PetscFunctionBeginUser;
    123:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    124:   PetscCall(SetFromOptions(&ctx));
    125:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    126:   PetscCall(TSSetType(ts, TSCN));
    127:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    128:   PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &ctx));
    129:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, -4, -4, PETSC_DECIDE, PETSC_DECIDE, N_SPECIES, 1, NULL, NULL, &da));
    130:   PetscCall(DMSetFromOptions(da));
    131:   PetscCall(DMSetUp(da));
    132:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
    133:   PetscCall(DMDASetFieldName(da, 0, "species A"));
    134:   PetscCall(DMDASetFieldName(da, 1, "species B"));
    135:   PetscCall(DMDASetFieldName(da, 2, "species C"));
    136:   PetscCall(DMSetApplicationContext(da, &ctx));
    137:   PetscCall(DMCreateGlobalVector(da, &x));
    138:   PetscCall(FormInitialGuess(da, &ctx, x));

    140:   PetscCall(TSSetDM(ts, da));
    141:   PetscCall(TSSetMaxTime(ts, 1000.0));
    142:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    143:   PetscCall(TSSetTimeStep(ts, 1.0));
    144:   PetscCall(TSSetSolution(ts, x));
    145:   PetscCall(TSSetFromOptions(ts));

    147:   PetscCall(TSGetSNES(ts, &snes));
    148:   PetscCall(SNESGetLineSearch(snes, &linesearch));
    149:   PetscCall(SNESLineSearchSetPostCheck(linesearch, ReactingFlowPostCheck, (void *)&ctx));
    150:   PetscCall(SNESSetFromOptions(snes));
    151:   PetscCall(TSSolve(ts, x));

    153:   PetscCall(VecDestroy(&x));
    154:   PetscCall(TSDestroy(&ts));
    155:   PetscCall(DMDestroy(&da));
    156:   PetscCall(PetscFinalize());
    157:   PetscFunctionReturn(PETSC_SUCCESS);
    158: }

    160: /* ------------------------------------------------------------------- */

    162: PetscErrorCode FormInitialGuess(DM da, AppCtx *ctx, Vec X)
    163: {
    164:   PetscInt i, j, l, Mx, My, xs, ys, xm, ym;
    165:   Field  **x;

    167:   PetscFunctionBeginUser;
    168:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
    169:   PetscCall(DMDAVecGetArray(da, X, &x));
    170:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

    172:   for (j = ys; j < ys + ym; j++) {
    173:     for (i = xs; i < xs + xm; i++) {
    174:       for (l = 0; l < N_SPECIES; l++) {
    175:         if (i == 0) {
    176:           if (l == 0) x[j][i].sp[l] = (ctx->x_inflow.sp[l] * ((PetscScalar)j) / (My - 1));
    177:           else if (l == 1) x[j][i].sp[l] = (ctx->x_inflow.sp[l] * (1. - ((PetscScalar)j) / (My - 1)));
    178:           else x[j][i].sp[l] = ctx->x_0.sp[l];
    179:         }
    180:       }
    181:     }
    182:   }
    183:   PetscCall(DMDAVecRestoreArray(da, X, &x));
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

    187: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscScalar ptime, Field **x, Field **xt, Field **f, AppCtx *ctx)
    188: {
    189:   PetscInt    i, j, l, m;
    190:   PetscReal   hx, hy, dhx, dhy, hxdhy, hydhx, scale;
    191:   PetscScalar u, uxx, uyy;
    192:   PetscScalar vx, vy, sxp, syp, sxm, sym, avx, vxp, vxm, avy, vyp, vym, f_advect;
    193:   PetscScalar rate;

    195:   PetscFunctionBeginUser;
    196:   hx = ctx->length[0] / ((PetscReal)(info->mx - 1));
    197:   hy = ctx->length[1] / ((PetscReal)(info->my - 1));

    199:   dhx   = 1. / hx;
    200:   dhy   = 1. / hy;
    201:   hxdhy = hx / hy;
    202:   hydhx = hy / hx;
    203:   scale = hx * hy;

    205:   for (j = info->ys; j < info->ys + info->ym; j++) {
    206:     for (i = info->xs; i < info->xs + info->xm; i++) {
    207:       vx  = ctx->gradq_inflow * ctx->porosity * ctx->saturation;
    208:       vy  = 0.0 * dhy;
    209:       avx = PetscAbsScalar(vx);
    210:       vxp = .5 * (vx + avx);
    211:       vxm = .5 * (vx - avx);
    212:       avy = PetscAbsScalar(vy);
    213:       vyp = .5 * (vy + avy);
    214:       vym = .5 * (vy - avy);
    215:       /* chemical reactions */
    216:       for (l = 0; l < N_SPECIES; l++) {
    217:         /* determine the velocites as the gradients of the pressure */
    218:         if (i == 0) {
    219:           sxp = (x[j][i + 1].sp[l] - x[j][i].sp[l]) * dhx;
    220:           sxm = sxp;
    221:         } else if (i == info->mx - 1) {
    222:           sxp = (x[j][i].sp[l] - x[j][i - 1].sp[l]) * dhx;
    223:           sxm = sxp;
    224:         } else {
    225:           sxm = (x[j][i + 1].sp[l] - x[j][i].sp[l]) * dhx;
    226:           sxp = (x[j][i].sp[l] - x[j][i - 1].sp[l]) * dhx;
    227:         }
    228:         if (j == 0) {
    229:           syp = (x[j + 1][i].sp[l] - x[j][i].sp[l]) * dhy;
    230:           sym = syp;
    231:         } else if (j == info->my - 1) {
    232:           syp = (x[j][i].sp[l] - x[j - 1][i].sp[l]) * dhy;
    233:           sym = syp;
    234:         } else {
    235:           sym = (x[j + 1][i].sp[l] - x[j][i].sp[l]) * dhy;
    236:           syp = (x[j][i].sp[l] - x[j - 1][i].sp[l]) * dhy;
    237:         } /* 4 flops per species*point */

    239:         if (i == 0) {
    240:           if (l == 0) f[j][i].sp[l] = (x[j][i].sp[l] - ctx->x_inflow.sp[l] * ((PetscScalar)j) / (info->my - 1));
    241:           else if (l == 1) f[j][i].sp[l] = (x[j][i].sp[l] - ctx->x_inflow.sp[l] * (1. - ((PetscScalar)j) / (info->my - 1)));
    242:           else f[j][i].sp[l] = x[j][i].sp[l];

    244:         } else {
    245:           f[j][i].sp[l] = xt[j][i].sp[l] * scale;
    246:           u             = x[j][i].sp[l];
    247:           if (j == 0) uyy = u - x[j + 1][i].sp[l];
    248:           else if (j == info->my - 1) uyy = u - x[j - 1][i].sp[l];
    249:           else uyy = (2.0 * u - x[j - 1][i].sp[l] - x[j + 1][i].sp[l]) * hxdhy;

    251:           if (i != info->mx - 1) uxx = (2.0 * u - x[j][i - 1].sp[l] - x[j][i + 1].sp[l]) * hydhx;
    252:           else uxx = u - x[j][i - 1].sp[l];
    253:           /* 10 flops per species*point */

    255:           f_advect = 0.;
    256:           f_advect += scale * (vxp * sxp + vxm * sxm);
    257:           f_advect += scale * (vyp * syp + vym * sym);
    258:           f[j][i].sp[l] += f_advect + ctx->dispersivity * (uxx + uyy);
    259:           /* 14 flops per species*point */
    260:         }
    261:       }
    262:       /* reaction */
    263:       if (i != 0) {
    264:         for (m = 0; m < N_REACTIONS; m++) {
    265:           rate = ctx->rate_constant[m];
    266:           for (l = 0; l < N_SPECIES; l++) {
    267:             if (stoich(m, l) < 0) {
    268:               /* assume an elementary reaction */
    269:               rate *= PetscPowScalar(x[j][i].sp[l], PetscAbsScalar(stoich(m, l)));
    270:               /* ~10 flops per reaction per species per point */
    271:             }
    272:           }
    273:           for (l = 0; l < N_SPECIES; l++) {
    274:             f[j][i].sp[l] += -scale * stoich(m, l) * rate; /* Reaction term */
    275:                                                            /* 3 flops per reaction per species per point */
    276:           }
    277:         }
    278:       }
    279:     }
    280:   }
    281:   PetscCall(PetscLogFlops((N_SPECIES * (28.0 + 13.0 * N_REACTIONS)) * info->ym * info->xm));
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }

    285: PetscErrorCode ReactingFlowPostCheck(SNESLineSearch linesearch, Vec X, Vec Y, Vec W, PetscBool *changed_y, PetscBool *changed_w, void *vctx)
    286: {
    287:   PetscInt    i, j, l, Mx, My, xs, ys, xm, ym;
    288:   Field     **x;
    289:   SNES        snes;
    290:   DM          da;
    291:   PetscScalar min;

    293:   PetscFunctionBeginUser;
    294:   *changed_w = PETSC_FALSE;
    295:   PetscCall(VecMin(X, NULL, &min));
    296:   if (min >= 0.) PetscFunctionReturn(PETSC_SUCCESS);

    298:   *changed_w = PETSC_TRUE;
    299:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
    300:   PetscCall(SNESGetDM(snes, &da));
    301:   PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
    302:   PetscCall(DMDAVecGetArray(da, W, &x));
    303:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
    304:   for (j = ys; j < ys + ym; j++) {
    305:     for (i = xs; i < xs + xm; i++) {
    306:       for (l = 0; l < N_SPECIES; l++) {
    307:         if (x[j][i].sp[l] < 0.) x[j][i].sp[l] = 0.;
    308:       }
    309:     }
    310:   }
    311:   PetscCall(DMDAVecRestoreArray(da, W, &x));
    312:   PetscFunctionReturn(PETSC_SUCCESS);
    313: }

    315: PetscErrorCode FormIFunction(TS ts, PetscReal ptime, Vec X, Vec Xdot, Vec F, void *user)
    316: {
    317:   DMDALocalInfo info;
    318:   Field       **u, **udot, **fu;
    319:   Vec           localX;
    320:   DM            da;

    322:   PetscFunctionBeginUser;
    323:   PetscCall(TSGetDM(ts, (DM *)&da));
    324:   PetscCall(DMGetLocalVector(da, &localX));
    325:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
    326:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
    327:   PetscCall(DMDAGetLocalInfo(da, &info));
    328:   PetscCall(DMDAVecGetArrayRead(da, localX, &u));
    329:   PetscCall(DMDAVecGetArrayRead(da, Xdot, &udot));
    330:   PetscCall(DMDAVecGetArray(da, F, &fu));
    331:   PetscCall(FormIFunctionLocal(&info, ptime, u, udot, fu, (AppCtx *)user));
    332:   PetscCall(DMDAVecRestoreArrayRead(da, localX, &u));
    333:   PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &udot));
    334:   PetscCall(DMDAVecRestoreArray(da, F, &fu));
    335:   PetscCall(DMRestoreLocalVector(da, &localX));
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }