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
main(int argc,char ** argv)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 */
FormDiffusionFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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 */
FormDiffusionJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void * ptr)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
FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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
FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void * ptr)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
FormInitialSolution(TS ts,Vec X,void * ctx)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
FormMoleFraction(UserLGCtx * ctx,Vec massf,Vec * molef)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
MonitorCellDestroy(UserLGCtx * uctx)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 */
MonitorCell(TS ts,User user,PetscInt cell)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