33template<
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];
241template<
class CompType,
class ThermoType>
249template<
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;
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;
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_);
#define R(A, B, C, D, E, F, K, M)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
void setSize(const label n)
Same as resize()
void append(const T &val)
Copy append an element to the end of this list.
A FIFO stack based on a singly-linked list.
void push(const T &elem)
Push an element onto the back of the stack.
T pop()
Pop the bottom element off the stack.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void setSize(const label n)
Alias for resize()
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
A templated (M x N) rectangular matrix of objects of <Type>, containing M*N elements,...
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,...
Extends StandardChemistryModel by adding the TDAC method.
PtrList< volScalarField > & Y()
void size(const label n)
Older name for setAddressableSize.
An abstract class for methods of chemical mechanism reduction.
TDACChemistryModel< CompType, ThermoType > & chemistry_
const label nSpecie_
Number of species.
const dictionary coeffsDict_
Dictionary that store the algorithm data.
The Dynamic Adaptive Chemistry (DAC) method [1] simplifies the chemistry using the matrix rAB defined...
virtual ~DAC()
Destructor.
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
wordList toc() const
Return the table of contents.
label nAtoms() const
Return the number of atoms of this element in the specie.
const word & name() const
Return the name of the element.
A class for handling words, derived from Foam::string.
BasicChemistryModel< psiReactionThermo > & chemistry
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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))
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)
errorManipArg< error, int > exit(error &err, const int errNo=1)
#define forAll(list, i)
Loop across all elements in list.