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 temperature 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.
35: These are other data sets for other possible runs
36: https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/n_heptane_v3.1_therm.dat
37: https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/nc7_ver3.1_mech.txt
39: */
40: typedef struct _User *User;
41: struct _User {
42: PetscReal pressure;
43: int Nspec;
44: int Nreac;
45: PetscReal Tini;
46: double *tchemwork;
47: double *Jdense; /* Dense array workspace where Tchem computes the Jacobian */
48: PetscInt *rows;
49: char **snames;
50: };
52: static PetscErrorCode PrintSpecies(User, Vec);
53: static PetscErrorCode MassFractionToMoleFraction(User, Vec, Vec *);
54: static PetscErrorCode MoleFractionToMassFraction(User, Vec, Vec *);
55: static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
56: static PetscErrorCode FormRHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
57: static PetscErrorCode FormInitialSolution(TS, Vec, void *);
58: static PetscErrorCode ComputeMassConservation(Vec, PetscReal *, void *);
59: static PetscErrorCode MonitorMassConservation(TS, PetscInt, PetscReal, Vec, void *);
60: static PetscErrorCode MonitorTempature(TS, PetscInt, PetscReal, Vec, void *);
62: #define PetscCallTC(ierr) \
63: do { \
65: } 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: PetscReal ftime, dt;
75: 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];
76: const char *periodic = "file://${PETSC_DIR}/${PETSC_ARCH}/share/periodictable.dat";
77: struct _User user; /* user-defined work context */
78: TSConvergedReason reason;
79: char **snames, *names;
80: PetscInt i;
81: TSTrajectory tj;
82: PetscBool flg = PETSC_FALSE, tflg = PETSC_FALSE, found;
85: PetscInitialize(&argc, &argv, (char *)0, help);
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);
90: PetscOptionsString("-thermo", "NASA thermo input file", "", thermofile, thermofile, sizeof(thermofile), NULL);
91: PetscFileRetrieve(PETSC_COMM_WORLD, thermofile, lthermofile, PETSC_MAX_PATH_LEN, &found);
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 temperature 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);
105: TC_initChem(lchemfile, lthermofile, 0, 1.0);
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) TSMonitorSet(ts, MonitorMassConservation, NULL, NULL);
140: if (tflg) TSMonitorSet(ts, MonitorTempature, &user, NULL);
142: ftime = 1.0;
143: TSSetMaxTime(ts, ftime);
144: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Set initial conditions
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: FormInitialSolution(ts, X, &user);
150: TSSetSolution(ts, X);
151: dt = 1e-10; /* Initial time step */
152: TSSetTimeStep(ts, dt);
153: TSGetAdapt(ts, &adapt);
154: TSAdaptSetStepLimits(adapt, 1e-12, 1e-4); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
155: TSSetMaxSNESFailures(ts, -1); /* Retry step an unlimited number of times */
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Set runtime options
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160: TSSetFromOptions(ts);
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Set final conditions for sensitivities
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: VecDuplicate(X, &lambda);
166: TSSetCostGradients(ts, 1, &lambda, NULL);
167: VecSetValue(lambda, 0, 1.0, INSERT_VALUES);
168: VecAssemblyBegin(lambda);
169: VecAssemblyEnd(lambda);
171: TSGetTrajectory(ts, &tj);
172: if (tj) {
173: TSTrajectorySetVariableNames(tj, (const char *const *)user.snames);
174: TSTrajectorySetTransform(tj, (PetscErrorCode(*)(void *, Vec, Vec *))MassFractionToMoleFraction, NULL, &user);
175: }
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Pass information to graphical monitoring routine
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: TSMonitorLGSetVariableNames(ts, (const char *const *)user.snames);
181: TSMonitorLGSetTransform(ts, (PetscErrorCode(*)(void *, Vec, Vec *))MassFractionToMoleFraction, NULL, &user);
183: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: Solve ODE
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: TSSolve(ts, X);
187: TSGetSolveTime(ts, &ftime);
188: TSGetStepNumber(ts, &steps);
189: TSGetConvergedReason(ts, &reason);
190: PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps);
192: /* {
193: Vec max;
194: PetscInt i;
195: const PetscReal *bmax;
197: TSMonitorEnvelopeGetBounds(ts,&max,NULL);
198: if (max) {
199: VecGetArrayRead(max,&bmax);
200: PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n");
201: for (i=1; i<user.Nspec; i++) {
202: if (bmax[i] > .01) PetscPrintf(PETSC_COMM_SELF,"%s %g\n",user.snames[i],(double)bmax[i]);
203: }
204: VecRestoreArrayRead(max,&bmax);
205: }
206: }
208: Vec y;
209: MassFractionToMoleFraction(&user,X,&y);
210: PrintSpecies(&user,y);
211: VecDestroy(&y); */
213: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: Free work space.
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: TC_reset();
217: PetscStrArrayDestroy(&user.snames);
218: MatDestroy(&J);
219: VecDestroy(&X);
220: VecDestroy(&lambda);
221: TSDestroy(&ts);
222: PetscFree3(user.tchemwork, user.Jdense, user.rows);
223: PetscFinalize();
224: return 0;
225: }
227: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
228: {
229: User user = (User)ptr;
230: PetscScalar *f;
231: const PetscScalar *x;
234: VecGetArrayRead(X, &x);
235: VecGetArray(F, &f);
237: PetscArraycpy(user->tchemwork, x, user->Nspec + 1);
238: user->tchemwork[0] *= user->Tini; /* Dimensionalize */
239: TC_getSrc(user->tchemwork, user->Nspec + 1, f);
240: f[0] /= user->Tini; /* Non-dimensionalize */
242: VecRestoreArrayRead(X, &x);
243: VecRestoreArray(F, &f);
244: return 0;
245: }
247: static PetscErrorCode FormRHSJacobian(TS ts, PetscReal t, Vec X, Mat Amat, Mat Pmat, void *ptr)
248: {
249: User user = (User)ptr;
250: const PetscScalar *x;
251: PetscInt M = user->Nspec + 1, i;
254: VecGetArrayRead(X, &x);
255: PetscArraycpy(user->tchemwork, x, user->Nspec + 1);
256: VecRestoreArrayRead(X, &x);
257: user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */
258: TC_getJacTYN(user->tchemwork, user->Nspec, user->Jdense, 1);
260: for (i = 0; i < M; i++) user->Jdense[i + 0 * M] /= user->Tini; /* Non-dimensionalize first column */
261: for (i = 0; i < M; i++) user->Jdense[0 + i * M] /= user->Tini; /* Non-dimensionalize first row */
262: for (i = 0; i < M; i++) user->rows[i] = i;
263: MatSetOption(Pmat, MAT_ROW_ORIENTED, PETSC_FALSE);
264: MatSetOption(Pmat, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
265: MatZeroEntries(Pmat);
266: MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES);
268: MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY);
269: MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY);
270: if (Amat != Pmat) {
271: MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
272: MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
273: }
274: return 0;
275: }
277: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
278: {
279: PetscScalar *x;
280: PetscInt i;
281: Vec y;
282: const PetscInt maxspecies = 10;
283: PetscInt smax = maxspecies, mmax = maxspecies;
284: char *names[maxspecies];
285: PetscReal molefracs[maxspecies], sum;
286: PetscBool flg;
289: VecZeroEntries(X);
290: VecGetArray(X, &x);
291: x[0] = 1.0; /* Non-dimensionalized by user->Tini */
293: PetscOptionsGetStringArray(NULL, NULL, "-initial_species", names, &smax, &flg);
295: PetscOptionsGetRealArray(NULL, NULL, "-initial_mole", molefracs, &mmax, &flg);
297: sum = 0;
298: for (i = 0; i < smax; i++) sum += molefracs[i];
299: for (i = 0; i < smax; i++) molefracs[i] = molefracs[i] / sum;
300: for (i = 0; i < smax; i++) {
301: int ispec = TC_getSpos(names[i], strlen(names[i]));
303: PetscPrintf(PETSC_COMM_SELF, "Species %" PetscInt_FMT ": %s %g\n", i, names[i], (double)molefracs[i]);
304: x[1 + ispec] = molefracs[i];
305: }
306: for (i = 0; i < smax; i++) PetscFree(names[i]);
307: VecRestoreArray(X, &x);
308: /* PrintSpecies((User)ctx,X); */
309: MoleFractionToMassFraction((User)ctx, X, &y);
310: VecCopy(y, X);
311: VecDestroy(&y);
312: return 0;
313: }
315: /*
316: Converts the input vector which is in mass fractions (used by tchem) to mole fractions
317: */
318: PetscErrorCode MassFractionToMoleFraction(User user, Vec massf, Vec *molef)
319: {
320: PetscScalar *mof;
321: const PetscScalar *maf;
324: VecDuplicate(massf, molef);
325: VecGetArrayRead(massf, &maf);
326: VecGetArray(*molef, &mof);
327: mof[0] = maf[0]; /* copy over temperature */
328: TC_getMs2Ml((double *)maf + 1, user->Nspec, mof + 1);
329: VecRestoreArray(*molef, &mof);
330: VecRestoreArrayRead(massf, &maf);
331: return 0;
332: }
334: /*
335: Converts the input vector which is in mole fractions to mass fractions (used by tchem)
336: */
337: PetscErrorCode MoleFractionToMassFraction(User user, Vec molef, Vec *massf)
338: {
339: const PetscScalar *mof;
340: PetscScalar *maf;
343: VecDuplicate(molef, massf);
344: VecGetArrayRead(molef, &mof);
345: VecGetArray(*massf, &maf);
346: maf[0] = mof[0]; /* copy over temperature */
347: TC_getMl2Ms((double *)mof + 1, user->Nspec, maf + 1);
348: VecRestoreArrayRead(molef, &mof);
349: VecRestoreArray(*massf, &maf);
350: return 0;
351: }
353: PetscErrorCode ComputeMassConservation(Vec x, PetscReal *mass, void *ctx)
354: {
356: VecSum(x, mass);
357: return 0;
358: }
360: PetscErrorCode MonitorMassConservation(TS ts, PetscInt step, PetscReal time, Vec x, void *ctx)
361: {
362: const PetscScalar *T;
363: PetscReal mass;
366: ComputeMassConservation(x, &mass, ctx);
367: VecGetArrayRead(x, &T);
368: mass -= PetscAbsScalar(T[0]);
369: VecRestoreArrayRead(x, &T);
370: PetscPrintf(PETSC_COMM_WORLD, "Timestep %" PetscInt_FMT " time %g percent mass lost or gained %g\n", step, (double)time, (double)(100. * (1.0 - mass)));
371: return 0;
372: }
374: PetscErrorCode MonitorTempature(TS ts, PetscInt step, PetscReal time, Vec x, void *ctx)
375: {
376: User user = (User)ctx;
377: const PetscScalar *T;
380: VecGetArrayRead(x, &T);
381: PetscPrintf(PETSC_COMM_WORLD, "Timestep %" PetscInt_FMT " time %g temperature %g\n", step, (double)time, (double)(T[0] * user->Tini));
382: VecRestoreArrayRead(x, &T);
383: return 0;
384: }
386: /*
387: Prints out each species with its name
388: */
389: PETSC_UNUSED PetscErrorCode PrintSpecies(User user, Vec molef)
390: {
391: const PetscScalar *mof;
392: PetscInt i, *idx, n = user->Nspec + 1;
395: PetscMalloc1(n, &idx);
396: for (i = 0; i < n; i++) idx[i] = i;
397: VecGetArrayRead(molef, &mof);
398: PetscSortRealWithPermutation(n, mof, idx);
399: for (i = 0; i < n; i++) PetscPrintf(PETSC_COMM_SELF, "%6s %g\n", user->snames[idx[n - i - 1]], (double)PetscRealPart(mof[idx[n - i - 1]]));
400: PetscFree(idx);
401: VecRestoreArrayRead(molef, &mof);
402: return 0;
403: }
405: /*TEST
406: build:
407: requires: tchem
409: test:
410: 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
411: requires: !single
412: filter: grep -v iterations
414: TEST*/