36template<
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())
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");
147template<
class ReactionThermo,
class ThermoType>
154template<
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;
208template<
class ReactionThermo,
class ThermoType>
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;
317template<
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;
401template<
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])
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)
586template<
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);
605template<
class ReactionThermo,
class ThermoType>
606template<
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_)
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();
684 phiq[this->
nSpecie() + 1] = pi;
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;
893template<
class ReactionThermo,
class ThermoType>
908template<
class ReactionThermo,
class ThermoType>
914 return this->solve<scalarField>(deltaT);
918template<
class ReactionThermo,
class ThermoType>
925 tabulationResults_[celli] = 0.0;
929template<
class ReactionThermo,
class ThermoType>
933 tabulationResults_[celli] = 1.0;
937template<
class ReactionThermo,
class ThermoType>
944 tabulationResults_[celli] = 2.0;
#define R(A, B, C, D, E, F, K, M)
ReactionThermo & thermo()
Return access to the thermo package.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
const Time & time() const
Return Time associated with the objectRegistry.
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.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
static void listCombineAllGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
After completion all processors have the same data.
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
Extends StandardChemistryModel by adding the TDAC method.
virtual ~TDACChemistryModel()
Destructor.
autoPtr< OFstream > logFile(const word &name) const
Create and return a TDAC log file of the given name.
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
void setTabulationResultsAdd(const label celli)
PtrList< volScalarField > & Y()
void setTabulationResultsRetrieve(const label celli)
void jacobian(const scalar t, const scalarField &c, scalarSquareMatrix &dfdc) const
Pure jacobian function for tabulation.
void setTabulationResultsGrow(const label celli)
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
void size(const label n)
Older name for setAddressableSize.
static bool & parRun() noexcept
Test if this a parallel run.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
T & ref()
Return reference to the managed object without nullptr checking.
const word & name() const
const fvMesh & mesh() const
Return const access to the mesh database.
void setActive(label speciei)
Set speciei active.
bool active(label speciei) const
Return true for active species.
const speciesTable & species() const
Return the table of species.
void setInactive(label speciei)
Set speciei inactive.
Specialization of basicMultiComponentMixture for a mixture consisting of a number for molecular speci...
Starts timing and returns elapsed time from start. Uses std::chrono::high_resolution_clock for better...
double timeIncrement() const
The time [seconds] since the last call to timeIncrement()
A wordList with hashed named lookup, which can be faster in some situations than using the normal lis...
const volScalarField & omega() const
Access functions.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
basicSpecieMixture & composition
PtrList< volScalarField > & Y
runTime controlDict().readEntry("adjustTimeStep"
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))
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar exp(const dimensionedScalar &ds)
const dimensionSet dimless
Dimensionless.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
static constexpr const zero Zero
Global zero (0)
To & dynamicCast(From &r)
const volScalarField & cp
#define forAll(list, i)
Loop across all elements in list.