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"),
64 this->time().timeName(),
79 dynamicCast<const reactingMixture<ThermoType>&>(this->
thermo())
84 specieComp_[i] = specComp[this->
Y()[i].member()];
95 if (mechRed_->active())
112 this->
Y()[i].writeOpt() = IOobject::NO_WRITE;
125 cpuReduceFile_ = logFile(
"cpu_reduce.out");
126 nActiveSpeciesFile_ = logFile(
"nActiveSpecies.out");
129 if (tabulation_->log())
131 cpuAddFile_ = logFile(
"cpu_add.out");
132 cpuGrowFile_ = logFile(
"cpu_grow.out");
133 cpuRetrieveFile_ = logFile(
"cpu_retrieve.out");
136 if (mechRed_->log() || tabulation_->log())
138 cpuSolveFile_ = logFile(
"cpu_solve.out");
145 template<
class ReactionThermo,
class ThermoType>
152 template<
class ReactionThermo,
class ThermoType>
161 const bool reduced = mechRed_->active();
163 scalar pf, cf, pr, cr;
168 forAll(this->reactions_, i)
170 if (!reactionsDisabled_[i])
174 scalar omegai = omega
176 R,
c,
T,
p, pf, cf, lRef, pr, cr, rRef
181 label si =
R.lhs()[
s].index;
184 si = completeToSimplifiedIndex_[si];
187 const scalar sl =
R.lhs()[
s].stoichCoeff;
188 dcdt[si] -= sl*omegai;
192 label si =
R.rhs()[
s].index;
195 si = completeToSimplifiedIndex_[si];
198 const scalar sr =
R.rhs()[
s].stoichCoeff;
199 dcdt[si] += sr*omegai;
206 template<
class ReactionThermo,
class ThermoType>
221 const scalar kf =
R.kf(
p,
T,
c);
222 const scalar kr =
R.kr(kf,
p,
T,
c);
224 const label Nl =
R.lhs().size();
225 const label Nr =
R.rhs().size();
228 lRef =
R.lhs()[slRef].index;
231 for (label
s=1;
s<Nl;
s++)
233 const label si =
R.lhs()[
s].index;
237 const scalar
exp =
R.lhs()[slRef].exponent;
244 const scalar
exp =
R.lhs()[
s].exponent;
248 cf =
max(
c[lRef], 0.0);
251 const scalar
exp =
R.lhs()[slRef].exponent;
270 rRef =
R.rhs()[srRef].index;
274 for (label
s=1;
s<Nr;
s++)
276 const label si =
R.rhs()[
s].index;
279 const scalar
exp =
R.rhs()[srRef].exponent;
286 const scalar
exp =
R.rhs()[
s].exponent;
290 cr =
max(
c[rRef], 0.0);
293 const scalar
exp =
R.rhs()[srRef].exponent;
311 return pf*cf - pr*cr;
315 template<
class ReactionThermo,
class ThermoType>
323 const bool reduced = mechRed_->active();
325 const scalar
T =
c[this->nSpecie_];
326 const scalar
p =
c[this->nSpecie_ + 1];
333 this->c_ = completeC_;
338 for (label i=0; i<NsDAC_; i++)
340 this->c_[simplifiedToCompleteIndex_[i]] =
max(
c[i], 0.0);
345 for (label i=0; i<this->
nSpecie(); i++)
347 this->c_[i] =
max(
c[i], 0.0);
351 omega(this->c_,
T,
p, dcdt);
356 for (label i=0; i<this->c_.size(); i++)
358 const scalar W = this->specieThermo_[i].W();
359 rho += W*this->c_[i];
363 for (label i=0; i<this->c_.size(); i++)
366 cp += this->c_[i]*this->specieThermo_[i].cp(
p,
T);
374 for (label i=0; i<this->nSpecie_; i++)
379 si = simplifiedToCompleteIndex_[i];
387 const scalar hi = this->specieThermo_[si].ha(
p,
T);
392 dcdt[this->nSpecie_] = -dT;
395 dcdt[this->nSpecie_ + 1] = 0;
399 template<
class ReactionThermo,
class ThermoType>
407 const bool reduced = mechRed_->active();
414 const label
nSpecie = this->nSpecie_;
421 this->c_ = completeC_;
422 for (label i=0; i<NsDAC_; ++i)
424 this->c_[simplifiedToCompleteIndex_[i]] =
max(
c[i], 0.0);
431 this->c_[i] =
max(
c[i], 0.0);
437 forAll(this->reactions_, ri)
439 if (!reactionsDisabled_[ri])
443 const scalar kf0 =
R.kf(
p,
T, this->c_);
444 const scalar kr0 =
R.kr(kf0,
p,
T, this->c_);
448 label sj =
R.lhs()[j].index;
451 sj = completeToSimplifiedIndex_[sj];
456 const label si =
R.lhs()[i].index;
457 const scalar el =
R.lhs()[i].exponent;
462 if (this->c_[si] > SMALL)
464 kf *= el*
pow(this->c_[si], el - 1);
473 kf *= el*
pow(this->c_[si], el - 1);
478 kf *=
pow(this->c_[si], el);
484 label si =
R.lhs()[i].index;
487 si = completeToSimplifiedIndex_[si];
489 const scalar sl =
R.lhs()[i].stoichCoeff;
490 dfdc(si, sj) -= sl*kf;
494 label si =
R.rhs()[i].index;
497 si = completeToSimplifiedIndex_[si];
499 const scalar sr =
R.rhs()[i].stoichCoeff;
500 dfdc(si, sj) += sr*kf;
506 label sj =
R.rhs()[j].index;
509 sj = completeToSimplifiedIndex_[sj];
514 const label si =
R.rhs()[i].index;
515 const scalar er =
R.rhs()[i].exponent;
520 if (this->c_[si] > SMALL)
522 kr *= er*
pow(this->c_[si], er - 1);
531 kr *= er*
pow(this->c_[si], er - 1);
536 kr *=
pow(this->c_[si], er);
542 label si =
R.lhs()[i].index;
545 si = completeToSimplifiedIndex_[si];
547 const scalar sl =
R.lhs()[i].stoichCoeff;
548 dfdc(si, sj) += sl*kr;
552 label si =
R.rhs()[i].index;
555 si = completeToSimplifiedIndex_[si];
557 const scalar sr =
R.rhs()[i].stoichCoeff;
558 dfdc(si, sj) -= sr*kr;
565 const scalar
delta = 1
e-3;
567 omega(this->c_,
T +
delta,
p, this->dcdt_);
568 for (label i=0; i<
nSpecie; ++i)
570 dfdc(i,
nSpecie) = this->dcdt_[i];
573 omega(this->c_,
T -
delta,
p, this->dcdt_);
574 for (label i=0; i<
nSpecie; ++i)
584 template<
class ReactionThermo,
class ThermoType>
593 jacobian(t,
c, dfdc);
595 const scalar
T =
c[this->nSpecie_];
596 const scalar
p =
c[this->nSpecie_ + 1];
599 omega(this->c_,
T,
p, dcdt);
603 template<
class ReactionThermo,
class ThermoType>
604 template<
class DeltaTType>
607 const DeltaTType& deltaT
613 const bool reduced = mechRed_->
active();
615 label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
622 scalar reduceMechCpuTime_ = 0;
623 scalar addNewLeafCpuTime_ = 0;
624 scalar growCpuTime_ = 0;
625 scalar solveChemistryCpuTime_ = 0;
626 scalar searchISATCpuTime_ = 0;
628 this->resetTabulationResults();
631 scalar nActiveSpecies = 0;
636 scalar deltaTMin = GREAT;
638 if (!this->chemistry_)
670 const scalar rhoi =
rho[celli];
671 scalar
pi =
p[celli];
672 scalar Ti =
T[celli];
674 for (label i=0; i<this->nSpecie_; i++)
677 c[i] = rhoi*Yi[celli]/this->specieThermo_[i].W();
683 if (tabulation_->variableTimeStep())
685 phiq[this->
nSpecie() + 2] = deltaT[celli];
690 scalar timeLeft = deltaT[celli];
700 if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
703 for (label i=0; i<this->
nSpecie(); ++i)
705 c[i] = rhoi*Rphiq[i]/this->specieThermo_[i].W();
722 mechRed_->reduceMechanism(
c, Ti,
pi);
723 nActiveSpecies += mechRed_->NsSimp();
726 reduceMechCpuTime_ += timeIncr;
731 while (timeLeft > SMALL)
733 scalar dt = timeLeft;
743 simplifiedC_, Ti,
pi, dt, this->deltaTChem_[celli]
746 for (label i=0; i<NsDAC_; ++i)
748 c[simplifiedToCompleteIndex_[i]] = simplifiedC_[i];
753 this->
solve(c, Ti,
pi, dt, this->deltaTChem_[celli]);
760 solveChemistryCpuTime_ += timeIncr;
766 if (tabulation_->active())
770 Rphiq[i] =
c[i]/rhoi*this->specieThermo_[i].W();
772 if (tabulation_->variableTimeStep())
774 Rphiq[Rphiq.size()-3] = Ti;
775 Rphiq[Rphiq.size()-2] =
pi;
776 Rphiq[Rphiq.size()-1] = deltaT[celli];
780 Rphiq[Rphiq.size()-2] = Ti;
781 Rphiq[Rphiq.size()-1] =
pi;
784 tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]);
788 this->setTabulationResultsAdd(celli);
793 this->setTabulationResultsGrow(celli);
803 this->nSpecie_ = mechRed_->nSpecie();
805 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
807 this->deltaTChem_[celli] =
808 min(this->deltaTChem_[celli], this->deltaTChemMax_);
812 for (label i=0; i<this->nSpecie_; ++i)
814 this->RR_[i][celli] =
815 (
c[i] - c0[i])*this->specieThermo_[i].W()/deltaT[celli];
819 if (mechRed_->log() || tabulation_->log())
822 << this->time().timeOutputValue()
823 <<
" " << solveChemistryCpuTime_ <<
endl;
829 << this->time().timeOutputValue()
830 <<
" " << reduceMechCpuTime_ <<
endl;
833 if (tabulation_->active())
836 tabulation_->update();
839 tabulation_->writePerformance();
841 if (tabulation_->log())
844 << this->time().timeOutputValue()
845 <<
" " << searchISATCpuTime_ <<
endl;
848 << this->time().timeOutputValue()
849 <<
" " << growCpuTime_ <<
endl;
852 << this->time().timeOutputValue()
853 <<
" " << addNewLeafCpuTime_ <<
endl;
857 if (reduced && nAvg && mechRed_->log())
860 nActiveSpeciesFile_()
861 << this->time().timeOutputValue()
862 <<
" " << nActiveSpecies/nAvg <<
endl;
865 if (Pstream::parRun())
869 Pstream::listCombineScatter(active);
884 this->
Y()[i].writeOpt() = IOobject::AUTO_WRITE;
892 template<
class ReactionThermo,
class ThermoType>
907 template<
class ReactionThermo,
class ThermoType>
913 return this->solve<scalarField>(deltaT);
917 template<
class ReactionThermo,
class ThermoType>
924 tabulationResults_[celli] = 0.0;
928 template<
class ReactionThermo,
class ThermoType>
932 tabulationResults_[celli] = 1.0;
936 template<
class ReactionThermo,
class ThermoType>
943 tabulationResults_[celli] = 2.0;