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_ = this->coeffsDict_.template lookup<scalar>(
"sortPart");
87 template<
class CompType,
class ThermoType>
94 template<
class CompType,
class ThermoType>
103 scalarField& completeC(this->chemistry_.completeC());
106 for (
label i=0; i<this->nSpecie_; i++)
112 c1[this->nSpecie_] =
T;
113 c1[this->nSpecie_+1] =
p;
125 scalar CFlux(0.0), HFlux(0.0), OFlux(0.0), NFlux(0.0);
131 scalar pf, cf, pr, cr;
133 forAll(this->chemistry_.reactions(), i)
138 this->chemistry_.omega
140 R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef
142 scalar fr =
mag(pf*cf)+
mag(pr*cr);
143 scalar NCi(0.0),NHi(0.0),NOi(0.0),NNi(0.0);
147 scalar stoicCoeff = R.
lhs()[
s].stoichCoeff;
148 NCi += sC_[curIndex]*stoicCoeff;
149 NHi += sH_[curIndex]*stoicCoeff;
150 NOi += sO_[curIndex]*stoicCoeff;
151 NNi += sN_[curIndex]*stoicCoeff;
163 scalar Acoeff = R.
lhs()[Ai].stoichCoeff;
168 scalar Bcoeff = R.
rhs()[Bi].stoichCoeff;
170 if (rABPos(B, A)==-1)
172 label otherS = rABPos(A, B);
175 rABPos(A, B) = NbrABInit[A]++;
176 otherS = rABPos(A, B);
178 rABOtherSpec(A, otherS) = B;
182 CFluxAB(A, otherS) +=
183 fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
187 HFluxAB(A, otherS) +=
188 fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
192 OFluxAB(A, otherS) +=
193 fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
197 NFluxAB(A, otherS) +=
198 fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
205 label otherS = rABPos(B, A);
208 CFluxAB(B, otherS) +=
209 fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
213 HFluxAB(B, otherS) +=
214 fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
218 OFluxAB(B, otherS) +=
219 fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
223 NFluxAB(B, otherS) +=
224 fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
229 CFlux += fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
233 HFlux += fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
237 OFlux += fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
241 NFlux += fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
249 label speciesNumber = 0;
250 for (
label i=0; i<this->nSpecie_; i++)
252 this->activeSpecies_[i] =
false;
261 for (
int A=0; A<this->nSpecie_ ; A++)
263 for (
int j=0; j<NbrABInit[A]; j++)
265 label pairIndex = nP++;
266 pairsFlux[pairIndex] = 0.0;
267 label B = rABOtherSpec(A, j);
268 pairsFlux[pairIndex] += CFluxAB(A, j);
269 source[pairIndex] = A;
279 scalar threshold((1-this->tolerance())*CFlux);
281 label nbToSort(static_cast<label> (nbPairs*sortPart_));
282 nbToSort =
max(nbToSort,1);
284 bool cumRespected(
false);
287 if (startPoint >= nbPairs)
290 <<
"startPoint outside number of pairs without reaching" 295 label nbi =
min(nbToSort, nbPairs-startPoint);
298 for (
int i=0; i<nbi; i++)
300 cumFlux += pairsFlux[idx[startPoint+i]];
301 if (!this->activeSpecies_[source[idx[startPoint+i]]])
303 this->activeSpecies_[source[idx[startPoint+i]]] =
true;
306 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
308 this->activeSpecies_[sink[idx[startPoint+i]]] =
true;
311 if (cumFlux >= threshold)
317 startPoint += nbToSort;
327 for (
int A=0; A<this->nSpecie_ ; A++)
329 for (
int j=0; j<NbrABInit[A]; j++)
331 label pairIndex = nP++;
332 pairsFlux[pairIndex] = 0.0;
333 label B = rABOtherSpec(A, j);
334 pairsFlux[pairIndex] += HFluxAB(A, j);
335 source[pairIndex] = A;
341 scalar threshold((1-this->tolerance())*HFlux);
343 label nbToSort(static_cast<label> (nbPairs*sortPart_));
344 nbToSort =
max(nbToSort,1);
346 bool cumRespected(
false);
349 if (startPoint >= nbPairs)
352 <<
"startPoint outside number of pairs without reaching" 357 label nbi =
min(nbToSort, nbPairs-startPoint);
360 for (
int i=0; i<nbi; i++)
362 cumFlux += pairsFlux[idx[startPoint+i]];
364 if (!this->activeSpecies_[source[idx[startPoint+i]]])
366 this->activeSpecies_[source[idx[startPoint+i]]] =
true;
369 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
371 this->activeSpecies_[sink[idx[startPoint+i]]] =
true;
374 if (cumFlux >= threshold)
380 startPoint += nbToSort;
390 for (
int A=0; A<this->nSpecie_ ; A++)
392 for (
int j=0; j<NbrABInit[A]; j++)
394 label pairIndex = nP++;
395 pairsFlux[pairIndex] = 0.0;
396 label B = rABOtherSpec(A, j);
397 pairsFlux[pairIndex] += OFluxAB(A, j);
398 source[pairIndex] = A;
403 scalar threshold((1-this->tolerance())*OFlux);
405 label nbToSort(static_cast<label> (nbPairs*sortPart_));
406 nbToSort =
max(nbToSort,1);
408 bool cumRespected(
false);
411 if (startPoint >= nbPairs)
414 <<
"startPoint outside number of pairs without reaching" 419 label nbi =
min(nbToSort, nbPairs-startPoint);
422 for (
int i=0; i<nbi; i++)
424 cumFlux += pairsFlux[idx[startPoint+i]];
426 if (!this->activeSpecies_[source[idx[startPoint+i]]])
428 this->activeSpecies_[source[idx[startPoint+i]]] =
true;
431 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
433 this->activeSpecies_[sink[idx[startPoint+i]]] =
true;
436 if (cumFlux >= threshold)
442 startPoint += nbToSort;
452 for (
int A=0; A<this->nSpecie_ ; A++)
454 for (
int j=0; j<NbrABInit[A]; j++)
456 label pairIndex = nP++;
457 pairsFlux[pairIndex] = 0.0;
458 label B = rABOtherSpec(A, j);
459 pairsFlux[pairIndex] += NFluxAB(A, j);
460 source[pairIndex] = A;
465 scalar threshold((1-this->tolerance())*NFlux);
467 label nbToSort(static_cast<label> (nbPairs*sortPart_));
468 nbToSort =
max(nbToSort,1);
470 bool cumRespected(
false);
473 if (startPoint >= nbPairs)
476 <<
"startPoint outside number of pairs without reaching" 481 label nbi =
min(nbToSort, nbPairs-startPoint);
484 for (
int i=0; i<nbi; i++)
486 cumFlux += pairsFlux[idx[startPoint+i]];
488 if (!this->activeSpecies_[source[idx[startPoint+i]]])
490 this->activeSpecies_[source[idx[startPoint+i]]] =
true;
493 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
495 this->activeSpecies_[sink[idx[startPoint+i]]] =
true;
498 if (cumFlux >= threshold)
504 startPoint += nbToSort;
509 forAll(this->chemistry_.reactions(), i)
512 this->chemistry_.reactionsDisabled()[i] =
false;
516 if (!this->activeSpecies_[ss])
518 this->chemistry_.reactionsDisabled()[i] =
true;
522 if (!this->chemistry_.reactionsDisabled()[i])
527 if (!this->activeSpecies_[ss])
529 this->chemistry_.reactionsDisabled()[i] =
true;
536 this->NsSimp_ = speciesNumber;
537 scalarField& simplifiedC(this->chemistry_.simplifiedC());
538 simplifiedC.
setSize(this->NsSimp_+2);
541 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
544 for (
label i=0; i<this->nSpecie_; i++)
546 if (this->activeSpecies_[i])
549 simplifiedC[j] = c[i];
551 if (!this->chemistry_.active(i))
553 this->chemistry_.setActive(i);
561 simplifiedC[this->NsSimp_] =
T;
562 simplifiedC[this->NsSimp_+1] =
p;
563 this->chemistry_.setNsDAC(this->NsSimp_);
566 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.
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.
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...
Ostream & endl(Ostream &os)
Add newline and flush stream.
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
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.
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, const label li)
Reduce the mechanism.
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)
dimensioned< scalar > mag(const dimensioned< Type > &)
void partialSort(int M, int start)
Partial sort the list (if changed after construction time)
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
An abstract class for methods of chemical mechanism reduction.