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].name()))
89 scaleFactor_[i] = otherScaleFactor;
96 scaleDict.lookup(this->chemistry_.Y()[i].name())
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 if (iter() == MRUList_.
last())
174 MRUList_.remove(iter);
175 MRUList_.insert(phi0);
180 <<
"wrong MRUList construction" 186 MRUList_.insert(phi0);
193 template<
class CompType,
class ThermoType>
201 label nEqns = this->chemistry_.nEqns();
202 bool mechRedActive = this->chemistry_.mechRed()->
active();
203 Rphiq = phi0->
Rphi();
210 for (
label i=0; i<nEqns-nAdditionalEqns_; i++)
214 label si = completeToSimplified[i];
218 for (
label j=0; j<nEqns-2; j++)
220 label sj = completeToSimplified[j];
223 Rphiq[i] += gradientsMatrix(si, sj)*dphi[j];
232 if (this->variableTimeStep())
241 Rphiq[i] =
max(0.0,Rphiq[i]);
247 Rphiq[i] =
max(0.0,Rphiq[i]);
252 for (
label j=0; j<nEqns; j++)
254 Rphiq[i] += gradientsMatrix(i, j)*dphi[j];
258 Rphiq[i] =
max(0.0,Rphiq[i]);
264 template<
class CompType,
class ThermoType>
280 if (phi0->
nGrowth() > maxGrowth_)
282 cleaningRequired_ =
true;
292 return phi0->
grow(phiq);
302 template<
class CompType,
class ThermoType>
306 bool treeModified(
false);
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;
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 >::iterator iter = MRUList_.begin();
468 for ( ; iter != MRUList_.end(); ++iter)
471 if (phi0->
inEOA(phiq))
483 lastSearch_ =
nullptr;
489 scalar elapsedTimeSteps =
490 this->chemistry_.timeSteps() - phi0->
timeTag();
494 if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->
toRemove())
496 cleaningRequired_ =
true;
499 lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
501 calcNewC(phi0,phiq, Rphiq);
514 template<
class CompType,
class ThermoType>
523 label growthOrAddFlag = 1;
526 if (lastSearch_ && growPoints_)
528 if (grow(lastSearch_,phiq, Rphiq))
533 return growthOrAddFlag;
540 if (chemisTree().isFull())
545 if (!cleanAndBalance())
555 >::iterator iter = MRUList_.begin();
556 for ( ; iter != MRUList_.end(); ++iter)
564 chemisTree().clear();
574 chemisTree().insertNewLeaf
590 lastSearch_ =
nullptr;
594 label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
596 computeA(A, Rphiq, rho, deltaT);
598 chemisTree().insertNewLeaf
611 return growthOrAddFlag;
615 template<
class CompType,
class ThermoType>
622 << runTime_.timeOutputValue() <<
" " << nRetrieved_ <<
endl;
626 << runTime_.timeOutputValue() <<
" " << nGrowth_ <<
endl;
630 << runTime_.timeOutputValue() <<
" " << nAdd_ <<
endl;
634 << runTime_.timeOutputValue() <<
" " << this->size() <<
endl;
Extends chemistryModel by adding the TDAC method.
#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)
errorManipArg< error, int > exit(error &err, const int errNo=1)
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 > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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()
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.
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.
virtual label nEqns() const
Number of ODE's to solve.
const scalarField & phi() const
T & last()
Return the last entry added.
const scalarField & Rphi() const
An abstract class for chemistry tabulation.
psiChemistryModel & chemistry
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()