34 template<
class CompType,
class ThermoType>
46 chemisTree_(
chemistry, this->coeffsDict_),
47 scaleFactor_(
chemistry.nEqns() + ((this->variableTimeStep()) ? 1 : 0), 1),
51 this->coeffsDict_.lookupOrDefault(
"chPMaxLifeTime", INT_MAX)
53 maxGrowth_(this->coeffsDict_.lookupOrDefault(
"maxGrowth", INT_MAX)),
54 checkEntireTreeInterval_
56 this->coeffsDict_.lookupOrDefault(
"checkEntireTreeInterval", INT_MAX)
60 this->coeffsDict_.lookupOrDefault
63 (chemisTree_.maxNLeafs() - 1)
64 /(
log(scalar(chemisTree_.maxNLeafs()))/
log(2.0))
69 this->coeffsDict_.lookupOrDefault
71 "minBalanceThreshold",0.1*chemisTree_.maxNLeafs()
74 MRURetrieve_(this->coeffsDict_.lookupOrDefault(
"MRURetrieve",
false)),
75 maxMRUSize_(this->coeffsDict_.lookupOrDefault(
"maxMRUSize", 0)),
77 growPoints_(this->coeffsDict_.lookupOrDefault(
"growPoints",
true)),
81 cleaningRequired_(
false)
85 dictionary scaleDict(this->coeffsDict_.subDict(
"scaleFactor"));
86 const label Ysize = this->chemistry_.Y().size();
87 const scalar otherScaleFactor = scaleDict.
get<scalar>(
"otherSpecies");
88 for (
label i=0; i<Ysize; ++i)
90 if (!scaleDict.found(this->chemistry_.Y()[i].member()))
92 scaleFactor_[i] = otherScaleFactor;
97 scaleDict.get<scalar>(this->chemistry_.Y()[i].member());
101 scaleDict.readEntry(
"Temperature", scaleFactor_[Ysize]);
102 scaleDict.readEntry(
"Pressure", scaleFactor_[Ysize + 1]);
104 if (this->variableTimeStep())
106 scaleDict.readEntry(
"deltaT", scaleFactor_[Ysize + 2]);
110 if (this->variableTimeStep())
112 nAdditionalEqns_ = 3;
116 nAdditionalEqns_ = 2;
121 nRetrievedFile_ =
chemistry.logFile(
"found_isat.out");
122 nGrowthFile_ =
chemistry.logFile(
"growth_isat.out");
123 nAddFile_ =
chemistry.logFile(
"add_isat.out");
124 sizeFile_ =
chemistry.logFile(
"size_isat.out");
131 template<
class CompType,
class ThermoType>
138 template<
class CompType,
class ThermoType>
144 if (maxMRUSize_ > 0 && MRURetrieve_)
146 auto iter = MRUList_.begin();
149 bool isInList =
false;
150 for (; iter.good(); ++iter)
162 if (iter != MRUList_.begin())
164 MRUList_.remove(iter);
165 MRUList_.insert(
phi0);
171 if (MRUList_.size() == maxMRUSize_)
173 if (iter == MRUList_.end())
175 MRUList_.remove(iter);
176 MRUList_.insert(
phi0);
181 <<
"Error in MRUList construction"
187 MRUList_.insert(
phi0);
194 template<
class CompType,
class ThermoType>
197 chemPointISAT<CompType, ThermoType>*
phi0,
202 label nEqns = this->chemistry_.nEqns();
203 bool mechRedActive = this->chemistry_.mechRed()->active();
204 Rphiq =
phi0->Rphi();
207 List<label>& completeToSimplified(
phi0->completeToSimplifiedIndex());
211 for (
label i=0; i<nEqns-nAdditionalEqns_; ++i)
215 label si = completeToSimplified[i];
219 for (
label j=0; j<nEqns-2; j++)
221 label sj = completeToSimplified[j];
224 Rphiq[i] += gradientsMatrix(si, sj)*dphi[j];
228 gradientsMatrix(si,
phi0->nActiveSpecies())*dphi[nEqns - 2];
230 gradientsMatrix(si,
phi0->nActiveSpecies() + 1)
233 if (this->variableTimeStep())
236 gradientsMatrix(si,
phi0->nActiveSpecies() + 2)
242 Rphiq[i] =
max(0.0, Rphiq[i]);
248 Rphiq[i] =
max(0.0, Rphiq[i]);
253 for (
label j=0; j<nEqns; ++j)
255 Rphiq[i] += gradientsMatrix(i, j)*dphi[j];
259 Rphiq[i] =
max(0.0, Rphiq[i]);
265 template<
class CompType,
class ThermoType>
268 chemPointISAT<CompType, ThermoType>*
phi0,
281 if (
phi0->nGrowth() > maxGrowth_)
283 cleaningRequired_ =
true;
284 phi0->toRemove() =
true;
291 if (
phi0->checkSolution(phiq, Rphiq))
293 return phi0->grow(phiq);
301 template<
class CompType,
class ThermoType>
305 bool treeModified(
false);
309 chemPointISAT<CompType, ThermoType>*
x = chemisTree_.treeMin();
312 chemPointISAT<CompType, ThermoType>* xtmp =
313 chemisTree_.treeSuccessor(
x);
315 scalar elapsedTimeSteps = this->chemistry_.timeSteps() -
x->timeTag();
317 if ((elapsedTimeSteps > chPMaxLifeTime_) || (
x->nGrowth() > maxGrowth_))
319 chemisTree_.deleteLeaf(
x);
329 chemisTree_.size() > minBalanceThreshold_
330 && chemisTree_.depth() >
331 maxDepthFactor_*
log(scalar(chemisTree_.size()))/
log(2.0)
334 chemisTree_.balance();
341 return (treeModified && !chemisTree_.isFull());
345 template<
class CompType,
class ThermoType>
354 bool mechRedActive = this->chemistry_.mechRed()->
active();
355 label speciesNumber = this->chemistry_.nSpecie();
356 scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
357 for (
label i=0; i<speciesNumber; ++i)
362 s2c = this->chemistry_.simplifiedToCompleteIndex()[i];
364 Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
366 Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
367 Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
368 if (this->variableTimeStep())
370 Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
383 this->chemistry_.jacobian(runTime_.value(), Rcq,
A);
387 for (
label i=0; i<speciesNumber; ++i)
393 si = this->chemistry_.simplifiedToCompleteIndex()[i];
396 for (
label j=0; j<speciesNumber; ++j)
401 sj = this->chemistry_.simplifiedToCompleteIndex()[j];
404 -dt*this->chemistry_.specieThermo()[si].W()
405 /this->chemistry_.specieThermo()[sj].W();
410 A(i, speciesNumber) *=
411 -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
412 A(i, speciesNumber+1) *=
413 -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
417 A(speciesNumber, speciesNumber) = 1;
418 A(speciesNumber + 1, speciesNumber + 1) = 1;
419 if (this->variableTimeStep())
421 A[speciesNumber + 2][speciesNumber + 2] = 1;
425 LUscalarMatrix LUA(
A);
432 template<
class CompType,
class ThermoType>
439 bool retrieved(
false);
443 if (chemisTree_.size())
445 chemisTree_.binaryTreeSearch(phiq, chemisTree_.root(),
phi0);
450 if (
phi0->inEOA(phiq))
456 else if (chemisTree_.secondaryBTSearch(phiq,
phi0))
460 else if (MRURetrieve_)
465 if (
phi0->inEOA(phiq))
477 lastSearch_ =
nullptr;
482 phi0->increaseNumRetrieve();
483 scalar elapsedTimeSteps =
484 this->chemistry_.timeSteps() -
phi0->timeTag();
488 if (elapsedTimeSteps > chPMaxLifeTime_ && !
phi0->toRemove())
490 cleaningRequired_ =
true;
491 phi0->toRemove() =
true;
493 lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
495 calcNewC(
phi0,phiq, Rphiq);
507 template<
class CompType,
class ThermoType>
516 label growthOrAddFlag = 1;
520 if (lastSearch_ && growPoints_)
522 if (grow(lastSearch_,phiq, Rphiq))
528 return growthOrAddFlag;
535 if (chemisTree().isFull())
540 if (!cleanAndBalance())
555 chemisTree().clear();
563 for (
auto& t : tempList)
565 chemisTree().insertNewLeaf
581 lastSearch_ =
nullptr;
585 label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
587 computeA(
A, Rphiq,
rho, deltaT);
589 chemisTree().insertNewLeaf
602 return growthOrAddFlag;
606 template<
class CompType,
class ThermoType>
613 << runTime_.timeOutputValue() <<
" " << nRetrieved_ <<
endl;
617 << runTime_.timeOutputValue() <<
" " << nGrowth_ <<
endl;
621 << runTime_.timeOutputValue() <<
" " << nAdd_ <<
endl;
625 << runTime_.timeOutputValue() <<
" " << this->size() <<
endl;