35template<
class CompType,
class ThermoType>
47 chemisTree_(
chemistry, this->coeffsDict_),
48 scaleFactor_(
chemistry.nEqns() + ((this->variableTimeStep()) ? 1 : 0), 1),
52 this->coeffsDict_.getOrDefault(
"chPMaxLifeTime", INT_MAX)
54 maxGrowth_(this->coeffsDict_.getOrDefault(
"maxGrowth", INT_MAX)),
55 checkEntireTreeInterval_
57 this->coeffsDict_.getOrDefault(
"checkEntireTreeInterval", INT_MAX)
61 this->coeffsDict_.getOrDefault
64 (chemisTree_.maxNLeafs() - 1)
65 /(
log(scalar(chemisTree_.maxNLeafs()))/
log(2.0))
70 this->coeffsDict_.getOrDefault
72 "minBalanceThreshold",0.1*chemisTree_.maxNLeafs()
75 MRURetrieve_(this->coeffsDict_.getOrDefault(
"MRURetrieve", false)),
76 maxMRUSize_(this->coeffsDict_.getOrDefault(
"maxMRUSize", 0)),
78 growPoints_(this->coeffsDict_.getOrDefault(
"growPoints", true)),
82 cleaningRequired_(false)
88 const scalar otherScaleFactor = scaleDict.
get<scalar>(
"otherSpecies");
89 for (label i=0; i<Ysize; ++i)
91 if (!scaleDict.
found(this->chemistry_.Y()[i].member()))
93 scaleFactor_[i] = otherScaleFactor;
102 scaleDict.
readEntry(
"Temperature", scaleFactor_[Ysize]);
103 scaleDict.
readEntry(
"Pressure", scaleFactor_[Ysize + 1]);
107 scaleDict.
readEntry(
"deltaT", scaleFactor_[Ysize + 2]);
113 nAdditionalEqns_ = 3;
117 nAdditionalEqns_ = 2;
122 nRetrievedFile_ =
chemistry.logFile(
"found_isat.out");
123 nGrowthFile_ =
chemistry.logFile(
"growth_isat.out");
124 nAddFile_ =
chemistry.logFile(
"add_isat.out");
125 sizeFile_ =
chemistry.logFile(
"size_isat.out");
132template<
class CompType,
class ThermoType>
139template<
class CompType,
class ThermoType>
145 if (maxMRUSize_ > 0 && MRURetrieve_)
147 auto iter = MRUList_.begin();
150 bool isInList =
false;
151 for (; iter.good(); ++iter)
163 if (iter != MRUList_.begin())
165 MRUList_.remove(iter);
166 MRUList_.insert(phi0);
172 if (MRUList_.size() == maxMRUSize_)
174 if (iter == MRUList_.end())
176 MRUList_.remove(iter);
177 MRUList_.insert(phi0);
182 <<
"Error in MRUList construction"
188 MRUList_.insert(phi0);
195template<
class CompType,
class ThermoType>
198 chemPointISAT<CompType, ThermoType>* phi0,
203 label nEqns = this->chemistry_.nEqns();
204 bool mechRedActive = this->chemistry_.mechRed()->
active();
205 Rphiq =
phi0->Rphi();
208 List<label>& completeToSimplified(
phi0->completeToSimplifiedIndex());
212 for (label i=0; i<nEqns-nAdditionalEqns_; ++i)
216 label si = completeToSimplified[i];
220 for (label j=0; j<nEqns-2; j++)
222 label sj = completeToSimplified[j];
225 Rphiq[i] += gradientsMatrix(si, sj)*dphi[j];
229 gradientsMatrix(si,
phi0->nActiveSpecies())*dphi[nEqns - 2];
231 gradientsMatrix(si,
phi0->nActiveSpecies() + 1)
234 if (this->variableTimeStep())
237 gradientsMatrix(si,
phi0->nActiveSpecies() + 2)
243 Rphiq[i] =
max(0.0, Rphiq[i]);
249 Rphiq[i] =
max(0.0, Rphiq[i]);
254 for (label j=0; j<nEqns; ++j)
256 Rphiq[i] += gradientsMatrix(i, j)*dphi[j];
260 Rphiq[i] =
max(0.0, Rphiq[i]);
266template<
class CompType,
class ThermoType>
269 chemPointISAT<CompType, ThermoType>* phi0,
282 if (
phi0->nGrowth() > maxGrowth_)
284 cleaningRequired_ =
true;
285 phi0->toRemove() =
true;
292 if (
phi0->checkSolution(phiq, Rphiq))
294 return phi0->grow(phiq);
302template<
class CompType,
class ThermoType>
306 bool treeModified(
false);
310 chemPointISAT<CompType, ThermoType>*
x = chemisTree_.treeMin();
313 chemPointISAT<CompType, ThermoType>* xtmp =
314 chemisTree_.treeSuccessor(
x);
316 scalar elapsedTimeSteps = this->chemistry_.timeSteps() -
x->timeTag();
318 if ((elapsedTimeSteps > chPMaxLifeTime_) || (
x->nGrowth() > maxGrowth_))
320 chemisTree_.deleteLeaf(
x);
330 chemisTree_.size() > minBalanceThreshold_
331 && chemisTree_.depth() >
332 maxDepthFactor_*
log(scalar(chemisTree_.size()))/
log(2.0)
335 chemisTree_.balance();
342 return (treeModified && !chemisTree_.isFull());
346template<
class CompType,
class ThermoType>
355 bool mechRedActive = this->chemistry_.mechRed()->
active();
356 label speciesNumber = this->chemistry_.nSpecie();
357 scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
358 for (label i=0; i<speciesNumber; ++i)
363 s2c = this->chemistry_.simplifiedToCompleteIndex()[i];
365 Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
367 Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
368 Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
369 if (this->variableTimeStep())
371 Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
384 this->chemistry_.jacobian(runTime_.value(), Rcq,
A);
388 for (label i=0; i<speciesNumber; ++i)
394 si = this->chemistry_.simplifiedToCompleteIndex()[i];
397 for (label j=0; j<speciesNumber; ++j)
402 sj = this->chemistry_.simplifiedToCompleteIndex()[j];
405 -dt*this->chemistry_.specieThermo()[si].W()
406 /this->chemistry_.specieThermo()[sj].W();
411 A(i, speciesNumber) *=
412 -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
413 A(i, speciesNumber+1) *=
414 -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
418 A(speciesNumber, speciesNumber) = 1;
419 A(speciesNumber + 1, speciesNumber + 1) = 1;
420 if (this->variableTimeStep())
422 A[speciesNumber + 2][speciesNumber + 2] = 1;
426 LUscalarMatrix LUA(
A);
433template<
class CompType,
class ThermoType>
440 bool retrieved(
false);
444 if (chemisTree_.size())
446 chemisTree_.binaryTreeSearch(phiq, chemisTree_.root(), phi0);
451 if (phi0->inEOA(phiq))
457 else if (chemisTree_.secondaryBTSearch(phiq, phi0))
461 else if (MRURetrieve_)
466 if (phi0->inEOA(phiq))
478 lastSearch_ =
nullptr;
483 phi0->increaseNumRetrieve();
484 scalar elapsedTimeSteps =
485 this->chemistry_.timeSteps() - phi0->timeTag();
489 if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->toRemove())
491 cleaningRequired_ =
true;
492 phi0->toRemove() =
true;
494 lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
496 calcNewC(phi0,phiq, Rphiq);
508template<
class CompType,
class ThermoType>
517 label growthOrAddFlag = 1;
521 if (lastSearch_ && growPoints_)
523 if (grow(lastSearch_,phiq, Rphiq))
529 return growthOrAddFlag;
536 if (chemisTree().isFull())
541 if (!cleanAndBalance())
556 chemisTree().clear();
564 for (
auto& t : tempList)
566 chemisTree().insertNewLeaf
582 lastSearch_ =
nullptr;
586 label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
588 computeA(
A, Rphiq,
rho, deltaT);
590 chemisTree().insertNewLeaf
603 return growthOrAddFlag;
607template<
class CompType,
class ThermoType>
614 << runTime_.timeOutputValue() <<
" " << nRetrieved_ <<
endl;
618 << runTime_.timeOutputValue() <<
" " << nGrowth_ <<
endl;
622 << runTime_.timeOutputValue() <<
" " << nAdd_ <<
endl;
626 << runTime_.timeOutputValue() <<
" " << this->size() <<
endl;
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
Extends StandardChemistryModel by adding the TDAC method.
PtrList< volScalarField > & Y()
Leaf of the binary tree. The chemPoint stores the composition 'phi', the mapping of this composition ...
An abstract class for chemistry tabulation.
Switch active_
Is tabulation active?
TDACChemistryModel< CompType, ThermoType > & chemistry_
const dictionary coeffsDict_
Implementation of the ISAT (In-situ adaptive tabulation), for chemistry calculation.
virtual bool retrieve(const Foam::scalarField &phiq, scalarField &Rphiq)
virtual void writePerformance()
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
Sums a given list of (at least two or more) fields and outputs the result into a new field,...
BasicChemistryModel< psiReactionThermo > & chemistry
Template functions to aid in the implementation of demand driven data.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const dimensionedScalar phi0
Magnetic flux quantum: default SI units: [Wb].
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar log(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
SquareMatrix< scalar > scalarSquareMatrix
static constexpr const zero Zero
Global zero (0)
errorManipArg< error, int > exit(error &err, const int errNo=1)
void deleteDemandDrivenData(DataPtr &dataPtr)
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.