58const char* Foam::chemkinReader::reactionTypeNames[4] =
62 "nonEquilibriumReversible",
66const char* Foam::chemkinReader::reactionRateTypeNames[8] =
70 "unimolecularFallOff",
71 "chemicallyActivatedBimolecular",
75 "unknownReactionRateType"
78const char* Foam::chemkinReader::fallOffFunctionNames[4] =
83 "unknownFallOffFunctionType"
86void Foam::chemkinReader::initReactionKeywordTable()
88 reactionKeywordTable_.insert(
"M", thirdBodyReactionType);
89 reactionKeywordTable_.insert(
"LOW", unimolecularFallOffReactionType);
90 reactionKeywordTable_.insert
93 chemicallyActivatedBimolecularReactionType
95 reactionKeywordTable_.insert(
"TROE", TroeReactionType);
96 reactionKeywordTable_.insert(
"SRI", SRIReactionType);
97 reactionKeywordTable_.insert(
"LT", LandauTellerReactionType);
98 reactionKeywordTable_.insert(
"RLT", reverseLandauTellerReactionType);
99 reactionKeywordTable_.insert(
"JAN", JanevReactionType);
100 reactionKeywordTable_.insert(
"FIT1", powerSeriesReactionRateType);
101 reactionKeywordTable_.insert(
"HV", radiationActivatedReactionType);
102 reactionKeywordTable_.insert(
"TDEP", speciesTempReactionType);
103 reactionKeywordTable_.insert(
"EXCI", energyLossReactionType);
104 reactionKeywordTable_.insert(
"MOME", plasmaMomentumTransfer);
105 reactionKeywordTable_.insert(
"XSMI", collisionCrossSection);
106 reactionKeywordTable_.insert(
"REV", nonEquilibriumReversibleReactionType);
107 reactionKeywordTable_.insert(
"DUPLICATE", duplicateReactionType);
108 reactionKeywordTable_.insert(
"DUP", duplicateReactionType);
109 reactionKeywordTable_.insert(
"FORD", speciesOrderForward);
110 reactionKeywordTable_.insert(
"RORD", speciesOrderReverse);
111 reactionKeywordTable_.insert(
"UNITS", UnitsOfReaction);
112 reactionKeywordTable_.insert(
"END", end);
116Foam::scalar Foam::chemkinReader::molecularWeight
123 forAll(specieComposition, i)
125 label nAtoms = specieComposition[i].nAtoms();
126 const word& elementName = specieComposition[i].name();
128 if (isotopeAtomicWts_.found(elementName))
130 molWt += nAtoms*isotopeAtomicWts_[elementName];
139 <<
"Unknown element " << elementName
140 <<
" on line " << lineNo_-1 <<
nl
141 <<
" specieComposition: " << specieComposition
150void Foam::chemkinReader::checkCoeffs
153 const char* reactionRateName,
157 if (reactionCoeffs.size() != nCoeffs)
160 <<
"Wrong number of coefficients for the " << reactionRateName
161 <<
" rate expression on line "
162 << lineNo_-1 <<
", should be "
163 << nCoeffs <<
" but " << reactionCoeffs.size() <<
" supplied." <<
nl
164 <<
"Coefficients are "
165 << reactionCoeffs <<
nl
170template<
class ReactionRateType>
171void Foam::chemkinReader::addReactionType
173 const reactionType rType,
176 const ReactionRateType& rr
226 <<
"Unhandled reaction type " << reactionTypeNames[rType]
227 <<
" on line " << lineNo_-1 <<
nl
233 <<
"Unknown reaction type (" << int(rType)
234 <<
"), on line " << lineNo_-1 <<
nl
241template<
template<
class,
class>
class PressureDependencyType>
242void Foam::chemkinReader::addPressureDependentReaction
244 const reactionType rType,
245 const fallOffFunctionType fofType,
252 const scalar Afactor0,
253 const scalar AfactorInf,
257 checkCoeffs(k0Coeffs,
"k0", 3);
258 checkCoeffs(kInfCoeffs,
"kInf", 3);
268 PressureDependencyType
273 Afactor0*k0Coeffs[0],
279 AfactorInf*kInfCoeffs[0],
293 reactionCoeffsTable[fallOffFunctionNames[fofType]]
296 if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
299 <<
"Wrong number of coefficients for Troe rate expression"
300 " on line " << lineNo_-1 <<
", should be 3 or 4 but "
301 << TroeCoeffs.size() <<
" supplied." <<
nl
302 <<
"Coefficients are "
307 if (TroeCoeffs.size() == 3)
309 TroeCoeffs.setSize(4);
310 TroeCoeffs[3] = GREAT;
317 PressureDependencyType
322 Afactor0*k0Coeffs[0],
328 AfactorInf*kInfCoeffs[0],
348 reactionCoeffsTable[fallOffFunctionNames[fofType]]
351 if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
354 <<
"Wrong number of coefficients for SRI rate expression"
355 " on line " << lineNo_-1 <<
", should be 3 or 5 but "
356 << SRICoeffs.size() <<
" supplied." <<
nl
357 <<
"Coefficients are "
362 if (SRICoeffs.size() == 3)
364 SRICoeffs.setSize(5);
373 PressureDependencyType
378 Afactor0*k0Coeffs[0],
384 AfactorInf*kInfCoeffs[0],
404 <<
"Fall-off function type (" << int(fofType)
405 <<
") not implemented, line " << lineNo_-1
412void Foam::chemkinReader::addReaction
417 const reactionType rType,
418 const reactionRateType rrType,
419 const fallOffFunctionType fofType,
425 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
432 speciesComposition_[speciesTable_[lhs[i].index]];
434 forAll(specieComposition, j)
436 label elementi = elementIndices_[specieComposition[j].name()];
438 lhs[i].stoichCoeff*specieComposition[j].nAtoms();
445 speciesComposition_[speciesTable_[
rhs[i].index]];
447 forAll(specieComposition, j)
449 label elementi = elementIndices_[specieComposition[j].name()];
451 rhs[i].stoichCoeff*specieComposition[j].nAtoms();
458 const scalar concFactor = 0.001;
462 sumExp += lhs[i].exponent;
464 scalar Afactor =
pow(concFactor, sumExp - 1.0);
466 scalar AfactorRev = Afactor;
468 if (rType == nonEquilibriumReversible)
473 sumExp +=
rhs[i].exponent;
475 AfactorRev =
pow(concFactor, sumExp - 1.0);
482 if (rType == nonEquilibriumReversible)
485 reactionCoeffsTable[reactionTypeNames[rType]];
487 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
503 Afactor*ArrheniusCoeffs[0],
505 ArrheniusCoeffs[2]/
RR
509 AfactorRev*reverseArrheniusCoeffs[0],
510 reverseArrheniusCoeffs[1],
511 reverseArrheniusCoeffs[2]/
RR
524 Afactor*ArrheniusCoeffs[0],
526 ArrheniusCoeffs[2]/
RR
532 case thirdBodyArrhenius:
534 if (rType == nonEquilibriumReversible)
537 reactionCoeffsTable[reactionTypeNames[rType]];
539 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
559 Afactor*concFactor*ArrheniusCoeffs[0],
561 ArrheniusCoeffs[2]/
RR,
566 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
567 reverseArrheniusCoeffs[1],
568 reverseArrheniusCoeffs[2]/
RR,
582 Afactor*concFactor*ArrheniusCoeffs[0],
584 ArrheniusCoeffs[2]/
RR,
591 case unimolecularFallOff:
593 addPressureDependentReaction<FallOffReactionRate>
600 reactionCoeffsTable[reactionRateTypeNames[rrType]],
609 case chemicallyActivatedBimolecular:
611 addPressureDependentReaction<ChemicallyActivatedReactionRate>
619 reactionCoeffsTable[reactionRateTypeNames[rrType]],
630 reactionCoeffsTable[reactionRateTypeNames[rrType]];
631 checkCoeffs(LandauTellerCoeffs,
"Landau-Teller", 2);
633 if (rType == nonEquilibriumReversible)
636 reactionCoeffsTable[reactionTypeNames[rType]];
637 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
642 word(reactionTypeNames[rType])
643 + reactionRateTypeNames[rrType]
645 checkCoeffs(LandauTellerCoeffs,
"reverse Landau-Teller", 2);
661 Afactor*ArrheniusCoeffs[0],
663 ArrheniusCoeffs[2]/
RR,
664 LandauTellerCoeffs[0],
665 LandauTellerCoeffs[1]
669 AfactorRev*reverseArrheniusCoeffs[0],
670 reverseArrheniusCoeffs[1],
671 reverseArrheniusCoeffs[2]/
RR,
672 reverseLandauTellerCoeffs[0],
673 reverseLandauTellerCoeffs[1]
686 Afactor*ArrheniusCoeffs[0],
688 ArrheniusCoeffs[2]/
RR,
689 LandauTellerCoeffs[0],
690 LandauTellerCoeffs[1]
699 reactionCoeffsTable[reactionRateTypeNames[rrType]];
701 checkCoeffs(JanevCoeffs,
"Janev", 9);
709 Afactor*ArrheniusCoeffs[0],
711 ArrheniusCoeffs[2]/
RR,
720 reactionCoeffsTable[reactionRateTypeNames[rrType]];
722 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
730 Afactor*ArrheniusCoeffs[0],
732 ArrheniusCoeffs[2]/
RR,
738 case unknownReactionRateType:
741 <<
"Internal error on line " << lineNo_-1
742 <<
": reaction rate type has not been set"
749 <<
"Reaction rate type (" << int(rrType)
750 <<
") not implemented, on line " << lineNo_-1
758 if (
mag(nAtoms[i]) > imbalanceTol_)
761 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
762 <<
" in " << elementNames_[i]
763 <<
" in reaction" <<
nl
764 << reactions_.last() <<
nl
765 <<
" on line " << lineNo_-1
772 reactionCoeffsTable.clear();
776void Foam::chemkinReader::read
783 transportDict_.read(
IFstream(transportFileName)());
785 if (!thermoFileName.empty())
787 std::ifstream thermoStream(thermoFileName);
792 <<
"file " << thermoFileName <<
" not found"
796 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
797 yy_switch_to_buffer(bufferPtr);
802 yy_delete_buffer(bufferPtr);
807 std::ifstream CHEMKINStream(CHEMKINFileName);
812 <<
"file " << CHEMKINFileName <<
" not found"
816 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
817 yy_switch_to_buffer(bufferPtr);
819 initReactionKeywordTable();
824 yy_delete_buffer(bufferPtr);
830Foam::chemkinReader::chemkinReader
842 reactions_(speciesTable_, speciesThermo_),
843 newFormat_(newFormat),
844 imbalanceTol_(ROOTSMALL)
846 read(CHEMKINFileName, thermoFileName, transportFileName);
850Foam::chemkinReader::chemkinReader
859 reactions_(speciesTable_, speciesThermo_),
860 newFormat_(
thermoDict.getOrDefault(
"newFormat", false)),
861 imbalanceTol_(
thermoDict.getOrDefault(
"imbalanceTolerance", ROOTSMALL))
865 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
868 fileName chemkinFile(
thermoDict.get<fileName>(
"CHEMKINFile"));
869 chemkinFile.expand();
872 thermoDict.readIfPresent(
"CHEMKINThermoFile", thermoFile);
875 fileName transportFile(
thermoDict.get<fileName>(
"CHEMKINTransportFile"));
876 transportFile.expand();
882 if (!chemkinFile.isAbsolute())
884 chemkinFile = relPath/chemkinFile;
887 if (!thermoFile.empty() && !thermoFile.isAbsolute())
889 thermoFile = relPath/thermoFile;
892 if (!transportFile.isAbsolute())
894 transportFile = relPath/transportFile;
898 read(chemkinFile, thermoFile, transportFile);
Macros for easy insertion into run-time selection tables.
#define addChemistryReaderType(Reader, Thermo)
Arrhenius reaction rate given by:
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
A 1D vector of objects of type <T> with a fixed length <N>.
A HashTable similar to std::unordered_map.
Input from file stream as an ISstream, normally using std::ifstream for the actual input.
Simple extension of Reaction to handle reversible reactions using equilibrium thermodynamics.
Janev, Langer, Evans and Post reaction rate.
Landau-Teller reaction rate.
Lindemann fall-off function.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Simple extension of Reaction to handle reversible reactions using equilibrium thermodynamics.
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Simple extension of Reaction to handle reversible reactions using equilibrium thermodynamics.
The SRI fall-off function.
The Troe fall-off function.
const speciesTable & species() const
Table of species.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
A class for handling file names.
static bool isAbsolute(const std::string &str)
Return true if filename starts with a '/' or '\' or (windows-only) with a filesystem-root.
Power series reaction rate.
string & expand(const bool allowEmpty=false)
Inplace expand initial tags, tildes, and all occurrences of environment variables as per stringOps::e...
Arrhenius reaction rate enhanced by third-body interaction.
A class for handling words, derived from Foam::string.
const dictionary & thermoDict
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const scalar RR
Universal gas constant: default in [J/(kmol K)].
bool read(const char *buf, int32_t &val)
Same as readInt32.
sutherlandTransport< species::thermo< janafThermo< perfectGas< specie > >, sensibleEnthalpy > > gasHThermoPhysics
messageStream Info
Information stream (stdout output on master, null elsewhere).
hashedWordList speciesTable
A table of species as a hashedWordList.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void rhs(fvMatrix< typename Expr::value_type > &m, const Expr &expression)
static constexpr const zero Zero
Global zero (0).
atomicWeightTable atomicWeights
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
errorManipArg< error, int > exit(error &err, const int errNo=1)
List< scalar > scalarList
List of scalar.
constexpr char nl
The newline '\n' character (0x0a).
#define forAll(list, i)
Loop across all elements in list.