33 template<
class CompType,
class ThermoType>
43 mesh.time().controlDict().lookupOrDefault(
"adjustTimeStep",
false)
44 || fv::localEulerDdt::enabled(mesh)
47 NsDAC_(this->nSpecie_),
48 completeC_(this->nSpecie_, 0),
49 reactionsDisabled_(this->reactions_.size(),
false),
50 specieComp_(this->nSpecie_),
51 completeToSimplifiedIndex_(this->nSpecie_, -1),
52 simplifiedToCompleteIndex_(this->nSpecie_),
73 dynamicCast<const reactingMixture<ThermoType>&>(this->
thermo())
78 specieComp_[i] = specComp[this->
Y()[i].name()];
89 if (mechRed_->active())
96 mesh.time().timeName(),
106 this->
Y()[i].writeOpt() = IOobject::NO_WRITE;
119 cpuReduceFile_ = logFile(
"cpu_reduce.out");
120 nActiveSpeciesFile_ = logFile(
"nActiveSpecies.out");
123 if (tabulation_->log())
125 cpuAddFile_ = logFile(
"cpu_add.out");
126 cpuGrowFile_ = logFile(
"cpu_grow.out");
127 cpuRetrieveFile_ = logFile(
"cpu_retrieve.out");
130 if (mechRed_->log() || tabulation_->log())
132 cpuSolveFile_ = logFile(
"cpu_solve.out");
139 template<
class CompType,
class ThermoType>
146 template<
class CompType,
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 = omega
170 R, c, T, p, 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 CompType,
class ThermoType>
215 const scalar kf = R.
kf(p, T, c);
216 const scalar kr = R.
kr(kf, p, T, c);
222 lRef = R.
lhs()[slRef].index;
231 const scalar
exp = R.
lhs()[slRef].exponent;
232 pf *=
pow(
max(0.0, c[lRef]), exp);
238 const scalar
exp = R.
lhs()[
s].exponent;
239 pf *=
pow(
max(0.0, c[si]), exp);
242 cf =
max(0.0, c[lRef]);
245 const scalar
exp = R.
lhs()[slRef].exponent;
250 pf *=
pow(cf, exp - 1);
259 pf *=
pow(cf, exp - 1);
264 rRef = R.
rhs()[srRef].index;
273 const scalar
exp = R.
rhs()[srRef].exponent;
274 pr *=
pow(
max(0.0, c[rRef]), exp);
280 const scalar
exp = R.
rhs()[
s].exponent;
281 pr *=
pow(
max(0.0, c[si]), exp);
284 cr =
max(0.0, c[rRef]);
287 const scalar
exp = R.
rhs()[srRef].exponent;
292 pr *=
pow(cr, exp - 1);
301 pr *=
pow(cr, exp - 1);
305 return pf*cf - pr*cr;
309 template<
class CompType,
class ThermoType>
317 const bool reduced = mechRed_->active();
319 const scalar T = c[this->nSpecie_];
320 const scalar p = c[this->nSpecie_ + 1];
327 this->c_ = completeC_;
332 for (
label i=0; i<NsDAC_; i++)
334 this->c_[simplifiedToCompleteIndex_[i]] =
max(0.0, c[i]);
341 this->c_[i] =
max(0.0, c[i]);
345 omega(this->c_, T, p, dcdt);
350 for (
label i=0; i<this->c_.size(); i++)
352 const scalar W = this->specieThermo_[i].W();
353 rho += W*this->c_[i];
357 for (
label i=0; i<this->c_.size(); i++)
360 cp += this->c_[i]*this->specieThermo_[i].cp(p, T);
368 for (
label i=0; i<this->nSpecie_; i++)
373 si = simplifiedToCompleteIndex_[i];
381 const scalar hi = this->specieThermo_[si].ha(p, T);
386 dcdt[this->nSpecie_] = -dT;
389 dcdt[this->nSpecie_ + 1] = 0;
393 template<
class CompType,
class ThermoType>
401 const bool reduced = mechRed_->active();
408 const scalar T = c[this->nSpecie_];
409 const scalar p = c[this->nSpecie_ + 1];
413 this->c_ = completeC_;
414 for (
label i=0; i<NsDAC_; i++)
416 this->c_[simplifiedToCompleteIndex_[i]] =
max(0.0, c[i]);
423 this->c_[i] =
max(c[i], 0.0);
429 forAll(this->reactions_, ri)
431 if (!reactionsDisabled_[ri])
435 const scalar kf0 = R.
kf(p, T, this->c_);
436 const scalar kr0 = R.
kr(kf0, p, T, this->c_);
443 sj = completeToSimplifiedIndex_[sj];
449 const scalar el = R.
lhs()[i].exponent;
454 if (this->c_[si] > SMALL)
456 kf *= el*
pow(this->c_[si] + VSMALL, el - 1);
465 kf *= el*
pow(this->c_[si], el - 1);
470 kf *=
pow(this->c_[si], el);
479 si = completeToSimplifiedIndex_[si];
481 const scalar sl = R.
lhs()[i].stoichCoeff;
482 dfdc(si, sj) -= sl*kf;
489 si = completeToSimplifiedIndex_[si];
491 const scalar sr = R.
rhs()[i].stoichCoeff;
492 dfdc(si, sj) += sr*kf;
501 sj = completeToSimplifiedIndex_[sj];
507 const scalar er = R.
rhs()[i].exponent;
512 if (this->c_[si] > SMALL)
514 kr *= er*
pow(this->c_[si] + VSMALL, er - 1);
523 kr *= er*
pow(this->c_[si], er - 1);
528 kr *=
pow(this->c_[si], er);
537 si = completeToSimplifiedIndex_[si];
539 const scalar sl = R.
lhs()[i].stoichCoeff;
540 dfdc(si, sj) += sl*kr;
547 si = completeToSimplifiedIndex_[si];
549 const scalar sr = R.
rhs()[i].stoichCoeff;
550 dfdc(si, sj) -= sr*kr;
557 const scalar
delta = 1
e-3;
559 omega(this->c_, T + delta, p, this->dcdt_);
560 for (
label i=0; i<this->nSpecie_; i++)
562 dfdc(i, this->nSpecie_) = this->dcdt_[i];
565 omega(this->c_, T - delta, p, this->dcdt_);
566 for (
label i=0; i<this->nSpecie_; i++)
568 dfdc(i, this->nSpecie_) =
569 0.5*(dfdc(i, this->nSpecie_) - this->dcdt_[i])/delta;
572 dfdc(this->nSpecie_, this->nSpecie_) = 0;
573 dfdc(this->nSpecie_ + 1, this->nSpecie_) = 0;
577 template<
class CompType,
class ThermoType>
586 jacobian(t, c, dfdc);
588 const scalar T = c[this->nSpecie_];
589 const scalar p = c[this->nSpecie_ + 1];
592 omega(this->c_, T, p, dcdt);
596 template<
class CompType,
class ThermoType>
597 template<
class DeltaTType>
600 const DeltaTType& deltaT
606 const bool reduced = mechRed_->
active();
608 label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
615 scalar reduceMechCpuTime_ = 0;
616 scalar addNewLeafCpuTime_ = 0;
617 scalar growCpuTime_ = 0;
618 scalar solveChemistryCpuTime_ = 0;
619 scalar searchISATCpuTime_ = 0;
621 this->resetTabulationResults();
624 scalar nActiveSpecies = 0;
629 scalar deltaTMin = GREAT;
631 if (!this->chemistry_)
663 const scalar rhoi = rho[celli];
664 scalar
pi = p[celli];
665 scalar Ti = T[celli];
667 for (
label i=0; i<this->nSpecie_; i++)
669 c[i] = rhoi*this->Y_[i][celli]/this->specieThermo_[i].W();
671 phiq[i] = this->
Y()[i][celli];
675 if (tabulation_->variableTimeStep())
677 phiq[this->
nSpecie() + 2] = deltaT[celli];
682 scalar timeLeft = deltaT[celli];
692 if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
697 c[i] = rhoi*Rphiq[i]/this->specieThermo_[i].W();
714 mechRed_->reduceMechanism(c, Ti, pi);
715 nActiveSpecies += mechRed_->NsSimp();
718 reduceMechCpuTime_ += timeIncr;
723 while (timeLeft > SMALL)
725 scalar dt = timeLeft;
735 simplifiedC_, Ti, pi, dt, this->deltaTChem_[celli]
738 for (
label i=0; i<NsDAC_; i++)
740 c[simplifiedToCompleteIndex_[i]] = simplifiedC_[i];
745 this->
solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
752 solveChemistryCpuTime_ += timeIncr;
758 if (tabulation_->active())
762 Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W();
764 if (tabulation_->variableTimeStep())
766 Rphiq[Rphiq.
size()-3] = Ti;
767 Rphiq[Rphiq.
size()-2] =
pi;
768 Rphiq[Rphiq.
size()-1] = deltaT[celli];
772 Rphiq[Rphiq.
size()-2] = Ti;
773 Rphiq[Rphiq.
size()-1] =
pi;
776 tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]);
779 this->setTabulationResultsAdd(celli);
784 this->setTabulationResultsGrow(celli);
794 this->nSpecie_ = mechRed_->nSpecie();
796 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
800 for (
label i=0; i<this->nSpecie_; i++)
802 this->RR_[i][celli] =
803 (c[i] - c0[i])*this->specieThermo_[i].W()/deltaT[celli];
807 if (mechRed_->log() || tabulation_->log())
810 << this->time().timeOutputValue()
811 <<
" " << solveChemistryCpuTime_ <<
endl;
817 << this->time().timeOutputValue()
818 <<
" " << reduceMechCpuTime_ <<
endl;
821 if (tabulation_->active())
824 tabulation_->update();
827 tabulation_->writePerformance();
829 if (tabulation_->log())
832 << this->time().timeOutputValue()
833 <<
" " << searchISATCpuTime_ <<
endl;
836 << this->time().timeOutputValue()
837 <<
" " << growCpuTime_ <<
endl;
840 << this->time().timeOutputValue()
841 <<
" " << addNewLeafCpuTime_ <<
endl;
845 if (reduced && nAvg && mechRed_->log())
848 nActiveSpeciesFile_()
849 << this->time().timeOutputValue()
850 <<
" " << nActiveSpecies/nAvg <<
endl;
853 if (Pstream::parRun())
857 Pstream::listCombineScatter(active);
870 if (composition.
active(i))
872 this->
Y()[i].writeOpt() = IOobject::AUTO_WRITE;
880 template<
class CompType,
class ThermoType>
895 template<
class CompType,
class ThermoType>
901 return this->solve<scalarField>(deltaT);
905 template<
class CompType,
class ThermoType>
911 tabulationResults_[celli] = 0.0;
915 template<
class CompType,
class ThermoType>
921 tabulationResults_[celli] = 1.0;
925 template<
class CompType,
class ThermoType>
932 tabulationResults_[celli] = 2.0;
virtual scalar kr(const scalar kfwd, const scalar p, const scalar T, const scalarField &c) const
Reverse rate constant from the given forward rate constant.
virtual scalar kf(const scalar p, const scalar T, const scalarField &c) const
Forward rate constant.
Extends chemistryModel 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 active(label speciei) const
Return true for active species.
void setTabulationResultsGrow(const label celli)
const List< specieCoeffs > & lhs() const
void setTabulationResultsRetrieve(const label celli)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
basicMultiComponentMixture & composition
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.
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.
psiReactionThermo & thermo
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))
const dimensionedScalar e
Elementary charge.
dimensionedScalar exp(const dimensionedScalar &ds)
A class for handling words, derived from string.
An STL-conforming hash table.
void setTabulationResultsAdd(const label celli)
void jacobian(const scalar t, const scalarField &c, scalarSquareMatrix &dfdc) const
Pure jacobian function for tabulation.
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(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/ubuntu/OpenFOAM-5.0/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
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
A wordList with hashed indices for faster lookup by name.
PtrList< volScalarField > & Y
Mesh data needed to do the Finite Volume discretisation.
const dimensionedScalar c
Speed of light in a vacuum.
double timeIncrement() const
Return time (in seconds) since last call to timeIncrement()
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
const List< specieCoeffs > & rhs() const
Starts timing (using rtc) and returns elapsed time from start. Better resolution (2uSec instead of ~2...