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 -1 -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 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 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 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, (char *)0, 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 */ 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 */ 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 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 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 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; 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, 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 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 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 */ 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, (PetscErrorCode(*)(void *))MonitorCellDestroy, uctx)); 480 PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy)); 481 PetscFunctionReturn(PETSC_SUCCESS); 482 } 483