36 template<
class ReactionThermo,
class ThermoType>
50 ||
fv::localEulerDdt::enabled(this->
mesh())
53 NsDAC_(this->nSpecie_),
54 completeC_(this->nSpecie_, 0),
55 reactionsDisabled_(this->reactions_.size(), false),
56 specieComp_(this->nSpecie_),
57 completeToSimplifiedIndex_(this->nSpecie_, -1),
58 simplifiedToCompleteIndex_(this->nSpecie_),
63 thermo.phasePropertyName(
"TabulationResults"),
86 specieComp_[i] = (specCompPtr.
ref())[this->
Y()[i].member()];
97 if (mechRed_->active())
114 this->
Y()[i].writeOpt(IOobject::NO_WRITE);
127 cpuReduceFile_ = logFile(
"cpu_reduce.out");
128 nActiveSpeciesFile_ = logFile(
"nActiveSpecies.out");
131 if (tabulation_->log())
133 cpuAddFile_ = logFile(
"cpu_add.out");
134 cpuGrowFile_ = logFile(
"cpu_grow.out");
135 cpuRetrieveFile_ = logFile(
"cpu_retrieve.out");
138 if (mechRed_->log() || tabulation_->log())
140 cpuSolveFile_ = logFile(
"cpu_solve.out");
147 template<
class ReactionThermo,
class ThermoType>
154 template<
class ReactionThermo,
class ThermoType>
163 const bool reduced = mechRed_->active();
165 scalar pf, cf, pr, cr;
170 forAll(this->reactions_, i)
172 if (!reactionsDisabled_[i])
176 scalar omegai = omega
178 R,
c,
T,
p, pf, cf, lRef, pr, cr, rRef
183 label si =
R.lhs()[
s].index;
186 si = completeToSimplifiedIndex_[si];
189 const scalar sl =
R.lhs()[
s].stoichCoeff;
190 dcdt[si] -= sl*omegai;
194 label si =
R.rhs()[
s].index;
197 si = completeToSimplifiedIndex_[si];
200 const scalar sr =
R.rhs()[
s].stoichCoeff;
201 dcdt[si] += sr*omegai;
208 template<
class ReactionThermo,
class ThermoType>
211 const Reaction<ThermoType>&
R,
223 const scalar kf =
R.kf(
p,
T,
c);
224 const scalar kr =
R.kr(kf,
p,
T,
c);
226 const label Nl =
R.lhs().size();
227 const label Nr =
R.rhs().size();
230 lRef =
R.lhs()[slRef].index;
233 for (label
s=1;
s<Nl;
s++)
235 const label si =
R.lhs()[
s].index;
239 const scalar
exp =
R.lhs()[slRef].exponent;
246 const scalar
exp =
R.lhs()[
s].exponent;
250 cf =
max(
c[lRef], 0.0);
253 const scalar
exp =
R.lhs()[slRef].exponent;
272 rRef =
R.rhs()[srRef].index;
276 for (label
s=1;
s<Nr;
s++)
278 const label si =
R.rhs()[
s].index;
281 const scalar
exp =
R.rhs()[srRef].exponent;
288 const scalar
exp =
R.rhs()[
s].exponent;
292 cr =
max(
c[rRef], 0.0);
295 const scalar
exp =
R.rhs()[srRef].exponent;
313 return pf*cf - pr*cr;
317 template<
class ReactionThermo,
class ThermoType>
325 const bool reduced = mechRed_->active();
327 const scalar
T =
c[this->nSpecie_];
328 const scalar
p =
c[this->nSpecie_ + 1];
335 this->c_ = completeC_;
340 for (label i=0; i<NsDAC_; i++)
342 this->c_[simplifiedToCompleteIndex_[i]] =
max(
c[i], 0.0);
347 for (label i=0; i<this->
nSpecie(); i++)
349 this->c_[i] =
max(
c[i], 0.0);
353 omega(this->c_,
T,
p, dcdt);
358 for (label i=0; i<this->c_.size(); i++)
360 const scalar W = this->specieThermo_[i].W();
361 rho += W*this->c_[i];
365 for (label i=0; i<this->c_.size(); i++)
368 cp += this->c_[i]*this->specieThermo_[i].cp(
p,
T);
376 for (label i=0; i<this->nSpecie_; i++)
381 si = simplifiedToCompleteIndex_[i];
389 const scalar hi = this->specieThermo_[si].ha(
p,
T);
394 dcdt[this->nSpecie_] = -dT;
397 dcdt[this->nSpecie_ + 1] = 0;
401 template<
class ReactionThermo,
class ThermoType>
409 const bool reduced = mechRed_->active();
416 const label
nSpecie = this->nSpecie_;
423 this->c_ = completeC_;
424 for (label i=0; i<NsDAC_; ++i)
426 this->c_[simplifiedToCompleteIndex_[i]] =
max(
c[i], 0.0);
433 this->c_[i] =
max(
c[i], 0.0);
439 forAll(this->reactions_, ri)
441 if (!reactionsDisabled_[ri])
443 const Reaction<ThermoType>&
R = this->reactions_[ri];
445 const scalar kf0 =
R.kf(
p,
T, this->c_);
446 const scalar kr0 =
R.kr(kf0,
p,
T, this->c_);
450 label sj =
R.lhs()[j].index;
453 sj = completeToSimplifiedIndex_[sj];
458 const label si =
R.lhs()[i].index;
459 const scalar el =
R.lhs()[i].exponent;
464 if (this->c_[si] > SMALL)
466 kf *= el*
pow(this->c_[si], el - 1);
475 kf *= el*
pow(this->c_[si], el - 1);
480 kf *=
pow(this->c_[si], el);
486 label si =
R.lhs()[i].index;
489 si = completeToSimplifiedIndex_[si];
491 const scalar sl =
R.lhs()[i].stoichCoeff;
492 dfdc(si, sj) -= sl*kf;
496 label si =
R.rhs()[i].index;
499 si = completeToSimplifiedIndex_[si];
501 const scalar sr =
R.rhs()[i].stoichCoeff;
502 dfdc(si, sj) += sr*kf;
508 label sj =
R.rhs()[j].index;
511 sj = completeToSimplifiedIndex_[sj];
516 const label si =
R.rhs()[i].index;
517 const scalar er =
R.rhs()[i].exponent;
522 if (this->c_[si] > SMALL)
524 kr *= er*
pow(this->c_[si], er - 1);
533 kr *= er*
pow(this->c_[si], er - 1);
538 kr *=
pow(this->c_[si], er);
544 label si =
R.lhs()[i].index;
547 si = completeToSimplifiedIndex_[si];
549 const scalar sl =
R.lhs()[i].stoichCoeff;
550 dfdc(si, sj) += sl*kr;
554 label si =
R.rhs()[i].index;
557 si = completeToSimplifiedIndex_[si];
559 const scalar sr =
R.rhs()[i].stoichCoeff;
560 dfdc(si, sj) -= sr*kr;
567 const scalar
delta = 1
e-3;
569 omega(this->c_,
T +
delta,
p, this->dcdt_);
570 for (label i=0; i<
nSpecie; ++i)
572 dfdc(i,
nSpecie) = this->dcdt_[i];
575 omega(this->c_,
T -
delta,
p, this->dcdt_);
576 for (label i=0; i<
nSpecie; ++i)
586 template<
class ReactionThermo,
class ThermoType>
595 jacobian(t,
c, dfdc);
597 const scalar
T =
c[this->nSpecie_];
598 const scalar
p =
c[this->nSpecie_ + 1];
601 omega(this->c_,
T,
p, dcdt);
605 template<
class ReactionThermo,
class ThermoType>
606 template<
class DeltaTType>
609 const DeltaTType& deltaT
615 const bool reduced = mechRed_->
active();
617 label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
624 scalar reduceMechCpuTime_ = 0;
625 scalar addNewLeafCpuTime_ = 0;
626 scalar growCpuTime_ = 0;
627 scalar solveChemistryCpuTime_ = 0;
628 scalar searchISATCpuTime_ = 0;
630 this->resetTabulationResults();
633 scalar nActiveSpecies = 0;
638 scalar deltaTMin = GREAT;
640 if (!this->chemistry_)
654 IOobject::NO_REGISTER
672 const scalar rhoi =
rho[celli];
673 scalar
pi =
p[celli];
674 scalar Ti =
T[celli];
676 for (label i=0; i<this->nSpecie_; i++)
679 c[i] = rhoi*Yi[celli]/this->specieThermo_[i].W();
685 if (tabulation_->variableTimeStep())
687 phiq[this->
nSpecie() + 2] = deltaT[celli];
692 scalar timeLeft = deltaT[celli];
702 if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
705 for (label i=0; i<this->
nSpecie(); ++i)
707 c[i] = rhoi*Rphiq[i]/this->specieThermo_[i].W();
724 mechRed_->reduceMechanism(
c, Ti,
pi);
725 nActiveSpecies += mechRed_->NsSimp();
728 reduceMechCpuTime_ += timeIncr;
733 while (timeLeft > SMALL)
735 scalar dt = timeLeft;
745 simplifiedC_, Ti,
pi, dt, this->deltaTChem_[celli]
748 for (label i=0; i<NsDAC_; ++i)
750 c[simplifiedToCompleteIndex_[i]] = simplifiedC_[i];
755 this->
solve(c, Ti,
pi, dt, this->deltaTChem_[celli]);
762 solveChemistryCpuTime_ += timeIncr;
768 if (tabulation_->active())
772 Rphiq[i] =
c[i]/rhoi*this->specieThermo_[i].W();
774 if (tabulation_->variableTimeStep())
776 Rphiq[Rphiq.size()-3] = Ti;
777 Rphiq[Rphiq.size()-2] =
pi;
778 Rphiq[Rphiq.size()-1] = deltaT[celli];
782 Rphiq[Rphiq.size()-2] = Ti;
783 Rphiq[Rphiq.size()-1] =
pi;
786 tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]);
790 this->setTabulationResultsAdd(celli);
795 this->setTabulationResultsGrow(celli);
805 this->nSpecie_ = mechRed_->nSpecie();
807 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
809 this->deltaTChem_[celli] =
810 min(this->deltaTChem_[celli], this->deltaTChemMax_);
814 for (label i=0; i<this->nSpecie_; ++i)
816 this->RR_[i][celli] =
817 (
c[i] - c0[i])*this->specieThermo_[i].W()/deltaT[celli];
821 if (mechRed_->log() || tabulation_->log())
824 << this->time().timeOutputValue()
825 <<
" " << solveChemistryCpuTime_ <<
endl;
831 << this->time().timeOutputValue()
832 <<
" " << reduceMechCpuTime_ <<
endl;
835 if (tabulation_->active())
838 tabulation_->update();
841 tabulation_->writePerformance();
843 if (tabulation_->log())
846 << this->time().timeOutputValue()
847 <<
" " << searchISATCpuTime_ <<
endl;
850 << this->time().timeOutputValue()
851 <<
" " << growCpuTime_ <<
endl;
854 << this->time().timeOutputValue()
855 <<
" " << addNewLeafCpuTime_ <<
endl;
859 if (reduced && nAvg && mechRed_->log())
862 nActiveSpeciesFile_()
863 << this->time().timeOutputValue()
864 <<
" " << nActiveSpecies/nAvg <<
endl;
867 if (reduced && Pstream::parRun())
870 Pstream::listCombineReduce(active, orEqOp<bool>());
885 this->
Y()[i].writeOpt(IOobject::AUTO_WRITE);
893 template<
class ReactionThermo,
class ThermoType>
908 template<
class ReactionThermo,
class ThermoType>
914 return this->solve<scalarField>(deltaT);
918 template<
class ReactionThermo,
class ThermoType>
925 tabulationResults_[celli] = 0.0;
929 template<
class ReactionThermo,
class ThermoType>
933 tabulationResults_[celli] = 1.0;
937 template<
class ReactionThermo,
class ThermoType>
944 tabulationResults_[celli] = 2.0;
Extends StandardChemistryModel by adding the TDAC method.
void setTabulationResultsGrow(const label celli)
void setTabulationResultsRetrieve(const label celli)
basicSpecieMixture & composition
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
virtual ~TDACChemistryModel()
Destructor.
const dimensionSet dimless
Dimensionless.
Specialization of basicMultiComponentMixture for a mixture consisting of a number for molecular speci...
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
virtual void omega(const scalarField &c, const scalar T, const scalar p, scalarField &dcdt) const
dc/dt = omega, rate of change in concentration, for each species
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
psiReactionThermo & thermo
const dimensionedScalar e
Elementary charge.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
dimensionedScalar exp(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr scalar pi(M_PI)
void setTabulationResultsAdd(const label celli)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
void jacobian(const scalar t, const scalarField &c, scalarSquareMatrix &dfdc) const
Pure jacobian function for tabulation.
runTime controlDict().readEntry("adjustTimeStep"
The central control dictionary, the contents of which are either taken directly from the FOAM_CONTROL...
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
const volScalarField & cp
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/chef2/andy/OpenFOAM/release/v2306/OpenFOAM-v2306/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=clamp((rho - rholSat)/(rhovSat - rholSat), zero_one{});alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Type & dynamicCast(U &obj)
A dynamic_cast (for references) that generates FatalError on failed casts.
T & ref()
Return reference to the managed object without nullptr checking.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
bool active(const label i) const
#define R(A, B, C, D, E, F, K, M)
A wordList with hashed named lookup, which can be faster in some situations than using the normal lis...
PtrList< volScalarField > & Y
const dimensionedScalar c
Speed of light in a vacuum.
double timeIncrement() const
The time [seconds] since the last call to timeIncrement()
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
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))
Defines the attributes of an object for which implicit objectRegistry management is supported...
Starts timing and returns elapsed time from start. Uses std::chrono::high_resolution_clock for better...
static constexpr const zero Zero
Global zero (0)