33 template<
class ReactionThermo,
class ThermoType>
36 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(),
76 dynamicCast<const reactingMixture<ThermoType>&>(this->
thermo())
81 specieComp_[i] = specComp[this->
Y()[i].member()];
92 if (mechRed_->active())
109 this->
Y()[i].writeOpt() = IOobject::NO_WRITE;
122 cpuReduceFile_ = logFile(
"cpu_reduce.out");
123 nActiveSpeciesFile_ = logFile(
"nActiveSpecies.out");
126 if (tabulation_->log())
128 cpuAddFile_ = logFile(
"cpu_add.out");
129 cpuGrowFile_ = logFile(
"cpu_grow.out");
130 cpuRetrieveFile_ = logFile(
"cpu_retrieve.out");
133 if (mechRed_->log() || tabulation_->log())
135 cpuSolveFile_ = logFile(
"cpu_solve.out");
142 template<
class ReactionThermo,
class ThermoType>
149 template<
class ReactionThermo,
class ThermoType>
158 const bool reduced = mechRed_->active();
160 scalar pf, cf, pr, cr;
165 forAll(this->reactions_, i)
167 if (!reactionsDisabled_[i])
171 scalar omegai = R.
omega 173 p, T, c, pf, cf, lRef, pr, cr, rRef
181 si = completeToSimplifiedIndex_[si];
184 const scalar sl = R.
lhs()[
s].stoichCoeff;
185 dcdt[si] -= sl*omegai;
192 si = completeToSimplifiedIndex_[si];
195 const scalar sr = R.
rhs()[
s].stoichCoeff;
196 dcdt[si] += sr*omegai;
203 template<
class ReactionThermo,
class ThermoType>
218 const scalar kf = R.
kf(p, T, c);
219 const scalar kr = R.
kr(kf, p, T, c);
225 lRef = R.
lhs()[slRef].index;
234 const scalar
exp = R.
lhs()[slRef].exponent;
235 pf *=
pow(
max(c[lRef], 0), exp);
241 const scalar
exp = R.
lhs()[
s].exponent;
242 pf *=
pow(
max(c[si], 0), exp);
245 cf =
max(c[lRef], 0);
248 const scalar
exp = R.
lhs()[slRef].exponent;
253 pf *=
pow(cf, exp - 1);
262 pf *=
pow(cf, exp - 1);
267 rRef = R.
rhs()[srRef].index;
276 const scalar
exp = R.
rhs()[srRef].exponent;
277 pr *=
pow(
max(c[rRef], 0), exp);
283 const scalar
exp = R.
rhs()[
s].exponent;
284 pr *=
pow(
max(c[si], 0), exp);
287 cr =
max(c[rRef], 0);
290 const scalar
exp = R.
rhs()[srRef].exponent;
295 pr *=
pow(cr, exp - 1);
304 pr *=
pow(cr, exp - 1);
308 return pf*cf - pr*cr;
312 template<
class ReactionThermo,
class ThermoType>
320 const bool reduced = mechRed_->active();
322 const scalar T = c[this->nSpecie_];
323 const scalar p = c[this->nSpecie_ + 1];
330 this->c_ = completeC_;
335 for (
label i=0; i<NsDAC_; i++)
337 this->c_[simplifiedToCompleteIndex_[i]] =
max(c[i], 0);
344 this->c_[i] =
max(c[i], 0);
348 omega(this->c_, T, p, dcdt);
353 for (
label i=0; i<this->c_.size(); i++)
355 const scalar
W = this->specieThermo_[i].W();
356 rho += W*this->c_[i];
360 for (
label i=0; i<this->c_.size(); i++)
363 cp += this->c_[i]*this->specieThermo_[i].cp(p, T);
371 for (
label i=0; i<this->nSpecie_; i++)
376 si = simplifiedToCompleteIndex_[i];
384 const scalar hi = this->specieThermo_[si].ha(p, T);
389 dcdt[this->nSpecie_] = -dT;
392 dcdt[this->nSpecie_ + 1] = 0;
396 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->specieThermo_[i].ha(p, T);
438 cpi[i] = this->specieThermo_[i].cp(p, T);
443 forAll(this->reactions_, ri)
445 if (!reactionsDisabled_[ri])
460 completeToSimplifiedIndex_
472 completeToSimplifiedIndex_,
481 scalar dcpdTMean = 0;
484 cpMean += this->c_[i]*cpi[i];
486 dcpdTMean += this->c_[i]*this->specieThermo_[i].dcpdT(p, T);
494 const label si = completeToSimplifiedIndex_[i];
497 dTdt += hi[i]*dcdt[si];
502 dTdt += hi[i]*dcdt[i];
506 dcdt[this->nSpecie_] = dTdt;
508 for (
label i = 0; i < this->nSpecie_; i++)
510 J(this->nSpecie_, i) = 0;
511 for (
label j = 0; j < this->nSpecie_; j++)
513 const label sj = reduced ? simplifiedToCompleteIndex_[j] : j;
514 J(this->nSpecie_, i) += hi[sj]*J(j, i);
516 const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
517 J(this->nSpecie_, i) += cpi[si]*dTdt;
518 J(this->nSpecie_, i) /= -cpMean;
522 J(this->nSpecie_, this->nSpecie_) = 0;
523 for (
label i = 0; i < this->nSpecie_; i++)
525 const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
526 J(this->nSpecie_, this->nSpecie_) +=
528 + hi[si]*J(i, this->nSpecie_);
530 J(this->nSpecie_, this->nSpecie_) += dTdt*dcpdTMean;
531 J(this->nSpecie_, this->nSpecie_) /= -cpMean;
532 J(this->nSpecie_, this->nSpecie_) += dTdt/
T;
536 template<
class ReactionThermo,
class ThermoType>
537 template<
class DeltaTType>
540 const DeltaTType& deltaT
546 const bool reduced = mechRed_->
active();
548 label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
555 scalar reduceMechCpuTime_ = 0;
556 scalar addNewLeafCpuTime_ = 0;
557 scalar growCpuTime_ = 0;
558 scalar solveChemistryCpuTime_ = 0;
559 scalar searchISATCpuTime_ = 0;
561 this->resetTabulationResults();
564 scalar nActiveSpecies = 0;
569 scalar deltaTMin = great;
571 if (!this->chemistry_)
603 const scalar rhoi = rho[celli];
604 scalar
pi = p[celli];
605 scalar Ti = T[celli];
607 for (
label i=0; i<this->nSpecie_; i++)
609 c[i] = rhoi*this->Y_[i][celli]/this->specieThermo_[i].W();
611 phiq[i] = this->
Y()[i][celli];
615 if (tabulation_->variableTimeStep())
617 phiq[this->
nSpecie() + 2] = deltaT[celli];
622 scalar timeLeft = deltaT[celli];
632 if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
637 c[i] = rhoi*Rphiq[i]/this->specieThermo_[i].W();
654 mechRed_->reduceMechanism(c, Ti, pi);
655 nActiveSpecies += mechRed_->NsSimp();
658 reduceMechCpuTime_ += timeIncr;
663 while (timeLeft > small)
665 scalar dt = timeLeft;
675 simplifiedC_, Ti, pi, dt, this->deltaTChem_[celli]
678 for (
label i=0; i<NsDAC_; i++)
680 c[simplifiedToCompleteIndex_[i]] = simplifiedC_[i];
685 this->
solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
692 solveChemistryCpuTime_ += timeIncr;
698 if (tabulation_->active())
702 Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W();
704 if (tabulation_->variableTimeStep())
706 Rphiq[Rphiq.
size()-3] = Ti;
707 Rphiq[Rphiq.
size()-2] =
pi;
708 Rphiq[Rphiq.
size()-1] = deltaT[celli];
712 Rphiq[Rphiq.
size()-2] = Ti;
713 Rphiq[Rphiq.
size()-1] =
pi;
716 tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]);
719 this->setTabulationResultsAdd(celli);
724 this->setTabulationResultsGrow(celli);
734 this->nSpecie_ = mechRed_->nSpecie();
736 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
738 this->deltaTChem_[celli] =
739 min(this->deltaTChem_[celli], this->deltaTChemMax_);
743 for (
label i=0; i<this->nSpecie_; i++)
745 this->RR_[i][celli] =
746 (c[i] - c0[i])*this->specieThermo_[i].
W()/deltaT[celli];
750 if (mechRed_->log() || tabulation_->log())
753 << this->time().timeOutputValue()
754 <<
" " << solveChemistryCpuTime_ <<
endl;
760 << this->time().timeOutputValue()
761 <<
" " << reduceMechCpuTime_ <<
endl;
764 if (tabulation_->active())
767 tabulation_->update();
770 tabulation_->writePerformance();
772 if (tabulation_->log())
775 << this->time().timeOutputValue()
776 <<
" " << searchISATCpuTime_ <<
endl;
779 << this->time().timeOutputValue()
780 <<
" " << growCpuTime_ <<
endl;
783 << this->time().timeOutputValue()
784 <<
" " << addNewLeafCpuTime_ <<
endl;
788 if (reduced && nAvg && mechRed_->log())
791 nActiveSpeciesFile_()
792 << this->time().timeOutputValue()
793 <<
" " << nActiveSpecies/nAvg <<
endl;
796 if (Pstream::parRun())
800 Pstream::listCombineScatter(active);
813 if (composition.
active(i))
815 this->
Y()[i].writeOpt() = IOobject::AUTO_WRITE;
823 template<
class ReactionThermo,
class ThermoType>
838 template<
class ReactionThermo,
class ThermoType>
844 return this->solve<scalarField>(deltaT);
848 template<
class ReactionThermo,
class ThermoType>
855 tabulationResults_[celli] = 0;
859 template<
class ReactionThermo,
class ThermoType>
863 tabulationResults_[celli] = 1;
867 template<
class ReactionThermo,
class ThermoType>
874 tabulationResults_[celli] = 2;
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.
bool active(label speciei) const
Return true for active species.
void setTabulationResultsGrow(const label celli)
void setTabulationResultsRetrieve(const label celli)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
void dwdT(const scalar p, const scalar T, const scalarField &c, 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.
basicSpecieMixture & composition
TDACChemistryModel(ReactionThermo &thermo)
Construct from thermo.
void size(const label)
Override size to be inconsistent with allocated storage.
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual ~TDACChemistryModel()
Destructor.
const speciesTable & species() const
Return the table of species.
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
virtual scalar kr(const scalar kfwd, const scalar p, const scalar T, const scalarField &c) const =0
Reverse rate constant from the given forward rate constant.
rhoReactionThermo & thermo
void omega(const scalar p, const scalar T, const scalarField &c, scalarField &dcdt) const
Net reaction rate for individual species.
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...
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
void setInactive(label speciei)
Set speciei inactive.
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)
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &J) const
Calculate the Jacobian of the system.
void dwdc(const scalar p, const scalar T, const scalarField &c, 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.
An STL-conforming hash table.
void setTabulationResultsAdd(const label celli)
void setActive(label speciei)
Set speciei active.
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(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-7/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
A wordList with hashed indices for faster lookup by name.
const scalarList W(::W(thermo))
PtrList< volScalarField > & Y
const dimensionedScalar c
Speed of light in a vacuum.
double timeIncrement() const
Return time (in seconds) since last call to timeIncrement()
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
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 =0
Forward rate constant.
Starts timing (using rtc) and returns elapsed time from start. Better resolution (2uSec instead of ~2...