31 template<
class CompType,
class ThermoType>
39 sC_(this->nSpecie_,0),
40 sH_(this->nSpecie_,0),
41 sO_(this->nSpecie_,0),
42 sN_(this->nSpecie_,0),
46 this->chemistry_.specieComp();
47 for (
label i=0; i<this->nSpecie_; i++)
52 forAll(curSpecieComposition, j)
55 curSpecieComposition[j];
56 if (curElement.
name() ==
"C")
58 sC_[i] = curElement.
nAtoms();
60 else if (curElement.
name() ==
"H")
62 sH_[i] = curElement.
nAtoms();
64 else if (curElement.
name() ==
"O")
66 sO_[i] = curElement.
nAtoms();
68 else if (curElement.
name() ==
"N")
70 sN_[i] = curElement.
nAtoms();
74 Info<<
"element not considered"<<
endl;
78 if (this->coeffsDict_.found(
"sortPart"))
80 sortPart_ =
readScalar(this->coeffsDict_.lookup(
"sortPart"));
87 template<
class CompType,
class ThermoType>
94 template<
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);
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;
161 scalar Acoeff = R.
lhs()[Ai].stoichCoeff;
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;
514 if (!this->activeSpecies_[ss])
516 this->chemistry_.reactionsDisabled()[i] =
true;
520 if (!this->chemistry_.reactionsDisabled()[i])
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_);
EFA(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
Extends chemistryModel by adding the TDAC method.
#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.
const List< specieCoeffs > & lhs() const
errorManipArg< error, int > exit(error &err, const int errNo=1)
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...
Ostream & endl(Ostream &os)
Add newline and flush stream.
A list that is sorted upon construction or when explicitly requested with the sort() method...
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
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.
A templated 2D rectangular m x n matrix of objects of <Type>.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
virtual ~EFA()
Destructor.
#define R(A, B, C, D, E, F, K, M)
psiChemistryModel & chemistry
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
dimensioned< scalar > mag(const dimensioned< Type > &)
void partialSort(int M, int start)
Partial sort the list (if changed after construction time)
An abstract class for methods of chemical mechanism reduction.
const List< specieCoeffs > & rhs() const