32 template<
class CompType,
class ThermoType>
40 searchInitSet_(this->coeffsDict_.subDict(
"initialSet").size())
46 for (label i=0; i<
chemistry.nSpecie(); i++)
50 searchInitSet_[j++] = i;
54 if (j<searchInitSet_.
size())
57 << searchInitSet_.
size()-j
58 <<
" species in the initial set is not in the mechanism " 67 template<
class CompType,
class ThermoType>
74 template<
class CompType,
class ThermoType>
82 scalarField& completeC(this->chemistry_.completeC());
85 for (label i=0; i<this->nSpecie_; i++)
91 c1[this->nSpecie_] =
T;
92 c1[this->nSpecie_+1] =
p;
107 scalar pf, cf, pr, cr;
109 forAll(this->chemistry_.reactions(), i)
113 scalar omegai = this->chemistry_.omega
115 R,
c1,
T,
p, pf, cf, lRef, pr, cr, rRef
127 label ss =
R.lhs()[
s].index;
128 scalar sl = -
R.lhs()[
s].stoichCoeff;
141 wA.append(sl*omegai);
147 label ss =
R.rhs()[
s].index;
148 scalar sl =
R.rhs()[
s].stoichCoeff;
161 wA.append(sl*omegai);
169 label curID = wAID[id];
170 scalar curwA = wA[id];
171 List<bool> deltaBi(this->nSpecie_,
false);
172 FIFOStack<label> usedIndex;
175 label sj =
R.lhs()[j].index;
181 label sj =
R.rhs()[j].index;
186 deltaBi[curID] =
false;
188 while(!usedIndex.empty())
190 label curIndex = usedIndex.pop();
192 if (deltaBi[curIndex])
194 deltaBi[curIndex] =
false;
195 if (rABPos(curID, curIndex)==-1)
197 rABPos(curID, curIndex) = NbrABInit[curID];
198 rABOtherSpec(curID, NbrABInit[curID]) = curIndex;
201 PAB(curID, NbrABInit[curID]) = curwA;
205 CAB(curID, NbrABInit[curID]) = -curwA;
213 PAB(curID, rABPos(curID, curIndex)) += curwA;
217 CAB(curID, rABPos(curID, curIndex)) += -curwA;
227 if (PA[curID] == 0.0)
238 if (CA[curID] == 0.0)
258 RectangularMatrix<scalar> PAB2nd(this->nSpecie_, this->nSpecie_,
Zero);
259 RectangularMatrix<scalar> CAB2nd(this->nSpecie_, this->nSpecie_,
Zero);
262 Field<label> NbrABInit2nd(this->nSpecie_,
Zero);
265 RectangularMatrix<label> rABPos2nd(this->nSpecie_, this->nSpecie_, -1);
268 RectangularMatrix<label> rABOtherSpec2nd
270 this->nSpecie_, this->nSpecie_, -1
275 for (
int i=0; i<NbrABInit[
A]; i++)
277 label ri = rABOtherSpec(
A, i);
278 scalar maxPACA =
max(PA[ri],CA[ri]);
279 if (maxPACA > VSMALL)
281 for (
int j=0; j<NbrABInit[ri]; j++)
283 label
B = rABOtherSpec(ri, j);
286 if (rABPos2nd(
A,
B)==-1)
288 rABPos2nd(
A,
B) = NbrABInit2nd[
A]++;
289 rABOtherSpec2nd(
A, rABPos2nd(
A,
B)) =
B;
290 PAB2nd(
A, rABPos2nd(
A,
B)) =
291 PAB(
A, i)*PAB(ri, j)/maxPACA;
292 CAB2nd(
A, rABPos2nd(
A,
B)) =
293 CAB(
A, i)*CAB(ri, j)/maxPACA;
297 PAB2nd(
A, rABPos2nd(
A,
B)) +=
298 PAB(
A, i)*PAB(ri, j)/maxPACA;
299 CAB2nd(
A, rABPos2nd(
A,
B)) +=
300 CAB(
A, i)*CAB(ri, j)/maxPACA;
309 label speciesNumber = 0;
313 for (label i=0; i<this->nSpecie_; i++)
315 this->activeSpecies_[i] =
false;
319 const labelList& SIS(this->searchInitSet_);
322 for (label i=0; i<SIS.size(); i++)
325 this->activeSpecies_[q] =
true;
334 scalar Den =
max(PA[u],CA[u]);
339 for (label v=0; v<NbrABInit[u]; v++)
341 label otherSpec = rABOtherSpec(u, v);
342 scalar rAB = (PAB(u, v)+CAB(u, v))/Den;
343 label id2nd = rABPos2nd(u, otherSpec);
346 rAB += (PAB2nd(u, id2nd)+CAB2nd(u, id2nd))/Den;
351 rAB >= this->tolerance()
352 && !this->activeSpecies_[otherSpec]
356 this->activeSpecies_[otherSpec] =
true;
362 for (label v=0; v<NbrABInit2nd[u]; v++)
364 label otherSpec = rABOtherSpec2nd(u, v);
365 scalar rAB = (PAB2nd(u, v)+CAB2nd(u, v))/Den;
369 rAB >= this->tolerance()
370 && !this->activeSpecies_[otherSpec]
374 this->activeSpecies_[otherSpec] =
true;
382 forAll(this->chemistry_.reactions(), i)
384 const Reaction<ThermoType>&
R = this->chemistry_.reactions()[i];
385 this->chemistry_.reactionsDisabled()[i] =
false;
389 label ss =
R.lhs()[
s].index;
390 if (!this->activeSpecies_[ss])
392 this->chemistry_.reactionsDisabled()[i] =
true;
396 if (!this->chemistry_.reactionsDisabled()[i])
400 label ss =
R.rhs()[
s].index;
401 if (!this->activeSpecies_[ss])
403 this->chemistry_.reactionsDisabled()[i] =
true;
410 this->NsSimp_ = speciesNumber;
411 scalarField& simplifiedC(this->chemistry_.simplifiedC());
412 simplifiedC.
setSize(this->NsSimp_+2);
413 DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
414 s2c.setSize(this->NsSimp_);
415 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
418 for (label i=0; i<this->nSpecie_; i++)
420 if (this->activeSpecies_[i])
423 simplifiedC[j] =
c[i];
425 if (!this->chemistry_.active(i))
427 this->chemistry_.setActive(i);
435 simplifiedC[this->NsSimp_] =
T;
436 simplifiedC[this->NsSimp_+1] =
p;
437 this->chemistry_.setNsDAC(this->NsSimp_);
440 this->chemistry_.setNSpecie(this->NsSimp_);
void size(const label n)
Older name for setAddressableSize.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
BasicChemistryModel< psiReactionThermo > & chemistry
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry (const access) with the given keyword.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
#define forAll(list, i)
Loop across all elements in list.
void setSize(const label n)
Alias for resize()
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
errorManip< error > abort(error &err)
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
virtual ~PFA()
Destructor.
const dictionary coeffsDict_
Dictionary that store the algorithm data.
#define R(A, B, C, D, E, F, K, M)
const dimensionedScalar c
Speed of light in a vacuum.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
An abstract class for methods of chemical mechanism reduction.
List< label > labelList
A List of labels.
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))
PFA(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
static constexpr const zero Zero
Global zero (0)