33template<
class CompType,
class ThermoType>
41 sC_(this->nSpecie_, Zero),
42 sH_(this->nSpecie_, Zero),
43 sO_(this->nSpecie_, Zero),
44 sN_(this->nSpecie_, Zero),
49 for (label i=0; i<this->
nSpecie_; i++)
54 forAll(curSpecieComposition, j)
57 curSpecieComposition[j];
58 if (curElement.
name() ==
"C")
60 sC_[i] = curElement.
nAtoms();
62 else if (curElement.
name() ==
"H")
64 sH_[i] = curElement.
nAtoms();
66 else if (curElement.
name() ==
"O")
68 sO_[i] = curElement.
nAtoms();
70 else if (curElement.
name() ==
"N")
72 sN_[i] = curElement.
nAtoms();
76 Info<<
"element not considered"<<
endl;
87template<
class CompType,
class ThermoType>
94template<
class CompType,
class ThermoType>
102 scalarField& completeC(this->chemistry_.completeC());
105 for (label i=0; i<this->nSpecie_; i++)
111 c1[this->nSpecie_] =
T;
112 c1[this->nSpecie_+1] =
p;
124 scalar CFlux(0.0), HFlux(0.0), OFlux(0.0), NFlux(0.0);
130 scalar pf, cf, pr, cr;
132 forAll(this->chemistry_.reactions(), i)
136 this->chemistry_.omega
138 R, c1,
T,
p, pf, cf, lRef, pr, cr, rRef
140 scalar fr =
mag(pf*cf)+
mag(pr*cr);
141 scalar NCi(0.0),NHi(0.0),NOi(0.0),NNi(0.0);
144 label curIndex =
R.lhs()[
s].index;
145 scalar stoicCoeff =
R.lhs()[
s].stoichCoeff;
146 NCi += sC_[curIndex]*stoicCoeff;
147 NHi += sH_[curIndex]*stoicCoeff;
148 NOi += sO_[curIndex]*stoicCoeff;
149 NNi += sN_[curIndex]*stoicCoeff;
160 label
A =
R.lhs()[Ai].index;
161 scalar Acoeff =
R.lhs()[Ai].stoichCoeff;
165 label
B =
R.rhs()[Bi].index;
166 scalar Bcoeff =
R.rhs()[Bi].stoichCoeff;
168 if (rABPos(
B,
A)==-1)
170 label otherS = rABPos(
A,
B);
173 rABPos(
A,
B) = NbrABInit[
A]++;
174 otherS = rABPos(
A,
B);
176 rABOtherSpec(
A, otherS) =
B;
180 CFluxAB(
A, otherS) +=
181 fr*Acoeff*sC_[
A]*Bcoeff*sC_[
B]/NCi;
185 HFluxAB(
A, otherS) +=
186 fr*Acoeff*sH_[
A]*Bcoeff*sH_[
B]/NHi;
190 OFluxAB(
A, otherS) +=
191 fr*Acoeff*sO_[
A]*Bcoeff*sO_[
B]/NOi;
195 NFluxAB(
A, otherS) +=
196 fr*Acoeff*sN_[
A]*Bcoeff*sN_[
B]/NNi;
203 label otherS = rABPos(
B,
A);
206 CFluxAB(
B, otherS) +=
207 fr*Acoeff*sC_[
A]*Bcoeff*sC_[
B]/NCi;
211 HFluxAB(
B, otherS) +=
212 fr*Acoeff*sH_[
A]*Bcoeff*sH_[
B]/NHi;
216 OFluxAB(
B, otherS) +=
217 fr*Acoeff*sO_[
A]*Bcoeff*sO_[
B]/NOi;
221 NFluxAB(
B, otherS) +=
222 fr*Acoeff*sN_[
A]*Bcoeff*sN_[
B]/NNi;
227 CFlux += fr*Acoeff*sC_[
A]*Bcoeff*sC_[
B]/NCi;
231 HFlux += fr*Acoeff*sH_[
A]*Bcoeff*sH_[
B]/NHi;
235 OFlux += fr*Acoeff*sO_[
A]*Bcoeff*sO_[
B]/NOi;
239 NFlux += fr*Acoeff*sN_[
A]*Bcoeff*sN_[
B]/NNi;
247 label speciesNumber = 0;
248 for (label i=0; i<this->nSpecie_; i++)
250 this->activeSpecies_[i] =
false;
259 for (
int A=0;
A<this->nSpecie_ ;
A++)
261 for (
int j=0; j<NbrABInit[
A]; j++)
263 label pairIndex = nP++;
264 pairsFlux[pairIndex] = 0.0;
265 label
B = rABOtherSpec(
A, j);
266 pairsFlux[pairIndex] += CFluxAB(
A, j);
267 source[pairIndex] =
A;
277 scalar threshold((1-this->tolerance())*CFlux);
279 label nbToSort(
static_cast<label
> (nbPairs*sortPart_));
280 nbToSort =
max(nbToSort,1);
282 bool cumRespected(
false);
285 if (startPoint >= nbPairs)
288 <<
"startPoint outside number of pairs without reaching"
293 label nbi =
min(nbToSort, nbPairs-startPoint);
296 for (
int i=0; i<nbi; i++)
298 cumFlux += pairsFlux[idx[startPoint+i]];
299 if (!this->activeSpecies_[source[idx[startPoint+i]]])
301 this->activeSpecies_[source[idx[startPoint+i]]] =
true;
304 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
306 this->activeSpecies_[sink[idx[startPoint+i]]] =
true;
309 if (cumFlux >= threshold)
315 startPoint += nbToSort;
325 for (
int A=0;
A<this->nSpecie_ ;
A++)
327 for (
int j=0; j<NbrABInit[
A]; j++)
329 label pairIndex = nP++;
330 pairsFlux[pairIndex] = 0.0;
331 label
B = rABOtherSpec(
A, j);
332 pairsFlux[pairIndex] += HFluxAB(
A, j);
333 source[pairIndex] =
A;
339 scalar threshold((1-this->tolerance())*HFlux);
341 label nbToSort(
static_cast<label
> (nbPairs*sortPart_));
342 nbToSort =
max(nbToSort,1);
344 bool cumRespected(
false);
347 if (startPoint >= nbPairs)
350 <<
"startPoint outside number of pairs without reaching"
355 label nbi =
min(nbToSort, nbPairs-startPoint);
358 for (
int i=0; i<nbi; i++)
360 cumFlux += pairsFlux[idx[startPoint+i]];
362 if (!this->activeSpecies_[source[idx[startPoint+i]]])
364 this->activeSpecies_[source[idx[startPoint+i]]] =
true;
367 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
369 this->activeSpecies_[sink[idx[startPoint+i]]] =
true;
372 if (cumFlux >= threshold)
378 startPoint += nbToSort;
388 for (
int A=0;
A<this->nSpecie_ ;
A++)
390 for (
int j=0; j<NbrABInit[
A]; j++)
392 label pairIndex = nP++;
393 pairsFlux[pairIndex] = 0.0;
394 label
B = rABOtherSpec(
A, j);
395 pairsFlux[pairIndex] += OFluxAB(
A, j);
396 source[pairIndex] =
A;
401 scalar threshold((1-this->tolerance())*OFlux);
403 label nbToSort(
static_cast<label
> (nbPairs*sortPart_));
404 nbToSort =
max(nbToSort,1);
406 bool cumRespected(
false);
409 if (startPoint >= nbPairs)
412 <<
"startPoint outside number of pairs without reaching"
417 label nbi =
min(nbToSort, nbPairs-startPoint);
420 for (
int i=0; i<nbi; i++)
422 cumFlux += pairsFlux[idx[startPoint+i]];
424 if (!this->activeSpecies_[source[idx[startPoint+i]]])
426 this->activeSpecies_[source[idx[startPoint+i]]] =
true;
429 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
431 this->activeSpecies_[sink[idx[startPoint+i]]] =
true;
434 if (cumFlux >= threshold)
440 startPoint += nbToSort;
450 for (
int A=0;
A<this->nSpecie_ ;
A++)
452 for (
int j=0; j<NbrABInit[
A]; j++)
454 label pairIndex = nP++;
455 pairsFlux[pairIndex] = 0.0;
456 label
B = rABOtherSpec(
A, j);
457 pairsFlux[pairIndex] += NFluxAB(
A, j);
458 source[pairIndex] =
A;
463 scalar threshold((1-this->tolerance())*NFlux);
465 label nbToSort(
static_cast<label
> (nbPairs*sortPart_));
466 nbToSort =
max(nbToSort,1);
468 bool cumRespected(
false);
471 if (startPoint >= nbPairs)
474 <<
"startPoint outside number of pairs without reaching"
479 label nbi =
min(nbToSort, nbPairs-startPoint);
482 for (
int i=0; i<nbi; i++)
484 cumFlux += pairsFlux[idx[startPoint+i]];
486 if (!this->activeSpecies_[source[idx[startPoint+i]]])
488 this->activeSpecies_[source[idx[startPoint+i]]] =
true;
491 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
493 this->activeSpecies_[sink[idx[startPoint+i]]] =
true;
496 if (cumFlux >= threshold)
502 startPoint += nbToSort;
507 forAll(this->chemistry_.reactions(), i)
510 this->chemistry_.reactionsDisabled()[i] =
false;
513 label ss =
R.lhs()[
s].index;
514 if (!this->activeSpecies_[ss])
516 this->chemistry_.reactionsDisabled()[i] =
true;
520 if (!this->chemistry_.reactionsDisabled()[i])
524 label ss =
R.rhs()[
s].index;
525 if (!this->activeSpecies_[ss])
527 this->chemistry_.reactionsDisabled()[i] =
true;
534 this->NsSimp_ = speciesNumber;
535 scalarField& simplifiedC(this->chemistry_.simplifiedC());
536 simplifiedC.
setSize(this->NsSimp_+2);
539 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
542 for (label i=0; i<this->nSpecie_; i++)
544 if (this->activeSpecies_[i])
547 simplifiedC[j] = c[i];
549 if (!this->chemistry_.active(i))
551 this->chemistry_.setActive(i);
559 simplifiedC[this->NsSimp_] =
T;
560 simplifiedC[this->NsSimp_+1] =
p;
561 this->chemistry_.setNsDAC(this->NsSimp_);
564 this->chemistry_.setNSpecie(this->NsSimp_);
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
#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.
void setSize(const label n)
Same as resize()
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,...
A list that is sorted upon construction or when explicitly requested with the sort() method.
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
void partialSort(int M, int start)
Partial sort the list (if changed after construction time)
Extends StandardChemistryModel by adding the TDAC method.
List< List< specieElement > > & specieComp()
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.
virtual ~EFA()
Destructor.
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
label nAtoms() const
Return the number of atoms of this element in the specie.
const word & name() const
Return the name of the element.
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)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
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.