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_time_step 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: const PetscReal *bmax;
196: PetscCall(TSMonitorEnvelopeGetBounds(ts,&max,NULL));
197: if (max) {
198: PetscCall(VecGetArrayRead(max,&bmax));
199: PetscCall(PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n"));
200: for (PetscInt i=1; i<user.Nspec; i++) {
201: if (bmax[i] > .01) PetscCall(PetscPrintf(PETSC_COMM_SELF,"%s %g\n",user.snames[i],(double)bmax[i]));
202: }
203: PetscCall(VecRestoreArrayRead(max,&bmax));
204: }
205: }
207: Vec y;
208: PetscCall(MassFractionToMoleFraction(&user,X,&y));
209: PetscCall(PrintSpecies(&user,y));
210: PetscCall(VecDestroy(&y)); */
212: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213: Free work space.
214: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
215: TC_reset();
216: PetscCall(PetscStrArrayDestroy(&user.snames));
217: PetscCall(MatDestroy(&J));
218: PetscCall(VecDestroy(&X));
219: PetscCall(VecDestroy(&lambda));
220: PetscCall(TSDestroy(&ts));
221: PetscCall(PetscFree3(user.tchemwork, user.Jdense, user.rows));
222: PetscCall(PetscFinalize());
223: return 0;
224: }
226: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
227: {
228: User user = (User)ptr;
229: PetscScalar *f;
230: const PetscScalar *x;
232: PetscFunctionBeginUser;
233: PetscCall(VecGetArrayRead(X, &x));
234: PetscCall(VecGetArray(F, &f));
236: PetscCall(PetscArraycpy(user->tchemwork, x, user->Nspec + 1));
237: user->tchemwork[0] *= user->Tini; /* Dimensionalize */
238: PetscCallTC(TC_getSrc(user->tchemwork, user->Nspec + 1, f));
239: f[0] /= user->Tini; /* Non-dimensionalize */
241: PetscCall(VecRestoreArrayRead(X, &x));
242: PetscCall(VecRestoreArray(F, &f));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: static PetscErrorCode FormRHSJacobian(TS ts, PetscReal t, Vec X, Mat Amat, Mat Pmat, void *ptr)
247: {
248: User user = (User)ptr;
249: const PetscScalar *x;
250: PetscInt M = user->Nspec + 1, i;
252: PetscFunctionBeginUser;
253: PetscCall(VecGetArrayRead(X, &x));
254: PetscCall(PetscArraycpy(user->tchemwork, x, user->Nspec + 1));
255: PetscCall(VecRestoreArrayRead(X, &x));
256: user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */
257: PetscCall(TC_getJacTYN(user->tchemwork, user->Nspec, user->Jdense, 1));
259: for (i = 0; i < M; i++) user->Jdense[i + 0 * M] /= user->Tini; /* Non-dimensionalize first column */
260: for (i = 0; i < M; i++) user->Jdense[0 + i * M] /= user->Tini; /* Non-dimensionalize first row */
261: for (i = 0; i < M; i++) user->rows[i] = i;
262: PetscCall(MatSetOption(Pmat, MAT_ROW_ORIENTED, PETSC_FALSE));
263: PetscCall(MatSetOption(Pmat, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
264: PetscCall(MatZeroEntries(Pmat));
265: PetscCall(MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES));
267: PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
268: PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
269: if (Amat != Pmat) {
270: PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
271: PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
272: }
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx ctx)
277: {
278: PetscScalar *x;
279: PetscInt i;
280: Vec y;
281: const PetscInt maxspecies = 10;
282: PetscInt smax = maxspecies, mmax = maxspecies;
283: char *names[maxspecies];
284: PetscReal molefracs[maxspecies], sum;
285: PetscBool flg;
287: PetscFunctionBeginUser;
288: PetscCall(VecZeroEntries(X));
289: PetscCall(VecGetArray(X, &x));
290: x[0] = 1.0; /* Non-dimensionalized by user->Tini */
292: PetscCall(PetscOptionsGetStringArray(NULL, NULL, "-initial_species", names, &smax, &flg));
293: PetscCheck(smax >= 2, PETSC_COMM_SELF, PETSC_ERR_USER, "Must provide at least two initial species");
294: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-initial_mole", molefracs, &mmax, &flg));
295: 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);
296: sum = 0;
297: for (i = 0; i < smax; i++) sum += molefracs[i];
298: for (i = 0; i < smax; i++) molefracs[i] = molefracs[i] / sum;
299: for (i = 0; i < smax; i++) {
300: int ispec = TC_getSpos(names[i], (int)strlen(names[i]));
301: PetscCheck(ispec >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Could not find species %s", names[i]);
302: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species %" PetscInt_FMT ": %s %g\n", i, names[i], (double)molefracs[i]));
303: x[1 + ispec] = molefracs[i];
304: }
305: for (i = 0; i < smax; i++) PetscCall(PetscFree(names[i]));
306: PetscCall(VecRestoreArray(X, &x));
307: /* PetscCall(PrintSpecies((User)ctx,X)); */
308: PetscCall(MoleFractionToMassFraction((User)ctx, X, &y));
309: PetscCall(VecCopy(y, X));
310: PetscCall(VecDestroy(&y));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*
315: Converts the input vector which is in mass fractions (used by tchem) to mole fractions
316: */
317: PetscErrorCode MassFractionToMoleFraction(User user, Vec massf, Vec *molef)
318: {
319: PetscScalar *mof;
320: const PetscScalar *maf;
322: PetscFunctionBeginUser;
323: PetscCall(VecDuplicate(massf, molef));
324: PetscCall(VecGetArrayRead(massf, &maf));
325: PetscCall(VecGetArray(*molef, &mof));
326: mof[0] = maf[0]; /* copy over temperature */
327: TC_getMs2Ml((double *)maf + 1, user->Nspec, mof + 1);
328: PetscCall(VecRestoreArray(*molef, &mof));
329: PetscCall(VecRestoreArrayRead(massf, &maf));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: /*
334: Converts the input vector which is in mole fractions to mass fractions (used by tchem)
335: */
336: PetscErrorCode MoleFractionToMassFraction(User user, Vec molef, Vec *massf)
337: {
338: const PetscScalar *mof;
339: PetscScalar *maf;
341: PetscFunctionBeginUser;
342: PetscCall(VecDuplicate(molef, massf));
343: PetscCall(VecGetArrayRead(molef, &mof));
344: PetscCall(VecGetArray(*massf, &maf));
345: maf[0] = mof[0]; /* copy over temperature */
346: TC_getMl2Ms((double *)mof + 1, user->Nspec, maf + 1);
347: PetscCall(VecRestoreArrayRead(molef, &mof));
348: PetscCall(VecRestoreArray(*massf, &maf));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: PetscErrorCode ComputeMassConservation(Vec x, PetscReal *mass, PetscCtx ctx)
353: {
354: PetscFunctionBeginUser;
355: PetscCall(VecSum(x, mass));
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: PetscErrorCode MonitorMassConservation(TS ts, PetscInt step, PetscReal time, Vec x, PetscCtx ctx)
360: {
361: const PetscScalar *T;
362: PetscReal mass;
364: PetscFunctionBeginUser;
365: PetscCall(ComputeMassConservation(x, &mass, ctx));
366: PetscCall(VecGetArrayRead(x, &T));
367: mass -= PetscAbsScalar(T[0]);
368: PetscCall(VecRestoreArrayRead(x, &T));
369: 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))));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: PetscErrorCode MonitorTempature(TS ts, PetscInt step, PetscReal time, Vec x, PetscCtx ctx)
374: {
375: User user = (User)ctx;
376: const PetscScalar *T;
378: PetscFunctionBeginUser;
379: PetscCall(VecGetArrayRead(x, &T));
380: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep %" PetscInt_FMT " time %g temperature %g\n", step, (double)time, (double)(T[0] * user->Tini)));
381: PetscCall(VecRestoreArrayRead(x, &T));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*
386: Prints out each species with its name
387: */
388: PETSC_UNUSED PetscErrorCode PrintSpecies(User user, Vec molef)
389: {
390: const PetscScalar *mof;
391: PetscInt i, *idx, n = user->Nspec + 1;
393: PetscFunctionBeginUser;
394: PetscCall(PetscMalloc1(n, &idx));
395: for (i = 0; i < n; i++) idx[i] = i;
396: PetscCall(VecGetArrayRead(molef, &mof));
397: PetscCall(PetscSortRealWithPermutation(n, mof, idx));
398: 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]])));
399: PetscCall(PetscFree(idx));
400: PetscCall(VecRestoreArrayRead(molef, &mof));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*TEST
405: build:
406: requires: tchem
408: test:
409: 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
410: requires: !single
411: filter: grep -v iterations
413: TEST*/