33 template<
class CompType,
class ThermoType>
41 searchInitSet_(this->coeffsDict_.subDict(
"initialSet").size()),
43 nbCLarge_(this->coeffsDict_.template getOrDefault<label>(
"nbCLarge", 3)),
44 sC_(this->nSpecie_, 0),
45 sH_(this->nSpecie_, 0),
46 sO_(this->nSpecie_, 0),
47 sN_(this->nSpecie_, 0),
55 this->coeffsDict_.template getOrDefault<
Switch>
63 this->coeffsDict_.template getOrDefault<scalar>
65 "phiTol", this->tolerance()
70 this->coeffsDict_.template getOrDefault<scalar>
76 CO2Name_(this->coeffsDict_.template getOrDefault<
word>(
"CO2",
"CO2")),
77 COName_(this->coeffsDict_.template getOrDefault<
word>(
"CO",
"CO")),
78 HO2Name_(this->coeffsDict_.template getOrDefault<
word>(
"HO2",
"HO2")),
79 H2OName_(this->coeffsDict_.template getOrDefault<
word>(
"H2O",
"H2O")),
80 NOName_(this->coeffsDict_.template getOrDefault<
word>(
"NO",
"NO")),
83 this->coeffsDict_.template getOrDefault<
Switch>
93 for (label i = 0; i <
chemistry.nSpecie(); i++)
97 searchInitSet_[j++] = i;
101 if (j < searchInitSet_.
size())
104 << searchInitSet_.
size()-j
105 <<
" species in the initial set is not in the mechanism " 114 for (label i=0; i<this->
nSpecie_; i++)
117 specieComposition[i];
120 forAll(curSpecieComposition, j)
124 if (curElement.
name() ==
"C")
126 sC_[i] = curElement.
nAtoms();
128 else if (curElement.
name() ==
"H")
130 sH_[i] = curElement.
nAtoms();
132 else if (curElement.
name() ==
"O")
134 sO_[i] = curElement.
nAtoms();
136 else if (curElement.
name() ==
"N")
138 sN_[i] = curElement.
nAtoms();
142 Info<<
" element " << curElement.
name() <<
" not considered" 150 else if (this->
chemistry_.
Y()[i].member() == COName_)
154 else if (this->
chemistry_.
Y()[i].member() == HO2Name_)
158 else if (this->
chemistry_.
Y()[i].member() == H2OName_)
162 else if (this->
chemistry_.
Y()[i].member() == NOName_)
168 if ((CO2Id_==-1 || COId_==-1 || HO2Id_==-1 || H2OId_==-1) && automaticSIS_)
171 <<
"The name of the species used in automatic SIS are not found in " 172 <<
" the mechanism. You should either set the name for CO2, CO, H2O" 173 <<
" and HO2 properly or set automaticSIS to off " 185 fuelSpecies_ = fuelDict.
toc();
186 if (fuelSpecies_.
size() == 0)
189 <<
"With automatic SIS, the fuel species should be " 190 <<
"specified in the fuelSpecies subDict" 197 <<
"With automatic SIS, the fuel species should be " 198 <<
"specified in the fuelSpecies subDict" 209 fuelSpeciesProp_[i] = fuelDict.
get<scalar>(fuelSpecies_[i]);
210 for (label j=0; j<this->
nSpecie_; j++)
212 if (this->
chemistry_.
Y()[j].member() == fuelSpecies_[i])
214 fuelSpeciesID_[i] = j;
220 Mmtot += fuelSpeciesProp_[i]/curMm;
228 label curID = fuelSpeciesID_[i];
231 nbC += fuelSpeciesProp_[i]*Mmtot/curMm*sC_[curID];
232 nbO += fuelSpeciesProp_[i]*Mmtot/curMm*sO_[curID];
241 template<
class CompType,
class ThermoType>
249 template<
class CompType,
class ThermoType>
257 scalarField& completeC(this->chemistry_.completeC());
259 for (label i=0; i<this->nSpecie_; i++)
265 c1[this->nSpecie_] =
T;
266 c1[this->nSpecie_+1] =
p;
280 scalar pf, cf, pr, cr;
285 scalar omegai = this->chemistry_.omega
287 R,
c1,
T,
p, pf, cf, lRef, pr, cr, rRef
305 label ss =
R.lhs()[
s].index;
306 scalar sl = -
R.lhs()[
s].stoichCoeff;
311 label sj =
R.lhs()[j].index;
317 label sj =
R.rhs()[j].index;
325 while (!usedIndex.empty())
327 label curIndex = usedIndex.
pop();
328 if (deltaBi[curIndex])
331 deltaBi[curIndex] =
false;
333 if (rABPos(ss, curIndex) == -1)
336 rABPos(ss, curIndex) = NbrABInit[ss];
339 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
341 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
345 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
361 wA.append(sl*omegai);
369 label ss =
R.rhs()[
s].index;
370 scalar sl =
R.rhs()[
s].stoichCoeff;
375 label sj =
R.lhs()[j].index;
381 label sj =
R.rhs()[j].index;
389 while (!usedIndex.empty())
391 label curIndex = usedIndex.
pop();
392 if (deltaBi[curIndex])
395 deltaBi[curIndex] =
false;
398 if (rABPos(ss, curIndex) == -1)
401 rABPos(ss, curIndex) = NbrABInit[ss];
403 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
404 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
408 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
424 wA.append(sl*omegai);
437 if (PA[wAID[
id]] == 0.0)
439 PA[wAID[id]] = wA[id];
443 PA[wAID[id]] += wA[id];
448 if (CA[wAID[
id]] == 0.0)
450 CA[wAID[id]] = -wA[id];
454 CA[wAID[id]] += -wA[id];
461 scalar phiLarge(0.0);
462 scalar phiProgress(0.0);
473 for (label i=0; i<this->nSpecie_; i++)
478 this->chemistry_.Y()[i].member() ==
"CO2" 479 || this->chemistry_.Y()[i].member() ==
"H2O" 485 Na[0] += sC_[i]*
c[i];
486 Na[1] += sH_[i]*
c[i];
487 Na[2] += sO_[i]*
c[i];
488 if (sC_[i]>nbCLarge_ || this->chemistry_.Y()[i].member() ==
"O2")
490 Nal[0] += sC_[i]*
c[i];
491 Nal[1] += sH_[i]*
c[i];
492 Nal[2] += sO_[i]*
c[i];
501 phiProgress = (2*Na[0]+Na[1]/2-zprime_*Na[0])/(Na[2]-zprime_*Na[0]);
506 phiLarge = (2*Nal[0] + Nal[1]/2)/Nal[2];
512 label speciesNumber = 0;
516 for (label i=0; i<this->nSpecie_; ++i)
518 this->activeSpecies_[i] =
false;
530 if (phiLarge >= phiTol_ && phiProgress >= phiTol_)
536 this->activeSpecies_[COId_] =
true;
540 this->activeSpecies_[HO2Id_] =
true;
541 Rvalue[HO2Id_] = 1.0;
542 for (
const label fuelId : fuelSpeciesID_)
546 this->activeSpecies_[fuelId] =
true;
547 Rvalue[fuelId] = 1.0;
551 else if (phiLarge < phiTol_ && phiProgress >= phiTol_)
557 this->activeSpecies_[COId_] =
true;
561 this->activeSpecies_[HO2Id_] =
true;
562 Rvalue[HO2Id_] = 1.0;
564 if (forceFuelInclusion_)
566 for (
const label fuelId : fuelSpeciesID_)
570 this->activeSpecies_[fuelId] =
true;
571 Rvalue[fuelId] = 1.0;
581 this->activeSpecies_[CO2Id_] =
true;
582 Rvalue[CO2Id_] = 1.0;
586 this->activeSpecies_[H2OId_] =
true;
587 Rvalue[H2OId_] = 1.0;
588 if (forceFuelInclusion_)
590 for (
const label fuelId : fuelSpeciesID_)
594 this->activeSpecies_[fuelId] =
true;
595 Rvalue[fuelId] = 1.0;
600 if (
T > NOxThreshold_ && NOId_ != -1)
604 this->activeSpecies_[NOId_] =
true;
610 for (label i=0; i<SIS.
size(); i++)
613 this->activeSpecies_[q] =
true;
624 scalar Den =
max(PA[u], CA[u]);
627 for (label v=0; v<NbrABInit[u]; ++v)
629 label otherSpec = rABOtherSpec(u, v);
630 scalar rAB =
mag(rABNum(u, v))/Den;
637 if (rAB >= this->tolerance())
639 scalar Rtemp = Rvalue[u]*rAB;
643 if ((Rvalue[otherSpec]<Rtemp) && (Rtemp>=this->tolerance()))
646 Rvalue[otherSpec] = Rtemp;
647 if (!this->activeSpecies_[otherSpec])
649 this->activeSpecies_[otherSpec] =
true;
659 forAll(this->chemistry_.reactions(), i)
662 this->chemistry_.reactionsDisabled()[i] =
false;
666 label ss =
R.lhs()[
s].index;
667 if (!this->activeSpecies_[ss])
670 this->chemistry_.reactionsDisabled()[i] =
true;
674 if (!this->chemistry_.reactionsDisabled()[i])
678 label ss =
R.rhs()[
s].index;
679 if (!this->activeSpecies_[ss])
682 this->chemistry_.reactionsDisabled()[i] =
true;
689 this->NsSimp_ = speciesNumber;
690 scalarField& simplifiedC(this->chemistry_.simplifiedC());
691 simplifiedC.
setSize(this->NsSimp_ + 2);
694 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
697 for (label i=0; i<this->nSpecie_; i++)
699 if (this->activeSpecies_[i])
702 simplifiedC[j] =
c[i];
704 if (!this->chemistry_.active(i))
706 this->chemistry_.setActive(i);
715 simplifiedC[this->NsSimp_] =
T;
716 simplifiedC[this->NsSimp_+1] =
p;
717 this->chemistry_.setNsDAC(this->NsSimp_);
721 this->chemistry_.setNSpecie(this->NsSimp_);
A FIFO stack based on a singly-linked list.
void size(const label n)
Older name for setAddressableSize.
PtrList< volScalarField > & Y()
TDACChemistryModel< CompType, ThermoType > & chemistry_
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const PtrList< ThermoType > & specieThermo() const
Thermodynamic data of the species.
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, any/none. Also accepts 0/1 as a string and shortcuts t/f, y/n.
void setSize(const label n)
Same as resize()
BasicChemistryModel< psiReactionThermo > & chemistry
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T. FatalIOError if not found, or if the number of tokens is incorrect.
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
wordList toc() const
Return the table of contents.
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry (const access) with the given keyword.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
#define forAll(list, i)
Loop across all elements in list.
void push(const T &elem)
Push an element onto the back of the stack.
void setSize(const label n)
Alias for resize()
A class for handling words, derived from Foam::string.
DAC(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
const label nSpecie_
Number of species.
const word & name() const
Return the name of the element.
label nAtoms() const
Return the number of atoms of this element in the specie.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
virtual ~DAC()
Destructor.
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
const dictionary coeffsDict_
Dictionary that store the algorithm data.
#define R(A, B, C, D, E, F, K, M)
const dimensionedScalar c
Speed of light in a vacuum.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
messageStream Info
Information stream (stdout output on master, null elsewhere)
T pop()
Pop the bottom element off the stack.
An abstract class for methods of chemical mechanism reduction.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static constexpr const zero Zero
Global zero (0)