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 = scaleDict.
lookup<scalar>(
"otherSpecies");
85 for (
label i=0; i<Ysize; i++)
87 if (!scaleDict.found(this->chemistry_.Y()[i].member()))
89 scaleFactor_[i] = otherScaleFactor;
94 scaleDict.lookup<scalar>(this->chemistry_.Y()[i].member());
97 scaleFactor_[Ysize] = scaleDict.lookup<scalar>(
"Temperature");
98 scaleFactor_[Ysize + 1] = scaleDict.lookup<scalar>(
"Pressure");
99 if (this->variableTimeStep())
101 scaleFactor_[Ysize + 2] = scaleDict.lookup<scalar>(
"deltaT");
105 if (this->variableTimeStep())
107 nAdditionalEqns_ = 3;
111 nAdditionalEqns_ = 2;
116 nRetrievedFile_ = chemistry.
logFile(
"found_isat.out");
117 nGrowthFile_ = chemistry.
logFile(
"growth_isat.out");
118 nAddFile_ = chemistry.
logFile(
"add_isat.out");
119 sizeFile_ = chemistry.
logFile(
"size_isat.out");
126 template<
class CompType,
class ThermoType>
133 template<
class CompType,
class ThermoType>
139 if (maxMRUSize_ > 0 && MRURetrieve_)
142 bool isInList =
false;
145 for ( ; iter != MRUList_.
end(); ++iter)
156 if (iter() != MRUList_.
first())
159 MRUList_.remove(iter);
162 MRUList_.insert(phi0);
167 if (MRUList_.size() == maxMRUSize_)
169 MRUList_.remove(iter);
170 MRUList_.insert(phi0);
174 MRUList_.insert(phi0);
181 template<
class CompType,
class ThermoType>
189 label nEqns = this->chemistry_.nEqns();
190 bool mechRedActive = this->chemistry_.mechRed()->
active();
191 Rphiq = phi0->
Rphi();
198 for (
label i=0; i<nEqns-nAdditionalEqns_; i++)
202 label si = completeToSimplified[i];
206 for (
label j=0; j<nEqns-2; j++)
208 label sj = completeToSimplified[j];
211 Rphiq[i] += gradientsMatrix(si, sj)*dphi[j];
220 if (this->variableTimeStep())
229 Rphiq[i] =
max(0, Rphiq[i]);
235 Rphiq[i] =
max(0, Rphiq[i]);
240 for (
label j=0; j<nEqns; j++)
242 Rphiq[i] += gradientsMatrix(i, j)*dphi[j];
246 Rphiq[i] =
max(0, Rphiq[i]);
252 template<
class CompType,
class ThermoType>
268 if (phi0->
nGrowth() > maxGrowth_)
270 cleaningRequired_ =
true;
280 return phi0->
grow(phiq);
290 template<
class CompType,
class ThermoType>
294 bool treeModified(
false);
302 chemisTree_.treeSuccessor(x);
304 scalar elapsedTimeSteps = this->chemistry_.timeSteps() - x->
timeTag();
306 if ((elapsedTimeSteps > chPMaxLifeTime_) || (x->
nGrowth() > maxGrowth_))
308 chemisTree_.deleteLeaf(x);
321 chemisTree_.size() > minBalanceThreshold_
322 && chemisTree_.depth() >
323 maxDepthFactor_*
log(scalar(chemisTree_.size()))/
log(2.0)
326 chemisTree_.balance();
332 return (treeModified && !chemisTree_.isFull());
336 template<
class CompType,
class ThermoType>
346 bool mechRedActive = this->chemistry_.mechRed()->
active();
347 label speciesNumber = this->chemistry_.nSpecie();
348 scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
349 for (
label i=0; i<speciesNumber; i++)
354 s2c = this->chemistry_.simplifiedToCompleteIndex()[i];
356 Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermos()[s2c].W();
358 Rcq[speciesNumber] = Rphiq[Rphiq.
size() - nAdditionalEqns_];
359 Rcq[speciesNumber + 1] = Rphiq[Rphiq.
size() - nAdditionalEqns_ + 1];
360 if (this->variableTimeStep())
362 Rcq[speciesNumber + 2] = Rphiq[Rphiq.
size() - nAdditionalEqns_ + 2];
375 this->chemistry_.jacobian(runTime_.value(), Rcq, li, dcdt, A);
379 for (
label i=0; i<speciesNumber; i++)
385 si = this->chemistry_.simplifiedToCompleteIndex()[i];
388 for (
label j=0; j<speciesNumber; j++)
393 sj = this->chemistry_.simplifiedToCompleteIndex()[j];
396 -dt*this->chemistry_.specieThermos()[si].W()
397 /this->chemistry_.specieThermos()[sj].W();
402 A(i, speciesNumber) *=
403 -dt*this->chemistry_.specieThermos()[si].W()/rhoi;
404 A(i, speciesNumber + 1) *=
405 -dt*this->chemistry_.specieThermos()[si].W()/rhoi;
410 for (
label i=0; i<speciesNumber; i++)
415 si = this->chemistry_.simplifiedToCompleteIndex()[i];
418 A(speciesNumber, i) *=
419 -dt*rhoi/this->chemistry_.specieThermos()[si].W();
420 A(speciesNumber + 1, i) *=
421 -dt*rhoi/this->chemistry_.specieThermos()[si].W();
424 A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
426 A(speciesNumber + 1, speciesNumber + 1) =
427 -dt*A(speciesNumber + 1, speciesNumber + 1) + 1;
429 if (this->variableTimeStep())
431 A(speciesNumber + 2, speciesNumber + 2) = 1;
441 for (
label i=0; i<speciesNumber; i++)
443 A(speciesNumber, i) = 0;
444 A(speciesNumber + 1, i) = 0;
451 template<
class CompType,
class ThermoType>
458 bool retrieved(
false);
462 if (chemisTree_.size())
464 chemisTree_.binaryTreeSearch(phiq, chemisTree_.root(),
phi0);
469 if (phi0->
inEOA(phiq))
475 else if (chemisTree_.secondaryBTSearch(phiq, phi0))
479 else if (MRURetrieve_)
484 >::iterator iter = MRUList_.begin();
486 for ( ; iter != MRUList_.end(); ++iter)
489 if (phi0->
inEOA(phiq))
501 lastSearch_ =
nullptr;
507 scalar elapsedTimeSteps =
508 this->chemistry_.timeSteps() - phi0->
timeTag();
512 if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->
toRemove())
514 cleaningRequired_ =
true;
517 lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
519 calcNewC(phi0, phiq, Rphiq);
532 template<
class CompType,
class ThermoType>
542 label growthOrAddFlag = 1;
545 if (lastSearch_ && growPoints_)
547 if (grow(lastSearch_, phiq, Rphiq))
551 addToMRU(lastSearch_);
553 return growthOrAddFlag;
560 if (chemisTree().isFull())
565 if (!cleanAndBalance())
575 >::iterator iter = MRUList_.begin();
576 for ( ; iter != MRUList_.end(); ++iter)
584 chemisTree().clear();
594 chemisTree().insertNewLeaf
610 lastSearch_ =
nullptr;
614 label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
616 computeA(A, Rphiq, li, rho, deltaT);
618 chemisTree().insertNewLeaf
628 if (lastSearch_ !=
nullptr)
630 addToMRU(lastSearch_);
634 return growthOrAddFlag;
638 template<
class CompType,
class ThermoType>
645 << runTime_.timeOutputValue() <<
" " << nRetrieved_ <<
endl;
649 << runTime_.timeOutputValue() <<
" " << nGrowth_ <<
endl;
653 << runTime_.timeOutputValue() <<
" " << nAdd_ <<
endl;
657 << 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.
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...
virtual label add(const scalarField &phiq, const scalarField &Rphiq, const label li, const scalar rho, const scalar deltaT)
Add information to the tabulation.
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 dimensionedScalar & phi0
Magnetic flux quantum: default SI units: [Wb].
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)
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
void inv(scalarSquareMatrix &M) const
Set M to the inverse of this square matrix.
virtual void writePerformance()