Actual source code: extchemfield.c

      1: static const char help[] = "Integrate chemistry using TChem.\n";

      3: #include <petscts.h>
      4: #include <petscdmda.h>

      6: #if defined(PETSC_HAVE_TCHEM)
      7:   #if defined(MAX)
      8:     #undef MAX
      9:   #endif
     10:   #if defined(MIN)
     11:     #undef MIN
     12:   #endif
     13:   #include <TC_params.h>
     14:   #include <TC_interface.h>
     15: #else
     16:   #error TChem is required for this example.  Reconfigure PETSc using --download-tchem.
     17: #endif
     18: /*

     20:     This is an extension of extchem.c to solve the reaction equations independently in each cell of a one dimensional field

     22:     Obtain the three files into this directory

     24:        curl http://combustion.berkeley.edu/gri_mech/version30/files30/grimech30.dat > chem.inp
     25:        curl http://combustion.berkeley.edu/gri_mech/version30/files30/thermo30.dat > therm.dat
     26:        cp $PETSC_DIR/$PETSC_ARCH/externalpackages/tchem/data/periodictable.dat .

     28:     Run with
     29:      ./extchemfield  -ts_arkimex_fully_implicit -ts_max_snes_failures unlimited -ts_adapt_monitor -ts_adapt_dt_max 1e-4 -ts_arkimex_type 4 -ts_max_time .005

     31:      Options for visualizing the solution:
     32:         Watch certain variables in each cell evolve with time
     33:         -draw_solution 1 -ts_monitor_lg_solution_variables Temp,H2,O2,H2O,CH4,CO,CO2,C2H2,N2 -lg_use_markers false  -draw_pause -2

     35:         Watch certain variables in all cells evolve with time
     36:         -da_refine 4 -ts_monitor_draw_solution -draw_fields_by_name Temp,H2 -draw_vec_mark_points  -draw_pause -2

     38:         Keep the initial temperature distribution as one monitors the current temperature distribution
     39:         -ts_monitor_draw_solution_initial -draw_bounds .9,1.7 -draw_fields_by_name Temp

     41:         Save the images in a .gif (movie) file
     42:         -draw_save -draw_save_single_file

     44:         Compute the sensitivities of the solution of the first temperature on the initial conditions
     45:         -ts_adjoint_solve  -ts_dt 1.e-5 -ts_type cn -ts_adjoint_view_solution draw

     47:         Turn off diffusion
     48:         -diffusion no

     50:         Turn off reactions
     51:         -reactions no

     53:     The solution for component i = 0 is the temperature.

     55:     The solution, i > 0, is the mass fraction, massf[i], of species i, i.e. mass of species i/ total mass of all species

     57:     The mole fraction molef[i], i > 0, is the number of moles of a species/ total number of moles of all species
     58:         Define M[i] = mass per mole of species i then
     59:         molef[i] = massf[i]/(M[i]*(sum_j massf[j]/M[j]))

     61:     FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species.

     63: */
     64: typedef struct _User *User;
     65: struct _User {
     66:   PetscReal pressure;
     67:   int       Nspec;
     68:   int       Nreac;
     69:   PetscReal Tini, dx;
     70:   PetscReal diffus;
     71:   DM        dm;
     72:   PetscBool diffusion, reactions;
     73:   double   *tchemwork;
     74:   double   *Jdense; /* Dense array workspace where Tchem computes the Jacobian */
     75:   PetscInt *rows;
     76: };

     78: static PetscErrorCode MonitorCell(TS, User, PetscInt);
     79: static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
     80: static PetscErrorCode FormRHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
     81: static PetscErrorCode FormInitialSolution(TS, Vec, void *);

     83: #define PetscCallTC(ierr) \
     84:   do { \
     85:     PetscCheck(!ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in TChem library, return code %d", ierr); \
     86:   } while (0)

     88: int main(int argc, char **argv)
     89: {
     90:   TS                ts; /* time integrator */
     91:   TSAdapt           adapt;
     92:   Vec               X; /* solution vector */
     93:   Mat               J; /* Jacobian matrix */
     94:   PetscInt          steps, ncells, xs, xm, i;
     95:   PetscReal         ftime, dt;
     96:   char              chemfile[PETSC_MAX_PATH_LEN] = "chem.inp", thermofile[PETSC_MAX_PATH_LEN] = "therm.dat";
     97:   struct _User      user;
     98:   TSConvergedReason reason;
     99:   PetscBool         showsolutions = PETSC_FALSE;
    100:   char            **snames, *names;
    101:   Vec               lambda; /* used with TSAdjoint for sensitivities */

    103:   PetscFunctionBeginUser;
    104:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    105:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Chemistry solver options", "");
    106:   PetscCall(PetscOptionsString("-chem", "CHEMKIN input file", "", chemfile, chemfile, sizeof(chemfile), NULL));
    107:   PetscCall(PetscOptionsString("-thermo", "NASA thermo input file", "", thermofile, thermofile, sizeof(thermofile), NULL));
    108:   user.pressure = 1.01325e5; /* Pascal */
    109:   PetscCall(PetscOptionsReal("-pressure", "Pressure of reaction [Pa]", "", user.pressure, &user.pressure, NULL));
    110:   user.Tini = 1550;
    111:   PetscCall(PetscOptionsReal("-Tini", "Initial temperature [K]", "", user.Tini, &user.Tini, NULL));
    112:   user.diffus = 100;
    113:   PetscCall(PetscOptionsReal("-diffus", "Diffusion constant", "", user.diffus, &user.diffus, NULL));
    114:   PetscCall(PetscOptionsBool("-draw_solution", "Plot the solution for each cell", "", showsolutions, &showsolutions, NULL));
    115:   user.diffusion = PETSC_TRUE;
    116:   PetscCall(PetscOptionsBool("-diffusion", "Have diffusion", "", user.diffusion, &user.diffusion, NULL));
    117:   user.reactions = PETSC_TRUE;
    118:   PetscCall(PetscOptionsBool("-reactions", "Have reactions", "", user.reactions, &user.reactions, NULL));
    119:   PetscOptionsEnd();

    121:   PetscCallTC(TC_initChem(chemfile, thermofile, 0, 1.0));
    122:   user.Nspec = TC_getNspec();
    123:   user.Nreac = TC_getNreac();

    125:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, user.Nspec + 1, 1, NULL, &user.dm));
    126:   PetscCall(DMSetFromOptions(user.dm));
    127:   PetscCall(DMSetUp(user.dm));
    128:   PetscCall(DMDAGetInfo(user.dm, NULL, &ncells, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    129:   user.dx = 1.0 / ncells; /* Set the coordinates of the cell centers; note final ghost cell is at x coordinate 1.0 */
    130:   PetscCall(DMDASetUniformCoordinates(user.dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));

    132:   /* set the names of each field in the DMDA based on the species name */
    133:   PetscCall(PetscMalloc1((user.Nspec + 1) * LENGTHOFSPECNAME, &names));
    134:   PetscCall(PetscStrncpy(names, "Temp", (user.Nspec + 1) * LENGTHOFSPECNAME);
    135:   TC_getSnames(user.Nspec, names + LENGTHOFSPECNAME);
    136:   PetscCall(PetscMalloc1( user.Nspec + 2 , &snames));
    137:   for (i = 0; i < user.Nspec + 1; i++) snames[i] = names + i * LENGTHOFSPECNAME;
    138:   snames[user.Nspec + 1] = NULL;
    139:   PetscCall(DMDASetFieldNames(user.dm, (const char *const *)snames));
    140:   PetscCall(PetscFree(snames));
    141:   PetscCall(PetscFree(names));

    143:   PetscCall(DMCreateMatrix(user.dm, &J));
    144:   PetscCall(DMCreateGlobalVector(user.dm, &X));

    146:   PetscCall(PetscMalloc3(user.Nspec + 1, &user.tchemwork, PetscSqr(user.Nspec + 1), &user.Jdense, user.Nspec + 1, &user.rows));

    148:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    149:      Create timestepping solver context
    150:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    151:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    152:   PetscCall(TSSetDM(ts, user.dm));
    153:   PetscCall(TSSetType(ts, TSARKIMEX));
    154:   PetscCall(TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE));
    155:   PetscCall(TSARKIMEXSetType(ts, TSARKIMEX4));
    156:   PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
    157:   PetscCall(TSSetRHSJacobian(ts, J, J, FormRHSJacobian, &user));

    159:   ftime = 1.0;
    160:   PetscCall(TSSetMaxTime(ts, ftime));
    161:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    164:      Set initial conditions
    165:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    166:   PetscCall(FormInitialSolution(ts, X, &user));
    167:   PetscCall(TSSetSolution(ts, X));
    168:   dt = 1e-10; /* Initial time step */
    169:   PetscCall(TSSetTimeStep(ts, dt));
    170:   PetscCall(TSGetAdapt(ts, &adapt));
    171:   PetscCall(TSAdaptSetStepLimits(adapt, 1e-12, 1e-4)); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
    172:   PetscCall(TSSetMaxSNESFailures(ts, -1));             /* Retry step an unlimited number of times */

    174:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    175:      Pass information to graphical monitoring routine
    176:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    177:   if (showsolutions) {
    178:     PetscCall(DMDAGetCorners(user.dm, &xs, NULL, NULL, &xm, NULL, NULL));
    179:     for (i = xs; i < xs + xm; i++) PetscCall(MonitorCell(ts, &user, i));
    180:   }

    182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    183:      Set runtime options
    184:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    185:   PetscCall(TSSetFromOptions(ts));

    187:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    188:      Set final conditions for sensitivities
    189:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    190:   PetscCall(DMCreateGlobalVector(user.dm, &lambda));
    191:   PetscCall(TSSetCostGradients(ts, 1, &lambda, NULL));
    192:   PetscCall(VecSetValue(lambda, 0, 1.0, INSERT_VALUES));
    193:   PetscCall(VecAssemblyBegin(lambda));
    194:   PetscCall(VecAssemblyEnd(lambda));

    196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    197:      Solve ODE
    198:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    199:   PetscCall(TSSolve(ts, X));
    200:   PetscCall(TSGetSolveTime(ts, &ftime));
    201:   PetscCall(TSGetStepNumber(ts, &steps));
    202:   PetscCall(TSGetConvergedReason(ts, &reason));
    203:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));

    205:   {
    206:     Vec                max;
    207:     const char *const *names;
    208:     PetscInt           i;
    209:     const PetscReal   *bmax;

    211:     PetscCall(TSMonitorEnvelopeGetBounds(ts, &max, NULL));
    212:     if (max) {
    213:       PetscCall(TSMonitorLGGetVariableNames(ts, &names));
    214:       if (names) {
    215:         PetscCall(VecGetArrayRead(max, &bmax));
    216:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species - maximum mass fraction\n"));
    217:         for (i = 1; i < user.Nspec; i++) {
    218:           if (bmax[i] > .01) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s %g\n", names[i], (double)bmax[i]));
    219:         }
    220:         PetscCall(VecRestoreArrayRead(max, &bmax));
    221:       }
    222:     }
    223:   }

    225:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    226:      Free work space.
    227:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    228:   TC_reset();
    229:   PetscCall(DMDestroy(&user.dm));
    230:   PetscCall(MatDestroy(&J));
    231:   PetscCall(VecDestroy(&X));
    232:   PetscCall(VecDestroy(&lambda));
    233:   PetscCall(TSDestroy(&ts));
    234:   PetscCall(PetscFree3(user.tchemwork, user.Jdense, user.rows));
    235:   PetscCall(PetscFinalize());
    236:   return 0;
    237: }

    239: /*
    240:    Applies the second order centered difference diffusion operator on a one dimensional periodic domain
    241: */
    242: static PetscErrorCode FormDiffusionFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
    243: {
    244:   User                user = (User)ptr;
    245:   PetscScalar       **f;
    246:   const PetscScalar **x;
    247:   DM                  dm;
    248:   PetscInt            i, xs, xm, j, dof;
    249:   Vec                 Xlocal;
    250:   PetscReal           idx;

    252:   PetscFunctionBeginUser;
    253:   PetscCall(TSGetDM(ts, &dm));
    254:   PetscCall(DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
    255:   PetscCall(DMGetLocalVector(dm, &Xlocal));
    256:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xlocal));
    257:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xlocal));
    258:   PetscCall(DMDAVecGetArrayDOFRead(dm, Xlocal, &x));
    259:   PetscCall(DMDAVecGetArrayDOF(dm, F, &f));
    260:   PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

    262:   idx = 1.0 * user->diffus / user->dx;
    263:   for (i = xs; i < xs + xm; i++) {
    264:     for (j = 0; j < dof; j++) f[i][j] += idx * (x[i + 1][j] - 2.0 * x[i][j] + x[i - 1][j]);
    265:   }
    266:   PetscCall(DMDAVecRestoreArrayDOFRead(dm, Xlocal, &x));
    267:   PetscCall(DMDAVecRestoreArrayDOF(dm, F, &f));
    268:   PetscCall(DMRestoreLocalVector(dm, &Xlocal));
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: }

    272: /*
    273:    Produces the second order centered difference diffusion operator on a one dimensional periodic domain
    274: */
    275: static PetscErrorCode FormDiffusionJacobian(TS ts, PetscReal t, Vec X, Mat Amat, Mat Pmat, void *ptr)
    276: {
    277:   User       user = (User)ptr;
    278:   DM         dm;
    279:   PetscInt   i, xs, xm, j, dof;
    280:   PetscReal  idx, values[3];
    281:   MatStencil row, col[3];

    283:   PetscFunctionBeginUser;
    284:   PetscCall(TSGetDM(ts, &dm));
    285:   PetscCall(DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
    286:   PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

    288:   idx       = 1.0 * user->diffus / user->dx;
    289:   values[0] = idx;
    290:   values[1] = -2.0 * idx;
    291:   values[2] = idx;
    292:   for (i = xs; i < xs + xm; i++) {
    293:     for (j = 0; j < dof; j++) {
    294:       row.i    = i;
    295:       row.c    = j;
    296:       col[0].i = i - 1;
    297:       col[0].c = j;
    298:       col[1].i = i;
    299:       col[1].c = j;
    300:       col[2].i = i + 1;
    301:       col[2].c = j;
    302:       PetscCall(MatSetValuesStencil(Pmat, 1, &row, 3, col, values, ADD_VALUES));
    303:     }
    304:   }
    305:   PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
    306:   PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }

    310: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
    311: {
    312:   User                user = (User)ptr;
    313:   PetscScalar       **f;
    314:   const PetscScalar **x;
    315:   DM                  dm;
    316:   PetscInt            i, xs, xm;

    318:   PetscFunctionBeginUser;
    319:   if (user->reactions) {
    320:     PetscCall(TSGetDM(ts, &dm));
    321:     PetscCall(DMDAVecGetArrayDOFRead(dm, X, &x));
    322:     PetscCall(DMDAVecGetArrayDOF(dm, F, &f));
    323:     PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

    325:     for (i = xs; i < xs + xm; i++) {
    326:       PetscCall(PetscArraycpy(user->tchemwork, x[i], user->Nspec + 1));
    327:       user->tchemwork[0] *= user->Tini; /* Dimensionalize */
    328:       PetscCallTC(TC_getSrc(user->tchemwork, user->Nspec + 1, f[i]));
    329:       f[i][0] /= user->Tini; /* Non-dimensionalize */
    330:     }

    332:     PetscCall(DMDAVecRestoreArrayDOFRead(dm, X, &x));
    333:     PetscCall(DMDAVecRestoreArrayDOF(dm, F, &f));
    334:   } else {
    335:     PetscCall(VecZeroEntries(F));
    336:   }
    337:   if (user->diffusion) PetscCall(FormDiffusionFunction(ts, t, X, F, ptr));
    338:   PetscFunctionReturn(PETSC_SUCCESS);
    339: }

    341: static PetscErrorCode FormRHSJacobian(TS ts, PetscReal t, Vec X, Mat Amat, Mat Pmat, void *ptr)
    342: {
    343:   User                user = (User)ptr;
    344:   const PetscScalar **x;
    345:   PetscInt            M = user->Nspec + 1, i, j, xs, xm;
    346:   DM                  dm;

    348:   PetscFunctionBeginUser;
    349:   if (user->reactions) {
    350:     PetscCall(TSGetDM(ts, &dm));
    351:     PetscCall(MatZeroEntries(Pmat));
    352:     PetscCall(MatSetOption(Pmat, MAT_ROW_ORIENTED, PETSC_FALSE));
    353:     PetscCall(MatSetOption(Pmat, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
    354:     PetscCall(DMDAVecGetArrayDOFRead(dm, X, &x));
    355:     PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

    357:     for (i = xs; i < xs + xm; i++) {
    358:       PetscCall(PetscArraycpy(user->tchemwork, x[i], user->Nspec + 1));
    359:       user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */
    360:       PetscCall(TC_getJacTYN(user->tchemwork, user->Nspec, user->Jdense, 1));

    362:       for (j = 0; j < M; j++) user->Jdense[j + 0 * M] /= user->Tini; /* Non-dimensionalize first column */
    363:       for (j = 0; j < M; j++) user->Jdense[0 + j * M] /= user->Tini; /* Non-dimensionalize first row */
    364:       for (j = 0; j < M; j++) user->rows[j] = i * M + j;
    365:       PetscCall(MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES));
    366:     }
    367:     PetscCall(DMDAVecRestoreArrayDOFRead(dm, X, &x));
    368:     PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
    369:     PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
    370:   } else {
    371:     PetscCall(MatZeroEntries(Pmat));
    372:   }
    373:   if (user->diffusion) PetscCall(FormDiffusionJacobian(ts, t, X, Amat, Pmat, ptr));
    374:   if (Amat != Pmat) {
    375:     PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
    376:     PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
    377:   }
    378:   PetscFunctionReturn(PETSC_SUCCESS);
    379: }

    381: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
    382: {
    383:   PetscScalar **x, *xc;
    384:   struct {
    385:     const char *name;
    386:     PetscReal   massfrac;
    387:   } initial[] = {
    388:     {"CH4", 0.0948178320887 },
    389:     {"O2",  0.189635664177  },
    390:     {"N2",  0.706766236705  },
    391:     {"AR",  0.00878026702874}
    392:   };
    393:   PetscInt i, j, xs, xm;
    394:   DM       dm;

    396:   PetscFunctionBeginUser;
    397:   PetscCall(VecZeroEntries(X));
    398:   PetscCall(TSGetDM(ts, &dm));
    399:   PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

    401:   PetscCall(DMDAGetCoordinateArray(dm, &xc));
    402:   PetscCall(DMDAVecGetArrayDOF(dm, X, &x));
    403:   for (i = xs; i < xs + xm; i++) {
    404:     x[i][0] = 1.0 + .05 * PetscSinScalar(2. * PETSC_PI * xc[i]); /* Non-dimensionalized by user->Tini */
    405:     for (j = 0; j < PETSC_STATIC_ARRAY_LENGTH(initial); j++) {
    406:       int ispec = TC_getSpos(initial[j].name, strlen(initial[j].name));
    407:       PetscCheck(ispec >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Could not find species %s", initial[j].name);
    408:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species %d: %s %g\n", j, initial[j].name, (double)initial[j].massfrac));
    409:       x[i][1 + ispec] = initial[j].massfrac;
    410:     }
    411:   }
    412:   PetscCall(DMDAVecRestoreArrayDOF(dm, X, &x));
    413:   PetscCall(DMDARestoreCoordinateArray(dm, &xc));
    414:   PetscFunctionReturn(PETSC_SUCCESS);
    415: }

    417: /*
    418:     Routines for displaying the solutions
    419: */
    420: typedef struct {
    421:   PetscInt cell;
    422:   User     user;
    423: } UserLGCtx;

    425: static PetscErrorCode FormMoleFraction(UserLGCtx *ctx, Vec massf, Vec *molef)
    426: {
    427:   User                user = ctx->user;
    428:   PetscReal          *M, tM = 0;
    429:   PetscInt            i, n  = user->Nspec + 1;
    430:   PetscScalar        *mof;
    431:   const PetscScalar **maf;

    433:   PetscFunctionBeginUser;
    434:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, molef));
    435:   PetscCall(PetscMalloc1(user->Nspec, &M));
    436:   TC_getSmass(user->Nspec, M);
    437:   PetscCall(DMDAVecGetArrayDOFRead(user->dm, massf, &maf));
    438:   PetscCall(VecGetArray(*molef, &mof));
    439:   mof[0] = maf[ctx->cell][0]; /* copy over temperature */
    440:   for (i = 1; i < n; i++) tM += maf[ctx->cell][i] / M[i - 1];
    441:   for (i = 1; i < n; i++) mof[i] = maf[ctx->cell][i] / (M[i - 1] * tM);
    442:   PetscCall(DMDAVecRestoreArrayDOFRead(user->dm, massf, &maf));
    443:   PetscCall(VecRestoreArray(*molef, &mof));
    444:   PetscCall(PetscFree(M));
    445:   PetscFunctionReturn(PETSC_SUCCESS);
    446: }

    448: static PetscErrorCode MonitorCellDestroy(UserLGCtx *uctx)
    449: {
    450:   PetscFunctionBeginUser;
    451:   PetscCall(PetscFree(uctx));
    452:   PetscFunctionReturn(PETSC_SUCCESS);
    453: }

    455: /*
    456:    Use TSMonitorLG to monitor the reactions in a particular cell
    457: */
    458: static PetscErrorCode MonitorCell(TS ts, User user, PetscInt cell)
    459: {
    460:   TSMonitorLGCtx ctx;
    461:   char         **snames;
    462:   UserLGCtx     *uctx;
    463:   char           label[128];
    464:   PetscReal      temp, *xc;
    465:   PetscMPIInt    rank;

    467:   PetscFunctionBeginUser;
    468:   PetscCall(DMDAGetCoordinateArray(user->dm, &xc));
    469:   temp = 1.0 + .05 * PetscSinScalar(2. * PETSC_PI * xc[cell]); /* Non-dimensionalized by user->Tini */
    470:   PetscCall(DMDARestoreCoordinateArray(user->dm, &xc));
    471:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    472:   PetscCall(PetscSNPrintf(label, sizeof(label), "Initial Temperature %g Cell %d Rank %d", (double)user->Tini * temp, (int)cell, rank));
    473:   PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, label, PETSC_DECIDE, PETSC_DECIDE, 600, 400, 1, &ctx));
    474:   PetscCall(DMDAGetFieldNames(user->dm, (const char *const **)&snames));
    475:   PetscCall(TSMonitorLGCtxSetVariableNames(ctx, (const char *const *)snames));
    476:   PetscCall(PetscNew(&uctx));
    477:   uctx->cell = cell;
    478:   uctx->user = user;
    479:   PetscCall(TSMonitorLGCtxSetTransform(ctx, (PetscErrorCode (*)(void *, Vec, Vec *))FormMoleFraction, (PetscErrorCode (*)(void *))MonitorCellDestroy, uctx));
    480:   PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
    481:   PetscFunctionReturn(PETSC_SUCCESS);
    482: }