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) \ 63 do { \ 64 PetscCheck(!ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in TChem library, return code %d", ierr); \ 65 } while (0) 66 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; 83 84 PetscFunctionBeginUser; 85 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 86 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Chemistry solver options", ""); 87 PetscCall(PetscOptionsString("-chem", "CHEMKIN input file", "", chemfile, chemfile, sizeof(chemfile), NULL)); 88 PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD, chemfile, lchemfile, PETSC_MAX_PATH_LEN, &found)); 89 PetscCheck(found, PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "Cannot download %s and no local version %s", chemfile, lchemfile); 90 PetscCall(PetscOptionsString("-thermo", "NASA thermo input file", "", thermofile, thermofile, sizeof(thermofile), NULL)); 91 PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD, thermofile, lthermofile, PETSC_MAX_PATH_LEN, &found)); 92 PetscCheck(found, PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "Cannot download %s and no local version %s", thermofile, lthermofile); 93 user.pressure = 1.01325e5; /* Pascal */ 94 PetscCall(PetscOptionsReal("-pressure", "Pressure of reaction [Pa]", "", user.pressure, &user.pressure, NULL)); 95 user.Tini = 1000; /* Kelvin */ 96 PetscCall(PetscOptionsReal("-Tini", "Initial temperature [K]", "", user.Tini, &user.Tini, NULL)); 97 PetscCall(PetscOptionsBool("-monitor_mass", "Monitor the total mass at each timestep", "", flg, &flg, NULL)); 98 PetscCall(PetscOptionsBool("-monitor_temp", "Monitor the temperature each timestep", "", tflg, &tflg, NULL)); 99 PetscOptionsEnd(); 100 101 /* tchem requires periodic table in current directory */ 102 PetscCall(PetscFileRetrieve(PETSC_COMM_WORLD, periodic, lperiodic, PETSC_MAX_PATH_LEN, &found)); 103 PetscCheck(found, PETSC_COMM_WORLD, PETSC_ERR_FILE_OPEN, "Cannot located required periodic table %s or local version %s", periodic, lperiodic); 104 105 PetscCallTC(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 PetscCall(PetscMalloc1((user.Nspec + 1) * LENGTHOFSPECNAME, &names)); 113 PetscCall(PetscStrncpy(names, "Temp", (user.Nspec + 1) * LENGTHOFSPECNAME)); 114 TC_getSnames(user.Nspec, names + LENGTHOFSPECNAME); 115 PetscCall(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 PetscCall(PetscStrArrayallocpy((const char *const *)snames, &user.snames)); 119 PetscCall(PetscFree(snames)); 120 PetscCall(PetscFree(names)); 121 122 PetscCall(PetscMalloc3(user.Nspec + 1, &user.tchemwork, PetscSqr(user.Nspec + 1), &user.Jdense, user.Nspec + 1, &user.rows)); 123 PetscCall(VecCreateSeq(PETSC_COMM_SELF, user.Nspec + 1, &X)); 124 125 /* PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,PETSC_DECIDE,NULL,&J)); */ 126 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, user.Nspec + 1, user.Nspec + 1, NULL, &J)); 127 PetscCall(MatSetFromOptions(J)); 128 129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130 Create timestepping solver context 131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 132 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 133 PetscCall(TSSetType(ts, TSARKIMEX)); 134 PetscCall(TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE)); 135 PetscCall(TSARKIMEXSetType(ts, TSARKIMEX4)); 136 PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user)); 137 PetscCall(TSSetRHSJacobian(ts, J, J, FormRHSJacobian, &user)); 138 139 if (flg) PetscCall(TSMonitorSet(ts, MonitorMassConservation, NULL, NULL)); 140 if (tflg) PetscCall(TSMonitorSet(ts, MonitorTempature, &user, NULL)); 141 142 ftime = 1.0; 143 PetscCall(TSSetMaxTime(ts, ftime)); 144 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 145 146 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147 Set initial conditions 148 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 149 PetscCall(FormInitialSolution(ts, X, &user)); 150 PetscCall(TSSetSolution(ts, X)); 151 dt = 1e-10; /* Initial time step */ 152 PetscCall(TSSetTimeStep(ts, dt)); 153 PetscCall(TSGetAdapt(ts, &adapt)); 154 PetscCall(TSAdaptSetStepLimits(adapt, 1e-12, 1e-4)); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */ 155 PetscCall(TSSetMaxSNESFailures(ts, -1)); /* Retry step an unlimited number of times */ 156 157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 158 Set runtime options 159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 160 PetscCall(TSSetFromOptions(ts)); 161 162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 163 Set final conditions for sensitivities 164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 165 PetscCall(VecDuplicate(X, &lambda)); 166 PetscCall(TSSetCostGradients(ts, 1, &lambda, NULL)); 167 PetscCall(VecSetValue(lambda, 0, 1.0, INSERT_VALUES)); 168 PetscCall(VecAssemblyBegin(lambda)); 169 PetscCall(VecAssemblyEnd(lambda)); 170 171 PetscCall(TSGetTrajectory(ts, &tj)); 172 if (tj) { 173 PetscCall(TSTrajectorySetVariableNames(tj, (const char *const *)user.snames)); 174 PetscCall(TSTrajectorySetTransform(tj, (PetscErrorCode(*)(void *, Vec, Vec *))MassFractionToMoleFraction, NULL, &user)); 175 } 176 177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178 Pass information to graphical monitoring routine 179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 180 PetscCall(TSMonitorLGSetVariableNames(ts, (const char *const *)user.snames)); 181 PetscCall(TSMonitorLGSetTransform(ts, (PetscErrorCode(*)(void *, Vec, Vec *))MassFractionToMoleFraction, NULL, &user)); 182 183 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184 Solve ODE 185 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 186 PetscCall(TSSolve(ts, X)); 187 PetscCall(TSGetSolveTime(ts, &ftime)); 188 PetscCall(TSGetStepNumber(ts, &steps)); 189 PetscCall(TSGetConvergedReason(ts, &reason)); 190 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps)); 191 192 /* { 193 Vec max; 194 PetscInt i; 195 const PetscReal *bmax; 196 197 PetscCall(TSMonitorEnvelopeGetBounds(ts,&max,NULL)); 198 if (max) { 199 PetscCall(VecGetArrayRead(max,&bmax)); 200 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n")); 201 for (i=1; i<user.Nspec; i++) { 202 if (bmax[i] > .01) PetscCall(PetscPrintf(PETSC_COMM_SELF,"%s %g\n",user.snames[i],(double)bmax[i])); 203 } 204 PetscCall(VecRestoreArrayRead(max,&bmax)); 205 } 206 } 207 208 Vec y; 209 PetscCall(MassFractionToMoleFraction(&user,X,&y)); 210 PetscCall(PrintSpecies(&user,y)); 211 PetscCall(VecDestroy(&y)); */ 212 213 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 214 Free work space. 215 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 216 TC_reset(); 217 PetscCall(PetscStrArrayDestroy(&user.snames)); 218 PetscCall(MatDestroy(&J)); 219 PetscCall(VecDestroy(&X)); 220 PetscCall(VecDestroy(&lambda)); 221 PetscCall(TSDestroy(&ts)); 222 PetscCall(PetscFree3(user.tchemwork, user.Jdense, user.rows)); 223 PetscCall(PetscFinalize()); 224 return 0; 225 } 226 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; 232 233 PetscFunctionBeginUser; 234 PetscCall(VecGetArrayRead(X, &x)); 235 PetscCall(VecGetArray(F, &f)); 236 237 PetscCall(PetscArraycpy(user->tchemwork, x, user->Nspec + 1)); 238 user->tchemwork[0] *= user->Tini; /* Dimensionalize */ 239 PetscCallTC(TC_getSrc(user->tchemwork, user->Nspec + 1, f)); 240 f[0] /= user->Tini; /* Non-dimensionalize */ 241 242 PetscCall(VecRestoreArrayRead(X, &x)); 243 PetscCall(VecRestoreArray(F, &f)); 244 PetscFunctionReturn(PETSC_SUCCESS); 245 } 246 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; 252 253 PetscFunctionBeginUser; 254 PetscCall(VecGetArrayRead(X, &x)); 255 PetscCall(PetscArraycpy(user->tchemwork, x, user->Nspec + 1)); 256 PetscCall(VecRestoreArrayRead(X, &x)); 257 user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */ 258 PetscCall(TC_getJacTYN(user->tchemwork, user->Nspec, user->Jdense, 1)); 259 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 PetscCall(MatSetOption(Pmat, MAT_ROW_ORIENTED, PETSC_FALSE)); 264 PetscCall(MatSetOption(Pmat, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); 265 PetscCall(MatZeroEntries(Pmat)); 266 PetscCall(MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES)); 267 268 PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY)); 269 PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY)); 270 if (Amat != Pmat) { 271 PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY)); 272 PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY)); 273 } 274 PetscFunctionReturn(PETSC_SUCCESS); 275 } 276 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; 287 288 PetscFunctionBeginUser; 289 PetscCall(VecZeroEntries(X)); 290 PetscCall(VecGetArray(X, &x)); 291 x[0] = 1.0; /* Non-dimensionalized by user->Tini */ 292 293 PetscCall(PetscOptionsGetStringArray(NULL, NULL, "-initial_species", names, &smax, &flg)); 294 PetscCheck(smax >= 2, PETSC_COMM_SELF, PETSC_ERR_USER, "Must provide at least two initial species"); 295 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-initial_mole", molefracs, &mmax, &flg)); 296 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); 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])); 302 PetscCheck(ispec >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Could not find species %s", names[i]); 303 PetscCall(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++) PetscCall(PetscFree(names[i])); 307 PetscCall(VecRestoreArray(X, &x)); 308 /* PetscCall(PrintSpecies((User)ctx,X)); */ 309 PetscCall(MoleFractionToMassFraction((User)ctx, X, &y)); 310 PetscCall(VecCopy(y, X)); 311 PetscCall(VecDestroy(&y)); 312 PetscFunctionReturn(PETSC_SUCCESS); 313 } 314 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; 322 323 PetscFunctionBeginUser; 324 PetscCall(VecDuplicate(massf, molef)); 325 PetscCall(VecGetArrayRead(massf, &maf)); 326 PetscCall(VecGetArray(*molef, &mof)); 327 mof[0] = maf[0]; /* copy over temperature */ 328 TC_getMs2Ml((double *)maf + 1, user->Nspec, mof + 1); 329 PetscCall(VecRestoreArray(*molef, &mof)); 330 PetscCall(VecRestoreArrayRead(massf, &maf)); 331 PetscFunctionReturn(PETSC_SUCCESS); 332 } 333 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; 341 342 PetscFunctionBeginUser; 343 PetscCall(VecDuplicate(molef, massf)); 344 PetscCall(VecGetArrayRead(molef, &mof)); 345 PetscCall(VecGetArray(*massf, &maf)); 346 maf[0] = mof[0]; /* copy over temperature */ 347 TC_getMl2Ms((double *)mof + 1, user->Nspec, maf + 1); 348 PetscCall(VecRestoreArrayRead(molef, &mof)); 349 PetscCall(VecRestoreArray(*massf, &maf)); 350 PetscFunctionReturn(PETSC_SUCCESS); 351 } 352 353 PetscErrorCode ComputeMassConservation(Vec x, PetscReal *mass, void *ctx) 354 { 355 PetscFunctionBeginUser; 356 PetscCall(VecSum(x, mass)); 357 PetscFunctionReturn(PETSC_SUCCESS); 358 } 359 360 PetscErrorCode MonitorMassConservation(TS ts, PetscInt step, PetscReal time, Vec x, void *ctx) 361 { 362 const PetscScalar *T; 363 PetscReal mass; 364 365 PetscFunctionBeginUser; 366 PetscCall(ComputeMassConservation(x, &mass, ctx)); 367 PetscCall(VecGetArrayRead(x, &T)); 368 mass -= PetscAbsScalar(T[0]); 369 PetscCall(VecRestoreArrayRead(x, &T)); 370 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)))); 371 PetscFunctionReturn(PETSC_SUCCESS); 372 } 373 374 PetscErrorCode MonitorTempature(TS ts, PetscInt step, PetscReal time, Vec x, void *ctx) 375 { 376 User user = (User)ctx; 377 const PetscScalar *T; 378 379 PetscFunctionBeginUser; 380 PetscCall(VecGetArrayRead(x, &T)); 381 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep %" PetscInt_FMT " time %g temperature %g\n", step, (double)time, (double)(T[0] * user->Tini))); 382 PetscCall(VecRestoreArrayRead(x, &T)); 383 PetscFunctionReturn(PETSC_SUCCESS); 384 } 385 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; 393 394 PetscFunctionBeginUser; 395 PetscCall(PetscMalloc1(n, &idx)); 396 for (i = 0; i < n; i++) idx[i] = i; 397 PetscCall(VecGetArrayRead(molef, &mof)); 398 PetscCall(PetscSortRealWithPermutation(n, mof, idx)); 399 for (i = 0; i < n; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%6s %g\n", user->snames[idx[n - i - 1]], (double)PetscRealPart(mof[idx[n - i - 1]]))); 400 PetscCall(PetscFree(idx)); 401 PetscCall(VecRestoreArrayRead(molef, &mof)); 402 PetscFunctionReturn(PETSC_SUCCESS); 403 } 404 405 /*TEST 406 build: 407 requires: tchem 408 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 413 414 TEST*/ 415