31 template<
class CompType,
class ThermoType>
39 searchInitSet_(this->coeffsDict_.subDict(
"initialSet").size()),
40 sC_(this->nSpecie_,0),
41 sH_(this->nSpecie_,0),
42 sO_(this->nSpecie_,0),
43 sN_(this->nSpecie_,0),
50 if (initSet.
found(chemistry.
Y()[i].member()))
52 searchInitSet_[j++]=i;
55 if (j<searchInitSet_.size())
58 << searchInitSet_.size()-j
59 <<
" species in the initial set is not in the mechanism " 64 if (this->coeffsDict_.found(
"NGroupBased"))
66 NGroupBased_ =
readLabel(this->coeffsDict_.lookup(
"NGroupBased"));
70 this->chemistry_.specieComp();
71 for (
label i=0; i<this->nSpecie_; i++)
76 forAll(curSpecieComposition, j)
79 curSpecieComposition[j];
80 if (curElement.
name() ==
"C")
82 sC_[i] = curElement.
nAtoms();
84 else if (curElement.
name() ==
"H")
86 sH_[i] = curElement.
nAtoms();
88 else if (curElement.
name() ==
"O")
90 sO_[i] = curElement.
nAtoms();
92 else if (curElement.
name() ==
"N")
94 sN_[i] = curElement.
nAtoms();
98 Info<<
"element not considered"<<
endl;
107 template<
class CompType,
class ThermoType>
114 template<
class CompType,
class ThermoType>
123 scalarField& completeC(this->chemistry_.completeC());
126 for (
label i=0; i<this->nSpecie_; i++)
132 c1[this->nSpecie_] =
T;
133 c1[this->nSpecie_+1] =
p;
147 scalar pf, cf, pr, cr;
149 scalarField omegaV(this->chemistry_.reactions().size());
150 forAll(this->chemistry_.reactions(), i)
154 scalar omegai = this->chemistry_.omega
156 R, c1, T, p, pf, cf, lRef, pr, cr, rRef
169 scalar sl = -R.
lhs()[
s].stoichCoeff;
188 while(!usedIndex.empty())
191 if (deltaBi[curIndex])
194 deltaBi[curIndex] =
false;
196 if (rABPos(ss, curIndex)==-1)
198 rABPos(ss, curIndex) = NbrABInit[ss];
200 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
201 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
205 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
220 wA.append(sl*omegai);
229 scalar sl = R.
rhs()[
s].stoichCoeff;
248 while(!usedIndex.empty())
251 if (deltaBi[curIndex])
254 deltaBi[curIndex] =
false;
256 if (rABPos(ss, curIndex)==-1)
258 rABPos(ss, curIndex) = NbrABInit[ss];
260 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
261 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
265 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
281 wA.append(sl*omegai);
294 if (PA[wAID[
id]] == 0.0)
296 PA[wAID[id]] = wA[id];
300 PA[wAID[id]] += wA[id];
305 if (CA[wAID[
id]] == 0.0)
307 CA[wAID[id]] = -wA[id];
311 CA[wAID[id]] += -wA[id];
324 for (
label i=0; i<this->nSpecie_; i++)
326 Pa[0] += sC_[i]*
max(0.0,(PA[i]-CA[i]));
327 Ca[0] += sC_[i]*
max(0.0,-(PA[i]-CA[i]));
328 Pa[1] += sH_[i]*
max(0.0,(PA[i]-CA[i]));
329 Ca[1] += sH_[i]*
max(0.0,-(PA[i]-CA[i]));
330 Pa[2] += sO_[i]*
max(0.0,(PA[i]-CA[i]));
331 Ca[2] += sO_[i]*
max(0.0,-(PA[i]-CA[i]));
332 Pa[3] += sN_[i]*
max(0.0,(PA[i]-CA[i]));
333 Ca[3] += sN_[i]*
max(0.0,-(PA[i]-CA[i]));
339 label speciesNumber = 0;
340 List<bool> disabledSpecies(this->nSpecie_,
false);
344 for (
label i=0; i<this->nSpecie_; i++)
346 this->activeSpecies_[i] =
false;
350 const labelList& SIS(this->searchInitSet_);
364 scalar alphaTmp = (sC_[q]*
mag(PA[q]-CA[q])/Pa[0]);
365 if (alphaTmp > alphaA)
373 scalar alphaTmp = (sH_[q]*
mag(PA[q]-CA[q])/Pa[1]);
374 if (alphaTmp > alphaA)
382 scalar alphaTmp = (sO_[q]*
mag(PA[q]-CA[q])/Pa[2]);
383 if (alphaTmp > alphaA)
391 scalar alphaTmp = (sN_[q]*
mag(PA[q]-CA[q])/Pa[3]);
392 if (alphaTmp > alphaA)
397 if (alphaA > this->tolerance())
399 this->activeSpecies_[q] =
true;
420 if (Rvalue[SIS[i]] > Rmax)
422 Rmax = Rvalue[SIS[i]];
427 QStart.append(specID);
430 Rvalue[specID] = 1.0;
431 this->activeSpecies_[specID] =
true;
438 scalar Den =
max(PA[u],CA[u]);
441 for (
label v=0; v<NbrABInit[u]; v++)
443 label otherSpec = rABOtherSpec(u, v);
444 scalar rAB =
mag(rABNum(u, v))/Den;
450 scalar Rtemp = Rvalue[u]*rAB;
452 if (Rvalue[otherSpec] < Rtemp)
454 Rvalue[otherSpec] = Rtemp;
456 if (Rtemp >= this->tolerance())
459 if (!this->activeSpecies_[otherSpec])
461 this->activeSpecies_[otherSpec] =
true;
472 label NDisabledSpecies(this->nSpecie_-speciesNumber);
479 while(NDisabledSpecies > NGroupBased_)
486 forAll(disabledSpecies, i)
489 if (!this->activeSpecies_[i] && !disabledSpecies[i])
492 Rdisabled[nD] = Rvalue[i];
501 for (
label i=0; i<NGroupBased_; i++)
503 disabledSpecies[Rindex[tmpIndex[i]]] =
true;
505 NDisabledSpecies -= NGroupBased_;
511 for (
label v=0; v<NbrABInit[i]; v++)
516 forAll(this->chemistry_.reactions(), i)
520 scalar omegai = omegaV[i];
525 scalar sl = -R.
lhs()[
s].stoichCoeff;
527 bool alreadyDisabled(
false);
534 if (disabledSpecies[sj])
536 alreadyDisabled=
true;
544 if (disabledSpecies[sj])
546 alreadyDisabled=
true;
556 for (
label v=0; v<NbrABInit[ss]; v++)
558 rABNum(ss, v) += sl*omegai;
563 while(!usedIndex.empty())
566 if (deltaBi[curIndex])
569 deltaBi[curIndex] =
false;
570 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
579 scalar sl = R.
rhs()[
s].stoichCoeff;
581 bool alreadyDisabled(
false);
588 if (disabledSpecies[sj])
590 alreadyDisabled=
true;
598 if (disabledSpecies[sj])
600 alreadyDisabled=
true;
610 for (
label v=0; v<NbrABInit[ss]; v++)
612 rABNum(ss, v) += sl*omegai;
617 while(!usedIndex.empty())
620 if (deltaBi[curIndex])
622 deltaBi[curIndex] =
false;
623 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
632 label u = QStart[qi];
639 scalar Den =
max(PA[u],CA[u]);
642 for (
label v=0; v<NbrABInit[u]; v++)
644 label otherSpec = rABOtherSpec(u, v);
645 if (!disabledSpecies[otherSpec])
647 scalar rAB =
mag(rABNum(u, v))/Den;
653 scalar Rtemp = Rvalue[u]*rAB;
655 if (Rvalue[otherSpec] < Rtemp)
657 Rvalue[otherSpec] = Rtemp;
658 if (Rtemp >= this->tolerance())
661 if (!this->activeSpecies_[otherSpec])
663 this->activeSpecies_[otherSpec] =
true;
678 forAll(this->chemistry_.reactions(), i)
681 this->chemistry_.reactionsDisabled()[i] =
false;
685 if (!this->activeSpecies_[ss])
687 this->chemistry_.reactionsDisabled()[i] =
true;
691 if (!this->chemistry_.reactionsDisabled()[i])
696 if (!this->activeSpecies_[ss])
698 this->chemistry_.reactionsDisabled()[i] =
true;
705 this->NsSimp_ = speciesNumber;
706 scalarField& simplifiedC(this->chemistry_.simplifiedC());
707 simplifiedC.
setSize(this->NsSimp_+2);
710 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
713 for (
label i=0; i<this->nSpecie_; i++)
715 if (this->activeSpecies_[i])
718 simplifiedC[j] = c[i];
720 if (!this->chemistry_.active(i))
722 this->chemistry_.setActive(i);
730 simplifiedC[this->NsSimp_] =
T;
731 simplifiedC[this->NsSimp_+1] =
p;
732 this->chemistry_.setNsDAC(this->NsSimp_);
735 this->chemistry_.setNSpecie(this->NsSimp_);
A FIFO stack based on a singly-linked list.
virtual label nSpecie() const
The number of species.
DRGEP(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
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.
PtrList< volScalarField > & Y()
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
errorManipArg< error, int > exit(error &err, const int errNo=1)
virtual ~DRGEP()
Destructor.
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.
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.
A list that is sorted upon construction or when explicitly requested with the sort() method...
BasicChemistryModel< rhoReactionThermo > & chemistry
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
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 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
const word & name() const
Return the name of the element.
label nAtoms() const
Return the number of atoms of this element in the specie.
label readLabel(Istream &is)
void setSize(const label)
Reset size of List.
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
#define R(A, B, C, D, E, F, K, M)
void push(const T &a)
Push an element onto the stack.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
dimensioned< scalar > mag(const dimensioned< Type > &)
T pop()
Pop the bottom element off the stack.
An abstract class for methods of chemical mechanism reduction.
void partialSort(int M)
Partial sort the list (if changed after construction time)
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.