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, (char *)0, 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 -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*/