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