30 template<
class CompType,
class ThermoType>
38 searchInitSet_(this->coeffsDict_.subDict(
"initialSet").size()),
41 sC_(this->nSpecie_,0),
42 sH_(this->nSpecie_,0),
43 sO_(this->nSpecie_,0),
44 sN_(this->nSpecie_,0),
51 phiTol_(this->tolerance()),
55 dict.subDict(
"reduction").lookupOrDefault<
word>
62 dict.subDict(
"reduction").lookupOrDefault<
word>
69 dict.subDict(
"reduction").lookupOrDefault<
word>
76 dict.subDict(
"reduction").lookupOrDefault<
word>
83 dict.subDict(
"reduction").lookupOrDefault<
word>
88 forceFuelInclusion_(
false)
94 if (initSet.
found(chemistry.
Y()[i].member()))
96 searchInitSet_[j++] = i;
99 if (j<searchInitSet_.size())
102 << searchInitSet_.size()-j
103 <<
" species in the initial set is not in the mechanism " 108 if (this->coeffsDict_.found(
"automaticSIS"))
110 automaticSIS_.readIfPresent(
"automaticSIS", this->coeffsDict_);
113 if (this->coeffsDict_.found(
"forceFuelInclusion"))
115 forceFuelInclusion_.readIfPresent
117 "forceFuelInclusion",this->coeffsDict_
121 if (this->coeffsDict_.found(
"phiTol"))
123 phiTol_ = this->coeffsDict_.template lookup<scalar>(
"phiTol");
126 if (this->coeffsDict_.found(
"NOxThreshold"))
129 this->coeffsDict_.template lookup<scalar>(
"NOxThreshold");
134 for (
label i=0; i<this->nSpecie_; i++)
137 specieComposition[i];
139 forAll(curSpecieComposition, j)
142 curSpecieComposition[j];
143 if (curElement.
name() ==
"C")
145 sC_[i] = curElement.
nAtoms();
147 else if (curElement.
name() ==
"H")
149 sH_[i] = curElement.
nAtoms();
151 else if (curElement.
name() ==
"O")
153 sO_[i] = curElement.
nAtoms();
155 else if (curElement.
name() ==
"N")
157 sN_[i] = curElement.
nAtoms();
161 Info<<
"element not considered"<<
endl;
164 if (this->chemistry_.Y()[i].member() == CO2Name_)
168 else if (this->chemistry_.Y()[i].member() == COName_)
172 else if (this->chemistry_.Y()[i].member() == HO2Name_)
176 else if (this->chemistry_.Y()[i].member() == H2OName_)
180 else if (this->chemistry_.Y()[i].member() == NOName_)
186 if ((CO2Id_==-1 || COId_==-1 || HO2Id_==-1 || H2OId_==-1) && automaticSIS_)
189 <<
"The name of the species used in automatic SIS are not found in " 190 <<
" the mechanism. You should either set the name for CO2, CO, H2O" 191 <<
" and HO2 properly or set automaticSIS to off " 200 if (this->coeffsDict_.found(
"fuelSpecies"))
202 fuelDict = this->coeffsDict_.
subDict(
"fuelSpecies");
203 fuelSpecies_ = fuelDict.
toc();
204 if (fuelSpecies_.size() == 0)
207 <<
"With automatic SIS, the fuel species should be " 208 <<
"specified in the fuelSpecies subDict" 215 <<
"With automatic SIS, the fuel species should be " 216 <<
"specified in the fuelSpecies subDict" 220 if (this->coeffsDict_.found(
"nbCLarge"))
225 fuelSpeciesID_.
setSize(fuelSpecies_.size());
226 fuelSpeciesProp_.setSize(fuelSpecies_.size());
231 fuelSpeciesProp_[i] = fuelDict.
lookup<scalar>(fuelSpecies_[i]);
232 for (
label j=0; j<this->nSpecie_; j++)
234 if (this->chemistry_.Y()[j].member() == fuelSpecies_[i])
236 fuelSpeciesID_[i] = j;
241 this->chemistry_.specieThermos()[fuelSpeciesID_[i]].W();
242 Mmtot += fuelSpeciesProp_[i]/curMm;
250 label curID = fuelSpeciesID_[i];
251 scalar curMm = this->chemistry_.specieThermos()[curID].W();
253 nbC += fuelSpeciesProp_[i]*Mmtot/curMm*sC_[curID];
254 nbO += fuelSpeciesProp_[i]*Mmtot/curMm*sO_[curID];
263 template<
class CompType,
class ThermoType>
271 template<
class CompType,
class ThermoType>
280 scalarField& completeC(this->chemistry_.completeC());
282 for(
label i=0; i<this->nSpecie_; i++)
288 c1[this->nSpecie_] =
T;
289 c1[this->nSpecie_+1] =
p;
303 scalar pf, cf, pr, cr;
305 forAll(this->chemistry_.reactions(), i)
310 scalar omegai = this->chemistry_.omega
312 R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef
331 scalar sl = -R.
lhs()[
s].stoichCoeff;
350 while(!usedIndex.empty())
353 if (deltaBi[curIndex])
356 deltaBi[curIndex] =
false;
358 if (rABPos(ss, curIndex)==-1)
361 rABPos(ss, curIndex) = NbrABInit[ss];
364 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
366 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
370 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
386 wA.append(sl*omegai);
394 scalar sl = R.
rhs()[
s].stoichCoeff;
413 while(!usedIndex.empty())
416 if (deltaBi[curIndex])
419 deltaBi[curIndex] =
false;
422 if (rABPos(ss, curIndex) == -1)
425 rABPos(ss, curIndex) = NbrABInit[ss];
427 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
428 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
432 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
448 wA.append(sl*omegai);
461 if (PA[wAID[
id]] == 0.0)
463 PA[wAID[id]] = wA[id];
467 PA[wAID[id]] += wA[id];
472 if (CA[wAID[
id]] == 0.0)
474 CA[wAID[id]] = -wA[id];
478 CA[wAID[id]] += -wA[id];
485 scalar phiLarge(0.0);
486 scalar phiProgress(0.0);
497 for (
label i=0; i<this->nSpecie_; i++)
502 this->chemistry_.Y()[i].member() ==
"CO2" 503 || this->chemistry_.Y()[i].member() ==
"H2O" 508 Na[0] += sC_[i]*c[i];
509 Na[1] += sH_[i]*c[i];
510 Na[2] += sO_[i]*c[i];
511 if (sC_[i]>nbCLarge_ || this->chemistry_.Y()[i].member() ==
"O2")
513 Nal[0] += sC_[i]*c[i];
514 Nal[1] += sH_[i]*c[i];
515 Nal[2] += sO_[i]*c[i];
524 phiProgress = (2*Na[0]+Na[1]/2-zprime_*Na[0])/(Na[2]-zprime_*Na[0]);
529 phiLarge = (2*Nal[0]+Nal[1]/2)/Nal[2];
535 label speciesNumber = 0;
539 for (
label i=0; i<this->nSpecie_; i++)
541 this->activeSpecies_[i] =
false;
553 if (phiLarge >= phiTol_ && phiProgress >= phiTol_)
559 this->activeSpecies_[COId_] =
true;
563 this->activeSpecies_[HO2Id_] =
true;
564 Rvalue[HO2Id_] = 1.0;
567 Q.
push(fuelSpeciesID_[i]);
569 this->activeSpecies_[fuelSpeciesID_[i]] =
true;
570 Rvalue[fuelSpeciesID_[i]] = 1.0;
574 else if (phiLarge < phiTol_ && phiProgress >= phiTol_)
580 this->activeSpecies_[COId_] =
true;
584 this->activeSpecies_[HO2Id_] =
true;
585 Rvalue[HO2Id_] = 1.0;
587 if (forceFuelInclusion_)
591 Q.
push(fuelSpeciesID_[i]);
593 this->activeSpecies_[fuelSpeciesID_[i]] =
true;
594 Rvalue[fuelSpeciesID_[i]] = 1.0;
604 this->activeSpecies_[CO2Id_] =
true;
605 Rvalue[CO2Id_] = 1.0;
609 this->activeSpecies_[H2OId_] =
true;
610 Rvalue[H2OId_] = 1.0;
611 if (forceFuelInclusion_)
615 Q.
push(fuelSpeciesID_[i]);
617 this->activeSpecies_[fuelSpeciesID_[i]] =
true;
618 Rvalue[fuelSpeciesID_[i]] = 1.0;
623 if (T>NOxThreshold_ && NOId_!=-1)
627 this->activeSpecies_[NOId_] =
true;
636 this->activeSpecies_[q] =
true;
647 scalar Den =
max(PA[u],CA[u]);
650 for (
label v=0; v<NbrABInit[u]; v++)
652 label otherSpec = rABOtherSpec(u, v);
653 scalar rAB =
mag(rABNum(u, v))/Den;
659 if (rAB >= this->tolerance())
661 scalar Rtemp = Rvalue[u]*rAB;
665 if ((Rvalue[otherSpec]<Rtemp) && (Rtemp>=this->tolerance()))
668 Rvalue[otherSpec] = Rtemp;
669 if (!this->activeSpecies_[otherSpec])
671 this->activeSpecies_[otherSpec] =
true;
681 forAll(this->chemistry_.reactions(), i)
684 this->chemistry_.reactionsDisabled()[i] =
false;
689 if (!this->activeSpecies_[ss])
692 this->chemistry_.reactionsDisabled()[i] =
true;
696 if (!this->chemistry_.reactionsDisabled()[i])
701 if (!this->activeSpecies_[ss])
704 this->chemistry_.reactionsDisabled()[i] =
true;
711 this->NsSimp_ = speciesNumber;
712 scalarField& simplifiedC(this->chemistry_.simplifiedC());
713 simplifiedC.
setSize(this->NsSimp_+2);
716 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
719 for (
label i=0; i<this->nSpecie_; i++)
721 if (this->activeSpecies_[i])
724 simplifiedC[j] = c[i];
726 if (!this->chemistry_.active(i))
728 this->chemistry_.setActive(i);
737 simplifiedC[this->NsSimp_] =
T;
738 simplifiedC[this->NsSimp_+1] =
p;
739 this->chemistry_.setNsDAC(this->NsSimp_);
743 this->chemistry_.setNSpecie(this->NsSimp_);
A FIFO stack based on a singly-linked list.
virtual label nSpecie() const
The number of species.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A list of keyword definitions, which are a keyword followed by any number of values (e...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, const label li)
Reduce the mechanism.
const dimensionedScalar & c1
First radiation constant: default SI units: [W/m^2].
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void size(const label)
Override size to be inconsistent with allocated storage.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const PtrList< volScalarField > & Y()
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
BasicChemistryModel< rhoReactionThermo > & chemistry
wordList toc() const
Return the table of contents.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
void setSize(const label)
Alter the addressed list size.
A class for handling words, derived from string.
DAC(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
const word & name() const
Return the name of the element.
label nAtoms() const
Return the number of atoms of this element in the specie.
List< List< specieElement > > & specieComp()
virtual ~DAC()
Destructor.
void setSize(const label)
Reset size of List.
#define R(A, B, C, D, E, F, K, M)
void push(const T &a)
Push an element onto the stack.
dimensioned< scalar > mag(const dimensioned< Type > &)
T pop()
Pop the bottom element off the stack.
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
An abstract class for methods of chemical mechanism reduction.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.