1 static const char help[] = "Integrate chemistry using TChem.\n";
2
3 #include <petscts.h>
4
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
19
20 See also h2_10sp.inp for another example
21
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
24
25 The solution for component i = 0 is the temperature.
26
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
28
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]))
32
33 FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species.
34
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
38
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 };
51
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 *);
61
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)
66
main(int argc,char ** argv)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;
83
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();
100
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);
104
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));
121
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));
124
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));
128
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));
138
139 if (flg) PetscCall(TSMonitorSet(ts, MonitorMassConservation, NULL, NULL));
140 if (tflg) PetscCall(TSMonitorSet(ts, MonitorTempature, &user, NULL));
141
142 ftime = 1.0;
143 PetscCall(TSSetMaxTime(ts, ftime));
144 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
145
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 */
156
157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158 Set runtime options
159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160 PetscCall(TSSetFromOptions(ts));
161
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));
170
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 }
176
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));
182
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));
191
192 /* {
193 Vec max;
194 PetscInt i;
195 const PetscReal *bmax;
196
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 }
207
208 Vec y;
209 PetscCall(MassFractionToMoleFraction(&user,X,&y));
210 PetscCall(PrintSpecies(&user,y));
211 PetscCall(VecDestroy(&y)); */
212
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 }
226
FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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;
232
233 PetscFunctionBeginUser;
234 PetscCall(VecGetArrayRead(X, &x));
235 PetscCall(VecGetArray(F, &f));
236
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 */
241
242 PetscCall(VecRestoreArrayRead(X, &x));
243 PetscCall(VecRestoreArray(F, &f));
244 PetscFunctionReturn(PETSC_SUCCESS);
245 }
246
FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void * ptr)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;
252
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));
259
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));
267
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 }
276
FormInitialSolution(TS ts,Vec X,PetscCtx ctx)277 PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx 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;
287
288 PetscFunctionBeginUser;
289 PetscCall(VecZeroEntries(X));
290 PetscCall(VecGetArray(X, &x));
291 x[0] = 1.0; /* Non-dimensionalized by user->Tini */
292
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], (int)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 }
314
315 /*
316 Converts the input vector which is in mass fractions (used by tchem) to mole fractions
317 */
MassFractionToMoleFraction(User user,Vec massf,Vec * molef)318 PetscErrorCode MassFractionToMoleFraction(User user, Vec massf, Vec *molef)
319 {
320 PetscScalar *mof;
321 const PetscScalar *maf;
322
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 }
333
334 /*
335 Converts the input vector which is in mole fractions to mass fractions (used by tchem)
336 */
MoleFractionToMassFraction(User user,Vec molef,Vec * massf)337 PetscErrorCode MoleFractionToMassFraction(User user, Vec molef, Vec *massf)
338 {
339 const PetscScalar *mof;
340 PetscScalar *maf;
341
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 }
352
ComputeMassConservation(Vec x,PetscReal * mass,PetscCtx ctx)353 PetscErrorCode ComputeMassConservation(Vec x, PetscReal *mass, PetscCtx ctx)
354 {
355 PetscFunctionBeginUser;
356 PetscCall(VecSum(x, mass));
357 PetscFunctionReturn(PETSC_SUCCESS);
358 }
359
MonitorMassConservation(TS ts,PetscInt step,PetscReal time,Vec x,PetscCtx ctx)360 PetscErrorCode MonitorMassConservation(TS ts, PetscInt step, PetscReal time, Vec x, PetscCtx ctx)
361 {
362 const PetscScalar *T;
363 PetscReal mass;
364
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 }
373
MonitorTempature(TS ts,PetscInt step,PetscReal time,Vec x,PetscCtx ctx)374 PetscErrorCode MonitorTempature(TS ts, PetscInt step, PetscReal time, Vec x, PetscCtx ctx)
375 {
376 User user = (User)ctx;
377 const PetscScalar *T;
378
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 }
385
386 /*
387 Prints out each species with its name
388 */
PrintSpecies(User user,Vec molef)389 PETSC_UNUSED PetscErrorCode PrintSpecies(User user, Vec molef)
390 {
391 const PetscScalar *mof;
392 PetscInt i, *idx, n = user->Nspec + 1;
393
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 }
404
405 /*TEST
406 build:
407 requires: tchem
408
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
413
414 TEST*/
415