Actual source code: extchem.c
petsc-3.7.5 2017-01-01
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: }