Actual source code: extchem.c
1: static const char help[] = "Integrate chemistry using TChem.\n";
3: #include <petscts.h>
5: #if defined(PETSC_HAVE_TCHEM)
6: #if defined(MAX)
7: #undef MAX
8: #endif
9: #if defined(MIN)
10: #undef MIN
11: #endif
12: # include <TC_params.h>
13: # include <TC_interface.h>
14: #else
15: # error TChem is required for this example. Reconfigure PETSc using --download-tchem.
16: #endif
17: /*
18: See extchem.example.1 for how to run an example
20: See also h2_10sp.inp for another example
22: Determine sensitivity of final tempature on each variables initial conditions
23: -ts_dt 1.e-5 -ts_type cn -ts_adjoint_solve -ts_adjoint_view_solution draw
25: The solution for component i = 0 is the temperature.
27: The solution, i > 0, is the mass fraction, massf[i], of species i, i.e. mass of species i/ total mass of all species
29: The mole fraction molef[i], i > 0, is the number of moles of a species/ total number of moles of all species
30: Define M[i] = mass per mole of species i then
31: molef[i] = massf[i]/(M[i]*(sum_j massf[j]/M[j]))
33: FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species.
36: These are other data sets for other possible runs
37: https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/n_heptane_v3.1_therm.dat
38: https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/nc7_ver3.1_mech.txt
41: */
42: typedef struct _User *User;
43: struct _User {
44: PetscReal pressure;
45: int Nspec;
46: int Nreac;
47: PetscReal Tini;
48: double *tchemwork;
49: double *Jdense; /* Dense array workspace where Tchem computes the Jacobian */
50: PetscInt *rows;
51: char **snames;
52: };
55: static PetscErrorCode PrintSpecies(User,Vec);
56: static PetscErrorCode MassFractionToMoleFraction(User,Vec,Vec*);
57: static PetscErrorCode MoleFractionToMassFraction(User,Vec,Vec*);
58: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
59: static PetscErrorCode FormRHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
60: static PetscErrorCode FormInitialSolution(TS,Vec,void*);
61: static PetscErrorCode ComputeMassConservation(Vec,PetscReal*,void*);
62: static PetscErrorCode MonitorMassConservation(TS,PetscInt,PetscReal,Vec,void*);
63: static PetscErrorCode MonitorTempature(TS,PetscInt,PetscReal,Vec,void*);
65: #define TCCHKERRQ(ierr) do {if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in TChem library, return code %d",ierr);} while (0)
67: int main(int argc,char **argv)
68: {
69: TS ts; /* time integrator */
70: TSAdapt adapt;
71: Vec X,lambda; /* solution vector */
72: Mat J; /* Jacobian matrix */
73: PetscInt steps;
74: PetscErrorCode ierr;
75: PetscReal ftime,dt;
76: char chemfile[PETSC_MAX_PATH_LEN],thermofile[PETSC_MAX_PATH_LEN],lchemfile[PETSC_MAX_PATH_LEN],lthermofile[PETSC_MAX_PATH_LEN],lperiodic[PETSC_MAX_PATH_LEN];
77: const char *periodic = "file://${PETSC_DIR}/${PETSC_ARCH}/share/periodictable.dat";
78: struct _User user; /* user-defined work context */
79: TSConvergedReason reason;
80: char **snames,*names;
81: PetscInt i;
82: TSTrajectory tj;
83: PetscBool flg = PETSC_FALSE,tflg = PETSC_FALSE,found;
85: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
86: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Chemistry solver options","");
87: PetscOptionsString("-chem","CHEMKIN input file","",chemfile,chemfile,sizeof(chemfile),NULL);
88: PetscFileRetrieve(PETSC_COMM_WORLD,chemfile,lchemfile,PETSC_MAX_PATH_LEN,&found);
89: if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot download %s and no local version %s",chemfile,lchemfile);
90: PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL);
91: PetscFileRetrieve(PETSC_COMM_WORLD,thermofile,lthermofile,PETSC_MAX_PATH_LEN,&found);
92: if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot download %s and no local version %s",thermofile,lthermofile);
93: user.pressure = 1.01325e5; /* Pascal */
94: PetscOptionsReal("-pressure","Pressure of reaction [Pa]","",user.pressure,&user.pressure,NULL);
95: user.Tini = 1000; /* Kelvin */
96: PetscOptionsReal("-Tini","Initial temperature [K]","",user.Tini,&user.Tini,NULL);
97: PetscOptionsBool("-monitor_mass","Monitor the total mass at each timestep","",flg,&flg,NULL);
98: PetscOptionsBool("-monitor_temp","Monitor the tempature each timestep","",tflg,&tflg,NULL);
99: PetscOptionsEnd();
101: /* tchem requires periodic table in current directory */
102: PetscFileRetrieve(PETSC_COMM_WORLD,periodic,lperiodic,PETSC_MAX_PATH_LEN,&found);
103: if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot located required periodic table %s or local version %s",periodic,lperiodic);
105: TC_initChem(lchemfile, lthermofile, 0, 1.0);TC
106: TC_setThermoPres(user.pressure);
107: user.Nspec = TC_getNspec();
108: user.Nreac = TC_getNreac();
109: /*
110: Get names of all species in easy to use array
111: */
112: PetscMalloc1((user.Nspec+1)*LENGTHOFSPECNAME,&names);
113: PetscStrcpy(names,"Temp");
114: TC_getSnames(user.Nspec,names+LENGTHOFSPECNAME);
115: PetscMalloc1((user.Nspec+2),&snames);
116: for (i=0; i<user.Nspec+1; i++) snames[i] = names+i*LENGTHOFSPECNAME;
117: snames[user.Nspec+1] = NULL;
118: PetscStrArrayallocpy((const char *const *)snames,&user.snames);
119: PetscFree(snames);
120: PetscFree(names);
122: PetscMalloc3(user.Nspec+1,&user.tchemwork,PetscSqr(user.Nspec+1),&user.Jdense,user.Nspec+1,&user.rows);
123: VecCreateSeq(PETSC_COMM_SELF,user.Nspec+1,&X);
125: /* MatCreateSeqAIJ(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,PETSC_DECIDE,NULL,&J); */
126: MatCreateSeqDense(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,NULL,&J);
127: MatSetFromOptions(J);
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Create timestepping solver context
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: TSCreate(PETSC_COMM_WORLD,&ts);
133: TSSetType(ts,TSARKIMEX);
134: TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);
135: TSARKIMEXSetType(ts,TSARKIMEX4);
136: TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);
137: TSSetRHSJacobian(ts,J,J,FormRHSJacobian,&user);
139: if (flg) {
140: TSMonitorSet(ts,MonitorMassConservation,NULL,NULL);
141: }
142: if (tflg) {
143: TSMonitorSet(ts,MonitorTempature,&user,NULL);
144: }
146: ftime = 1.0;
147: TSSetMaxTime(ts,ftime);
148: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Set initial conditions
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: FormInitialSolution(ts,X,&user);
154: TSSetSolution(ts,X);
155: dt = 1e-10; /* Initial time step */
156: TSSetTimeStep(ts,dt);
157: TSGetAdapt(ts,&adapt);
158: TSAdaptSetStepLimits(adapt,1e-12,1e-4); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
159: TSSetMaxSNESFailures(ts,-1); /* Retry step an unlimited number of times */
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Set runtime options
163: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164: TSSetFromOptions(ts);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Set final conditions for sensitivities
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: VecDuplicate(X,&lambda);
170: TSSetCostGradients(ts,1,&lambda,NULL);
171: VecSetValue(lambda,0,1.0,INSERT_VALUES);
172: VecAssemblyBegin(lambda);
173: VecAssemblyEnd(lambda);
175: TSGetTrajectory(ts,&tj);
176: if (tj) {
177: TSTrajectorySetVariableNames(tj,(const char * const *)user.snames);
178: TSTrajectorySetTransform(tj,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user);
179: }
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: Pass information to graphical monitoring routine
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185: TSMonitorLGSetVariableNames(ts,(const char * const *)user.snames);
186: TSMonitorLGSetTransform(ts,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user);
188: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: Solve ODE
190: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191: TSSolve(ts,X);
192: TSGetSolveTime(ts,&ftime);
193: TSGetStepNumber(ts,&steps);
194: TSGetConvergedReason(ts,&reason);
195: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
197: /* {
198: Vec max;
199: PetscInt i;
200: const PetscReal *bmax;
202: TSMonitorEnvelopeGetBounds(ts,&max,NULL);
203: if (max) {
204: VecGetArrayRead(max,&bmax);
205: PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n");
206: for (i=1; i<user.Nspec; i++) {
207: if (bmax[i] > .01) {PetscPrintf(PETSC_COMM_SELF,"%s %g\n",user.snames[i],(double)bmax[i]);}
208: }
209: VecRestoreArrayRead(max,&bmax);
210: }
211: }
213: Vec y;
214: MassFractionToMoleFraction(&user,X,&y);
215: PrintSpecies(&user,y);
216: VecDestroy(&y); */
218: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219: Free work space.
220: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221: TC_reset();
222: PetscStrArrayDestroy(&user.snames);
223: MatDestroy(&J);
224: VecDestroy(&X);
225: VecDestroy(&lambda);
226: TSDestroy(&ts);
227: PetscFree3(user.tchemwork,user.Jdense,user.rows);
228: PetscFinalize();
229: return ierr;
230: }
232: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
233: {
234: User user = (User)ptr;
235: PetscErrorCode ierr;
236: PetscScalar *f;
237: const PetscScalar *x;
240: VecGetArrayRead(X,&x);
241: VecGetArray(F,&f);
243: PetscArraycpy(user->tchemwork,x,user->Nspec+1);
244: user->tchemwork[0] *= user->Tini; /* Dimensionalize */
245: TC_getSrc(user->tchemwork,user->Nspec+1,f);TC
246: f[0] /= user->Tini; /* Non-dimensionalize */
248: VecRestoreArrayRead(X,&x);
249: VecRestoreArray(F,&f);
250: return(0);
251: }
253: static PetscErrorCode FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr)
254: {
255: User user = (User)ptr;
256: PetscErrorCode ierr;
257: const PetscScalar *x;
258: PetscInt M = user->Nspec+1,i;
261: VecGetArrayRead(X,&x);
262: PetscArraycpy(user->tchemwork,x,user->Nspec+1);
263: VecRestoreArrayRead(X,&x);
264: user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */
265: TC_getJacTYN(user->tchemwork,user->Nspec,user->Jdense,1);
267: for (i=0; i<M; i++) user->Jdense[i + 0*M] /= user->Tini; /* Non-dimensionalize first column */
268: for (i=0; i<M; i++) user->Jdense[0 + i*M] /= user->Tini; /* Non-dimensionalize first row */
269: for (i=0; i<M; i++) user->rows[i] = i;
270: MatSetOption(Pmat,MAT_ROW_ORIENTED,PETSC_FALSE);
271: MatSetOption(Pmat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
272: MatZeroEntries(Pmat);
273: MatSetValues(Pmat,M,user->rows,M,user->rows,user->Jdense,INSERT_VALUES);
275: MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
276: MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
277: if (Amat != Pmat) {
278: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
279: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
280: }
281: return(0);
282: }
284: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
285: {
286: PetscScalar *x;
288: PetscInt i;
289: Vec y;
290: const PetscInt maxspecies = 10;
291: PetscInt smax = maxspecies,mmax = maxspecies;
292: char *names[maxspecies];
293: PetscReal molefracs[maxspecies],sum;
294: PetscBool flg;
297: VecZeroEntries(X);
298: VecGetArray(X,&x);
299: x[0] = 1.0; /* Non-dimensionalized by user->Tini */
301: PetscOptionsGetStringArray(NULL,NULL,"-initial_species",names,&smax,&flg);
302: if (smax < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Must provide at least two initial species");
303: PetscOptionsGetRealArray(NULL,NULL,"-initial_mole",molefracs,&mmax,&flg);
304: if (smax != mmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Must provide same number of initial species %D as initial moles %D",smax,mmax);
305: sum = 0;
306: for (i=0; i<smax; i++) sum += molefracs[i];
307: for (i=0; i<smax; i++) molefracs[i] = molefracs[i]/sum;
308: for (i=0; i<smax; i++) {
309: int ispec = TC_getSpos(names[i], strlen(names[i]));
310: if (ispec < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Could not find species %s",names[i]);
311: PetscPrintf(PETSC_COMM_SELF,"Species %d: %s %g\n",i,names[i],molefracs[i]);
312: x[1+ispec] = molefracs[i];
313: }
314: for (i=0; i<smax; i++) {
315: PetscFree(names[i]);
316: }
317: VecRestoreArray(X,&x);
318: /* PrintSpecies((User)ctx,X); */
319: MoleFractionToMassFraction((User)ctx,X,&y);
320: VecCopy(y,X);
321: VecDestroy(&y);
322: return(0);
323: }
325: /*
326: Converts the input vector which is in mass fractions (used by tchem) to mole fractions
327: */
328: PetscErrorCode MassFractionToMoleFraction(User user,Vec massf,Vec *molef)
329: {
330: PetscErrorCode ierr;
331: PetscScalar *mof;
332: const PetscScalar *maf;
335: VecDuplicate(massf,molef);
336: VecGetArrayRead(massf,&maf);
337: VecGetArray(*molef,&mof);
338: mof[0] = maf[0]; /* copy over temperature */
339: TC_getMs2Ml((double*)maf+1,user->Nspec,mof+1);
340: VecRestoreArray(*molef,&mof);
341: VecRestoreArrayRead(massf,&maf);
342: return(0);
343: }
345: /*
346: Converts the input vector which is in mole fractions to mass fractions (used by tchem)
347: */
348: PetscErrorCode MoleFractionToMassFraction(User user,Vec molef,Vec *massf)
349: {
350: PetscErrorCode ierr;
351: const PetscScalar *mof;
352: PetscScalar *maf;
355: VecDuplicate(molef,massf);
356: VecGetArrayRead(molef,&mof);
357: VecGetArray(*massf,&maf);
358: maf[0] = mof[0]; /* copy over temperature */
359: TC_getMl2Ms((double*)mof+1,user->Nspec,maf+1);
360: VecRestoreArrayRead(molef,&mof);
361: VecRestoreArray(*massf,&maf);
362: return(0);
363: }
365: PetscErrorCode ComputeMassConservation(Vec x,PetscReal *mass,void* ctx)
366: {
370: VecSum(x,mass);
371: return(0);
372: }
374: PetscErrorCode MonitorMassConservation(TS ts,PetscInt step,PetscReal time,Vec x,void* ctx)
375: {
376: const PetscScalar *T;
377: PetscReal mass;
378: PetscErrorCode ierr;
381: ComputeMassConservation(x,&mass,ctx);
382: VecGetArrayRead(x,&T);
383: mass -= PetscAbsScalar(T[0]);
384: VecRestoreArrayRead(x,&T);
385: PetscPrintf(PETSC_COMM_WORLD,"Timestep %D time %g percent mass lost or gained %g\n",step,(double)time,(double)100.*(1.0 - mass));
386: return(0);
387: }
389: PetscErrorCode MonitorTempature(TS ts,PetscInt step,PetscReal time,Vec x,void* ctx)
390: {
391: User user = (User) ctx;
392: const PetscScalar *T;
393: PetscErrorCode ierr;
396: VecGetArrayRead(x,&T);
397: PetscPrintf(PETSC_COMM_WORLD,"Timestep %D time %g tempature %g\n",step,(double)time,(double)T[0]*user->Tini);
398: VecRestoreArrayRead(x,&T);
399: return(0);
400: }
402: /*
403: Prints out each species with its name
404: */
405: PETSC_UNUSED PetscErrorCode PrintSpecies(User user,Vec molef)
406: {
407: PetscErrorCode ierr;
408: const PetscScalar *mof;
409: PetscInt i,*idx,n = user->Nspec+1;
412: PetscMalloc1(n,&idx);
413: for (i=0; i<n;i++) idx[i] = i;
414: VecGetArrayRead(molef,&mof);
415: PetscSortRealWithPermutation(n,mof,idx);
416: for (i=0; i<n; i++) {
417: PetscPrintf(PETSC_COMM_SELF,"%6s %g\n",user->snames[idx[n-i-1]],mof[idx[n-i-1]]);
418: }
419: PetscFree(idx);
420: VecRestoreArrayRead(molef,&mof);
421: return(0);
422: }
424: /*TEST
425: build:
426: requires: tchem
428: test:
429: args: -chem http://combustion.berkeley.edu/gri_mech/version30/files30/grimech30.dat -thermo http://combustion.berkeley.edu/gri_mech/version30/files30/thermo30.dat -initial_species CH4,O2,N2,AR -initial_mole 0.0948178320887,0.189635664177,0.706766236705,0.00878026702874 -Tini 1500 -Tini 1500 -ts_arkimex_fully_implicit -ts_max_snes_failures -1 -ts_adapt_dt_max 1e-4 -ts_arkimex_type 4 -ts_max_time .005
430: requires: !single
431: filter: grep -v iterations
433: TEST*/