35 template<
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");
132 template<
class CompType,
class ThermoType>
139 template<
class CompType,
class ThermoType>
142 chemPointISAT<CompType, ThermoType>*
phi0 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);
195 template<
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]);
266 template<
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);
302 template<
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());
346 template<
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);
433 template<
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);
508 template<
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();
563 chemPointISAT<CompType, ThermoType>* nulPhi =
nullptr;
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;
607 template<
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;
Switch active_
Is tabulation active?
PtrList< volScalarField > & Y()
dimensionedScalar log(const dimensionedScalar &ds)
errorManipArg< error, int > exit(error &err, const int errNo=1)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Leaf of the binary tree. The chemPoint stores the composition 'phi', the mapping of this composition ...
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type get(const label i) const
Return bool value at specified position, always false for out-of-range access.
const dimensionedScalar phi0
Magnetic flux quantum: default SI units: [Wb].
BasicChemistryModel< psiReactionThermo > & chemistry
virtual label add(const scalarField &phiq, const scalarField &Rphiq, const scalar rho, const scalar deltaT)
Add information to the tabulation.
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
virtual bool retrieve(const Foam::scalarField &phiq, scalarField &Rphiq)
Find the closest stored leaf of phiQ and store the result in RphiQ or return false.
TDACChemistryModel< CompType, ThermoType > & chemistry_
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void append(const T &val)
Copy append an element to the end of this list.
Template functions to aid in the implementation of demand driven data.
An abstract class for chemistry tabulation.
Implementation of the ISAT (In-situ adaptive tabulation), for chemistry calculation.
const dictionary coeffsDict_
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
void deleteDemandDrivenData(DataPtr &dataPtr)
SquareMatrix< scalar > scalarSquareMatrix
forAllConstIters(mixture.phases(), phase)
virtual void writePerformance()
static constexpr const zero Zero
Global zero (0)