Actual source code: extchemfield.c
petsc-3.14.6 2021-03-30
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 tempature 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
54: The solution for component i = 0 is the temperature.
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
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]))
62: FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species.
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: };
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*);
84: #define TCCHKERRQ(ierr) do {if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in TChem library, return code %d",ierr);} while (0)
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 */
102: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
103: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Chemistry solver options","");
104: PetscOptionsString("-chem","CHEMKIN input file","",chemfile,chemfile,sizeof(chemfile),NULL);
105: PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL);
106: user.pressure = 1.01325e5; /* Pascal */
107: PetscOptionsReal("-pressure","Pressure of reaction [Pa]","",user.pressure,&user.pressure,NULL);
108: user.Tini = 1550;
109: PetscOptionsReal("-Tini","Initial temperature [K]","",user.Tini,&user.Tini,NULL);
110: user.diffus = 100;
111: PetscOptionsReal("-diffus","Diffusion constant","",user.diffus,&user.diffus,NULL);
112: PetscOptionsBool("-draw_solution","Plot the solution for each cell","",showsolutions,&showsolutions,NULL);
113: user.diffusion = PETSC_TRUE;
114: PetscOptionsBool("-diffusion","Have diffusion","",user.diffusion,&user.diffusion,NULL);
115: user.reactions = PETSC_TRUE;
116: PetscOptionsBool("-reactions","Have reactions","",user.reactions,&user.reactions,NULL);
117: PetscOptionsEnd();
119: TC_initChem(chemfile, thermofile, 0, 1.0);TC
120: user.Nspec = TC_getNspec();
121: user.Nreac = TC_getNreac();
123: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,10,user.Nspec+1,1,NULL,&user.dm);
124: DMSetFromOptions(user.dm);
125: DMSetUp(user.dm);
126: DMDAGetInfo(user.dm,NULL,&ncells,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
127: user.dx = 1.0/ncells; /* Set the coordinates of the cell centers; note final ghost cell is at x coordinate 1.0 */
128: DMDASetUniformCoordinates(user.dm,0.0,1.0,0.0,1.0,0.0,1.0);
130: /* set the names of each field in the DMDA based on the species name */
131: PetscMalloc1((user.Nspec+1)*LENGTHOFSPECNAME,&names);
132: PetscStrcpy(names,"Temp");
133: TC_getSnames(user.Nspec,names+LENGTHOFSPECNAME);
134: PetscMalloc1((user.Nspec+2),&snames);
135: for (i=0; i<user.Nspec+1; i++) snames[i] = names+i*LENGTHOFSPECNAME;
136: snames[user.Nspec+1] = NULL;
137: DMDASetFieldNames(user.dm,(const char * const *)snames);
138: PetscFree(snames);
139: PetscFree(names);
142: DMCreateMatrix(user.dm,&J);
143: DMCreateGlobalVector(user.dm,&X);
145: PetscMalloc3(user.Nspec+1,&user.tchemwork,PetscSqr(user.Nspec+1),&user.Jdense,user.Nspec+1,&user.rows);
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Create timestepping solver context
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: TSCreate(PETSC_COMM_WORLD,&ts);
151: TSSetDM(ts,user.dm);
152: TSSetType(ts,TSARKIMEX);
153: TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);
154: TSARKIMEXSetType(ts,TSARKIMEX4);
155: TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);
156: TSSetRHSJacobian(ts,J,J,FormRHSJacobian,&user);
158: ftime = 1.0;
159: TSSetMaxTime(ts,ftime);
160: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Set initial conditions
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: FormInitialSolution(ts,X,&user);
166: TSSetSolution(ts,X);
167: dt = 1e-10; /* Initial time step */
168: TSSetTimeStep(ts,dt);
169: TSGetAdapt(ts,&adapt);
170: TSAdaptSetStepLimits(adapt,1e-12,1e-4); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
171: TSSetMaxSNESFailures(ts,-1); /* Retry step an unlimited number of times */
174: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: Pass information to graphical monitoring routine
176: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: if (showsolutions) {
178: DMDAGetCorners(user.dm,&xs,NULL,NULL,&xm,NULL,NULL);
179: for (i=xs;i<xs+xm;i++) {
180: MonitorCell(ts,&user,i);
181: }
182: }
184: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: Set runtime options
186: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187: TSSetFromOptions(ts);
189: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: Set final conditions for sensitivities
191: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192: DMCreateGlobalVector(user.dm,&lambda);
193: TSSetCostGradients(ts,1,&lambda,NULL);
194: VecSetValue(lambda,0,1.0,INSERT_VALUES);
195: VecAssemblyBegin(lambda);
196: VecAssemblyEnd(lambda);
198: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: Solve ODE
200: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201: TSSolve(ts,X);
202: TSGetSolveTime(ts,&ftime);
203: TSGetStepNumber(ts,&steps);
204: TSGetConvergedReason(ts,&reason);
205: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
207: {
208: Vec max;
209: const char * const *names;
210: PetscInt i;
211: const PetscReal *bmax;
213: TSMonitorEnvelopeGetBounds(ts,&max,NULL);
214: if (max) {
215: TSMonitorLGGetVariableNames(ts,&names);
216: if (names) {
217: VecGetArrayRead(max,&bmax);
218: PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n");
219: for (i=1; i<user.Nspec; i++) {
220: if (bmax[i] > .01) {PetscPrintf(PETSC_COMM_SELF,"%s %g\n",names[i],bmax[i]);}
221: }
222: VecRestoreArrayRead(max,&bmax);
223: }
224: }
225: }
227: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228: Free work space.
229: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
230: TC_reset();
231: DMDestroy(&user.dm);
232: MatDestroy(&J);
233: VecDestroy(&X);
234: VecDestroy(&lambda);
235: TSDestroy(&ts);
236: PetscFree3(user.tchemwork,user.Jdense,user.rows);
237: PetscFinalize();
238: return ierr;
239: }
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;
256: TSGetDM(ts,&dm);
257: DMDAGetInfo(dm,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
258: DMGetLocalVector(dm,&Xlocal);
259: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xlocal);
260: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xlocal);
261: DMDAVecGetArrayDOFRead(dm,Xlocal,&x);
262: DMDAVecGetArrayDOF(dm,F,&f);
263: DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);
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: DMDAVecRestoreArrayDOFRead(dm,Xlocal,&x);
272: DMDAVecRestoreArrayDOF(dm,F,&f);
273: DMRestoreLocalVector(dm,&Xlocal);
274: return(0);
275: }
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];
290: TSGetDM(ts,&dm);
291: DMDAGetInfo(dm,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
292: DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);
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: MatSetValuesStencil(Pmat,1,&row,3,col,values,ADD_VALUES);
305: }
306: }
307: MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
308: MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
309: return(0);
310: }
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;
322: if (user->reactions) {
323: TSGetDM(ts,&dm);
324: DMDAVecGetArrayDOFRead(dm,X,&x);
325: DMDAVecGetArrayDOF(dm,F,&f);
326: DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);
328: for (i=xs; i<xs+xm; i++) {
329: PetscArraycpy(user->tchemwork,x[i],user->Nspec+1);
330: user->tchemwork[0] *= user->Tini; /* Dimensionalize */
331: TC_getSrc(user->tchemwork,user->Nspec+1,f[i]);TC
332: f[i][0] /= user->Tini; /* Non-dimensionalize */
333: }
335: DMDAVecRestoreArrayDOFRead(dm,X,&x);
336: DMDAVecRestoreArrayDOF(dm,F,&f);
337: } else {
338: VecZeroEntries(F);
339: }
340: if (user->diffusion) {
341: FormDiffusionFunction(ts,t,X,F,ptr);
342: }
343: return(0);
344: }
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;
355: if (user->reactions) {
356: TSGetDM(ts,&dm);
357: MatZeroEntries(Pmat);
358: MatSetOption(Pmat,MAT_ROW_ORIENTED,PETSC_FALSE);
359: MatSetOption(Pmat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
360: DMDAVecGetArrayDOFRead(dm,X,&x);
361: DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);
363: for (i=xs; i<xs+xm; i++) {
364: PetscArraycpy(user->tchemwork,x[i],user->Nspec+1);
365: user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */
366: TC_getJacTYN(user->tchemwork,user->Nspec,user->Jdense,1);
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: MatSetValues(Pmat,M,user->rows,M,user->rows,user->Jdense,INSERT_VALUES);
372: }
373: DMDAVecRestoreArrayDOFRead(dm,X,&x);
374: MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
375: MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
376: } else {
377: MatZeroEntries(Pmat);
378: }
379: if (user->diffusion) {
380: FormDiffusionJacobian(ts,t,X,Amat,Pmat,ptr);
381: }
382: if (Amat != Pmat) {
383: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
384: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
385: }
386: return(0);
387: }
389: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
390: {
391: PetscScalar **x,*xc;
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;
403: VecZeroEntries(X);
404: TSGetDM(ts,&dm);
405: DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);
407: DMDAGetCoordinateArray(dm,&xc);
408: DMDAVecGetArrayDOF(dm,X,&x);
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: PetscPrintf(PETSC_COMM_SELF,"Species %d: %s %g\n",j,initial[j].name,initial[j].massfrac);
415: x[i][1+ispec] = initial[j].massfrac;
416: }
417: }
418: DMDAVecRestoreArrayDOF(dm,X,&x);
419: DMDARestoreCoordinateArray(dm,&xc);
420: return(0);
421: }
423: /*
424: Routines for displaying the solutions
425: */
426: typedef struct {
427: PetscInt cell;
428: User user;
429: } UserLGCtx;
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;
441: VecCreateSeq(PETSC_COMM_SELF,n,molef);
442: PetscMalloc1(user->Nspec,&M);
443: TC_getSmass(user->Nspec, M);
444: DMDAVecGetArrayDOFRead(user->dm,massf,&maf);
445: VecGetArray(*molef,&mof);
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: DMDAVecRestoreArrayDOFRead(user->dm,massf,&maf);
452: VecRestoreArray(*molef,&mof);
453: PetscFree(M);
454: return(0);
455: }
457: static PetscErrorCode MonitorCellDestroy(UserLGCtx *uctx)
458: {
462: PetscFree(uctx);
463: return(0);
464: }
466: /*
467: Use TSMonitorLG to monitor the reactions in a particular cell
468: */
469: static PetscErrorCode MonitorCell(TS ts,User user,PetscInt cell)
470: {
472: TSMonitorLGCtx ctx;
473: char **snames;
474: UserLGCtx *uctx;
475: char label[128];
476: PetscReal temp,*xc;
477: PetscMPIInt rank;
480: DMDAGetCoordinateArray(user->dm,&xc);
481: temp = 1.0 + .05*PetscSinScalar(2.*PETSC_PI*xc[cell]); /* Non-dimensionalized by user->Tini */
482: DMDARestoreCoordinateArray(user->dm,&xc);
483: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
484: PetscSNPrintf(label,sizeof(label),"Initial Temperature %g Cell %d Rank %d",(double)user->Tini*temp,(int)cell,rank);
485: TSMonitorLGCtxCreate(PETSC_COMM_SELF,NULL,label,PETSC_DECIDE,PETSC_DECIDE,600,400,1,&ctx);
486: DMDAGetFieldNames(user->dm,(const char * const **)&snames);
487: TSMonitorLGCtxSetVariableNames(ctx,(const char * const *)snames);
488: PetscNew(&uctx);
489: uctx->cell = cell;
490: uctx->user = user;
491: TSMonitorLGCtxSetTransform(ctx,(PetscErrorCode (*)(void*,Vec,Vec*))FormMoleFraction,(PetscErrorCode (*)(void*))MonitorCellDestroy,uctx);
492: TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
493: return(0);
494: }