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.
#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
Return the components of the left hand side.
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...
BasicChemistryModel< rhoReactionThermo > & chemistry
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)
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
Return the components of the right hand side.