1 static const char help[] = "Integrate chemistry using TChem.\n";
2
3 #include <petscts.h>
4 #include <petscdmda.h>
5
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 /*
19
20 This is an extension of extchem.c to solve the reaction equations independently in each cell of a one dimensional field
21
22 Obtain the three files into this directory
23
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 .
27
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
30
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
34
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
37
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
40
41 Save the images in a .gif (movie) file
42 -draw_save -draw_save_single_file
43
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
46
47 Turn off diffusion
48 -diffusion no
49
50 Turn off reactions
51 -reactions no
52
53 The solution for component i = 0 is the temperature.
54
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
56
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]))
60
61 FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species.
62
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 };
77
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 *);
82
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)
87
main(int argc,char ** argv)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 */
102
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();
120
121 PetscCallTC(TC_initChem(chemfile, thermofile, 0, 1.0));
122 user.Nspec = TC_getNspec();
123 user.Nreac = TC_getNreac();
124
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));
131
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));
142
143 PetscCall(DMCreateMatrix(user.dm, &J));
144 PetscCall(DMCreateGlobalVector(user.dm, &X));
145
146 PetscCall(PetscMalloc3(user.Nspec + 1, &user.tchemwork, PetscSqr(user.Nspec + 1), &user.Jdense, user.Nspec + 1, &user.rows));
147
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));
158
159 ftime = 1.0;
160 PetscCall(TSSetMaxTime(ts, ftime));
161 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
162
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 */
173
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 }
181
182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183 Set runtime options
184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185 PetscCall(TSSetFromOptions(ts));
186
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));
195
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));
204
205 {
206 Vec max;
207 const char *const *names;
208 PetscInt i;
209 const PetscReal *bmax;
210
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 }
224
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 }
238
239 /*
240 Applies the second order centered difference diffusion operator on a one dimensional periodic domain
241 */
FormDiffusionFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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;
251
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));
261
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 }
271
272 /*
273 Produces the second order centered difference diffusion operator on a one dimensional periodic domain
274 */
FormDiffusionJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void * ptr)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];
282
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));
287
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 }
309
FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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;
317
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));
324
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 }
331
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 }
340
FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void * ptr)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;
347
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));
356
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));
361
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 }
380
FormInitialSolution(TS ts,Vec X,PetscCtx ctx)381 PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx 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;
395
396 PetscFunctionBeginUser;
397 PetscCall(VecZeroEntries(X));
398 PetscCall(TSGetDM(ts, &dm));
399 PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));
400
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, (int)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 }
416
417 /*
418 Routines for displaying the solutions
419 */
420 typedef struct {
421 PetscInt cell;
422 User user;
423 } UserLGCtx;
424
FormMoleFraction(UserLGCtx * ctx,Vec massf,Vec * molef)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;
432
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 }
447
MonitorCellDestroy(UserLGCtx ** uctx)448 static PetscErrorCode MonitorCellDestroy(UserLGCtx **uctx)
449 {
450 PetscFunctionBeginUser;
451 PetscCall(PetscFree(*uctx));
452 PetscFunctionReturn(PETSC_SUCCESS);
453 }
454
455 /*
456 Use TSMonitorLG to monitor the reactions in a particular cell
457 */
MonitorCell(TS ts,User user,PetscInt cell)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;
466
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, (PetscCtxDestroyFn *)MonitorCellDestroy, uctx));
480 PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
481 PetscFunctionReturn(PETSC_SUCCESS);
482 }
483