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_ = this->coeffsDict_.template lookup<label>(
"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>
124 scalarField& completeC(this->chemistry_.completeC());
127 for (
label i=0; i<this->nSpecie_; i++)
133 c1[this->nSpecie_] =
T;
134 c1[this->nSpecie_+1] =
p;
148 scalar pf, cf, pr, cr;
150 scalarField omegaV(this->chemistry_.reactions().size());
151 forAll(this->chemistry_.reactions(), i)
155 scalar omegai = this->chemistry_.omega
157 R, p,T, c1, li, pf, cf, lRef, pr, cr, rRef
170 scalar sl = -R.
lhs()[
s].stoichCoeff;
189 while(!usedIndex.empty())
192 if (deltaBi[curIndex])
195 deltaBi[curIndex] =
false;
197 if (rABPos(ss, curIndex)==-1)
199 rABPos(ss, curIndex) = NbrABInit[ss];
201 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
202 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
206 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
221 wA.append(sl*omegai);
230 scalar sl = R.
rhs()[
s].stoichCoeff;
249 while(!usedIndex.empty())
252 if (deltaBi[curIndex])
255 deltaBi[curIndex] =
false;
257 if (rABPos(ss, curIndex)==-1)
259 rABPos(ss, curIndex) = NbrABInit[ss];
261 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
262 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
266 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
282 wA.append(sl*omegai);
295 if (PA[wAID[
id]] == 0.0)
297 PA[wAID[id]] = wA[id];
301 PA[wAID[id]] += wA[id];
306 if (CA[wAID[
id]] == 0.0)
308 CA[wAID[id]] = -wA[id];
312 CA[wAID[id]] += -wA[id];
325 for (
label i=0; i<this->nSpecie_; i++)
327 Pa[0] += sC_[i]*
max(0.0,(PA[i]-CA[i]));
328 Ca[0] += sC_[i]*
max(0.0,-(PA[i]-CA[i]));
329 Pa[1] += sH_[i]*
max(0.0,(PA[i]-CA[i]));
330 Ca[1] += sH_[i]*
max(0.0,-(PA[i]-CA[i]));
331 Pa[2] += sO_[i]*
max(0.0,(PA[i]-CA[i]));
332 Ca[2] += sO_[i]*
max(0.0,-(PA[i]-CA[i]));
333 Pa[3] += sN_[i]*
max(0.0,(PA[i]-CA[i]));
334 Ca[3] += sN_[i]*
max(0.0,-(PA[i]-CA[i]));
340 label speciesNumber = 0;
341 List<bool> disabledSpecies(this->nSpecie_,
false);
345 for (
label i=0; i<this->nSpecie_; i++)
347 this->activeSpecies_[i] =
false;
351 const labelList& SIS(this->searchInitSet_);
365 scalar alphaTmp = (sC_[q]*
mag(PA[q]-CA[q])/Pa[0]);
366 if (alphaTmp > alphaA)
374 scalar alphaTmp = (sH_[q]*
mag(PA[q]-CA[q])/Pa[1]);
375 if (alphaTmp > alphaA)
383 scalar alphaTmp = (sO_[q]*
mag(PA[q]-CA[q])/Pa[2]);
384 if (alphaTmp > alphaA)
392 scalar alphaTmp = (sN_[q]*
mag(PA[q]-CA[q])/Pa[3]);
393 if (alphaTmp > alphaA)
398 if (alphaA > this->tolerance())
400 this->activeSpecies_[q] =
true;
421 if (Rvalue[SIS[i]] > Rmax)
423 Rmax = Rvalue[SIS[i]];
428 QStart.append(specID);
431 Rvalue[specID] = 1.0;
432 this->activeSpecies_[specID] =
true;
439 scalar Den =
max(PA[u],CA[u]);
442 for (
label v=0; v<NbrABInit[u]; v++)
444 label otherSpec = rABOtherSpec(u, v);
445 scalar rAB =
mag(rABNum(u, v))/Den;
451 scalar Rtemp = Rvalue[u]*rAB;
453 if (Rvalue[otherSpec] < Rtemp)
455 Rvalue[otherSpec] = Rtemp;
457 if (Rtemp >= this->tolerance())
460 if (!this->activeSpecies_[otherSpec])
462 this->activeSpecies_[otherSpec] =
true;
473 label NDisabledSpecies(this->nSpecie_-speciesNumber);
480 while(NDisabledSpecies > NGroupBased_)
487 forAll(disabledSpecies, i)
490 if (!this->activeSpecies_[i] && !disabledSpecies[i])
493 Rdisabled[nD] = Rvalue[i];
502 for (
label i=0; i<NGroupBased_; i++)
504 disabledSpecies[Rindex[tmpIndex[i]]] =
true;
506 NDisabledSpecies -= NGroupBased_;
512 for (
label v=0; v<NbrABInit[i]; v++)
517 forAll(this->chemistry_.reactions(), i)
521 scalar omegai = omegaV[i];
526 scalar sl = -R.
lhs()[
s].stoichCoeff;
528 bool alreadyDisabled(
false);
535 if (disabledSpecies[sj])
537 alreadyDisabled=
true;
545 if (disabledSpecies[sj])
547 alreadyDisabled=
true;
557 for (
label v=0; v<NbrABInit[ss]; v++)
559 rABNum(ss, v) += sl*omegai;
564 while(!usedIndex.empty())
567 if (deltaBi[curIndex])
570 deltaBi[curIndex] =
false;
571 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
580 scalar sl = R.
rhs()[
s].stoichCoeff;
582 bool alreadyDisabled(
false);
589 if (disabledSpecies[sj])
591 alreadyDisabled=
true;
599 if (disabledSpecies[sj])
601 alreadyDisabled=
true;
611 for (
label v=0; v<NbrABInit[ss]; v++)
613 rABNum(ss, v) += sl*omegai;
618 while(!usedIndex.empty())
621 if (deltaBi[curIndex])
623 deltaBi[curIndex] =
false;
624 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
633 label u = QStart[qi];
640 scalar Den =
max(PA[u],CA[u]);
643 for (
label v=0; v<NbrABInit[u]; v++)
645 label otherSpec = rABOtherSpec(u, v);
646 if (!disabledSpecies[otherSpec])
648 scalar rAB =
mag(rABNum(u, v))/Den;
654 scalar Rtemp = Rvalue[u]*rAB;
656 if (Rvalue[otherSpec] < Rtemp)
658 Rvalue[otherSpec] = Rtemp;
659 if (Rtemp >= this->tolerance())
662 if (!this->activeSpecies_[otherSpec])
664 this->activeSpecies_[otherSpec] =
true;
679 forAll(this->chemistry_.reactions(), i)
682 this->chemistry_.reactionsDisabled()[i] =
false;
686 if (!this->activeSpecies_[ss])
688 this->chemistry_.reactionsDisabled()[i] =
true;
692 if (!this->chemistry_.reactionsDisabled()[i])
697 if (!this->activeSpecies_[ss])
699 this->chemistry_.reactionsDisabled()[i] =
true;
706 this->NsSimp_ = speciesNumber;
707 scalarField& simplifiedC(this->chemistry_.simplifiedC());
708 simplifiedC.
setSize(this->NsSimp_+2);
711 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
714 for (
label i=0; i<this->nSpecie_; i++)
716 if (this->activeSpecies_[i])
719 simplifiedC[j] = c[i];
721 if (!this->chemistry_.active(i))
723 this->chemistry_.setActive(i);
731 simplifiedC[this->NsSimp_] =
T;
732 simplifiedC[this->NsSimp_+1] =
p;
733 this->chemistry_.setNsDAC(this->NsSimp_);
736 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.
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.
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.
A list that is sorted upon construction or when explicitly requested with the sort() method...
BasicChemistryModel< rhoReactionThermo > & chemistry
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.
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 setSize(const label)
Reset size of List.
#define R(A, B, C, D, E, F, K, M)
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, const label li)
Reduce the mechanism.
void push(const T &a)
Push an element onto the stack.
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.
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
An abstract class for methods of chemical mechanism reduction.
void partialSort(int M)
Partial sort the list (if changed after construction time)