Actual source code: extchem.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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 tempature 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.


 36:     These are other data sets for other possible runs
 37:        https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/n_heptane_v3.1_therm.dat
 38:        https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/nc7_ver3.1_mech.txt


 41: */
 42: typedef struct _User *User;
 43: struct _User {
 44:   PetscReal pressure;
 45:   int       Nspec;
 46:   int       Nreac;
 47:   PetscReal Tini;
 48:   double    *tchemwork;
 49:   double    *Jdense;        /* Dense array workspace where Tchem computes the Jacobian */ 
 50:   PetscInt  *rows;
 51:   char      **snames;
 52: };


 55: static PetscErrorCode PrintSpecies(User,Vec);
 56: static PetscErrorCode MassFractionToMoleFraction(User,Vec,Vec*);
 57: static PetscErrorCode MoleFractionToMassFraction(User,Vec,Vec*);
 58: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
 59: static PetscErrorCode FormRHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
 60: static PetscErrorCode FormInitialSolution(TS,Vec,void*);
 61: static PetscErrorCode ComputeMassConservation(Vec,PetscReal*,void*);
 62: static PetscErrorCode MonitorMassConservation(TS,PetscInt,PetscReal,Vec,void*);
 63: static PetscErrorCode MonitorTempature(TS,PetscInt,PetscReal,Vec,void*);

 65: #define TCCHKERRQ(ierr) do {if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in TChem library, return code %d",ierr);} 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:   PetscErrorCode    ierr;
 75:   PetscReal         ftime,dt;
 76:   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];
 77:   const char        *periodic = "file://${PETSC_DIR}/${PETSC_ARCH}/share/periodictable.dat";
 78:   struct _User      user;       /* user-defined work context */
 79:   TSConvergedReason reason;
 80:   char              **snames,*names;
 81:   PetscInt          i;
 82:   TSTrajectory      tj;
 83:   PetscBool         flg = PETSC_FALSE,tflg = PETSC_FALSE,found;

 85:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 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);
 89:   if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot download %s and no local version %s",chemfile,lchemfile);
 90:   PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL);
 91:   PetscFileRetrieve(PETSC_COMM_WORLD,thermofile,lthermofile,PETSC_MAX_PATH_LEN,&found);
 92:   if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot download %s and no local version %s",thermofile,lthermofile);
 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 tempature 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);
103:   if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot located required periodic table %s or local version %s",periodic,lperiodic);

105:   TC_initChem(lchemfile, lthermofile, 0, 1.0);TC
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) {
140:     TSMonitorSet(ts,MonitorMassConservation,NULL,NULL);
141:   }
142:   if (tflg) {
143:     TSMonitorSet(ts,MonitorTempature,&user,NULL);
144:   }

146:   ftime = 1.0;
147:   TSSetMaxTime(ts,ftime);
148:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:      Set initial conditions
152:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153:   FormInitialSolution(ts,X,&user);
154:   TSSetSolution(ts,X);
155:   dt   = 1e-10;                 /* Initial time step */
156:   TSSetTimeStep(ts,dt);
157:   TSGetAdapt(ts,&adapt);
158:   TSAdaptSetStepLimits(adapt,1e-12,1e-4); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
159:   TSSetMaxSNESFailures(ts,-1);            /* Retry step an unlimited number of times */

161:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162:      Set runtime options
163:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164:   TSSetFromOptions(ts);

166:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167:      Set final conditions for sensitivities
168:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169:   VecDuplicate(X,&lambda);
170:   TSSetCostGradients(ts,1,&lambda,NULL);
171:   VecSetValue(lambda,0,1.0,INSERT_VALUES);
172:   VecAssemblyBegin(lambda);
173:   VecAssemblyEnd(lambda);

175:   TSGetTrajectory(ts,&tj);
176:   if (tj) {
177:     TSTrajectorySetVariableNames(tj,(const char * const *)user.snames);
178:     TSTrajectorySetTransform(tj,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user);
179:   }


182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:      Pass information to graphical monitoring routine
184:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185:   TSMonitorLGSetVariableNames(ts,(const char * const *)user.snames);
186:   TSMonitorLGSetTransform(ts,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user);

188:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189:      Solve ODE
190:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191:   TSSolve(ts,X);
192:   TSGetSolveTime(ts,&ftime);
193:   TSGetStepNumber(ts,&steps);
194:   TSGetConvergedReason(ts,&reason);
195:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);

