31 template<
class CompType,
class ThermoType>
43 chemisTree_(chemistry, this->coeffsDict_),
44 scaleFactor_(chemistry.
nEqns() + ((this->variableTimeStep()) ? 1 : 0), 1),
45 runTime_(chemistry.
time()),
48 this->coeffsDict_.lookupOrDefault(
"chPMaxLifeTime", INT_MAX)
50 maxGrowth_(this->coeffsDict_.lookupOrDefault(
"maxGrowth", INT_MAX)),
51 checkEntireTreeInterval_
53 this->coeffsDict_.lookupOrDefault(
"checkEntireTreeInterval", INT_MAX)
57 this->coeffsDict_.lookupOrDefault
60 (chemisTree_.maxNLeafs() - 1)
61 /(
log(scalar(chemisTree_.maxNLeafs()))/
log(2.0))
66 this->coeffsDict_.lookupOrDefault
68 "minBalanceThreshold", 0.1*chemisTree_.maxNLeafs()
71 MRURetrieve_(this->coeffsDict_.lookupOrDefault(
"MRURetrieve",
false)),
72 maxMRUSize_(this->coeffsDict_.lookupOrDefault(
"maxMRUSize", 0)),
74 growPoints_(this->coeffsDict_.lookupOrDefault(
"growPoints",
true)),
78 cleaningRequired_(
false)
82 dictionary scaleDict(this->coeffsDict_.subDict(
"scaleFactor"));
83 label Ysize = this->chemistry_.Y().size();
84 scalar otherScaleFactor =
readScalar(scaleDict.lookup(
"otherSpecies"));
85 for (
label i=0; i<Ysize; i++)
87 if (!scaleDict.found(this->chemistry_.Y()[i].member()))
89 scaleFactor_[i] = otherScaleFactor;
96 scaleDict.lookup(this->chemistry_.Y()[i].member())
100 scaleFactor_[Ysize] =
readScalar(scaleDict.lookup(
"Temperature"));
101 scaleFactor_[Ysize + 1] =
readScalar(scaleDict.lookup(
"Pressure"));
102 if (this->variableTimeStep())
104 scaleFactor_[Ysize + 2] =
readScalar(scaleDict.lookup(
"deltaT"));
108 if (this->variableTimeStep())
110 nAdditionalEqns_ = 3;
114 nAdditionalEqns_ = 2;
119 nRetrievedFile_ = chemistry.
logFile(
"found_isat.out");
120 nGrowthFile_ = chemistry.
logFile(
"growth_isat.out");
121 nAddFile_ = chemistry.
logFile(
"add_isat.out");
122 sizeFile_ = chemistry.
logFile(
"size_isat.out");
129 template<
class CompType,
class ThermoType>
136 template<
class CompType,
class ThermoType>
142 if (maxMRUSize_ > 0 && MRURetrieve_)
145 bool isInList =
false;
148 for ( ; iter != MRUList_.
end(); ++iter)
159 if (iter() != MRUList_.
first())
162 MRUList_.remove(iter);
165 MRUList_.insert(phi0);
170 if (MRUList_.size() == maxMRUSize_)
172 MRUList_.remove(iter);
173 MRUList_.insert(phi0);
177 MRUList_.insert(phi0);
184 template<
class CompType,
class ThermoType>
192 label nEqns = this->chemistry_.nEqns();
193 bool mechRedActive = this->chemistry_.mechRed()->
active();
194 Rphiq = phi0->
Rphi();
201 for (
label i=0; i<nEqns-nAdditionalEqns_; i++)
205 label si = completeToSimplified[i];
209 for (
label j=0; j<nEqns-2; j++)
211 label sj = completeToSimplified[j];
214 Rphiq[i] += gradientsMatrix(si, sj)*dphi[j];
223 if (this->variableTimeStep())
232 Rphiq[i] =
max(0, Rphiq[i]);
238 Rphiq[i] =
max(0, Rphiq[i]);
243 for (
label j=0; j<nEqns; j++)
245 Rphiq[i] += gradientsMatrix(i, j)*dphi[j];
249 Rphiq[i] =
max(0, Rphiq[i]);
255 template<
class CompType,
class ThermoType>
271 if (phi0->
nGrowth() > maxGrowth_)
273 cleaningRequired_ =
true;
283 return phi0->
grow(phiq);
293 template<
class CompType,
class ThermoType>
297 bool treeModified(
false);
305 chemisTree_.treeSuccessor(x);
307 scalar elapsedTimeSteps = this->chemistry_.timeSteps() - x->
timeTag();
309 if ((elapsedTimeSteps > chPMaxLifeTime_) || (x->
nGrowth() > maxGrowth_))
311 chemisTree_.deleteLeaf(x);
324 chemisTree_.size() > minBalanceThreshold_
325 && chemisTree_.depth() >
326 maxDepthFactor_*
log(scalar(chemisTree_.size()))/
log(2.0)
329 chemisTree_.balance();
335 return (treeModified && !chemisTree_.isFull());
339 template<
class CompType,
class ThermoType>
348 bool mechRedActive = this->chemistry_.mechRed()->
active();
349 label speciesNumber = this->chemistry_.nSpecie();
350 scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
351 for (
label i=0; i<speciesNumber; i++)
356 s2c = this->chemistry_.simplifiedToCompleteIndex()[i];
358 Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
360 Rcq[speciesNumber] = Rphiq[Rphiq.
size() - nAdditionalEqns_];
361 Rcq[speciesNumber + 1] = Rphiq[Rphiq.
size() - nAdditionalEqns_ + 1];
362 if (this->variableTimeStep())
364 Rcq[speciesNumber + 2] = Rphiq[Rphiq.
size() - nAdditionalEqns_ + 2];
377 this->chemistry_.jacobian(runTime_.value(), Rcq, dcdt, A);
381 for (
label i=0; i<speciesNumber; i++)
387 si = this->chemistry_.simplifiedToCompleteIndex()[i];
390 for (
label j=0; j<speciesNumber; j++)
395 sj = this->chemistry_.simplifiedToCompleteIndex()[j];
398 -dt*this->chemistry_.specieThermo()[si].W()
399 /this->chemistry_.specieThermo()[sj].W();
404 A(i, speciesNumber) *=
405 -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
406 A(i, speciesNumber + 1) *=
407 -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
412 for (
label i=0; i<speciesNumber; i++)
417 si = this->chemistry_.simplifiedToCompleteIndex()[i];
420 A(speciesNumber, i) *=
421 -dt*rhoi/this->chemistry_.specieThermo()[si].W();
422 A(speciesNumber + 1, i) *=
423 -dt*rhoi/this->chemistry_.specieThermo()[si].W();
426 A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
428 A(speciesNumber + 1, speciesNumber + 1) =
429 -dt*A(speciesNumber + 1, speciesNumber + 1) + 1;
431 if (this->variableTimeStep())
433 A(speciesNumber + 2, speciesNumber + 2) = 1;
443 for (
label i=0; i<speciesNumber; i++)
445 A(speciesNumber, i) = 0;
446 A(speciesNumber + 1, i) = 0;
453 template<
class CompType,
class ThermoType>
460 bool retrieved(
false);
464 if (chemisTree_.size())
466 chemisTree_.binaryTreeSearch(phiq, chemisTree_.root(),
phi0);
471 if (phi0->
inEOA(phiq))
477 else if (chemisTree_.secondaryBTSearch(phiq, phi0))
481 else if (MRURetrieve_)
486 >::iterator iter = MRUList_.begin();
488 for ( ; iter != MRUList_.end(); ++iter)
491 if (phi0->
inEOA(phiq))
503 lastSearch_ =
nullptr;
509 scalar elapsedTimeSteps =
510 this->chemistry_.timeSteps() - phi0->
timeTag();
514 if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->
toRemove())
516 cleaningRequired_ =
true;
519 lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
521 calcNewC(phi0, phiq, Rphiq);
534 template<
class CompType,
class ThermoType>
543 label growthOrAddFlag = 1;
546 if (lastSearch_ && growPoints_)
548 if (grow(lastSearch_, phiq, Rphiq))
552 addToMRU(lastSearch_);
554 return growthOrAddFlag;
561 if (chemisTree().isFull())
566 if (!cleanAndBalance())
576 >::iterator iter = MRUList_.begin();
577 for ( ; iter != MRUList_.end(); ++iter)
585 chemisTree().clear();
595 chemisTree().insertNewLeaf
611 lastSearch_ =
nullptr;
615 label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
617 computeA(A, Rphiq, rho, deltaT);
619 chemisTree().insertNewLeaf
629 if (lastSearch_ !=
nullptr)
631 addToMRU(lastSearch_);
635 return growthOrAddFlag;
639 template<
class CompType,
class ThermoType>
646 << runTime_.timeOutputValue() <<
" " << nRetrieved_ <<
endl;
650 << runTime_.timeOutputValue() <<
" " << nGrowth_ <<
endl;
654 << runTime_.timeOutputValue() <<
" " << nAdd_ <<
endl;
658 << runTime_.timeOutputValue() <<
" " << this->size() <<
endl;
#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.
const scalarSquareMatrix & A() const
dimensionedScalar log(const dimensionedScalar &ds)
A list of keyword definitions, which are a keyword followed by any number of values (e...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Class to perform the LU decomposition on a symmetric matrix.
Template class for non-intrusive linked lists.
T & first()
Return the first entry added.
void size(const label)
Override size to be inconsistent with allocated storage.
bool checkSolution(const scalarField &phiq, const scalarField &Rphiq)
If phiq is not in the EOA, then the mapping is computed.
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 ...
List< label > & completeToSimplifiedIndex()
ISAT(const dictionary &chemistryProperties, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from dictionary.
const dimensionedScalar phi0
Magnetic flux quantum: default SI units: [Wb].
virtual label add(const scalarField &phiq, const scalarField &Rphiq, const scalar rho, const scalar deltaT)
Add information to the tabulation.
BasicChemistryModel< rhoReactionThermo > & chemistry
virtual bool retrieve(const Foam::scalarField &phiq, scalarField &Rphiq)
Find the closest stored leaf of phiQ and store the result in.
label size() const
Return the number of elements in matrix (m*n)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
bool grow(const scalarField &phiq)
More details about the minimum-volume ellipsoid covering an.
const scalarField & phi() const
const scalarField & Rphi() const
An abstract class for chemistry tabulation.
virtual label nEqns() const
Number of ODE's to solve.
const Time & time() const
Return time.
void increaseNumRetrieve()
Increases the number of retrieves the chempoint has generated.
bool inEOA(const scalarField &phiq)
To RETRIEVE the mapping from the stored chemPoint phi, the query.
autoPtr< OFstream > logFile(const word &name) const
Create and return a TDAC log file of the given name.
Implementation of the ISAT (In-situ adaptive tabulation), for chemistry calculation.
void deleteDemandDrivenData(DataPtr &dataPtr)
void inv(scalarSquareMatrix &M) const
Set M to the inverse of this square matrix.
virtual void writePerformance()