33 template<
class CompType,
class ThermoType>
41 searchInitSet_(this->coeffsDict_.subDict(
"initialSet").size()),
42 sC_(this->nSpecie_, 0),
43 sH_(this->nSpecie_, 0),
44 sO_(this->nSpecie_, 0),
45 sN_(this->nSpecie_, 0),
50 for (label i=0; i<
chemistry.nSpecie(); ++i)
54 searchInitSet_[j++] = i;
57 if (j<searchInitSet_.
size())
60 << searchInitSet_.
size() - j
61 <<
" species in the initial set is not in the mechanism " 71 for (label i=0; i<this->
nSpecie_; ++i)
78 if (curElement.name() ==
"C")
80 sC_[i] = curElement.nAtoms();
82 else if (curElement.name() ==
"H")
84 sH_[i] = curElement.nAtoms();
86 else if (curElement.name() ==
"O")
88 sO_[i] = curElement.nAtoms();
90 else if (curElement.name() ==
"N")
92 sN_[i] = curElement.nAtoms();
96 Info<<
"element not considered"<<
endl;
105 template<
class CompType,
class ThermoType>
112 template<
class CompType,
class ThermoType>
121 scalarField& completeC(this->chemistry_.completeC());
124 for (label i=0; i<this->nSpecie_; ++i)
130 c1[this->nSpecie_] =
T;
131 c1[this->nSpecie_+1] =
p;
145 scalar pf, cf, pr, cr;
147 scalarField omegaV(this->chemistry_.reactions().size());
148 forAll(this->chemistry_.reactions(), i)
153 scalar omegai = this->chemistry_.omega
155 R,
c1,
T,
p, pf, cf, lRef, pr, cr, rRef
167 label ss =
R.lhs()[
s].index;
168 scalar sl = -
R.lhs()[
s].stoichCoeff;
173 label sj =
R.lhs()[j].index;
179 label sj =
R.rhs()[j].index;
187 while (!usedIndex.empty())
189 label curIndex = usedIndex.
pop();
190 if (deltaBi[curIndex])
193 deltaBi[curIndex] =
false;
195 if (rABPos(ss, curIndex) == -1)
197 rABPos(ss, curIndex) = NbrABInit[ss];
199 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
200 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
204 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
219 wA.append(sl*omegai);
227 label ss =
R.rhs()[
s].index;
228 scalar sl =
R.rhs()[
s].stoichCoeff;
233 label sj =
R.lhs()[j].index;
239 label sj =
R.rhs()[j].index;
247 while (!usedIndex.empty())
249 label curIndex = usedIndex.
pop();
250 if (deltaBi[curIndex])
253 deltaBi[curIndex] =
false;
255 if (rABPos(ss, curIndex) == -1)
257 rABPos(ss, curIndex) = NbrABInit[ss];
259 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
260 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
264 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
280 wA.append(sl*omegai);
293 if (PA[wAID[
id]] == 0.0)
295 PA[wAID[id]] = wA[id];
299 PA[wAID[id]] += wA[id];
304 if (CA[wAID[
id]] == 0.0)
306 CA[wAID[id]] = -wA[id];
310 CA[wAID[id]] += -wA[id];
323 for (label i=0; i<this->nSpecie_; ++i)
325 Pa[0] += sC_[i]*
max(0.0, (PA[i]-CA[i]));
326 Ca[0] += sC_[i]*
max(0.0,-(PA[i]-CA[i]));
327 Pa[1] += sH_[i]*
max(0.0, (PA[i]-CA[i]));
328 Ca[1] += sH_[i]*
max(0.0,-(PA[i]-CA[i]));
329 Pa[2] += sO_[i]*
max(0.0, (PA[i]-CA[i]));
330 Ca[2] += sO_[i]*
max(0.0,-(PA[i]-CA[i]));
331 Pa[3] += sN_[i]*
max(0.0, (PA[i]-CA[i]));
332 Ca[3] += sN_[i]*
max(0.0,-(PA[i]-CA[i]));
338 label speciesNumber = 0;
339 List<bool> disabledSpecies(this->nSpecie_,
false);
343 for (label i=0; i<this->nSpecie_; ++i)
345 this->activeSpecies_[i] =
false;
349 const labelList& SIS(this->searchInitSet_);
355 for (label i=0; i<SIS.
size(); ++i)
363 scalar alphaTmp = (sC_[q]*
mag(PA[q]-CA[q])/Pa[0]);
364 if (alphaTmp > alphaA)
372 scalar alphaTmp = (sH_[q]*
mag(PA[q]-CA[q])/Pa[1]);
373 if (alphaTmp > alphaA)
381 scalar alphaTmp = (sO_[q]*
mag(PA[q]-CA[q])/Pa[2]);
382 if (alphaTmp > alphaA)
390 scalar alphaTmp = (sN_[q]*
mag(PA[q]-CA[q])/Pa[3]);
391 if (alphaTmp > alphaA)
396 if (alphaA > this->tolerance())
398 this->activeSpecies_[q] =
true;
417 for (
const label sis : SIS)
419 if (Rvalue[sis] > Rmax)
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];
524 label ss =
R.lhs()[
s].index;
525 scalar sl = -
R.lhs()[
s].stoichCoeff;
527 bool alreadyDisabled(
false);
531 label sj =
R.lhs()[j].index;
534 if (disabledSpecies[sj])
536 alreadyDisabled=
true;
541 label sj =
R.rhs()[j].index;
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())
565 label curIndex = usedIndex.
pop();
566 if (deltaBi[curIndex])
569 deltaBi[curIndex] =
false;
570 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
578 label ss =
R.rhs()[
s].index;
579 scalar sl =
R.rhs()[
s].stoichCoeff;
581 bool alreadyDisabled(
false);
585 label sj =
R.lhs()[j].index;
588 if (disabledSpecies[sj])
590 alreadyDisabled =
true;
595 label sj =
R.rhs()[j].index;
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())
619 label curIndex = usedIndex.
pop();
620 if (deltaBi[curIndex])
622 deltaBi[curIndex] =
false;
623 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
630 for (
const label q : QStart)
638 scalar Den =
max(PA[u],CA[u]);
641 for (label v=0; v<NbrABInit[u]; ++v)
643 label otherSpec = rABOtherSpec(u, v);
644 if (!disabledSpecies[otherSpec])
646 scalar rAB =
mag(rABNum(u, v))/Den;
652 scalar Rtemp = Rvalue[u]*rAB;
654 if (Rvalue[otherSpec] < Rtemp)
656 Rvalue[otherSpec] = Rtemp;
657 if (Rtemp >= this->tolerance())
660 if (!this->activeSpecies_[otherSpec])
662 this->activeSpecies_[otherSpec] =
true;
677 forAll(this->chemistry_.reactions(), i)
680 this->chemistry_.reactionsDisabled()[i] =
false;
683 label ss =
R.lhs()[
s].index;
684 if (!this->activeSpecies_[ss])
686 this->chemistry_.reactionsDisabled()[i] =
true;
690 if (!this->chemistry_.reactionsDisabled()[i])
694 label ss =
R.rhs()[
s].index;
695 if (!this->activeSpecies_[ss])
697 this->chemistry_.reactionsDisabled()[i] =
true;
704 this->NsSimp_ = speciesNumber;
705 scalarField& simplifiedC(this->chemistry_.simplifiedC());
706 simplifiedC.
setSize(this->NsSimp_+2);
709 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
712 for (label i=0; i<this->nSpecie_; ++i)
714 if (this->activeSpecies_[i])
717 simplifiedC[j] =
c[i];
719 if (!this->chemistry_.active(i))
721 this->chemistry_.setActive(i);
729 simplifiedC[this->NsSimp_] =
T;
730 simplifiedC[this->NsSimp_+1] =
p;
731 this->chemistry_.setNsDAC(this->NsSimp_);
734 this->chemistry_.setNSpecie(this->NsSimp_);
A FIFO stack based on a singly-linked list.
DRGEP(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
void size(const label n)
Older name for setAddressableSize.
TDACChemistryModel< CompType, ThermoType > & chemistry_
errorManipArg< error, int > exit(error &err, const int errNo=1)
virtual ~DRGEP()
Destructor.
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.
A list that is sorted upon construction or when explicitly requested with the sort() method...
void setSize(const label n)
Same as resize()
BasicChemistryModel< psiReactionThermo > & chemistry
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...
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()
const label nSpecie_
Number of species.
List< List< specieElement > > & specieComp()
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry if present, and assign to T val. FatalIOError if it is found and the number of tokens i...
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
const dictionary coeffsDict_
Dictionary that store the algorithm data.
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)
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)
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
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))
void partialSort(int M)
Partial sort the list (if changed after construction time)
static constexpr const zero Zero
Global zero (0)