197:   /* {
198:     Vec                max;
199:     PetscInt           i;
200:     const PetscReal    *bmax;

202:     TSMonitorEnvelopeGetBounds(ts,&max,NULL);
203:     if (max) {
204:       VecGetArrayRead(max,&bmax);
205:       PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n");
206:       for (i=1; i<user.Nspec; i++) {
207:         if (bmax[i] > .01) {PetscPrintf(PETSC_COMM_SELF,"%s %g\n",user.snames[i],(double)bmax[i]);}
208:       }
209:       VecRestoreArrayRead(max,&bmax);
210:     }
211:   }

213:   Vec y;
214:   MassFractionToMoleFraction(&user,X,&y);
215:   PrintSpecies(&user,y);
216:   VecDestroy(&y); */

218:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219:      Free work space.
220:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221:   TC_reset();
222:   PetscStrArrayDestroy(&user.snames);
223:   MatDestroy(&J);
224:   VecDestroy(&X);
225:   VecDestroy(&lambda);
226:   TSDestroy(&ts);
227:   PetscFree3(user.tchemwork,user.Jdense,user.rows);
228:   PetscFinalize();
229:   return ierr;
230: }

232: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
233: {
234:   User              user = (User)ptr;
235:   PetscErrorCode    ierr;
236:   PetscScalar       *f;
237:   const PetscScalar *x;

240:   VecGetArrayRead(X,&x);
241:   VecGetArray(F,&f);

243:   PetscArraycpy(user->tchemwork,x,user->Nspec+1);
244:   user->tchemwork[0] *= user->Tini; /* Dimensionalize */
245:   TC_getSrc(user->tchemwork,user->Nspec+1,f);TC
246:   f[0] /= user->Tini;           /* Non-dimensionalize */

248:   VecRestoreArrayRead(X,&x);
249:   VecRestoreArray(F,&f);
250:   return(0);
251: }

253: static PetscErrorCode FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr)
254: {
255:   User              user = (User)ptr;
256:   PetscErrorCode    ierr;
257:   const PetscScalar *x;
258:   PetscInt          M = user->Nspec+1,i;

261:   VecGetArrayRead(X,&x);
262:   PetscArraycpy(user->tchemwork,x,user->Nspec+1);
263:   VecRestoreArrayRead(X,&x);
264:   user->tchemwork[0] *= user->Tini;  /* Dimensionalize temperature (first row) because that is what Tchem wants */
265:   TC_getJacTYN(user->tchemwork,user->Nspec,user->Jdense,1);

267:   for (i=0; i<M; i++) user->Jdense[i + 0*M] /= user->Tini; /* Non-dimensionalize first column */
268:   for (i=0; i<M; i++) user->Jdense[0 + i*M] /= user->Tini; /* Non-dimensionalize first row */
269:   for (i=0; i<M; i++) user->rows[i] = i;
270:   MatSetOption(Pmat,MAT_ROW_ORIENTED,PETSC_FALSE);
271:   MatSetOption(Pmat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
272:   MatZeroEntries(Pmat);
273:   MatSetValues(Pmat,M,user->rows,M,user->rows,user->Jdense,INSERT_VALUES);

275:   MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
276:   MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
277:   if (Amat != Pmat) {
278:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
279:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
280:   }
281:   return(0);
282: }

284: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
285: {
286:   PetscScalar    *x;
288:   PetscInt       i;
289:   Vec            y;
290:   const PetscInt maxspecies = 10;
291:   PetscInt       smax = maxspecies,mmax = maxspecies;
292:   char           *names[maxspecies];
293:   PetscReal      molefracs[maxspecies],sum;
294:   PetscBool      flg;

297:   VecZeroEntries(X);
298:   VecGetArray(X,&x);
299:   x[0] = 1.0;  /* Non-dimensionalized by user->Tini */

301:   PetscOptionsGetStringArray(NULL,NULL,"-initial_species",names,&smax,&flg);
302:   if (smax < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Must provide at least two initial species");
303:   PetscOptionsGetRealArray(NULL,NULL,"-initial_mole",molefracs,&mmax,&flg);
304:   if (smax != mmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Must provide same number of initial species %D as initial moles %D",smax,mmax);
305:   sum = 0;
306:   for (i=0; i<smax; i++) sum += molefracs[i];
307:   for (i=0; i<smax; i++) molefracs[i] = molefracs[i]/sum;
308:   for (i=0; i<smax; i++) {
309:     int ispec = TC_getSpos(names[i], strlen(names[i]));
310:     if (ispec < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Could not find species %s",names[i]);
311:     PetscPrintf(PETSC_COMM_SELF,"Species %d: %s %g\n",i,names[i],molefracs[i]);
312:     x[1+ispec] = molefracs[i];
313:   }
314:   for (i=0; i<smax; i++) {
315:     PetscFree(names[i]);
316:   }
317:   VecRestoreArray(X,&x);
318:   /* PrintSpecies((User)ctx,X); */
319:   MoleFractionToMassFraction((User)ctx,X,&y);
320:   VecCopy(y,X);
321:   VecDestroy(&y);
322:   return(0);
323: }

325: /*
326:    Converts the input vector which is in mass fractions (used by tchem) to mole fractions
327: */
328: PetscErrorCode MassFractionToMoleFraction(User user,Vec massf,Vec *molef)
329: {
330:   PetscErrorCode    ierr;
331:   PetscScalar       *mof;
332:   const PetscScalar *maf;

335:   VecDuplicate(massf,molef);
336:   VecGetArrayRead(massf,&maf);
337:   VecGetArray(*molef,&mof);
338:   mof[0] = maf[0]; /* copy over temperature */
339:   TC_getMs2Ml((double*)maf+1,user->Nspec,mof+1);
340:   VecRestoreArray(*molef,&mof);
341:   VecRestoreArrayRead(massf,&maf);
342:   return(0);
343: }

345: /*
346:    Converts the input vector which is in mole fractions to mass fractions (used by tchem)
347: */
348: PetscErrorCode MoleFractionToMassFraction(User user,Vec molef,Vec *massf)
349: {
350:   PetscErrorCode    ierr;
351:   const PetscScalar *mof;
352:   PetscScalar       *maf;

355:   VecDuplicate(molef,massf);
356:   VecGetArrayRead(molef,&mof);
357:   VecGetArray(*massf,&maf);
358:   maf[0] = mof[0]; /* copy over temperature */
359:   TC_getMl2Ms((double*)mof+1,user->Nspec,maf+1);
360:   VecRestoreArrayRead(molef,&mof);
361:   VecRestoreArray(*massf,&maf);
362:   return(0);
363: }

365: PetscErrorCode ComputeMassConservation(Vec x,PetscReal *mass,void* ctx)
366: {

370:   VecSum(x,mass);
371:   return(0);
372: }

374: PetscErrorCode MonitorMassConservation(TS ts,PetscInt step,PetscReal time,Vec x,void* ctx)
375: {
376:   const PetscScalar  *T;
377:   PetscReal          mass;
378:   PetscErrorCode     ierr;

381:   ComputeMassConservation(x,&mass,ctx);
382:   VecGetArrayRead(x,&T);
383:   mass -= PetscAbsScalar(T[0]);
384:   VecRestoreArrayRead(x,&T);
385:   PetscPrintf(PETSC_COMM_WORLD,"Timestep %D time %g percent mass lost or gained %g\n",step,(double)time,(double)100.*(1.0 - mass));
386:   return(0);
387: }

389: PetscErrorCode MonitorTempature(TS ts,PetscInt step,PetscReal time,Vec x,void* ctx)
390: {
391:   User               user = (User) ctx;
392:   const PetscScalar  *T;
393:   PetscErrorCode     ierr;

396:   VecGetArrayRead(x,&T);
397:   PetscPrintf(PETSC_COMM_WORLD,"Timestep %D time %g tempature %g\n",step,(double)time,(double)T[0]*user->Tini);
398:   VecRestoreArrayRead(x,&T);
399:   return(0);
400: }

402: /*
403:    Prints out each species with its name
404: */
405: PETSC_UNUSED PetscErrorCode PrintSpecies(User user,Vec molef)
406: {
407:   PetscErrorCode    ierr;
408:   const PetscScalar *mof;
409:   PetscInt          i,*idx,n = user->Nspec+1;

412:   PetscMalloc1(n,&idx);
413:   for (i=0; i<n;i++) idx[i] = i;
414:   VecGetArrayRead(molef,&mof);
415:   PetscSortRealWithPermutation(n,mof,idx);
416:   for (i=0; i<n; i++) {
417:     PetscPrintf(PETSC_COMM_SELF,"%6s %g\n",user->snames[idx[n-i-1]],mof[idx[n-i-1]]);
418:   }
419:   PetscFree(idx);
420:   VecRestoreArrayRead(molef,&mof);
421:   return(0);
422: }

424: /*TEST
425:     build:
426:       requires: tchem

428:     test:
429:       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
430:       requires: !single
431:       filter: grep -v iterations

433: TEST*/