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