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)
86 dictionary scaleDict(this->coeffsDict_.subDict(
"scaleFactor"));
87 const label Ysize = this->chemistry_.Y().size();
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;
98 scaleDict.get<scalar>(this->chemistry_.Y()[i].member());
102 scaleDict.readEntry(
"Temperature", scaleFactor_[Ysize]);
103 scaleDict.readEntry(
"Pressure", scaleFactor_[Ysize + 1]);
105 if (this->variableTimeStep())
107 scaleDict.readEntry(
"deltaT", scaleFactor_[Ysize + 2]);
111 if (this->variableTimeStep())
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>
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();
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;