1 static const char help[] = "Integrate chemistry using TChem.\n"; 2 3 #include <petscts.h> 4 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 19 20 See also h2_10sp.inp for another example 21 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 24 25 The solution for component i = 0 is the temperature. 26 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 28 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])) 32 33 FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species. 34 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 38 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 }; 51 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*); 61 62 #define PetscCallTC(ierr) do {PetscCheck(!ierr,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in TChem library, return code %d",ierr);} while (0) 63 64 int main(int argc,char **argv) 65 { 66 TS ts; /* time integrator */ 67 TSAdapt adapt; 68 Vec X,lambda; /* solution vector */ 69 Mat J; /* Jacobian matrix */ 70 PetscInt steps; 71 PetscReal ftime,dt; 72 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]; 73 const char *periodic = "file://${PETSC_DIR}/${PETSC_ARCH}/share/periodictable.dat"; 74 struct _User user; /* user-defined work context */ 75 TSConvergedReason reason; 76 char **snames,*names; 77 PetscInt i; 78 TSTrajectory tj; 79 PetscBool flg = PETSC_FALSE,tflg = PETSC_FALSE,found; 80 81 PetscFunctionBeginUser; 82 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 83 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Chemistry solver options",""); 84 PetscCall(PetscOptionsString("-chem","CHEMKIN input file","",chemfile,chemfile,sizeof(chemfile),NULL)); 85 PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD,chemfile,lchemfile,PETSC_MAX_PATH_LEN,&found)); 86 PetscCheck(found,PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot download %s and no local version %s",chemfile,lchemfile); 87 PetscCall(PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL)); 88 PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD,thermofile,lthermofile,PETSC_MAX_PATH_LEN,&found)); 89 PetscCheck(found,PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot download %s and no local version %s",thermofile,lthermofile); 90 user.pressure = 1.01325e5; /* Pascal */ 91 PetscCall(PetscOptionsReal("-pressure","Pressure of reaction [Pa]","",user.pressure,&user.pressure,NULL)); 92 user.Tini = 1000; /* Kelvin */ 93 PetscCall(PetscOptionsReal("-Tini","Initial temperature [K]","",user.Tini,&user.Tini,NULL)); 94 PetscCall(PetscOptionsBool("-monitor_mass","Monitor the total mass at each timestep","",flg,&flg,NULL)); 95 PetscCall(PetscOptionsBool("-monitor_temp","Monitor the temperature each timestep","",tflg,&tflg,NULL)); 96 PetscOptionsEnd(); 97 98 /* tchem requires periodic table in current directory */ 99 PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD,periodic,lperiodic,PETSC_MAX_PATH_LEN,&found)); 100 PetscCheck(found,PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot located required periodic table %s or local version %s",periodic,lperiodic); 101 102 PetscCallTC(TC_initChem(lchemfile, lthermofile, 0, 1.0)); 103 TC_setThermoPres(user.pressure); 104 user.Nspec = TC_getNspec(); 105 user.Nreac = TC_getNreac(); 106 /* 107 Get names of all species in easy to use array 108 */ 109 PetscCall(PetscMalloc1((user.Nspec+1)*LENGTHOFSPECNAME,&names)); 110 PetscCall(PetscStrcpy(names,"Temp")); 111 TC_getSnames(user.Nspec,names+LENGTHOFSPECNAME); 112 PetscCall(PetscMalloc1((user.Nspec+2),&snames)); 113 for (i=0; i<user.Nspec+1; i++) snames[i] = names+i*LENGTHOFSPECNAME; 114 snames[user.Nspec+1] = NULL; 115 PetscCall(PetscStrArrayallocpy((const char *const *)snames,&user.snames)); 116 PetscCall(PetscFree(snames)); 117 PetscCall(PetscFree(names)); 118 119 PetscCall(PetscMalloc3(user.Nspec+1,&user.tchemwork,PetscSqr(user.Nspec+1),&user.Jdense,user.Nspec+1,&user.rows)); 120 PetscCall(VecCreateSeq(PETSC_COMM_SELF,user.Nspec+1,&X)); 121 122 /* PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,PETSC_DECIDE,NULL,&J)); */ 123 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,NULL,&J)); 124 PetscCall(MatSetFromOptions(J)); 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 Create timestepping solver context 128 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 129 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 130 PetscCall(TSSetType(ts,TSARKIMEX)); 131 PetscCall(TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE)); 132 PetscCall(TSARKIMEXSetType(ts,TSARKIMEX4)); 133 PetscCall(TSSetRHSFunction(ts,NULL,FormRHSFunction,&user)); 134 PetscCall(TSSetRHSJacobian(ts,J,J,FormRHSJacobian,&user)); 135 136 if (flg) PetscCall(TSMonitorSet(ts,MonitorMassConservation,NULL,NULL)); 137 if (tflg) { 138 PetscCall(TSMonitorSet(ts,MonitorTempature,&user,NULL)); 139 } 140 141 ftime = 1.0; 142 PetscCall(TSSetMaxTime(ts,ftime)); 143 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 144 145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 Set initial conditions 147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 148 PetscCall(FormInitialSolution(ts,X,&user)); 149 PetscCall(TSSetSolution(ts,X)); 150 dt = 1e-10; /* Initial time step */ 151 PetscCall(TSSetTimeStep(ts,dt)); 152 PetscCall(TSGetAdapt(ts,&adapt)); 153 PetscCall(TSAdaptSetStepLimits(adapt,1e-12,1e-4)); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */ 154 PetscCall(TSSetMaxSNESFailures(ts,-1)); /* Retry step an unlimited number of times */ 155 156 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 157 Set runtime options 158 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 159 PetscCall(TSSetFromOptions(ts)); 160 161 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 162 Set final conditions for sensitivities 163 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 164 PetscCall(VecDuplicate(X,&lambda)); 165 PetscCall(TSSetCostGradients(ts,1,&lambda,NULL)); 166 PetscCall(VecSetValue(lambda,0,1.0,INSERT_VALUES)); 167 PetscCall(VecAssemblyBegin(lambda)); 168 PetscCall(VecAssemblyEnd(lambda)); 169 170 PetscCall(TSGetTrajectory(ts,&tj)); 171 if (tj) { 172 PetscCall(TSTrajectorySetVariableNames(tj,(const char * const *)user.snames)); 173 PetscCall(TSTrajectorySetTransform(tj,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user)); 174 } 175 176 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 177 Pass information to graphical monitoring routine 178 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 179 PetscCall(TSMonitorLGSetVariableNames(ts,(const char * const *)user.snames)); 180 PetscCall(TSMonitorLGSetTransform(ts,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user)); 181 182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 183 Solve ODE 184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 185 PetscCall(TSSolve(ts,X)); 186 PetscCall(TSGetSolveTime(ts,&ftime)); 187 PetscCall(TSGetStepNumber(ts,&steps)); 188 PetscCall(TSGetConvergedReason(ts,&reason)); 189 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %" PetscInt_FMT " steps\n",TSConvergedReasons[reason],(double)ftime,steps)); 190 191 /* { 192 Vec max; 193 PetscInt i; 194 const PetscReal *bmax; 195 196 PetscCall(TSMonitorEnvelopeGetBounds(ts,&max,NULL)); 197 if (max) { 198 PetscCall(VecGetArrayRead(max,&bmax)); 199 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n")); 200 for (i=1; i<user.Nspec; i++) { 201 if (bmax[i] > .01) PetscCall(PetscPrintf(PETSC_COMM_SELF,"%s %g\n",user.snames[i],(double)bmax[i])); 202 } 203 PetscCall(VecRestoreArrayRead(max,&bmax)); 204 } 205 } 206 207 Vec y; 208 PetscCall(MassFractionToMoleFraction(&user,X,&y)); 209 PetscCall(PrintSpecies(&user,y)); 210 PetscCall(VecDestroy(&y)); */ 211 212 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 213 Free work space. 214 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 215 TC_reset(); 216 PetscCall(PetscStrArrayDestroy(&user.snames)); 217 PetscCall(MatDestroy(&J)); 218 PetscCall(VecDestroy(&X)); 219 PetscCall(VecDestroy(&lambda)); 220 PetscCall(TSDestroy(&ts)); 221 PetscCall(PetscFree3(user.tchemwork,user.Jdense,user.rows)); 222 PetscCall(PetscFinalize()); 223 return 0; 224 } 225 226 static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 227 { 228 User user = (User)ptr; 229 PetscScalar *f; 230 const PetscScalar *x; 231 232 PetscFunctionBeginUser; 233 PetscCall(VecGetArrayRead(X,&x)); 234 PetscCall(VecGetArray(F,&f)); 235 236 PetscCall(PetscArraycpy(user->tchemwork,x,user->Nspec+1)); 237 user->tchemwork[0] *= user->Tini; /* Dimensionalize */ 238 PetscCallTC(TC_getSrc(user->tchemwork,user->Nspec+1,f)); 239 f[0] /= user->Tini; /* Non-dimensionalize */ 240 241 PetscCall(VecRestoreArrayRead(X,&x)); 242 PetscCall(VecRestoreArray(F,&f)); 243 PetscFunctionReturn(0); 244 } 245 246 static PetscErrorCode FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr) 247 { 248 User user = (User)ptr; 249 const PetscScalar *x; 250 PetscInt M = user->Nspec+1,i; 251 252 PetscFunctionBeginUser; 253 PetscCall(VecGetArrayRead(X,&x)); 254 PetscCall(PetscArraycpy(user->tchemwork,x,user->Nspec+1)); 255 PetscCall(VecRestoreArrayRead(X,&x)); 256 user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */ 257 PetscCall(TC_getJacTYN(user->tchemwork,user->Nspec,user->Jdense,1)); 258 259 for (i=0; i<M; i++) user->Jdense[i + 0*M] /= user->Tini; /* Non-dimensionalize first column */ 260 for (i=0; i<M; i++) user->Jdense[0 + i*M] /= user->Tini; /* Non-dimensionalize first row */ 261 for (i=0; i<M; i++) user->rows[i] = i; 262 PetscCall(MatSetOption(Pmat,MAT_ROW_ORIENTED,PETSC_FALSE)); 263 PetscCall(MatSetOption(Pmat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE)); 264 PetscCall(MatZeroEntries(Pmat)); 265 PetscCall(MatSetValues(Pmat,M,user->rows,M,user->rows,user->Jdense,INSERT_VALUES)); 266 267 PetscCall(MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY)); 268 PetscCall(MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY)); 269 if (Amat != Pmat) { 270 PetscCall(MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY)); 271 PetscCall(MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY)); 272 } 273 PetscFunctionReturn(0); 274 } 275 276 PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx) 277 { 278 PetscScalar *x; 279 PetscInt i; 280 Vec y; 281 const PetscInt maxspecies = 10; 282 PetscInt smax = maxspecies,mmax = maxspecies; 283 char *names[maxspecies]; 284 PetscReal molefracs[maxspecies],sum; 285 PetscBool flg; 286 287 PetscFunctionBeginUser; 288 PetscCall(VecZeroEntries(X)); 289 PetscCall(VecGetArray(X,&x)); 290 x[0] = 1.0; /* Non-dimensionalized by user->Tini */ 291 292 PetscCall(PetscOptionsGetStringArray(NULL,NULL,"-initial_species",names,&smax,&flg)); 293 PetscCheck(smax >= 2,PETSC_COMM_SELF,PETSC_ERR_USER,"Must provide at least two initial species"); 294 PetscCall(PetscOptionsGetRealArray(NULL,NULL,"-initial_mole",molefracs,&mmax,&flg)); 295 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); 296 sum = 0; 297 for (i=0; i<smax; i++) sum += molefracs[i]; 298 for (i=0; i<smax; i++) molefracs[i] = molefracs[i]/sum; 299 for (i=0; i<smax; i++) { 300 int ispec = TC_getSpos(names[i], strlen(names[i])); 301 PetscCheck(ispec >= 0,PETSC_COMM_SELF,PETSC_ERR_USER,"Could not find species %s",names[i]); 302 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Species %" PetscInt_FMT ": %s %g\n",i,names[i],(double)molefracs[i])); 303 x[1+ispec] = molefracs[i]; 304 } 305 for (i=0; i<smax; i++) { 306 PetscCall(PetscFree(names[i])); 307 } 308 PetscCall(VecRestoreArray(X,&x)); 309 /* PetscCall(PrintSpecies((User)ctx,X)); */ 310 PetscCall(MoleFractionToMassFraction((User)ctx,X,&y)); 311 PetscCall(VecCopy(y,X)); 312 PetscCall(VecDestroy(&y)); 313 PetscFunctionReturn(0); 314 } 315 316 /* 317 Converts the input vector which is in mass fractions (used by tchem) to mole fractions 318 */ 319 PetscErrorCode MassFractionToMoleFraction(User user,Vec massf,Vec *molef) 320 { 321 PetscScalar *mof; 322 const PetscScalar *maf; 323 324 PetscFunctionBeginUser; 325 PetscCall(VecDuplicate(massf,molef)); 326 PetscCall(VecGetArrayRead(massf,&maf)); 327 PetscCall(VecGetArray(*molef,&mof)); 328 mof[0] = maf[0]; /* copy over temperature */ 329 TC_getMs2Ml((double*)maf+1,user->Nspec,mof+1); 330 PetscCall(VecRestoreArray(*molef,&mof)); 331 PetscCall(VecRestoreArrayRead(massf,&maf)); 332 PetscFunctionReturn(0); 333 } 334 335 /* 336 Converts the input vector which is in mole fractions to mass fractions (used by tchem) 337 */ 338 PetscErrorCode MoleFractionToMassFraction(User user,Vec molef,Vec *massf) 339 { 340 const PetscScalar *mof; 341 PetscScalar *maf; 342 343 PetscFunctionBeginUser; 344 PetscCall(VecDuplicate(molef,massf)); 345 PetscCall(VecGetArrayRead(molef,&mof)); 346 PetscCall(VecGetArray(*massf,&maf)); 347 maf[0] = mof[0]; /* copy over temperature */ 348 TC_getMl2Ms((double*)mof+1,user->Nspec,maf+1); 349 PetscCall(VecRestoreArrayRead(molef,&mof)); 350 PetscCall(VecRestoreArray(*massf,&maf)); 351 PetscFunctionReturn(0); 352 } 353 354 PetscErrorCode ComputeMassConservation(Vec x,PetscReal *mass,void* ctx) 355 { 356 PetscFunctionBeginUser; 357 PetscCall(VecSum(x,mass)); 358 PetscFunctionReturn(0); 359 } 360 361 PetscErrorCode MonitorMassConservation(TS ts,PetscInt step,PetscReal time,Vec x,void* ctx) 362 { 363 const PetscScalar *T; 364 PetscReal mass; 365 366 PetscFunctionBeginUser; 367 PetscCall(ComputeMassConservation(x,&mass,ctx)); 368 PetscCall(VecGetArrayRead(x,&T)); 369 mass -= PetscAbsScalar(T[0]); 370 PetscCall(VecRestoreArrayRead(x,&T)); 371 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)))); 372 PetscFunctionReturn(0); 373 } 374 375 PetscErrorCode MonitorTempature(TS ts,PetscInt step,PetscReal time,Vec x,void* ctx) 376 { 377 User user = (User) ctx; 378 const PetscScalar *T; 379 380 PetscFunctionBeginUser; 381 PetscCall(VecGetArrayRead(x,&T)); 382 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Timestep %" PetscInt_FMT " time %g temperature %g\n",step,(double)time,(double)(T[0]*user->Tini))); 383 PetscCall(VecRestoreArrayRead(x,&T)); 384 PetscFunctionReturn(0); 385 } 386 387 /* 388 Prints out each species with its name 389 */ 390 PETSC_UNUSED PetscErrorCode PrintSpecies(User user,Vec molef) 391 { 392 const PetscScalar *mof; 393 PetscInt i,*idx,n = user->Nspec+1; 394 395 PetscFunctionBeginUser; 396 PetscCall(PetscMalloc1(n,&idx)); 397 for (i=0; i<n;i++) idx[i] = i; 398 PetscCall(VecGetArrayRead(molef,&mof)); 399 PetscCall(PetscSortRealWithPermutation(n,mof,idx)); 400 for (i=0; i<n; i++) { 401 PetscCall(PetscPrintf(PETSC_COMM_SELF,"%6s %g\n",user->snames[idx[n-i-1]],(double)PetscRealPart(mof[idx[n-i-1]]))); 402 } 403 PetscCall(PetscFree(idx)); 404 PetscCall(VecRestoreArrayRead(molef,&mof)); 405 PetscFunctionReturn(0); 406 } 407 408 /*TEST 409 build: 410 requires: tchem 411 412 test: 413 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 414 requires: !single 415 filter: grep -v iterations 416 417 TEST*/ 418