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 { PetscCheck(!ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in TChem library, return code %d", ierr); } while (0) 85 86 int main(int argc, char **argv) { 87 TS ts; /* time integrator */ 88 TSAdapt adapt; 89 Vec X; /* solution vector */ 90 Mat J; /* Jacobian matrix */ 91 PetscInt steps, ncells, xs, xm, i; 92 PetscReal ftime, dt; 93 char chemfile[PETSC_MAX_PATH_LEN] = "chem.inp", thermofile[PETSC_MAX_PATH_LEN] = "therm.dat"; 94 struct _User user; 95 TSConvergedReason reason; 96 PetscBool showsolutions = PETSC_FALSE; 97 char **snames, *names; 98 Vec lambda; /* used with TSAdjoint for sensitivities */ 99 100 PetscFunctionBeginUser; 101 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 102 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Chemistry solver options", ""); 103 PetscCall(PetscOptionsString("-chem", "CHEMKIN input file", "", chemfile, chemfile, sizeof(chemfile), NULL)); 104 PetscCall(PetscOptionsString("-thermo", "NASA thermo input file", "", thermofile, thermofile, sizeof(thermofile), NULL)); 105 user.pressure = 1.01325e5; /* Pascal */ 106 PetscCall(PetscOptionsReal("-pressure", "Pressure of reaction [Pa]", "", user.pressure, &user.pressure, NULL)); 107 user.Tini = 1550; 108 PetscCall(PetscOptionsReal("-Tini", "Initial temperature [K]", "", user.Tini, &user.Tini, NULL)); 109 user.diffus = 100; 110 PetscCall(PetscOptionsReal("-diffus", "Diffusion constant", "", user.diffus, &user.diffus, NULL)); 111 PetscCall(PetscOptionsBool("-draw_solution", "Plot the solution for each cell", "", showsolutions, &showsolutions, NULL)); 112 user.diffusion = PETSC_TRUE; 113 PetscCall(PetscOptionsBool("-diffusion", "Have diffusion", "", user.diffusion, &user.diffusion, NULL)); 114 user.reactions = PETSC_TRUE; 115 PetscCall(PetscOptionsBool("-reactions", "Have reactions", "", user.reactions, &user.reactions, NULL)); 116 PetscOptionsEnd(); 117 118 PetscCallTC(TC_initChem(chemfile, thermofile, 0, 1.0)); 119 user.Nspec = TC_getNspec(); 120 user.Nreac = TC_getNreac(); 121 122 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, user.Nspec + 1, 1, NULL, &user.dm)); 123 PetscCall(DMSetFromOptions(user.dm)); 124 PetscCall(DMSetUp(user.dm)); 125 PetscCall(DMDAGetInfo(user.dm, NULL, &ncells, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 126 user.dx = 1.0 / ncells; /* Set the coordinates of the cell centers; note final ghost cell is at x coordinate 1.0 */ 127 PetscCall(DMDASetUniformCoordinates(user.dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 128 129 /* set the names of each field in the DMDA based on the species name */ 130 PetscCall(PetscMalloc1((user.Nspec + 1) * LENGTHOFSPECNAME, &names)); 131 PetscCall(PetscStrcpy(names, "Temp")); 132 TC_getSnames(user.Nspec, names + LENGTHOFSPECNAME); 133 PetscCall(PetscMalloc1((user.Nspec + 2), &snames)); 134 for (i = 0; i < user.Nspec + 1; i++) snames[i] = names + i * LENGTHOFSPECNAME; 135 snames[user.Nspec + 1] = NULL; 136 PetscCall(DMDASetFieldNames(user.dm, (const char *const *)snames)); 137 PetscCall(PetscFree(snames)); 138 PetscCall(PetscFree(names)); 139 140 PetscCall(DMCreateMatrix(user.dm, &J)); 141 PetscCall(DMCreateGlobalVector(user.dm, &X)); 142 143 PetscCall(PetscMalloc3(user.Nspec + 1, &user.tchemwork, PetscSqr(user.Nspec + 1), &user.Jdense, user.Nspec + 1, &user.rows)); 144 145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 Create timestepping solver context 147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 148 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 149 PetscCall(TSSetDM(ts, user.dm)); 150 PetscCall(TSSetType(ts, TSARKIMEX)); 151 PetscCall(TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE)); 152 PetscCall(TSARKIMEXSetType(ts, TSARKIMEX4)); 153 PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user)); 154 PetscCall(TSSetRHSJacobian(ts, J, J, FormRHSJacobian, &user)); 155 156 ftime = 1.0; 157 PetscCall(TSSetMaxTime(ts, ftime)); 158 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 159 160 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 161 Set initial conditions 162 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 163 PetscCall(FormInitialSolution(ts, X, &user)); 164 PetscCall(TSSetSolution(ts, X)); 165 dt = 1e-10; /* Initial time step */ 166 PetscCall(TSSetTimeStep(ts, dt)); 167 PetscCall(TSGetAdapt(ts, &adapt)); 168 PetscCall(TSAdaptSetStepLimits(adapt, 1e-12, 1e-4)); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */ 169 PetscCall(TSSetMaxSNESFailures(ts, -1)); /* Retry step an unlimited number of times */ 170 171 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 172 Pass information to graphical monitoring routine 173 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 174 if (showsolutions) { 175 PetscCall(DMDAGetCorners(user.dm, &xs, NULL, NULL, &xm, NULL, NULL)); 176 for (i = xs; i < xs + xm; i++) PetscCall(MonitorCell(ts, &user, i)); 177 } 178 179 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 180 Set runtime options 181 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 182 PetscCall(TSSetFromOptions(ts)); 183 184 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185 Set final conditions for sensitivities 186 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 187 PetscCall(DMCreateGlobalVector(user.dm, &lambda)); 188 PetscCall(TSSetCostGradients(ts, 1, &lambda, NULL)); 189 PetscCall(VecSetValue(lambda, 0, 1.0, INSERT_VALUES)); 190 PetscCall(VecAssemblyBegin(lambda)); 191 PetscCall(VecAssemblyEnd(lambda)); 192 193 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194 Solve ODE 195 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 196 PetscCall(TSSolve(ts, X)); 197 PetscCall(TSGetSolveTime(ts, &ftime)); 198 PetscCall(TSGetStepNumber(ts, &steps)); 199 PetscCall(TSGetConvergedReason(ts, &reason)); 200 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps)); 201 202 { 203 Vec max; 204 const char *const *names; 205 PetscInt i; 206 const PetscReal *bmax; 207 208 PetscCall(TSMonitorEnvelopeGetBounds(ts, &max, NULL)); 209 if (max) { 210 PetscCall(TSMonitorLGGetVariableNames(ts, &names)); 211 if (names) { 212 PetscCall(VecGetArrayRead(max, &bmax)); 213 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species - maximum mass fraction\n")); 214 for (i = 1; i < user.Nspec; i++) { 215 if (bmax[i] > .01) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s %g\n", names[i], (double)bmax[i])); 216 } 217 PetscCall(VecRestoreArrayRead(max, &bmax)); 218 } 219 } 220 } 221 222 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 223 Free work space. 224 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 225 TC_reset(); 226 PetscCall(DMDestroy(&user.dm)); 227 PetscCall(MatDestroy(&J)); 228 PetscCall(VecDestroy(&X)); 229 PetscCall(VecDestroy(&lambda)); 230 PetscCall(TSDestroy(&ts)); 231 PetscCall(PetscFree3(user.tchemwork, user.Jdense, user.rows)); 232 PetscCall(PetscFinalize()); 233 return 0; 234 } 235 236 /* 237 Applies the second order centered difference diffusion operator on a one dimensional periodic domain 238 */ 239 static PetscErrorCode FormDiffusionFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 240 User user = (User)ptr; 241 PetscScalar **f; 242 const PetscScalar **x; 243 DM dm; 244 PetscInt i, xs, xm, j, dof; 245 Vec Xlocal; 246 PetscReal idx; 247 248 PetscFunctionBeginUser; 249 PetscCall(TSGetDM(ts, &dm)); 250 PetscCall(DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 251 PetscCall(DMGetLocalVector(dm, &Xlocal)); 252 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xlocal)); 253 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xlocal)); 254 PetscCall(DMDAVecGetArrayDOFRead(dm, Xlocal, &x)); 255 PetscCall(DMDAVecGetArrayDOF(dm, F, &f)); 256 PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL)); 257 258 idx = 1.0 * user->diffus / user->dx; 259 for (i = xs; i < xs + xm; i++) { 260 for (j = 0; j < dof; j++) f[i][j] += idx * (x[i + 1][j] - 2.0 * x[i][j] + x[i - 1][j]); 261 } 262 PetscCall(DMDAVecRestoreArrayDOFRead(dm, Xlocal, &x)); 263 PetscCall(DMDAVecRestoreArrayDOF(dm, F, &f)); 264 PetscCall(DMRestoreLocalVector(dm, &Xlocal)); 265 PetscFunctionReturn(0); 266 } 267 268 /* 269 Produces the second order centered difference diffusion operator on a one dimensional periodic domain 270 */ 271 static PetscErrorCode FormDiffusionJacobian(TS ts, PetscReal t, Vec X, Mat Amat, Mat Pmat, void *ptr) { 272 User user = (User)ptr; 273 DM dm; 274 PetscInt i, xs, xm, j, dof; 275 PetscReal idx, values[3]; 276 MatStencil row, col[3]; 277 278 PetscFunctionBeginUser; 279 PetscCall(TSGetDM(ts, &dm)); 280 PetscCall(DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 281 PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL)); 282 283 idx = 1.0 * user->diffus / user->dx; 284 values[0] = idx; 285 values[1] = -2.0 * idx; 286 values[2] = idx; 287 for (i = xs; i < xs + xm; i++) { 288 for (j = 0; j < dof; j++) { 289 row.i = i; 290 row.c = j; 291 col[0].i = i - 1; 292 col[0].c = j; 293 col[1].i = i; 294 col[1].c = j; 295 col[2].i = i + 1; 296 col[2].c = j; 297 PetscCall(MatSetValuesStencil(Pmat, 1, &row, 3, col, values, ADD_VALUES)); 298 } 299 } 300 PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY)); 301 PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY)); 302 PetscFunctionReturn(0); 303 } 304 305 static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 306 User user = (User)ptr; 307 PetscScalar **f; 308 const PetscScalar **x; 309 DM dm; 310 PetscInt i, xs, xm; 311 312 PetscFunctionBeginUser; 313 if (user->reactions) { 314 PetscCall(TSGetDM(ts, &dm)); 315 PetscCall(DMDAVecGetArrayDOFRead(dm, X, &x)); 316 PetscCall(DMDAVecGetArrayDOF(dm, F, &f)); 317 PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL)); 318 319 for (i = xs; i < xs + xm; i++) { 320 PetscCall(PetscArraycpy(user->tchemwork, x[i], user->Nspec + 1)); 321 user->tchemwork[0] *= user->Tini; /* Dimensionalize */ 322 PetscCallTC(TC_getSrc(user->tchemwork, user->Nspec + 1, f[i])); 323 f[i][0] /= user->Tini; /* Non-dimensionalize */ 324 } 325 326 PetscCall(DMDAVecRestoreArrayDOFRead(dm, X, &x)); 327 PetscCall(DMDAVecRestoreArrayDOF(dm, F, &f)); 328 } else { 329 PetscCall(VecZeroEntries(F)); 330 } 331 if (user->diffusion) PetscCall(FormDiffusionFunction(ts, t, X, F, ptr)); 332 PetscFunctionReturn(0); 333 } 334 335 static PetscErrorCode FormRHSJacobian(TS ts, PetscReal t, Vec X, Mat Amat, Mat Pmat, void *ptr) { 336 User user = (User)ptr; 337 const PetscScalar **x; 338 PetscInt M = user->Nspec + 1, i, j, xs, xm; 339 DM dm; 340 341 PetscFunctionBeginUser; 342 if (user->reactions) { 343 PetscCall(TSGetDM(ts, &dm)); 344 PetscCall(MatZeroEntries(Pmat)); 345 PetscCall(MatSetOption(Pmat, MAT_ROW_ORIENTED, PETSC_FALSE)); 346 PetscCall(MatSetOption(Pmat, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); 347 PetscCall(DMDAVecGetArrayDOFRead(dm, X, &x)); 348 PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL)); 349 350 for (i = xs; i < xs + xm; i++) { 351 PetscCall(PetscArraycpy(user->tchemwork, x[i], user->Nspec + 1)); 352 user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */ 353 PetscCall(TC_getJacTYN(user->tchemwork, user->Nspec, user->Jdense, 1)); 354 355 for (j = 0; j < M; j++) user->Jdense[j + 0 * M] /= user->Tini; /* Non-dimensionalize first column */ 356 for (j = 0; j < M; j++) user->Jdense[0 + j * M] /= user->Tini; /* Non-dimensionalize first row */ 357 for (j = 0; j < M; j++) user->rows[j] = i * M + j; 358 PetscCall(MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES)); 359 } 360 PetscCall(DMDAVecRestoreArrayDOFRead(dm, X, &x)); 361 PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY)); 362 PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY)); 363 } else { 364 PetscCall(MatZeroEntries(Pmat)); 365 } 366 if (user->diffusion) PetscCall(FormDiffusionJacobian(ts, t, X, Amat, Pmat, ptr)); 367 if (Amat != Pmat) { 368 PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY)); 369 PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY)); 370 } 371 PetscFunctionReturn(0); 372 } 373 374 PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx) { 375 PetscScalar **x, *xc; 376 struct { 377 const char *name; 378 PetscReal massfrac; 379 } initial[] = { 380 {"CH4", 0.0948178320887 }, 381 {"O2", 0.189635664177 }, 382 {"N2", 0.706766236705 }, 383 {"AR", 0.00878026702874} 384 }; 385 PetscInt i, j, xs, xm; 386 DM dm; 387 388 PetscFunctionBeginUser; 389 PetscCall(VecZeroEntries(X)); 390 PetscCall(TSGetDM(ts, &dm)); 391 PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL)); 392 393 PetscCall(DMDAGetCoordinateArray(dm, &xc)); 394 PetscCall(DMDAVecGetArrayDOF(dm, X, &x)); 395 for (i = xs; i < xs + xm; i++) { 396 x[i][0] = 1.0 + .05 * PetscSinScalar(2. * PETSC_PI * xc[i]); /* Non-dimensionalized by user->Tini */ 397 for (j = 0; j < PETSC_STATIC_ARRAY_LENGTH(initial); j++) { 398 int ispec = TC_getSpos(initial[j].name, strlen(initial[j].name)); 399 PetscCheck(ispec >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Could not find species %s", initial[j].name); 400 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species %d: %s %g\n", j, initial[j].name, (double)initial[j].massfrac)); 401 x[i][1 + ispec] = initial[j].massfrac; 402 } 403 } 404 PetscCall(DMDAVecRestoreArrayDOF(dm, X, &x)); 405 PetscCall(DMDARestoreCoordinateArray(dm, &xc)); 406 PetscFunctionReturn(0); 407 } 408 409 /* 410 Routines for displaying the solutions 411 */ 412 typedef struct { 413 PetscInt cell; 414 User user; 415 } UserLGCtx; 416 417 static PetscErrorCode FormMoleFraction(UserLGCtx *ctx, Vec massf, Vec *molef) { 418 User user = ctx->user; 419 PetscReal *M, tM = 0; 420 PetscInt i, n = user->Nspec + 1; 421 PetscScalar *mof; 422 const PetscScalar **maf; 423 424 PetscFunctionBeginUser; 425 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, molef)); 426 PetscCall(PetscMalloc1(user->Nspec, &M)); 427 TC_getSmass(user->Nspec, M); 428 PetscCall(DMDAVecGetArrayDOFRead(user->dm, massf, &maf)); 429 PetscCall(VecGetArray(*molef, &mof)); 430 mof[0] = maf[ctx->cell][0]; /* copy over temperature */ 431 for (i = 1; i < n; i++) tM += maf[ctx->cell][i] / M[i - 1]; 432 for (i = 1; i < n; i++) mof[i] = maf[ctx->cell][i] / (M[i - 1] * tM); 433 PetscCall(DMDAVecRestoreArrayDOFRead(user->dm, massf, &maf)); 434 PetscCall(VecRestoreArray(*molef, &mof)); 435 PetscCall(PetscFree(M)); 436 PetscFunctionReturn(0); 437 } 438 439 static PetscErrorCode MonitorCellDestroy(UserLGCtx *uctx) { 440 PetscFunctionBeginUser; 441 PetscCall(PetscFree(uctx)); 442 PetscFunctionReturn(0); 443 } 444 445 /* 446 Use TSMonitorLG to monitor the reactions in a particular cell 447 */ 448 static PetscErrorCode MonitorCell(TS ts, User user, PetscInt cell) { 449 TSMonitorLGCtx ctx; 450 char **snames; 451 UserLGCtx *uctx; 452 char label[128]; 453 PetscReal temp, *xc; 454 PetscMPIInt rank; 455 456 PetscFunctionBeginUser; 457 PetscCall(DMDAGetCoordinateArray(user->dm, &xc)); 458 temp = 1.0 + .05 * PetscSinScalar(2. * PETSC_PI * xc[cell]); /* Non-dimensionalized by user->Tini */ 459 PetscCall(DMDARestoreCoordinateArray(user->dm, &xc)); 460 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 461 PetscCall(PetscSNPrintf(label, sizeof(label), "Initial Temperature %g Cell %d Rank %d", (double)user->Tini * temp, (int)cell, rank)); 462 PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, label, PETSC_DECIDE, PETSC_DECIDE, 600, 400, 1, &ctx)); 463 PetscCall(DMDAGetFieldNames(user->dm, (const char *const **)&snames)); 464 PetscCall(TSMonitorLGCtxSetVariableNames(ctx, (const char *const *)snames)); 465 PetscCall(PetscNew(&uctx)); 466 uctx->cell = cell; 467 uctx->user = user; 468 PetscCall(TSMonitorLGCtxSetTransform(ctx, (PetscErrorCode(*)(void *, Vec, Vec *))FormMoleFraction, (PetscErrorCode(*)(void *))MonitorCellDestroy, uctx)); 469 PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy)); 470 PetscFunctionReturn(0); 471 } 472