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 -1 -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 sensitivies 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 { \
 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 */

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

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

125:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, user.Nspec + 1, 1, NULL, &user.dm);
126:   DMSetFromOptions(user.dm);
127:   DMSetUp(user.dm);
128:   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:   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:   PetscMalloc1((user.Nspec + 1) * LENGTHOFSPECNAME, &names);
134:   PetscStrcpy(names, "Temp");
135:   TC_getSnames(user.Nspec, names + LENGTHOFSPECNAME);
136:   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:   DMDASetFieldNames(user.dm, (const char *const *)snames);
140:   PetscFree(snames);
141:   PetscFree(names);

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

146:   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:   TSCreate(PETSC_COMM_WORLD, &ts);
152:   TSSetDM(ts, user.dm);
153:   TSSetType(ts, TSARKIMEX);
154:   TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE);
155:   TSARKIMEXSetType(ts, TSARKIMEX4);
156:   TSSetRHSFunction(ts, NULL, FormRHSFunction, &user);
157:   TSSetRHSJacobian(ts, J, J, FormRHSJacobian, &user);

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

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

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

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

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

196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197:      Solve ODE
198:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199:   TSSolve(ts, X);
200:   TSGetSolveTime(ts, &ftime);
201:   TSGetStepNumber(ts, &steps);
202:   TSGetConvergedReason(ts, &reason);
203:   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:     TSMonitorEnvelopeGetBounds(ts, &max, NULL);
212:     if (max) {
213:       TSMonitorLGGetVariableNames(ts, &names);
214:       if (names) {
215:         VecGetArrayRead(max, &bmax);
216:         PetscPrintf(PETSC_COMM_SELF, "Species - maximum mass fraction\n");
217:         for (i = 1; i < user.Nspec; i++) {
218:           if (bmax[i] > .01) PetscPrintf(PETSC_COMM_SELF, "%s %g\n", names[i], (double)bmax[i]);
219:         }
220:         VecRestoreArrayRead(max, &bmax);
221:       }
222:     }
223:   }

225:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
226:      Free work space.
227:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
228:   TC_reset();
229:   DMDestroy(&user.dm);
230:   MatDestroy(&J);
231:   VecDestroy(&X);
232:   VecDestroy(&lambda);
233:   TSDestroy(&ts);
234:   PetscFree3(user.tchemwork, user.Jdense, user.rows);
235:   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;

253:   TSGetDM(ts, &dm);
254:   DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
255:   DMGetLocalVector(dm, &Xlocal);
256:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xlocal);
257:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xlocal);
258:   DMDAVecGetArrayDOFRead(dm, Xlocal, &x);
259:   DMDAVecGetArrayDOF(dm, F, &f);
260:   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:   DMDAVecRestoreArrayDOFRead(dm, Xlocal, &x);
267:   DMDAVecRestoreArrayDOF(dm, F, &f);
268:   DMRestoreLocalVector(dm, &Xlocal);
269:   return 0;
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];

284:   TSGetDM(ts, &dm);
285:   DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
286:   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:       MatSetValuesStencil(Pmat, 1, &row, 3, col, values, ADD_VALUES);
303:     }
304:   }
305:   MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY);
306:   MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY);
307:   return 0;
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;

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

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

332:     DMDAVecRestoreArrayDOFRead(dm, X, &x);
333:     DMDAVecRestoreArrayDOF(dm, F, &f);
334:   } else {
335:     VecZeroEntries(F);
336:   }
337:   if (user->diffusion) FormDiffusionFunction(ts, t, X, F, ptr);
338:   return 0;
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;

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

357:     for (i = xs; i < xs + xm; i++) {
358:       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:       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:       MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES);
366:     }
367:     DMDAVecRestoreArrayDOFRead(dm, X, &x);
368:     MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY);
369:     MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY);
370:   } else {
371:     MatZeroEntries(Pmat);
372:   }
373:   if (user->diffusion) FormDiffusionJacobian(ts, t, X, Amat, Pmat, ptr);
374:   if (Amat != Pmat) {
375:     MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
376:     MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
377:   }
378:   return 0;
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;

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

401:   DMDAGetCoordinateArray(dm, &xc);
402:   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));
408:       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:   DMDAVecRestoreArrayDOF(dm, X, &x);
413:   DMDARestoreCoordinateArray(dm, &xc);
414:   return 0;
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;

434:   VecCreateSeq(PETSC_COMM_SELF, n, molef);
435:   PetscMalloc1(user->Nspec, &M);
436:   TC_getSmass(user->Nspec, M);
437:   DMDAVecGetArrayDOFRead(user->dm, massf, &maf);
438:   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:   DMDAVecRestoreArrayDOFRead(user->dm, massf, &maf);
443:   VecRestoreArray(*molef, &mof);
444:   PetscFree(M);
445:   return 0;
446: }

448: static PetscErrorCode MonitorCellDestroy(UserLGCtx *uctx)
449: {
451:   PetscFree(uctx);
452:   return 0;
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;

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