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