33 template<
class ReactionThermo,
class ThermoType>
36 const ReactionThermo& thermo
42 this->
mesh().time().controlDict().lookupOrDefault
47 || fv::localEulerDdt::enabled(this->
mesh())
50 NsDAC_(this->nSpecie_),
51 completeC_(this->nSpecie_, 0),
52 reactionsDisabled_(this->reactions_.size(),
false),
53 specieComp_(this->nSpecie_),
54 completeToSimplifiedIndex_(this->nSpecie_, -1),
55 simplifiedToCompleteIndex_(this->nSpecie_),
60 thermo.phasePropertyName(
"TabulationResults"),
61 this->time().timeName(),
73 dynamicCast<const multiComponentMixture<ThermoType>&>(this->
thermo())
78 specieComp_[i] = specComp[this->
Y()[i].member()];
89 if (mechRed_->active())
118 cpuReduceFile_ = logFile(
"cpu_reduce.out");
119 nActiveSpeciesFile_ = logFile(
"nActiveSpecies.out");
122 if (tabulation_->log())
124 cpuAddFile_ = logFile(
"cpu_add.out");
125 cpuGrowFile_ = logFile(
"cpu_grow.out");
126 cpuRetrieveFile_ = logFile(
"cpu_retrieve.out");
129 if (mechRed_->log() || tabulation_->log())
131 cpuSolveFile_ = logFile(
"cpu_solve.out");
138 template<
class ReactionThermo,
class ThermoType>
145 template<
class ReactionThermo,
class ThermoType>
155 const bool reduced = mechRed_->active();
157 scalar pf, cf, pr, cr;
162 forAll(this->reactions_, i)
164 if (!reactionsDisabled_[i])
168 scalar omegai = R.
omega 170 p, T, c, li, pf, cf, lRef, pr, cr, rRef
178 si = completeToSimplifiedIndex_[si];
181 const scalar sl = R.
lhs()[
s].stoichCoeff;
182 dcdt[si] -= sl*omegai;
189 si = completeToSimplifiedIndex_[si];
192 const scalar sr = R.
rhs()[
s].stoichCoeff;
193 dcdt[si] += sr*omegai;
200 template<
class ReactionThermo,
class ThermoType>
216 const scalar kf = R.
kf(p, T, c, li);
217 const scalar kr = R.
kr(kf, p, T, c, li);
223 lRef = R.
lhs()[slRef].index;
232 const scalar
exp = R.
lhs()[slRef].exponent;
233 pf *=
pow(
max(c[lRef], 0), exp);
239 const scalar
exp = R.
lhs()[
s].exponent;
240 pf *=
pow(
max(c[si], 0), exp);
243 cf =
max(c[lRef], 0);
246 const scalar
exp = R.
lhs()[slRef].exponent;
251 pf *=
pow(cf, exp - 1);
260 pf *=
pow(cf, exp - 1);
265 rRef = R.
rhs()[srRef].index;
274 const scalar
exp = R.
rhs()[srRef].exponent;
275 pr *=
pow(
max(c[rRef], 0), exp);
281 const scalar
exp = R.
rhs()[
s].exponent;
282 pr *=
pow(
max(c[si], 0), exp);
285 cr =
max(c[rRef], 0);
288 const scalar
exp = R.
rhs()[srRef].exponent;
293 pr *=
pow(cr, exp - 1);
302 pr *=
pow(cr, exp - 1);
306 return pf*cf - pr*cr;
310 template<
class ReactionThermo,
class ThermoType>
319 const bool reduced = mechRed_->active();
321 const scalar T = c[this->nSpecie_];
322 const scalar p = c[this->nSpecie_ + 1];
329 this->c_ = completeC_;
334 for (
label i=0; i<NsDAC_; i++)
336 this->c_[simplifiedToCompleteIndex_[i]] =
max(c[i], 0);
343 this->c_[i] =
max(c[i], 0);
347 omega(p, T, this->c_, li, dcdt);
352 for (
label i=0; i<this->c_.size(); i++)
354 const scalar
W = this->specieThermos_[i].W();
355 rho += W*this->c_[i];
359 for (
label i=0; i<this->c_.size(); i++)
362 cp += this->c_[i]*this->specieThermos_[i].cp(p, T);
370 for (
label i=0; i<this->nSpecie_; i++)
375 si = simplifiedToCompleteIndex_[i];
383 const scalar hi = this->specieThermos_[si].ha(p, T);
388 dcdt[this->nSpecie_] = -dT;
391 dcdt[this->nSpecie_ + 1] = 0;
395 template<
class ReactionThermo,
class ThermoType>
405 const bool reduced = mechRed_->active();
412 const scalar T = c[this->nSpecie_];
413 const scalar p = c[this->nSpecie_ + 1];
417 this->c_ = completeC_;
418 for (
label i=0; i<NsDAC_; i++)
420 this->c_[simplifiedToCompleteIndex_[i]] =
max(c[i], 0);
427 this->c_[i] =
max(c[i], 0);
437 hi[i] = this->specieThermos_[i].ha(p, T);
438 cpi[i] = this->specieThermos_[i].cp(p, T);
443 forAll(this->reactions_, ri)
445 if (!reactionsDisabled_[ri])
461 completeToSimplifiedIndex_
474 completeToSimplifiedIndex_,
483 scalar dcpdTMean = 0;
486 cpMean += this->c_[i]*cpi[i];
488 dcpdTMean += this->c_[i]*this->specieThermos_[i].dcpdT(p, T);
496 const label si = completeToSimplifiedIndex_[i];
499 dTdt += hi[i]*dcdt[si];
504 dTdt += hi[i]*dcdt[i];
508 dcdt[this->nSpecie_] = dTdt;
510 for (
label i = 0; i < this->nSpecie_; i++)
512 J(this->nSpecie_, i) = 0;
513 for (
label j = 0; j < this->nSpecie_; j++)
515 const label sj = reduced ? simplifiedToCompleteIndex_[j] : j;
516 J(this->nSpecie_, i) += hi[sj]*J(j, i);
518 const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
519 J(this->nSpecie_, i) += cpi[si]*dTdt;
520 J(this->nSpecie_, i) /= -cpMean;
524 J(this->nSpecie_, this->nSpecie_) = 0;
525 for (
label i = 0; i < this->nSpecie_; i++)
527 const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
528 J(this->nSpecie_, this->nSpecie_) +=
530 + hi[si]*J(i, this->nSpecie_);
532 J(this->nSpecie_, this->nSpecie_) += dTdt*dcpdTMean;
533 J(this->nSpecie_, this->nSpecie_) /= -cpMean;
534 J(this->nSpecie_, this->nSpecie_) += dTdt/
T;
538 template<
class ReactionThermo,
class ThermoType>
539 template<
class DeltaTType>
542 const DeltaTType& deltaT
548 const bool reduced = mechRed_->
active();
550 label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
557 scalar reduceMechCpuTime_ = 0;
558 scalar addNewLeafCpuTime_ = 0;
559 scalar growCpuTime_ = 0;
560 scalar solveChemistryCpuTime_ = 0;
561 scalar searchISATCpuTime_ = 0;
563 this->resetTabulationResults();
566 scalar nActiveSpecies = 0;
571 scalar deltaTMin = great;
573 if (!this->chemistry_)
605 const scalar rhoi = rho[celli];
606 scalar
pi = p[celli];
607 scalar Ti = T[celli];
609 for (
label i=0; i<this->nSpecie_; i++)
611 c[i] = rhoi*this->Y_[i][celli]/this->specieThermos_[i].W();
613 phiq[i] = this->
Y()[i][celli];
617 if (tabulation_->variableTimeStep())
619 phiq[this->
nSpecie() + 2] = deltaT[celli];
624 scalar timeLeft = deltaT[celli];
634 if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
639 c[i] = rhoi*Rphiq[i]/this->specieThermos_[i].W();
656 mechRed_->reduceMechanism(pi, Ti, c, celli);
657 nActiveSpecies += mechRed_->NsSimp();
663 while (timeLeft > small)
665 scalar dt = timeLeft;
680 this->deltaTChem_[celli]
683 for (
label i=0; i<NsDAC_; i++)
685 c[simplifiedToCompleteIndex_[i]] = simplifiedC_[i];
690 this->
solve(pi, Ti, c, celli, dt, this->deltaTChem_[celli]);
701 if (tabulation_->active())
705 Rphiq[i] = c[i]/rhoi*this->specieThermos_[i].W();
707 if (tabulation_->variableTimeStep())
709 Rphiq[Rphiq.
size()-3] = Ti;
710 Rphiq[Rphiq.
size()-2] =
pi;
711 Rphiq[Rphiq.
size()-1] = deltaT[celli];
715 Rphiq[Rphiq.
size()-2] = Ti;
716 Rphiq[Rphiq.
size()-1] =
pi;
719 tabulation_->add(phiq, Rphiq, celli, rhoi, deltaT[celli]);
722 this->setTabulationResultsAdd(celli);
727 this->setTabulationResultsGrow(celli);
737 this->nSpecie_ = mechRed_->nSpecie();
739 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
741 this->deltaTChem_[celli] =
742 min(this->deltaTChem_[celli], this->deltaTChemMax_);
746 for (
label i=0; i<this->nSpecie_; i++)
748 this->RR_[i][celli] =
749 (c[i] - c0[i])*this->specieThermos_[i].
W()/deltaT[celli];
753 if (mechRed_->log() || tabulation_->log())
756 << this->time().timeOutputValue()
757 <<
" " << solveChemistryCpuTime_ <<
endl;
763 << this->time().timeOutputValue()
764 <<
" " << reduceMechCpuTime_ <<
endl;
767 if (tabulation_->active())
770 tabulation_->update();
773 tabulation_->writePerformance();
775 if (tabulation_->log())
778 << this->time().timeOutputValue()
779 <<
" " << searchISATCpuTime_ <<
endl;
782 << this->time().timeOutputValue()
783 <<
" " << growCpuTime_ <<
endl;
786 << this->time().timeOutputValue()
787 <<
" " << addNewLeafCpuTime_ <<
endl;
791 if (reduced && nAvg && mechRed_->log())
794 nActiveSpeciesFile_()
795 << this->time().timeOutputValue()
796 <<
" " << nActiveSpecies/nAvg <<
endl;
799 if (Pstream::parRun())
803 Pstream::listCombineScatter(active);
818 template<
class ReactionThermo,
class ThermoType>
833 template<
class ReactionThermo,
class ThermoType>
839 return this->solve<scalarField>(deltaT);
843 template<
class ReactionThermo,
class ThermoType>
850 tabulationResults_[celli] = 0;
854 template<
class ReactionThermo,
class ThermoType>
858 tabulationResults_[celli] = 1;
862 template<
class ReactionThermo,
class ThermoType>
869 tabulationResults_[celli] = 2;
void setInactive(label speciei) const
Set speciei inactive.
Extends StandardChemistryModel 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.
bool typeHeaderOk(const bool checkType=true)
Read header (uses typeFilePath to find file) and check header.
void setTabulationResultsGrow(const label celli)
void setTabulationResultsRetrieve(const label celli)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
bool active(label speciei) const
Return true for active species.
basicSpecieMixture & composition
void dwdc(const scalar p, const scalar T, const scalarField &c, const label li, scalarSquareMatrix &J, scalarField &dcdt, scalar &omegaI, scalar &kfwd, scalar &kbwd, const bool reduced, const List< label > &c2s) const
Derivative of the net reaction rate for each species involved.
void size(const label)
Override size to be inconsistent with allocated storage.
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool cp(const fileName &src, const fileName &dst, const bool followLink=true)
Copy, recursively if necessary, the source to the destination.
virtual ~TDACChemistryModel()
Destructor.
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
rhoReactionThermo & thermo
const dimensionedScalar & c
Speed of light in a vacuum.
Specialization of basicMixture for a mixture consisting of a number for molecular species...
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
void omega(const scalar p, const scalar T, const scalarField &c, const label li, scalarField &dcdt) const
Net reaction rate for individual species.
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))
dimensionedScalar exp(const dimensionedScalar &ds)
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
virtual scalar kr(const scalar kfwd, const scalar p, const scalar T, const scalarField &c, const label li) const =0
Reverse rate constant from the given forward rate constant.
virtual void omega(const scalar p, const scalar T, const scalarField &c, const label li, scalarField &dcdt) const
dc/dt = omega, rate of change in concentration, for each species
TDACChemistryModel(const ReactionThermo &thermo)
Construct from thermo.
An STL-conforming hash table.
void setTabulationResultsAdd(const label celli)
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-8/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
bool active(const label i) const
const scalarList W(::W(thermo))
PtrList< volScalarField > & Y
virtual void derivatives(const scalar t, const scalarField &c, const label li, scalarField &dcdt) const
Calculate the derivatives in dydx.
double timeIncrement() const
Return time (in seconds) since last call to timeIncrement()
void setActive(label speciei) const
Set speciei active.
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
void dwdT(const scalar p, const scalar T, const scalarField &c, const label li, const scalar omegaI, const scalar kfwd, const scalar kbwd, scalarSquareMatrix &J, const bool reduced, const List< label > &c2s, const label indexT) const
Derivative of the net reaction rate for each species involved.
virtual void jacobian(const scalar t, const scalarField &c, const label li, scalarField &dcdt, scalarSquareMatrix &J) const
Calculate the Jacobian of the system.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
virtual scalar kf(const scalar p, const scalar T, const scalarField &c, const label li) const =0
Forward rate constant.
Starts timing (using rtc) and returns elapsed time from start. Better resolution (2uSec instead of ~2...