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_time_step 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: const PetscReal *bmax;
210: PetscCall(TSMonitorEnvelopeGetBounds(ts, &max, NULL));
211: if (max) {
212: PetscCall(TSMonitorLGGetVariableNames(ts, &names));
213: if (names) {
214: PetscCall(VecGetArrayRead(max, &bmax));
215: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species - maximum mass fraction\n"));
216: for (PetscInt i = 1; i < user.Nspec; i++) {
217: if (bmax[i] > .01) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s %g\n", names[i], (double)bmax[i]));
218: }
219: PetscCall(VecRestoreArrayRead(max, &bmax));
220: }
221: }
222: }
224: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225: Free work space.
226: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227: TC_reset();
228: PetscCall(DMDestroy(&user.dm));
229: PetscCall(MatDestroy(&J));
230: PetscCall(VecDestroy(&X));
231: PetscCall(VecDestroy(&lambda));
232: PetscCall(TSDestroy(&ts));
233: PetscCall(PetscFree3(user.tchemwork, user.Jdense, user.rows));
234: PetscCall(PetscFinalize());
235: return 0;
236: }
238: /*
239: Applies the second order centered difference diffusion operator on a one dimensional periodic domain
240: */
241: static PetscErrorCode FormDiffusionFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
242: {
243: User user = (User)ptr;
244: PetscScalar **f;
245: const PetscScalar **x;
246: DM dm;
247: PetscInt i, xs, xm, j, dof;
248: Vec Xlocal;
249: PetscReal idx;
251: PetscFunctionBeginUser;
252: PetscCall(TSGetDM(ts, &dm));
253: PetscCall(DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
254: PetscCall(DMGetLocalVector(dm, &Xlocal));
255: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xlocal));
256: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xlocal));
257: PetscCall(DMDAVecGetArrayDOFRead(dm, Xlocal, &x));
258: PetscCall(DMDAVecGetArrayDOF(dm, F, &f));
259: PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));
261: idx = 1.0 * user->diffus / user->dx;
262: for (i = xs; i < xs + xm; i++) {
263: for (j = 0; j < dof; j++) f[i][j] += idx * (x[i + 1][j] - 2.0 * x[i][j] + x[i - 1][j]);
264: }
265: PetscCall(DMDAVecRestoreArrayDOFRead(dm, Xlocal, &x));
266: PetscCall(DMDAVecRestoreArrayDOF(dm, F, &f));
267: PetscCall(DMRestoreLocalVector(dm, &Xlocal));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /*
272: Produces the second order centered difference diffusion operator on a one dimensional periodic domain
273: */
274: static PetscErrorCode FormDiffusionJacobian(TS ts, PetscReal t, Vec X, Mat Amat, Mat Pmat, void *ptr)
275: {
276: User user = (User)ptr;
277: DM dm;
278: PetscInt i, xs, xm, j, dof;
279: PetscReal idx, values[3];
280: MatStencil row, col[3];
282: PetscFunctionBeginUser;
283: PetscCall(TSGetDM(ts, &dm));
284: PetscCall(DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
285: PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));
287: idx = 1.0 * user->diffus / user->dx;
288: values[0] = idx;
289: values[1] = -2.0 * idx;
290: values[2] = idx;
291: for (i = xs; i < xs + xm; i++) {
292: for (j = 0; j < dof; j++) {
293: row.i = i;
294: row.c = j;
295: col[0].i = i - 1;
296: col[0].c = j;
297: col[1].i = i;
298: col[1].c = j;
299: col[2].i = i + 1;
300: col[2].c = j;
301: PetscCall(MatSetValuesStencil(Pmat, 1, &row, 3, col, values, ADD_VALUES));
302: }
303: }
304: PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
305: PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
310: {
311: User user = (User)ptr;
312: PetscScalar **f;
313: const PetscScalar **x;
314: DM dm;
315: PetscInt i, xs, xm;
317: PetscFunctionBeginUser;
318: if (user->reactions) {
319: PetscCall(TSGetDM(ts, &dm));
320: PetscCall(DMDAVecGetArrayDOFRead(dm, X, &x));
321: PetscCall(DMDAVecGetArrayDOF(dm, F, &f));
322: PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));
324: for (i = xs; i < xs + xm; i++) {
325: PetscCall(PetscArraycpy(user->tchemwork, x[i], user->Nspec + 1));
326: user->tchemwork[0] *= user->Tini; /* Dimensionalize */
327: PetscCallTC(TC_getSrc(user->tchemwork, user->Nspec + 1, f[i]));
328: f[i][0] /= user->Tini; /* Non-dimensionalize */
329: }
331: PetscCall(DMDAVecRestoreArrayDOFRead(dm, X, &x));
332: PetscCall(DMDAVecRestoreArrayDOF(dm, F, &f));
333: } else {
334: PetscCall(VecZeroEntries(F));
335: }
336: if (user->diffusion) PetscCall(FormDiffusionFunction(ts, t, X, F, ptr));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode FormRHSJacobian(TS ts, PetscReal t, Vec X, Mat Amat, Mat Pmat, void *ptr)
341: {
342: User user = (User)ptr;
343: const PetscScalar **x;
344: PetscInt M = user->Nspec + 1, i, j, xs, xm;
345: DM dm;
347: PetscFunctionBeginUser;
348: if (user->reactions) {
349: PetscCall(TSGetDM(ts, &dm));
350: PetscCall(MatZeroEntries(Pmat));
351: PetscCall(MatSetOption(Pmat, MAT_ROW_ORIENTED, PETSC_FALSE));
352: PetscCall(MatSetOption(Pmat, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
353: PetscCall(DMDAVecGetArrayDOFRead(dm, X, &x));
354: PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));
356: for (i = xs; i < xs + xm; i++) {
357: PetscCall(PetscArraycpy(user->tchemwork, x[i], user->Nspec + 1));
358: user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */
359: PetscCall(TC_getJacTYN(user->tchemwork, user->Nspec, user->Jdense, 1));
361: for (j = 0; j < M; j++) user->Jdense[j + 0 * M] /= user->Tini; /* Non-dimensionalize first column */
362: for (j = 0; j < M; j++) user->Jdense[0 + j * M] /= user->Tini; /* Non-dimensionalize first row */
363: for (j = 0; j < M; j++) user->rows[j] = i * M + j;
364: PetscCall(MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES));
365: }
366: PetscCall(DMDAVecRestoreArrayDOFRead(dm, X, &x));
367: PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
368: PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
369: } else {
370: PetscCall(MatZeroEntries(Pmat));
371: }
372: if (user->diffusion) PetscCall(FormDiffusionJacobian(ts, t, X, Amat, Pmat, ptr));
373: if (Amat != Pmat) {
374: PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
375: PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
376: }
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx ctx)
381: {
382: PetscScalar **x, *xc;
383: struct {
384: const char *name;
385: PetscReal massfrac;
386: } initial[] = {
387: {"CH4", 0.0948178320887 },
388: {"O2", 0.189635664177 },
389: {"N2", 0.706766236705 },
390: {"AR", 0.00878026702874}
391: };
392: PetscInt i, j, xs, xm;
393: DM dm;
395: PetscFunctionBeginUser;
396: PetscCall(VecZeroEntries(X));
397: PetscCall(TSGetDM(ts, &dm));
398: PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));
400: PetscCall(DMDAGetCoordinateArray(dm, &xc));
401: PetscCall(DMDAVecGetArrayDOF(dm, X, &x));
402: for (i = xs; i < xs + xm; i++) {
403: x[i][0] = 1.0 + .05 * PetscSinScalar(2. * PETSC_PI * xc[i]); /* Non-dimensionalized by user->Tini */
404: for (j = 0; j < PETSC_STATIC_ARRAY_LENGTH(initial); j++) {
405: int ispec = TC_getSpos(initial[j].name, (int)strlen(initial[j].name));
406: PetscCheck(ispec >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Could not find species %s", initial[j].name);
407: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species %d: %s %g\n", j, initial[j].name, (double)initial[j].massfrac));
408: x[i][1 + ispec] = initial[j].massfrac;
409: }
410: }
411: PetscCall(DMDAVecRestoreArrayDOF(dm, X, &x));
412: PetscCall(DMDARestoreCoordinateArray(dm, &xc));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: /*
417: Routines for displaying the solutions
418: */
419: typedef struct {
420: PetscInt cell;
421: User user;
422: } UserLGCtx;
424: static PetscErrorCode FormMoleFraction(UserLGCtx *ctx, Vec massf, Vec *molef)
425: {
426: User user = ctx->user;
427: PetscReal *M, tM = 0;
428: PetscInt i, n = user->Nspec + 1;
429: PetscScalar *mof;
430: const PetscScalar **maf;
432: PetscFunctionBeginUser;
433: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, molef));
434: PetscCall(PetscMalloc1(user->Nspec, &M));
435: TC_getSmass(user->Nspec, M);
436: PetscCall(DMDAVecGetArrayDOFRead(user->dm, massf, &maf));
437: PetscCall(VecGetArray(*molef, &mof));
438: mof[0] = maf[ctx->cell][0]; /* copy over temperature */
439: for (i = 1; i < n; i++) tM += maf[ctx->cell][i] / M[i - 1];
440: for (i = 1; i < n; i++) mof[i] = maf[ctx->cell][i] / (M[i - 1] * tM);
441: PetscCall(DMDAVecRestoreArrayDOFRead(user->dm, massf, &maf));
442: PetscCall(VecRestoreArray(*molef, &mof));
443: PetscCall(PetscFree(M));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode MonitorCellDestroy(UserLGCtx **uctx)
448: {
449: PetscFunctionBeginUser;
450: PetscCall(PetscFree(*uctx));
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: /*
455: Use TSMonitorLG to monitor the reactions in a particular cell
456: */
457: static PetscErrorCode MonitorCell(TS ts, User user, PetscInt cell)
458: {
459: TSMonitorLGCtx ctx;
460: char **snames;
461: UserLGCtx *uctx;
462: char label[128];
463: PetscReal temp, *xc;
464: PetscMPIInt rank;
466: PetscFunctionBeginUser;
467: PetscCall(DMDAGetCoordinateArray(user->dm, &xc));
468: temp = 1.0 + .05 * PetscSinScalar(2. * PETSC_PI * xc[cell]); /* Non-dimensionalized by user->Tini */
469: PetscCall(DMDARestoreCoordinateArray(user->dm, &xc));
470: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
471: PetscCall(PetscSNPrintf(label, sizeof(label), "Initial Temperature %g Cell %d Rank %d", (double)user->Tini * temp, (int)cell, rank));
472: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, label, PETSC_DECIDE, PETSC_DECIDE, 600, 400, 1, &ctx));
473: PetscCall(DMDAGetFieldNames(user->dm, (const char *const **)&snames));
474: PetscCall(TSMonitorLGCtxSetVariableNames(ctx, (const char *const *)snames));
475: PetscCall(PetscNew(&uctx));
476: uctx->cell = cell;
477: uctx->user = user;
478: PetscCall(TSMonitorLGCtxSetTransform(ctx, (PetscErrorCode (*)(void *, Vec, Vec *))FormMoleFraction, (PetscCtxDestroyFn *)MonitorCellDestroy, uctx));
479: PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }