Actual source code: extchem.c

petsc-3.7.3 2016-08-01
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:     Obtain the three files into this directory

 20:        curl http://combustion.berkeley.edu/gri_mech/version30/files30/grimech30.dat > chem.inp
 21:        curl http://combustion.berkeley.edu/gri_mech/version30/files30/thermo30.dat > therm.dat

 23:        https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/n_heptane_v3.1_therm.dat
 24:        https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/nc7_ver3.1_mech.txt

 26:        cp $PETSC_DIR/$PETSC_ARCH/externalpackages/tchem/data/periodictable.dat .

 28:     Run with
 29:    ./extchem -Tini 1500 -ts_arkimex_fully_implicit -ts_max_snes_failures -1 -ts_adapt_monitor -ts_adapt_dt_max 1e-4 -ts_arkimex_type 4 -ts_monitor_lg_solution -ts_final_time .005 -draw_pause -2 -lg_use_markers false -ts_monitor_lg_solution_variables H2,O2,H2O,CH4,CO,CO2,C2H2,N2  -ts_monitor_envelope

 31:     Determine sensitivity of final tempature on each variables initial conditions
 32:     -ts_dt 1.e-5 -ts_type cn -ts_adjoint_solve -ts_adjoint_view_solution draw

 34:     The solution for component i = 0 is the temperature.

 36:     The solution, i > 0, is the mass fraction, massf[i], of species i, i.e. mass of species i/ total mass of all species

 38:     The mole fraction molef[i], i > 0, is the number of moles of a species/ total number of moles of all species
 39:         Define M[i] = mass per mole of species i then
 40:         molef[i] = massf[i]/(M[i]*(sum_j massf[j]/M[j]))

 42:     FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species.

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


 58: static PetscErrorCode PrintSpecies(User,Vec);
 59: static PetscErrorCode MassFractionToMoleFraction(User,Vec,Vec*);
 60: static PetscErrorCode MoleFractionToMassFraction(User,Vec,Vec*);
 61: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
 62: static PetscErrorCode FormRHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
 63: static PetscErrorCode FormInitialSolution(TS,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)

 69: int main(int argc,char **argv)
 70: {
 71:   TS                ts;         /* time integrator */
 72:   TSAdapt           adapt;
 73:   Vec               X,lambda;          /* solution vector */
 74:   Mat               J;          /* Jacobian matrix */
 75:   PetscInt          steps,maxsteps;
 76:   PetscErrorCode    ierr;
 77:   PetscReal         ftime,dt;
 78:   char              chemfile[PETSC_MAX_PATH_LEN] = "chem.inp",thermofile[PETSC_MAX_PATH_LEN] = "therm.dat";
 79:   struct _User      user;       /* user-defined work context */
 80:   TSConvergedReason reason;
 81:   char              **snames,*names;
 82:   PetscInt          i;

 84:   PetscInitialize(&argc,&argv,(char*)0,help);
 85:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Chemistry solver options","");
 86:   PetscOptionsString("-chem","CHEMKIN input file","",chemfile,chemfile,sizeof(chemfile),NULL);
 87:   PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL);
 88:   user.pressure = 1.01325e5;    /* Pascal */
 89:   PetscOptionsReal("-pressure","Pressure of reaction [Pa]","",user.pressure,&user.pressure,NULL);
 90:   user.Tini = 1000;             /* Kelvin */
 91:   PetscOptionsReal("-Tini","Initial temperature [K]","",user.Tini,&user.Tini,NULL);
 92:   PetscOptionsEnd();

 94:   TC_initChem(chemfile, thermofile, 0, 1.0);TC
 95:   user.Nspec = TC_getNspec();
 96:   user.Nreac = TC_getNreac();
 97:   /*
 98:       Get names of all species in easy to use array
 99:   */
100:   PetscMalloc1((user.Nspec+1)*LENGTHOFSPECNAME,&names);
101:   PetscStrcpy(names,"Temp");
102:   TC_getSnames(user.Nspec,names+LENGTHOFSPECNAME);
103:   PetscMalloc1((user.Nspec+2),&snames);
104:   for (i=0; i<user.Nspec+1; i++) snames[i] = names+i*LENGTHOFSPECNAME;
105:   snames[user.Nspec+1] = NULL;
106:   PetscStrArrayallocpy((const char *const *)snames,&user.snames);
107:   PetscFree(snames);
108:   PetscFree(names);

110:   PetscMalloc3(user.Nspec+1,&user.tchemwork,PetscSqr(user.Nspec+1),&user.Jdense,user.Nspec+1,&user.rows);
111:   VecCreateSeq(PETSC_COMM_SELF,user.Nspec+1,&X);

113:   MatCreateSeqAIJ(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,PETSC_DECIDE,NULL,&J);
114:   /*MatCreateSeqDense(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,NULL,&J);*/
115:   MatSetFromOptions(J);

117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:      Create timestepping solver context
119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120:   TSCreate(PETSC_COMM_WORLD,&ts);
121:   TSSetType(ts,TSARKIMEX);
122:   TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);
123:   TSARKIMEXSetType(ts,TSARKIMEX4);
124:   TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);
125:   TSSetRHSJacobian(ts,J,J,FormRHSJacobian,&user);

127:   ftime    = 1.0;
128:   maxsteps = 10000;
129:   TSSetDuration(ts,maxsteps,ftime);
130:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

132:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133:      Set initial conditions
134:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135:   FormInitialSolution(ts,X,&user);
136:   TSSetSolution(ts,X);
137:   dt   = 1e-10;                 /* Initial time step */
138:   TSSetInitialTimeStep(ts,0.0,dt);
139:   TSGetAdapt(ts,&adapt);
140:   TSAdaptSetStepLimits(adapt,1e-12,1e-4); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
141:   TSSetMaxSNESFailures(ts,-1);            /* Retry step an unlimited number of times */

143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:      Set runtime options
145:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146:   TSSetFromOptions(ts);

148:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149:      Set final conditions for sensitivities
150:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151:   VecDuplicate(X,&lambda);
152:   TSSetCostGradients(ts,1,&lambda,NULL);
153:   VecSetValue(lambda,0,1.0,INSERT_VALUES);
154:   VecAssemblyBegin(lambda);
155:   VecAssemblyEnd(lambda);



159:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160:      Pass information to graphical monitoring routine
161:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162:   TSMonitorLGSetVariableNames(ts,(const char * const *)user.snames);
163:   TSMonitorLGSetTransform(ts,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user);

165:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166:      Solve ODE
167:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168:   TSSolve(ts,X);
169:   TSGetSolveTime(ts,&ftime);
170:   TSGetTimeStepNumber(ts,&steps);
171:   TSGetConvergedReason(ts,&reason);
172:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);

174:   {
175:     Vec                max;
176:     const char * const *names;
177:     PetscInt           i;
178:     const PetscReal    *bmax;

180:     TSMonitorEnvelopeGetBounds(ts,&max,NULL);
181:     if (max) {
182:       TSMonitorLGGetVariableNames(ts,&names);
183:       VecGetArrayRead(max,&bmax);
184:       PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n");
185:       for (i=1; i<user.Nspec; i++) {
186:         if (bmax[i] > .01) {PetscPrintf(PETSC_COMM_SELF,"%s %g\n",names[i],bmax[i]);}
187:       }
188:       VecRestoreArrayRead(max,&bmax);
189:     }
190:   }

192:   Vec y;
193:   MassFractionToMoleFraction(&user,X,&y);
194:   PrintSpecies(&user,y);
195:   VecDestroy(&y);
196: 
197:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198:      Free work space.
199:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200:   TC_reset();
201:   PetscStrArrayDestroy(&user.snames);
202:   MatDestroy(&J);
203:   VecDestroy(&X);
204:   VecDestroy(&lambda);
205:   TSDestroy(&ts);
206:   PetscFree3(user.tchemwork,user.Jdense,user.rows);
207:   PetscFinalize();
208:   return 0;
209: }

213: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
214: {
215:   User              user = (User)ptr;
216:   PetscErrorCode    ierr;
217:   PetscScalar       *f;
218:   const PetscScalar *x;

221:   VecGetArrayRead(X,&x);
222:   VecGetArray(F,&f);

224:   PetscMemcpy(user->tchemwork,x,(user->Nspec+1)*sizeof(x[0]));
225:   user->tchemwork[0] *= user->Tini; /* Dimensionalize */
226:   TC_getSrc(user->tchemwork,user->Nspec+1,f);TC
227:   f[0] /= user->Tini;           /* Non-dimensionalize */

229:   VecRestoreArrayRead(X,&x);
230:   VecRestoreArray(F,&f);
231:   return(0);
232: }

236: static PetscErrorCode FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr)
237: {
238:   User              user = (User)ptr;
239:   PetscErrorCode    ierr;
240:   const PetscScalar *x;
241:   PetscInt          M = user->Nspec+1,i;

244:   VecGetArrayRead(X,&x);
245:   PetscMemcpy(user->tchemwork,x,(user->Nspec+1)*sizeof(x[0]));
246:   VecRestoreArrayRead(X,&x);
247:   user->tchemwork[0] *= user->Tini;  /* Dimensionalize temperature (first row) because that is what Tchem wants */
248:   TC_getJacTYN(user->tchemwork,user->Nspec,user->Jdense,1);

250:   for (i=0; i<M; i++) user->Jdense[i + 0*M] /= user->Tini; /* Non-dimensionalize first column */
251:   for (i=0; i<M; i++) user->Jdense[0 + i*M] /= user->Tini; /* Non-dimensionalize first row */
252:   for (i=0; i<M; i++) user->rows[i] = i;
253:   MatSetOption(Pmat,MAT_ROW_ORIENTED,PETSC_FALSE);
254:   MatSetOption(Pmat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
255:   MatZeroEntries(Pmat);
256:   MatSetValues(Pmat,M,user->rows,M,user->rows,user->Jdense,INSERT_VALUES);

258:   MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
259:   MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
260:   if (Amat != Pmat) {
261:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
262:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
263:   }
264:   return(0);
265: }

269: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
270: {
271:   PetscScalar    *x;
273:   struct {const char *name; PetscReal molefrac;} initial[] = {
274:     {"CH4", 0.0948178320887},
275:     {"O2", 0.189635664177},
276:     {"N2", 0.706766236705},
277:     {"AR", 0.00878026702874}
278:   };
279:   PetscInt i;
280:   Vec      y;

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

287:   for (i=0; i<sizeof(initial)/sizeof(initial[0]); i++) {
288:     int ispec = TC_getSpos(initial[i].name, strlen(initial[i].name));
289:     if (ispec < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Could not find species %s",initial[i].name);
290:     PetscPrintf(PETSC_COMM_SELF,"Species %d: %s %g\n",i,initial[i].name,initial[i].molefrac);
291:     x[1+ispec] = initial[i].molefrac;
292:   }
293:   VecRestoreArray(X,&x);
294:   PrintSpecies((User)ctx,X);
295:   MoleFractionToMassFraction((User)ctx,X,&y);
296:   VecCopy(y,X);
297:   VecDestroy(&y);
298:   return(0);
299: }

303: /*
304:    Converts the input vector which is in mass fractions (used by tchem) to mole fractions
305: */
306: PetscErrorCode MassFractionToMoleFraction(User user,Vec massf,Vec *molef)
307: {
308:   PetscErrorCode    ierr;
309:   PetscScalar       *mof;
310:   const PetscScalar *maf;

313:   VecDuplicate(massf,molef);
314:   VecGetArrayRead(massf,&maf);
315:   VecGetArray(*molef,&mof);
316:   mof[0] = maf[0]; /* copy over temperature */
317:   TC_getMs2Ml((double*)maf+1,user->Nspec,mof+1);
318:   VecRestoreArray(*molef,&mof);
319:   VecRestoreArrayRead(massf,&maf);
320:   return(0);
321: }

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

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

347: /*
348:    Prints out each species with its name
349: */
350: PetscErrorCode PrintSpecies(User user,Vec molef)
351: {
352:   PetscErrorCode    ierr;
353:   const PetscScalar *mof;
354:   PetscInt          i,*idx,n = user->Nspec+1;

357:   PetscMalloc1(n,&idx);
358:   for (i=0; i<n;i++) idx[i] = i;
359:   VecGetArrayRead(molef,&mof);
360:   PetscSortRealWithPermutation(n,mof,idx);
361:   for (i=0; i<n; i++) {
362:     PetscPrintf(PETSC_COMM_SELF,"%6s %g\n",user->snames[idx[n-i-1]],mof[idx[n-i-1]]);
363:   }
364:   PetscFree(idx);
365:   VecRestoreArrayRead(molef,&mof);
366:   return(0);
367: }