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 { \
64: PetscCheck(!ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in TChem library, return code %d", ierr); \
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;
84: PetscFunctionBeginUser;
85: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
86: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Chemistry solver options", "");
87: PetscCall(PetscOptionsString("-chem", "CHEMKIN input file", "", chemfile, chemfile, sizeof(chemfile), NULL));
88: PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD, chemfile, lchemfile, PETSC_MAX_PATH_LEN, &found));
89: PetscCheck(found, PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "Cannot download %s and no local version %s", chemfile, lchemfile);
90: PetscCall(PetscOptionsString("-thermo", "NASA thermo input file", "", thermofile, thermofile, sizeof(thermofile), NULL));
91: PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD, thermofile, lthermofile, PETSC_MAX_PATH_LEN, &found));
92: PetscCheck(found, PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "Cannot download %s and no local version %s", thermofile, lthermofile);
93: user.pressure = 1.01325e5; /* Pascal */
94: PetscCall(PetscOptionsReal("-pressure", "Pressure of reaction [Pa]", "", user.pressure, &user.pressure, NULL));
95: user.Tini = 1000; /* Kelvin */
96: PetscCall(PetscOptionsReal("-Tini", "Initial temperature [K]", "", user.Tini, &user.Tini, NULL));
97: PetscCall(PetscOptionsBool("-monitor_mass", "Monitor the total mass at each timestep", "", flg, &flg, NULL));
98: PetscCall(PetscOptionsBool("-monitor_temp", "Monitor the temperature each timestep", "", tflg, &tflg, NULL));
99: PetscOptionsEnd();
101: /* tchem requires periodic table in current directory */
102: PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD, periodic, lperiodic, PETSC_MAX_PATH_LEN, &found));
103: PetscCheck(found, PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "Cannot located required periodic table %s or local version %s", periodic, lperiodic);
105: PetscCallTC(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: PetscCall(PetscMalloc1((user.Nspec + 1) * LENGTHOFSPECNAME, &names));
113: PetscCall(PetscStrncpy(names, "Temp", (user.Nspec + 1) * LENGTHOFSPECNAME));
114: TC_getSnames(user.Nspec, names + LENGTHOFSPECNAME);
115: PetscCall(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: PetscCall(PetscStrArrayallocpy((const char *const *)snames, &user.snames));
119: PetscCall(PetscFree(snames));
120: PetscCall(PetscFree(names));
122: PetscCall(PetscMalloc3(user.Nspec + 1, &user.tchemwork, PetscSqr(user.Nspec + 1), &user.Jdense, user.Nspec + 1, &user.rows));
123: PetscCall(VecCreateSeq(PETSC_COMM_SELF, user.Nspec + 1, &X));
125: /* PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,PETSC_DECIDE,NULL,&J)); */
126: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, user.Nspec + 1, user.Nspec + 1, NULL, &J));
127: PetscCall(MatSetFromOptions(J));
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Create timestepping solver context
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
133: PetscCall(TSSetType(ts, TSARKIMEX));
134: PetscCall(TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE));
135: PetscCall(TSARKIMEXSetType(ts, TSARKIMEX4));
136: PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
137: PetscCall(TSSetRHSJacobian(ts, J, J, FormRHSJacobian, &user));
139: if (flg) PetscCall(TSMonitorSet(ts, MonitorMassConservation, NULL, NULL));
140: if (tflg) PetscCall(TSMonitorSet(ts, MonitorTempature, &user, NULL));
142: ftime = 1.0;
143: PetscCall(TSSetMaxTime(ts, ftime));
144: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Set initial conditions
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: PetscCall(FormInitialSolution(ts, X, &user));
150: PetscCall(TSSetSolution(ts, X));
151: dt = 1e-10; /* Initial time step */
152: PetscCall(TSSetTimeStep(ts, dt));
153: PetscCall(TSGetAdapt(ts, &adapt));
154: PetscCall(TSAdaptSetStepLimits(adapt, 1e-12, 1e-4)); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
155: PetscCall(TSSetMaxSNESFailures(ts, -1)); /* Retry step an unlimited number of times */
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Set runtime options
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160: PetscCall(TSSetFromOptions(ts));
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Set final conditions for sensitivities
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: PetscCall(VecDuplicate(X, &lambda));
166: PetscCall(TSSetCostGradients(ts, 1, &lambda, NULL));
167: PetscCall(VecSetValue(lambda, 0, 1.0, INSERT_VALUES));
168: PetscCall(VecAssemblyBegin(lambda));
169: PetscCall(VecAssemblyEnd(lambda));
171: PetscCall(TSGetTrajectory(ts, &tj));
172: if (tj) {
173: PetscCall(TSTrajectorySetVariableNames(tj, (const char *const *)user.snames));
174: PetscCall(TSTrajectorySetTransform(tj, (PetscErrorCode (*)(void *, Vec, Vec *))MassFractionToMoleFraction, NULL, &user));
175: }
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Pass information to graphical monitoring routine
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: PetscCall(TSMonitorLGSetVariableNames(ts, (const char *const *)user.snames));
181: PetscCall(TSMonitorLGSetTransform(ts, (PetscErrorCode (*)(void *, Vec, Vec *))MassFractionToMoleFraction, NULL, &user));
183: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: Solve ODE
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: PetscCall(TSSolve(ts, X));
187: PetscCall(TSGetSolveTime(ts, &ftime));
188: PetscCall(TSGetStepNumber(ts, &steps));
189: PetscCall(TSGetConvergedReason(ts, &reason));
190: PetscCall(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: PetscCall(TSMonitorEnvelopeGetBounds(ts,&max,NULL));
198: if (max) {
199: PetscCall(VecGetArrayRead(max,&bmax));
200: PetscCall(PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n"));
201: for (i=1; i<user.Nspec; i++) {
202: if (bmax[i] > .01) PetscCall(PetscPrintf(PETSC_COMM_SELF,"%s %g\n",user.snames[i],(double)bmax[i]));
203: }
204: PetscCall(VecRestoreArrayRead(max,&bmax));
205: }
206: }
208: Vec y;
209: PetscCall(MassFractionToMoleFraction(&user,X,&y));
210: PetscCall(PrintSpecies(&user,y));
211: PetscCall(VecDestroy(&y)); */
213: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: Free work space.
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: TC_reset();
217: PetscCall(PetscStrArrayDestroy(&user.snames));
218: PetscCall(MatDestroy(&J));
219: PetscCall(VecDestroy(&X));
220: PetscCall(VecDestroy(&lambda));
221: PetscCall(TSDestroy(&ts));
222: PetscCall(PetscFree3(user.tchemwork, user.Jdense, user.rows));
223: PetscCall(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;
233: PetscFunctionBeginUser;
234: PetscCall(VecGetArrayRead(X, &x));
235: PetscCall(VecGetArray(F, &f));
237: PetscCall(PetscArraycpy(user->tchemwork, x, user->Nspec + 1));
238: user->tchemwork[0] *= user->Tini; /* Dimensionalize */
239: PetscCallTC(TC_getSrc(user->tchemwork, user->Nspec + 1, f));
240: f[0] /= user->Tini; /* Non-dimensionalize */
242: PetscCall(VecRestoreArrayRead(X, &x));
243: PetscCall(VecRestoreArray(F, &f));
244: PetscFunctionReturn(PETSC_SUCCESS);
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;
253: PetscFunctionBeginUser;
254: PetscCall(VecGetArrayRead(X, &x));
255: PetscCall(PetscArraycpy(user->tchemwork, x, user->Nspec + 1));
256: PetscCall(VecRestoreArrayRead(X, &x));
257: user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */
258: PetscCall(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: PetscCall(MatSetOption(Pmat, MAT_ROW_ORIENTED, PETSC_FALSE));
264: PetscCall(MatSetOption(Pmat, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
265: PetscCall(MatZeroEntries(Pmat));
266: PetscCall(MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES));
268: PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
269: PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
270: if (Amat != Pmat) {
271: PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
272: PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
273: }
274: PetscFunctionReturn(PETSC_SUCCESS);
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;
288: PetscFunctionBeginUser;
289: PetscCall(VecZeroEntries(X));
290: PetscCall(VecGetArray(X, &x));
291: x[0] = 1.0; /* Non-dimensionalized by user->Tini */
293: PetscCall(PetscOptionsGetStringArray(NULL, NULL, "-initial_species", names, &smax, &flg));
294: PetscCheck(smax >= 2, PETSC_COMM_SELF, PETSC_ERR_USER, "Must provide at least two initial species");
295: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-initial_mole", molefracs, &mmax, &flg));
296: PetscCheck(smax == mmax, PETSC_COMM_SELF, PETSC_ERR_USER, "Must provide same number of initial species %" PetscInt_FMT " as initial moles %" PetscInt_FMT, smax, mmax);
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]));
302: PetscCheck(ispec >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Could not find species %s", names[i]);
303: PetscCall(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++) PetscCall(PetscFree(names[i]));
307: PetscCall(VecRestoreArray(X, &x));
308: /* PetscCall(PrintSpecies((User)ctx,X)); */
309: PetscCall(MoleFractionToMassFraction((User)ctx, X, &y));
310: PetscCall(VecCopy(y, X));
311: PetscCall(VecDestroy(&y));
312: PetscFunctionReturn(PETSC_SUCCESS);
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;
323: PetscFunctionBeginUser;
324: PetscCall(VecDuplicate(massf, molef));
325: PetscCall(VecGetArrayRead(massf, &maf));
326: PetscCall(VecGetArray(*molef, &mof));
327: mof[0] = maf[0]; /* copy over temperature */
328: TC_getMs2Ml((double *)maf + 1, user->Nspec, mof + 1);
329: PetscCall(VecRestoreArray(*molef, &mof));
330: PetscCall(VecRestoreArrayRead(massf, &maf));
331: PetscFunctionReturn(PETSC_SUCCESS);
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;
342: PetscFunctionBeginUser;
343: PetscCall(VecDuplicate(molef, massf));
344: PetscCall(VecGetArrayRead(molef, &mof));
345: PetscCall(VecGetArray(*massf, &maf));
346: maf[0] = mof[0]; /* copy over temperature */
347: TC_getMl2Ms((double *)mof + 1, user->Nspec, maf + 1);
348: PetscCall(VecRestoreArrayRead(molef, &mof));
349: PetscCall(VecRestoreArray(*massf, &maf));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: PetscErrorCode ComputeMassConservation(Vec x, PetscReal *mass, void *ctx)
354: {
355: PetscFunctionBeginUser;
356: PetscCall(VecSum(x, mass));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: PetscErrorCode MonitorMassConservation(TS ts, PetscInt step, PetscReal time, Vec x, void *ctx)
361: {
362: const PetscScalar *T;
363: PetscReal mass;
365: PetscFunctionBeginUser;
366: PetscCall(ComputeMassConservation(x, &mass, ctx));
367: PetscCall(VecGetArrayRead(x, &T));
368: mass -= PetscAbsScalar(T[0]);
369: PetscCall(VecRestoreArrayRead(x, &T));
370: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: PetscErrorCode MonitorTempature(TS ts, PetscInt step, PetscReal time, Vec x, void *ctx)
375: {
376: User user = (User)ctx;
377: const PetscScalar *T;
379: PetscFunctionBeginUser;
380: PetscCall(VecGetArrayRead(x, &T));
381: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep %" PetscInt_FMT " time %g temperature %g\n", step, (double)time, (double)(T[0] * user->Tini)));
382: PetscCall(VecRestoreArrayRead(x, &T));
383: PetscFunctionReturn(PETSC_SUCCESS);
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;
394: PetscFunctionBeginUser;
395: PetscCall(PetscMalloc1(n, &idx));
396: for (i = 0; i < n; i++) idx[i] = i;
397: PetscCall(VecGetArrayRead(molef, &mof));
398: PetscCall(PetscSortRealWithPermutation(n, mof, idx));
399: for (i = 0; i < n; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%6s %g\n", user->snames[idx[n - i - 1]], (double)PetscRealPart(mof[idx[n - i - 1]])));
400: PetscCall(PetscFree(idx));
401: PetscCall(VecRestoreArrayRead(molef, &mof));
402: PetscFunctionReturn(PETSC_SUCCESS);
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 unlimited -ts_adapt_dt_max 1e-4 -ts_arkimex_type 4 -ts_max_time .005
411: requires: !single
412: filter: grep -v iterations
414: TEST*/