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*/