zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

STYLE: Code clean-up

Browse Source
This commit is contained in:
Andrew Heather
2017-09-18 13:39:37 +01:00
parent f80334415c
commit af48c843f8
10 changed files with 389 additions and 359 deletions
Download Patch File Download Diff File Expand all files Collapse all files

39
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
View File

@ -405,13 +405,15 @@ void Foam::TDACChemistryModel<CompType, ThermoType>::jacobian
// but according to the informations of the complete set
// (i.e. for the third-body efficiencies)
const scalar T = c[this->nSpecie_];
const scalar p = c[this->nSpecie_ + 1];
const label nSpecie = this->nSpecie_;
const scalar T = c[nSpecie];
const scalar p = c[nSpecie + 1];
if (reduced)
{
this->c_ = completeC_;
for (label i=0; i<NsDAC_; i++)
for (label i=0; i<NsDAC_; ++i)
{
this->c_[simplifiedToCompleteIndex_[i]] = max(0.0, c[i]);
}
@ -557,20 +559,19 @@ void Foam::TDACChemistryModel<CompType, ThermoType>::jacobian
const scalar delta = 1e-3;
omega(this->c_, T + delta, p, this->dcdt_);
for (label i=0; i<this->nSpecie_; i++)
for (label i=0; i<nSpecie; ++i)
{
dfdc(i, this->nSpecie_) = this->dcdt_[i];
dfdc(i, nSpecie) = this->dcdt_[i];
}
omega(this->c_, T - delta, p, this->dcdt_);
for (label i=0; i<this->nSpecie_; i++)
for (label i=0; i<nSpecie; ++i)
{
dfdc(i, this->nSpecie_) =
0.5*(dfdc(i, this->nSpecie_) - this->dcdt_[i])/delta;
dfdc(i, nSpecie) = 0.5*(dfdc(i, nSpecie) - this->dcdt_[i])/delta;
}
dfdc(this->nSpecie_, this->nSpecie_) = 0;
dfdc(this->nSpecie_ + 1, this->nSpecie_) = 0;
dfdc(nSpecie, nSpecie) = 0;
dfdc(nSpecie + 1, nSpecie) = 0;
}
@ -666,12 +667,13 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
for (label i=0; i<this->nSpecie_; i++)
{
c[i] = rhoi*this->Y_[i][celli]/this->specieThermo_[i].W();
const volScalarField& Yi = this->Y_[i];
c[i] = rhoi*Yi[celli]/this->specieThermo_[i].W();
c0[i] = c[i];
phiq[i] = this->Y()[i][celli];
phiq[i] = Yi[celli];
}
phiq[this->nSpecie()]=Ti;
phiq[this->nSpecie() + 1]=pi;
phiq[this->nSpecie()] = Ti;
phiq[this->nSpecie() + 1] = pi;
if (tabulation_->variableTimeStep())
{
phiq[this->nSpecie() + 2] = deltaT[celli];
@ -692,7 +694,7 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
{
// Retrieved solution stored in Rphiq
for (label i=0; i<this->nSpecie(); i++)
for (label i=0; i<this->nSpecie(); ++i)
{
c[i] = rhoi*Rphiq[i]/this->specieThermo_[i].W();
}
@ -713,7 +715,7 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
// Reduce mechanism change the number of species (only active)
mechRed_->reduceMechanism(c, Ti, pi);
nActiveSpecies += mechRed_->NsSimp();
nAvg++;
++nAvg;
scalar timeIncr = clockTime_.timeIncrement();
reduceMechCpuTime_ += timeIncr;
timeTmp += timeIncr;
@ -735,7 +737,7 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
simplifiedC_, Ti, pi, dt, this->deltaTChem_[celli]
);
for (label i=0; i<NsDAC_; i++)
for (label i=0; i<NsDAC_; ++i)
{
c[simplifiedToCompleteIndex_[i]] = simplifiedC_[i];
}
@ -774,6 +776,7 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
}
label growOrAdd =
tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]);
if (growOrAdd)
{
this->setTabulationResultsAdd(celli);
@ -797,7 +800,7 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
}
// Set the RR vector (used in the solver)
for (label i=0; i<this->nSpecie_; i++)
for (label i=0; i<this->nSpecie_; ++i)
{
this->RR_[i][celli] =
(c[i] - c0[i])*this->specieThermo_[i].W()/deltaT[celli];

4
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H
View File

@ -168,9 +168,9 @@ Foam::TDACChemistryModel<CompType, ThermoType>::specieComp()
template<class CompType, class ThermoType>
void Foam::TDACChemistryModel<CompType, ThermoType>::resetTabulationResults()
{
forAll(tabulationResults_, tabi)
for (auto& res : tabulationResults_)
{
tabulationResults_[tabi] = 2;
res = 2;
}
}

66
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
View File

@ -253,7 +253,7 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
{
scalarField& completeC(this->chemistry_.completeC());
scalarField c1(this->chemistry_.nEqns(), 0.0);
for(label i=0; i<this->nSpecie_; i++)
for (label i=0; i<this->nSpecie_; i++)
{
c1[i] = c[i];
completeC[i] = c[i];
@ -276,9 +276,8 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
scalar pf, cf, pr, cr;
label lRef, rRef;
forAll(this->chemistry_.reactions(), i)
for (const Reaction<ThermoType>& R : this->chemistry_.reactions())
{
const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
// for each reaction compute omegai
scalar omegai = this->chemistry_.omega
(
@ -511,7 +510,7 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
// Set all species to inactive and activate them according
// to rAB and initial set
for (label i=0; i<this->nSpecie_; i++)
for (label i=0; i<this->nSpecie_; ++i)
{
this->activeSpecies_[i] = false;
}
@ -530,19 +529,19 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
// When phiLarge and phiProgress >= phiTol then
// CO, HO2 and fuel are in the SIS
Q.push(COId_);
speciesNumber++;
++speciesNumber;
this->activeSpecies_[COId_] = true;
Rvalue[COId_] = 1.0;
Q.push(HO2Id_);
speciesNumber++;
++speciesNumber;
this->activeSpecies_[HO2Id_] = true;
Rvalue[HO2Id_] = 1.0;
forAll(fuelSpeciesID_,i)
for (const label fuelId : fuelSpeciesID_)
{
Q.push(fuelSpeciesID_[i]);
speciesNumber++;
this->activeSpecies_[fuelSpeciesID_[i]] = true;
Rvalue[fuelSpeciesID_[i]] = 1.0;
Q.push(fuelId);
++speciesNumber;
this->activeSpecies_[fuelId] = true;
Rvalue[fuelId] = 1.0;
}
}
@ -551,22 +550,22 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
// When phiLarge < phiTol and phiProgress >= phiTol then
// CO, HO2 are in the SIS
Q.push(COId_);
speciesNumber++;
++speciesNumber;
this->activeSpecies_[COId_] = true;
Rvalue[COId_] = 1.0;
Q.push(HO2Id_);
speciesNumber++;
++speciesNumber;
this->activeSpecies_[HO2Id_] = true;
Rvalue[HO2Id_] = 1.0;
if (forceFuelInclusion_)
{
forAll(fuelSpeciesID_, i)
for (const label fuelId : fuelSpeciesID_)
{
Q.push(fuelSpeciesID_[i]);
speciesNumber++;
this->activeSpecies_[fuelSpeciesID_[i]] = true;
Rvalue[fuelSpeciesID_[i]] = 1.0;
Q.push(fuelId);
++speciesNumber;
this->activeSpecies_[fuelId] = true;
Rvalue[fuelId] = 1.0;
}
}
}
@ -585,12 +584,12 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
Rvalue[H2OId_] = 1.0;
if (forceFuelInclusion_)
{
forAll(fuelSpeciesID_, i)
for (const label fuelId : fuelSpeciesID_)
{
Q.push(fuelSpeciesID_[i]);
speciesNumber++;
this->activeSpecies_[fuelSpeciesID_[i]] = true;
Rvalue[fuelSpeciesID_[i]] = 1.0;
Q.push(fuelId);
++speciesNumber;
this->activeSpecies_[fuelId] = true;
Rvalue[fuelId] = 1.0;
}
}
}
@ -598,7 +597,7 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
if (T > NOxThreshold_ && NOId_ != -1)
{
Q.push(NOId_);
speciesNumber++;
++speciesNumber;
this->activeSpecies_[NOId_] = true;
Rvalue[NOId_] = 1.0;
}
@ -609,7 +608,7 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
{
label q = SIS[i];
this->activeSpecies_[q] = true;
speciesNumber++;
++speciesNumber;
Q.push(q);
Rvalue[q] = 1.0;
}
@ -620,18 +619,23 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
{
label u = Q.pop();
scalar Den = max(PA[u], CA[u]);
if (Den!=0.0)
if (Den != 0.0)
{
for (label v=0; v<NbrABInit[u]; v++)
for (label v=0; v<NbrABInit[u]; ++v)
{
label otherSpec = rABOtherSpec(u, v);
scalar rAB = mag(rABNum(u, v))/Den;
if (rAB>1)
if (rAB > 1)
{
Info<< "Badly Conditioned rAB : " << rAB
<< "species involved : "<<u << "," << otherSpec << endl;
rAB=1.0;
<< " for species : "
<< this->chemistry_.Y()[u].name() << ","
<< this->chemistry_.Y()[otherSpec].name()
<< endl;
rAB = 1.0;
}
// The direct link is weaker than the user-defined tolerance
if (rAB >= this->tolerance())
{
@ -646,7 +650,7 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
if (!this->activeSpecies_[otherSpec])
{
this->activeSpecies_[otherSpec] = true;
speciesNumber++;
++speciesNumber;
}
}
}

14
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
View File

@ -105,8 +105,8 @@ void Foam::chemistryReductionMethods::DRG<CompType, ThermoType>::reduceMechanism
// For each reaction compute omegai
scalar omegai = this->chemistry_.omega
(
R, c1, T, p, pf, cf, lRef, pr, cr, rRef
);
R, c1, T, p, pf, cf, lRef, pr, cr, rRef
);
// Then for each pair of species composing this reaction,
@ -253,12 +253,14 @@ void Foam::chemistryReductionMethods::DRG<CompType, ThermoType>::reduceMechanism
label otherSpec = rABOtherSpec(u, v);
scalar rAB = rABNum(u, v)/Den;
if (rAB>1)
if (rAB > 1)
{
Info<< "Badly Conditioned rAB : " << rAB
<< "species involved : " << u << "," << otherSpec
<< endl;
rAB = 1;
<< " for species : "
<< this->chemistry_.Y()[u].name() << ","
<< this->chemistry_.Y()[otherSpec].name()
<< endl;
rAB = 1.0;
}
// Include B only if rAB is above the tolerance and if the

142
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
View File

@ -37,25 +37,25 @@ Foam::chemistryReductionMethods::DRGEP<CompType, ThermoType>::DRGEP
:
chemistryReductionMethod<CompType, ThermoType>(dict, chemistry),
searchInitSet_(this->coeffsDict_.subDict("initialSet").size()),
sC_(this->nSpecie_,0),
sH_(this->nSpecie_,0),
sO_(this->nSpecie_,0),
sN_(this->nSpecie_,0),
sC_(this->nSpecie_, 0),
sH_(this->nSpecie_, 0),
sO_(this->nSpecie_, 0),
sN_(this->nSpecie_, 0),
NGroupBased_(50)
{
label j=0;
label j = 0;
dictionary initSet = this->coeffsDict_.subDict("initialSet");
for (label i=0; i<chemistry.nSpecie(); i++)
for (label i=0; i<chemistry.nSpecie(); ++i)
{
if (initSet.found(chemistry.Y()[i].name()))
{
searchInitSet_[j++]=i;
searchInitSet_[j++] = i;
}
}
if (j<searchInitSet_.size())
{
FatalErrorInFunction
<< searchInitSet_.size()-j
<< searchInitSet_.size() - j
<< " species in the intial set is not in the mechanism "
<< initSet
<< exit(FatalError);
@ -68,15 +68,14 @@ Foam::chemistryReductionMethods::DRGEP<CompType, ThermoType>::DRGEP
const List<List<specieElement>>& specieComposition =
this->chemistry_.specieComp();
for (label i=0; i<this->nSpecie_; i++)
for (label i=0; i<this->nSpecie_; ++i)
{
const List<specieElement>& curSpecieComposition =
specieComposition[i];
const List<specieElement>& curSpecieComposition = specieComposition[i];
// for all elements in the current species
forAll(curSpecieComposition, j)
for (const specieElement& curElement : curSpecieComposition)
{
const specieElement& curElement =
curSpecieComposition[j];
if (curElement.name() == "C")
{
sC_[i] = curElement.nAtoms();
@ -123,7 +122,7 @@ reduceMechanism
scalarField& completeC(this->chemistry_.completeC());
scalarField c1(this->chemistry_.nEqns(), 0.0);
for (label i=0; i<this->nSpecie_; i++)
for (label i=0; i<this->nSpecie_; ++i)
{
c1[i] = c[i];
completeC[i] = c[i];
@ -133,12 +132,12 @@ reduceMechanism
c1[this->nSpecie_+1] = p;
// Compute the rAB matrix
RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
scalarField PA(this->nSpecie_,0.0);
scalarField CA(this->nSpecie_,0.0);
RectangularMatrix<scalar> rABNum(this->nSpecie_, this->nSpecie_, 0.0);
scalarField PA(this->nSpecie_, 0.0);
scalarField CA(this->nSpecie_, 0.0);
// Number of initialized rAB for each lines
Field<label> NbrABInit(this->nSpecie_,0);
Field<label> NbrABInit(this->nSpecie_, 0);
// Position of the initialized rAB, -1 when not initialized
RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
// Index of the other species involved in the rABNum
@ -150,11 +149,12 @@ reduceMechanism
forAll(this->chemistry_.reactions(), i)
{
const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
// for each reaction compute omegai
scalar omegai = this->chemistry_.omega
(
R, c1, T, p, pf, cf, lRef, pr, cr, rRef
);
R, c1, T, p, pf, cf, lRef, pr, cr, rRef
);
omegaV[i] = omegai;
// then for each pair of species composing this reaction,
@ -185,7 +185,7 @@ reduceMechanism
// Disable for self reference (by definition rAA=0)
deltaBi[ss] = false;
while(!usedIndex.empty())
while (!usedIndex.empty())
{
label curIndex = usedIndex.pop();
if (deltaBi[curIndex])
@ -193,7 +193,7 @@ reduceMechanism
// disable to avoid counting it more than once
deltaBi[curIndex] = false;
// test if this rAB is not initialized
if (rABPos(ss, curIndex)==-1)
if (rABPos(ss, curIndex) == -1)
{
rABPos(ss, curIndex) = NbrABInit[ss];
NbrABInit[ss]++;
@ -209,7 +209,7 @@ reduceMechanism
bool found(false);
forAll(wAID, id)
{
if (ss==wAID[id])
if (ss == wAID[id])
{
wA[id] += sl*omegai;
found = true;
@ -245,7 +245,7 @@ reduceMechanism
// Disable for self reference (by definition rAA=0)
deltaBi[ss] = false;
while(!usedIndex.empty())
while (!usedIndex.empty())
{
label curIndex = usedIndex.pop();
if (deltaBi[curIndex])
@ -253,7 +253,7 @@ reduceMechanism
// disable to avoid counting it more than once
deltaBi[curIndex] = false;
// test if this rAB is not initialized
if (rABPos(ss, curIndex)==-1)
if (rABPos(ss, curIndex) == -1)
{
rABPos(ss, curIndex) = NbrABInit[ss];
NbrABInit[ss]++;
@ -270,7 +270,7 @@ reduceMechanism
bool found(false);
forAll(wAID, id)
{
if (ss==wAID[id])
if (ss == wAID[id])
{
wA[id] += sl*omegai;
found = true;
@ -320,28 +320,28 @@ reduceMechanism
scalarList Pa(nElements,0.0);
scalarList Ca(nElements,0.0);
// for (label q=0; q<SIS.size(); q++)
for (label i=0; i<this->nSpecie_; i++)
// for (label q=0; q<SIS.size(); ++q)
for (label i=0; i<this->nSpecie_; ++i)
{
Pa[0] += sC_[i]*max(0.0,(PA[i]-CA[i]));
Pa[0] += sC_[i]*max(0.0, (PA[i]-CA[i]));
Ca[0] += sC_[i]*max(0.0,-(PA[i]-CA[i]));
Pa[1] += sH_[i]*max(0.0,(PA[i]-CA[i]));
Pa[1] += sH_[i]*max(0.0, (PA[i]-CA[i]));
Ca[1] += sH_[i]*max(0.0,-(PA[i]-CA[i]));
Pa[2] += sO_[i]*max(0.0,(PA[i]-CA[i]));
Pa[2] += sO_[i]*max(0.0, (PA[i]-CA[i]));
Ca[2] += sO_[i]*max(0.0,-(PA[i]-CA[i]));
Pa[3] += sN_[i]*max(0.0,(PA[i]-CA[i]));
Pa[3] += sN_[i]*max(0.0, (PA[i]-CA[i]));
Ca[3] += sN_[i]*max(0.0,-(PA[i]-CA[i]));
}
// Using the rAB matrix (numerator and denominator separated)
// compute the R value according to the search initiating set
scalarField Rvalue(this->nSpecie_,0.0);
scalarField Rvalue(this->nSpecie_, 0.0);
label speciesNumber = 0;
List<bool> disabledSpecies(this->nSpecie_,false);
List<bool> disabledSpecies(this->nSpecie_, false);
// set all species to inactive and activate them according
// to rAB and initial set
for (label i=0; i<this->nSpecie_; i++)
for (label i=0; i<this->nSpecie_; ++i)
{
this->activeSpecies_[i] = false;
}
@ -353,7 +353,7 @@ reduceMechanism
// Compute the alpha coefficient and initialize the R value of the species
// in the SIS
for (label i=0; i<SIS.size(); i++)
for (label i=0; i<SIS.size(); ++i)
{
label q = SIS[i];
// compute alpha
@ -397,7 +397,7 @@ reduceMechanism
if (alphaA > this->tolerance())
{
this->activeSpecies_[q] = true;
speciesNumber++;
++speciesNumber;
Q.push(q);
QStart.append(q);
alphaQ.append(1.0);
@ -415,18 +415,19 @@ reduceMechanism
{
scalar Rmax=0.0;
label specID=-1;
forAll(SIS, i)
for (const label sis : SIS)
{
if (Rvalue[SIS[i]] > Rmax)
if (Rvalue[sis] > Rmax)
{
Rmax = Rvalue[SIS[i]];
specID=SIS[i];
Rmax = Rvalue[sis];
specID = sis;
}
}
Q.push(specID);
QStart.append(specID);
alphaQ.append(1.0);
speciesNumber++;
++speciesNumber;
Rvalue[specID] = 1.0;
this->activeSpecies_[specID] = true;
}
@ -435,18 +436,22 @@ reduceMechanism
while (!Q.empty())
{
label u = Q.pop();
scalar Den = max(PA[u],CA[u]);
scalar Den = max(PA[u], CA[u]);
if (Den > VSMALL)
{
for (label v=0; v<NbrABInit[u]; v++)
for (label v=0; v<NbrABInit[u]; ++v)
{
label otherSpec = rABOtherSpec(u, v);
scalar rAB = mag(rABNum(u, v))/Den;
if (rAB>1)
if (rAB > 1)
{
Info<< "Badly Conditioned rAB : " << rAB
<< "species involved : "<<u << "," << otherSpec << endl;
rAB=1.0;
<< " for species : "
<< this->chemistry_.Y()[u].name() << ","
<< this->chemistry_.Y()[otherSpec].name()
<< endl;
rAB = 1.0;
}
scalar Rtemp = Rvalue[u]*rAB;
@ -461,7 +466,7 @@ reduceMechanism
if (!this->activeSpecies_[otherSpec])
{
this->activeSpecies_[otherSpec] = true;
speciesNumber++;
++speciesNumber;
}
}
}
@ -471,14 +476,14 @@ reduceMechanism
// Group-based reduction
// number of species disabled in the first step
label NDisabledSpecies(this->nSpecie_-speciesNumber);
label NDisabledSpecies(this->nSpecie_ - speciesNumber);
// while the number of removed species is greater than NGroupBased, the rAB
// are reevaluated according to the group based definition for each loop the
// temporary disabled species (in the first reduction) are sorted to disable
// definitely the NGroupBased species with lower R then these R value a
// reevaluated taking into account these disabled species
while(NDisabledSpecies > NGroupBased_)
while (NDisabledSpecies > NGroupBased_)
{
// keep track of disabled species using sortablelist to extract only
// NGroupBased lower R value
@ -500,7 +505,7 @@ reduceMechanism
labelList tmpIndex(Rdisabled.indices());
// disable definitely NGroupBased species in this loop
for (label i=0; i<NGroupBased_; i++)
for (label i=0; i<NGroupBased_; ++i)
{
disabledSpecies[Rindex[tmpIndex[i]]] = true;
}
@ -510,7 +515,7 @@ reduceMechanism
// only update the numerator
forAll(NbrABInit, i)
{
for (label v=0; v<NbrABInit[i]; v++)
for (label v=0; v<NbrABInit[i]; ++v)
{
rABNum(i, v) = 0.0;
}
@ -555,7 +560,7 @@ reduceMechanism
{
// if one of the species in this reaction is disabled, all
// species connected to species ss are modified
for (label v=0; v<NbrABInit[ss]; v++)
for (label v=0; v<NbrABInit[ss]; ++v)
{
rABNum(ss, v) += sl*omegai;
}
@ -589,7 +594,7 @@ reduceMechanism
deltaBi[sj] = true;
if (disabledSpecies[sj])
{
alreadyDisabled=true;
alreadyDisabled = true;
}
}
forAll(R.rhs(), j)
@ -599,7 +604,7 @@ reduceMechanism
deltaBi[sj] = true;
if (disabledSpecies[sj])
{
alreadyDisabled=true;
alreadyDisabled = true;
}
}
@ -609,7 +614,7 @@ reduceMechanism
{
// if one of the species in this reaction is disabled, all
// species connected to species ss are modified
for (label v=0; v<NbrABInit[ss]; v++)
for (label v=0; v<NbrABInit[ss]; ++v)
{
rABNum(ss, v) += sl*omegai;
}
@ -629,32 +634,31 @@ reduceMechanism
}
}
forAll(QStart, qi)
for (const label q : QStart)
{
label u = QStart[qi];
Q.push(u);
Q.push(q);
}
while (!Q.empty())
{
label u = Q.pop();
scalar Den = max(PA[u],CA[u]);
if (Den!=0.0)
if (Den != 0.0)
{
for (label v=0; v<NbrABInit[u]; v++)
for (label v=0; v<NbrABInit[u]; ++v)
{
label otherSpec = rABOtherSpec(u, v);
if (!disabledSpecies[otherSpec])
{
scalar rAB = mag(rABNum(u, v))/Den;
if (rAB>1.0)
if (rAB > 1)
{
Info<< "Badly Conditioned rAB : " << rAB
<< "species involved : "
<<this->chemistry_.Y()[u].name() << ","
<< " for species : "
<< this->chemistry_.Y()[u].name() << ","
<< this->chemistry_.Y()[otherSpec].name()
<< endl;
rAB=1.0;
rAB = 1.0;
}
scalar Rtemp = Rvalue[u]*rAB;
@ -668,7 +672,7 @@ reduceMechanism
if (!this->activeSpecies_[otherSpec])
{
this->activeSpecies_[otherSpec] = true;
speciesNumber++;
++speciesNumber;
NDisabledSpecies--;
}
}
@ -717,13 +721,13 @@ reduceMechanism
Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
label j = 0;
for (label i=0; i<this->nSpecie_; i++)
for (label i=0; i<this->nSpecie_; ++i)
{
if (this->activeSpecies_[i])
{
s2c[j] = i;
simplifiedC[j] = c[i];
c2s[i] = j++;
c2s[i] = ++j;
if (!this->chemistry_.active(i))
{
this->chemistry_.setActive(i);

93
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
View File

@ -82,19 +82,16 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
dictionary scaleDict(this->coeffsDict_.subDict("scaleFactor"));
label Ysize = this->chemistry_.Y().size();
scalar otherScaleFactor = readScalar(scaleDict.lookup("otherSpecies"));
for (label i=0; i<Ysize; i++)
for (label i=0; i<Ysize; ++i)
{
if (!scaleDict.found(this->chemistry_.Y()[i].name()))
const word& yName = this->chemistry_.Y()[i].name();
if (!scaleDict.found(yName))
{
scaleFactor_[i] = otherScaleFactor;
}
else
{
scaleFactor_[i] =
readScalar
(
scaleDict.lookup(this->chemistry_.Y()[i].name())
);
scaleFactor_[i] = readScalar(scaleDict.lookup(yName));
}
}
scaleFactor_[Ysize] = readScalar(scaleDict.lookup("Temperature"));
@ -141,10 +138,11 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::addToMRU
{
if (maxMRUSize_ > 0 && MRURetrieve_)
{
typename SLList<chemPointISAT<CompType, ThermoType>*>::iterator iter =
MRUList_.begin();
// First search if the chemPoint is already in the list
bool isInList = false;
typename SLList <chemPointISAT<CompType, ThermoType>*>::iterator iter =
MRUList_.begin();
for ( ; iter != MRUList_.end(); ++iter)
{
if (iter() == phi0)
@ -153,10 +151,11 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::addToMRU
break;
}
}
// If it is in the list, then move it to front
if (isInList)
{
if (iter() != MRUList_.first())
// If it is in the list, then move it to front
if (iter != MRUList_.begin())
{
// iter hold the position of the element to move
MRUList_.remove(iter);
@ -165,11 +164,12 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::addToMRU
MRUList_.insert(phi0);
}
}
else // chemPoint not yet in the list, iter is last
else
{
// chemPoint not yet in the list, iter is last
if (MRUList_.size() == maxMRUSize_)
{
if (iter() == MRUList_.last())
if (iter == MRUList_.end())
{
MRUList_.remove(iter);
MRUList_.insert(phi0);
@ -177,7 +177,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::addToMRU
else
{
FatalErrorInFunction
<< "wrong MRUList construction"
<< "Error in MRUList construction"
<< exit(FatalError);
}
}
@ -207,7 +207,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
// Rphiq[i]=Rphi0[i]+A(i, j)dphi[j]
// where Aij is dRi/dphi_j
for (label i=0; i<nEqns-nAdditionalEqns_; i++)
for (label i=0; i<nEqns-nAdditionalEqns_; ++i)
{
if (mechRedActive)
{
@ -238,24 +238,24 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
// As we use an approximation of A, Rphiq should be checked for
// negative values
Rphiq[i] = max(0.0,Rphiq[i]);
Rphiq[i] = max(0.0, Rphiq[i]);
}
// The species is not active A(i, j) = I(i, j)
else
{
Rphiq[i] += dphi[i];
Rphiq[i] = max(0.0,Rphiq[i]);
Rphiq[i] = max(0.0, Rphiq[i]);
}
}
else // Mechanism reduction is not active
{
for (label j=0; j<nEqns; j++)
for (label j=0; j<nEqns; ++j)
{
Rphiq[i] += gradientsMatrix(i, j)*dphi[j];
}
// As we use a first order gradient matrix, Rphiq should be checked
// for negative values
Rphiq[i] = max(0.0,Rphiq[i]);
Rphiq[i] = max(0.0, Rphiq[i]);
}
}
}
@ -270,7 +270,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::grow
)
{
// If the pointer to the chemPoint is nullptr, the function stops
if (!phi0)
if (phi0 == nullptr)
{
return false;
}
@ -306,9 +306,9 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
bool treeModified(false);
// Check all chemPoints to see if we need to delete some of the chemPoints
// according to the ellapsed time and number of growths
// according to the elapsed time and number of growths
chemPointISAT<CompType, ThermoType>* x = chemisTree_.treeMin();
while(x != nullptr)
while (x != nullptr)
{
chemPointISAT<CompType, ThermoType>* xtmp =
chemisTree_.treeSuccessor(x);
@ -323,7 +323,7 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
x = xtmp;
}
// Check if the tree should be balanced according to criterion:
// -the depth of the tree bigger than a*log2(size), log2(size) being the
// - the depth of the tree bigger than a*log2(size), log2(size) being the
// ideal depth (e.g. 4 leafs can be stored in a tree of depth 2)
if
(
@ -355,7 +355,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
bool mechRedActive = this->chemistry_.mechRed()->active();
label speciesNumber = this->chemistry_.nSpecie();
scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
for (label i=0; i<speciesNumber; i++)
for (label i=0; i<speciesNumber; ++i)
{
label s2c = i;
if (mechRedActive)
@ -371,10 +371,10 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
}
// Aaa is computed implicitely,
// Aaa is computed implicitly,
// A is given by A = C(psi0, t0+dt), where C is obtained through solving
// d/dt C(psi0,t) = J(psi(t))C(psi0,t)
// If we solve it implicitely:
// If we solve it implicitly:
// (C(psi0, t0+dt) - C(psi0,t0))/dt = J(psi(t0+dt))C(psi0,t0+dt)
// The Jacobian is thus computed according to the mapping
// C(psi0,t0+dt)*(I-dt*J(psi(t0+dt))) = C(psi0, t0)
@ -385,7 +385,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
// The jacobian is computed according to the molar concentration
// the following conversion allows the code to use A with mass fraction
for (label i=0; i<speciesNumber; i++)
for (label i=0; i<speciesNumber; ++i)
{
label si = i;
@ -394,7 +394,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
si = this->chemistry_.simplifiedToCompleteIndex()[i];
}
for (label j=0; j<speciesNumber; j++)
for (label j=0; j<speciesNumber; ++j)
{
label sj = j;
if (mechRedActive)
@ -460,12 +460,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
}
else if (MRURetrieve_)
{
typename SLList
<
chemPointISAT<CompType, ThermoType>*
>::iterator iter = MRUList_.begin();
for ( ; iter != MRUList_.end(); ++iter)
forAllConstIters(MRUList_, iter)
{
phi0 = iter();
if (phi0->inEOA(phiq))
@ -499,7 +494,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
addToMRU(phi0);
calcNewC(phi0,phiq, Rphiq);
nRetrieved_++;
++nRetrieved_;
return true;
}
else
@ -521,15 +516,17 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
)
{
label growthOrAddFlag = 1;
// If lastSearch_ holds a valid pointer to a chemPoint AND the growPoints_
// option is on, the code first tries to grow the point hold by lastSearch_
if (lastSearch_ && growPoints_)
{
if (grow(lastSearch_,phiq, Rphiq))
{
nGrowth_++;
++nGrowth_;
growthOrAddFlag = 0;
//the structure of the tree is not modified, return false
// The structure of the tree is not modified, return false
return growthOrAddFlag;
}
}
@ -541,7 +538,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
{
// If cleanAndBalance operation do not result in a reduction of the tree
// size, the last possibility is to delete completely the tree.
// It can be partially rebuild with the MRU list if this is used.
// It can be partially rebuilt with the MRU list if this is used.
if (!cleanAndBalance())
{
DynamicList<chemPointISAT<CompType, ThermoType>*> tempList;
@ -549,11 +546,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
{
// Create a copy of each chemPointISAT of the MRUList_ before
// they are deleted
typename SLList
<
chemPointISAT<CompType, ThermoType>*
>::iterator iter = MRUList_.begin();
for ( ; iter != MRUList_.end(); ++iter)
forAllConstIters(MRUList_, iter)
{
tempList.append
(
@ -568,20 +561,20 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
// Construct the tree without giving a reference to attach to it
// since the structure has been completely discarded
chemPointISAT<CompType, ThermoType>* nulPhi = 0;
forAll(tempList, i)
chemPointISAT<CompType, ThermoType>* nulPhi = nullptr;
for (auto& t : tempList)
{
chemisTree().insertNewLeaf
(
tempList[i]->phi(),
tempList[i]->Rphi(),
tempList[i]->A(),
t->phi(),
t->Rphi(),
t->A(),
scaleFactor(),
this->tolerance(),
scaleFactor_.size(),
nulPhi
);
deleteDemandDrivenData(tempList[i]);
deleteDemandDrivenData(t);
}
}
@ -606,7 +599,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
lastSearch_ // lastSearch_ may be nullptr (handled by binaryTree)
);
nAdd_++;
++nAdd_;
return growthOrAddFlag;
}

190
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryTree/binaryTree.C
View File

@ -31,18 +31,20 @@ License
template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::insertNode
(
chP*& phi0,
bn*& newNode
chemPoint*& phi0,
node*& newNode
)
{
if (phi0 == phi0->node()->leafRight())// phi0 is on the right
if (phi0 == phi0->node()->leafRight())
{
// phi0 is on the right
phi0->node()->leafRight() = nullptr;
phi0->node()->nodeRight() = newNode;
return;
}
else if (phi0 == phi0->node()->leafLeft())// phi0 is on the left
else if (phi0 == phi0->node()->leafLeft())
{
// phi0 is on the left
phi0->node()->leafLeft() = nullptr;
phi0->node()->nodeLeft() = newNode;
return;
@ -60,34 +62,38 @@ template<class CompType, class ThermoType>
bool Foam::binaryTree<CompType, ThermoType>::inSubTree
(
const scalarField& phiq,
bn* y,
chP* x
node* y,
chemPoint* x
)
{
if ((n2ndSearch_ < max2ndSearch_) && (y!=nullptr))
if ((n2ndSearch_ < max2ndSearch_) && (y != nullptr))
{
scalar vPhi = 0;
const scalarField& v = y->v();
const scalar a = y->a();
// compute v*phi
for (label i=0; i<phiq.size(); i++)
for (label i=0; i<phiq.size(); ++i)
{
vPhi += phiq[i]*v[i];
}
if (vPhi <= a)// on the left side of the node
if (vPhi <= a)
{
if (y->nodeLeft() == nullptr)// left is a chemPoint
// on the left side of the node
if (y->nodeLeft() == nullptr)
{
n2ndSearch_++;
// left is a chemPoint
++n2ndSearch_;
if (y->leafLeft()->inEOA(phiq))
{
x = y->leafLeft();
return true;
}
}
else // the left side is a node
else
{
if (inSubTree(phiq, y->nodeLeft(),x))
// the left side is a node
if (inSubTree(phiq, y->nodeLeft(), x))
{
return true;
}
@ -96,7 +102,7 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
// not on the left side, try the right side
if ((n2ndSearch_ < max2ndSearch_) && y->nodeRight() == nullptr)
{
n2ndSearch_++;
++n2ndSearch_;
// we reach the end of the subTree we can return the result
if (y->leafRight()->inEOA(phiq))
{
@ -111,14 +117,16 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
}
else // test for n2ndSearch is done in the call of inSubTree
{
return inSubTree(phiq, y->nodeRight(),x);
return inSubTree(phiq, y->nodeRight(), x);
}
}
else // on right side (symetric of above)
else
{
// on right side (symetric of above)
if (y->nodeRight() == nullptr)
{
n2ndSearch_++;
++n2ndSearch_;
if (y->leafRight()->inEOA(phiq))
{
return true;
@ -126,17 +134,18 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
}
else // the right side is a node
{
if (inSubTree(phiq, y->nodeRight(),x))
if (inSubTree(phiq, y->nodeRight(), x))
{
x = y->leafRight();
return true;
}
}
// if we reach this point, the retrieve has
// failed on the right side, explore the left side
if ((n2ndSearch_ < max2ndSearch_) && y->nodeLeft() == nullptr)
{
n2ndSearch_++;
++n2ndSearch_;
if (y->leafLeft()->inEOA(phiq))
{
x = y->leafLeft();
@ -150,7 +159,7 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
}
else
{
return inSubTree(phiq, y->nodeLeft(),x);
return inSubTree(phiq, y->nodeLeft(), x);
}
}
}
@ -162,7 +171,7 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::deleteSubTree(bn* subTreeRoot)
void Foam::binaryTree<CompType, ThermoType>::deleteSubTree(node* subTreeRoot)
{
if (subTreeRoot != nullptr)
{
@ -176,7 +185,7 @@ void Foam::binaryTree<CompType, ThermoType>::deleteSubTree(bn* subTreeRoot)
template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::transplant(bn* u, bn* v)
void Foam::binaryTree<CompType, ThermoType>::transplant(node* u, node* v)
{
if (v != nullptr)
{
@ -214,7 +223,7 @@ void Foam::binaryTree<CompType, ThermoType>::transplant(bn* u, bn* v)
template<class CompType, class ThermoType>
Foam::chemPointISAT<CompType, ThermoType>*
Foam::binaryTree<CompType, ThermoType>::chemPSibling(bn* y)
Foam::binaryTree<CompType, ThermoType>::chemPSibling(node* y)
{
if (y->parent() != nullptr)
{
@ -243,9 +252,9 @@ Foam::binaryTree<CompType, ThermoType>::chemPSibling(bn* y)
template<class CompType, class ThermoType>
Foam::chemPointISAT<CompType, ThermoType>*
Foam::binaryTree<CompType, ThermoType>::chemPSibling(chP* x)
Foam::binaryTree<CompType, ThermoType>::chemPSibling(chemPoint* x)
{
if (size_>1)
if (size_ > 1)
{
if (x == x->node()->leafLeft())
{
@ -273,9 +282,9 @@ Foam::binaryTree<CompType, ThermoType>::chemPSibling(chP* x)
template<class CompType, class ThermoType>
Foam::binaryNode<CompType, ThermoType>*
Foam::binaryTree<CompType, ThermoType>::nodeSibling(bn* y)
Foam::binaryTree<CompType, ThermoType>::nodeSibling(node* y)
{
if (y->parent()!=nullptr)
if (y->parent() != nullptr)
{
if (y == y->parent()->nodeLeft())
{
@ -300,9 +309,9 @@ Foam::binaryTree<CompType, ThermoType>::nodeSibling(bn* y)
template<class CompType, class ThermoType>
Foam::binaryNode<CompType, ThermoType>*
Foam::binaryTree<CompType, ThermoType>::nodeSibling(chP* x)
Foam::binaryTree<CompType, ThermoType>::nodeSibling(chemPoint* x)
{
if (size_>1)
if (size_ > 1)
{
if (x == x->node()->leafLeft())
{
@ -327,7 +336,7 @@ Foam::binaryTree<CompType, ThermoType>::nodeSibling(chP* x)
template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::deleteAllNode(bn* subTreeRoot)
void Foam::binaryTree<CompType, ThermoType>::deleteAllNode(node* subTreeRoot)
{
if (subTreeRoot != nullptr)
{
@ -359,7 +368,7 @@ Foam::binaryTree<CompType, ThermoType>::binaryTree
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CompType, class ThermoType>
Foam::label Foam::binaryTree<CompType, ThermoType>::depth(bn* subTreeRoot)
Foam::label Foam::binaryTree<CompType, ThermoType>::depth(node* subTreeRoot)
{
// when we reach the leaf, we return 0
if (subTreeRoot == nullptr)
@ -388,17 +397,17 @@ void Foam::binaryTree<CompType, ThermoType>::insertNewLeaf
const scalarField& scaleFactor,
const scalar& epsTol,
const label nCols,
chP*& phi0
chemPoint*& phi0
)
{
if (size_ == 0) // no points are stored
{
// create an empty binary node and point root_ to it
root_ = new bn();
root_ = new node();
// create the new chemPoint which holds the composition point
// phiq and the data to initialize the EOA
chP* newChemPoint =
new chP
chemPoint* newChemPoint =
new chemPoint
(
chemistry_,
phiq,
@ -410,22 +419,22 @@ void Foam::binaryTree<CompType, ThermoType>::insertNewLeaf
coeffsDict_,
root_
);
root_->leafLeft()=newChemPoint;
root_->leafLeft() = newChemPoint;
}
else // at least one point stored
{
// no reference chemPoint, a BT search is required
if (phi0 == nullptr)
{
binaryTreeSearch(phiq, root_,phi0);
binaryTreeSearch(phiq, root_, phi0);
}
// access to the parent node of the chemPoint
bn* parentNode = phi0->node();
node* parentNode = phi0->node();
// create the new chemPoint which holds the composition point
// phiq and the data to initialize the EOA
chP* newChemPoint =
new chP
chemPoint* newChemPoint =
new chemPoint
(
chemistry_,
phiq,
@ -440,10 +449,10 @@ void Foam::binaryTree<CompType, ThermoType>::insertNewLeaf
// previously stored leaf (phi0)
// the new node contains phi0 on the left and phiq on the right
// the hyper plane is computed in the binaryNode constructor
bn* newNode;
if (size_>1)
node* newNode;
if (size_ > 1)
{
newNode = new bn(phi0, newChemPoint, parentNode);
newNode = new node(phi0, newChemPoint, parentNode);
// make the parent of phi0 point to the newly created node
insertNode(phi0, newNode);
}
@ -451,14 +460,15 @@ void Foam::binaryTree<CompType, ThermoType>::insertNewLeaf
{
// when size is 1, the binaryNode is without hyperplane
deleteDemandDrivenData(root_);
newNode = new bn(phi0, newChemPoint, nullptr);
newNode = new node(phi0, newChemPoint, nullptr);
root_ = newNode;
}
phi0->node() = newNode;
newChemPoint->node()=newNode;
}
size_++;
++size_;
}
@ -466,8 +476,8 @@ template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::binaryTreeSearch
(
const scalarField& phiq,
bn* node,
chP*& nearest
node* node,
chemPoint*& nearest
)
{
if (size_ > 1)
@ -476,12 +486,15 @@ void Foam::binaryTree<CompType, ThermoType>::binaryTreeSearch
const scalarField& v = node->v();
const scalar& a = node->a();
// compute v*phi
for (label i=0; i<phiq.size(); i++) vPhi += phiq[i]*v[i];
for (label i=0; i<phiq.size(); ++i)
{
vPhi += phiq[i]*v[i];
}
if (vPhi > a) // on right side (side of the newly added point)
{
if (node->nodeRight()!=nullptr)
if (node->nodeRight() != nullptr)
{
binaryTreeSearch(phiq, node->nodeRight(), nearest);
}
@ -493,7 +506,7 @@ void Foam::binaryTree<CompType, ThermoType>::binaryTreeSearch
}
else // on left side (side of the previously stored point)
{
if (node->nodeLeft()!=nullptr)
if (node->nodeLeft() != nullptr)
{
binaryTreeSearch(phiq, node->nodeLeft(), nearest);
}
@ -520,14 +533,14 @@ template<class CompType, class ThermoType>
bool Foam::binaryTree<CompType, ThermoType>::secondaryBTSearch
(
const scalarField& phiq,
chP*& x
chemPoint*& x
)
{
// initialize n2ndSearch_
n2ndSearch_ = 0;
if ((n2ndSearch_ < max2ndSearch_) && (size_ > 1))
{
chP* xS = chemPSibling(x);
chemPoint* xS = chemPSibling(x);
if (xS != nullptr)
{
n2ndSearch_++;
@ -537,13 +550,13 @@ bool Foam::binaryTree<CompType, ThermoType>::secondaryBTSearch
return true;
}
}
else if (inSubTree(phiq, nodeSibling(x),x))
else if (inSubTree(phiq, nodeSibling(x), x))
{
return true;
}
// if we reach this point, no leafs were found at this depth or lower
// we move upward in the tree
bn* y = x->node();
node* y = x->node();
while((y->parent()!= nullptr) && (n2ndSearch_ < max2ndSearch_))
{
xS = chemPSibling(y);
@ -575,7 +588,7 @@ bool Foam::binaryTree<CompType, ThermoType>::secondaryBTSearch
template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::deleteLeaf(chP*& phi0)
void Foam::binaryTree<CompType, ThermoType>::deleteLeaf(chemPoint*& phi0)
{
if (size_ == 1) // only one point is stored
{
@ -584,16 +597,16 @@ void Foam::binaryTree<CompType, ThermoType>::deleteLeaf(chP*& phi0)
}
else if (size_ > 1)
{
bn* z = phi0->node();
bn* x;
chP* siblingPhi0 = chemPSibling(phi0);
node* z = phi0->node();
node* x;
chemPoint* siblingPhi0 = chemPSibling(phi0);
if (siblingPhi0 != nullptr)// the sibling of phi0 is a chemPoint
{
// z was root (only two chemPoints in the tree)
if (z->parent() == nullptr)
{
root_ = new bn();
root_ = new node();
root_->leafLeft()=siblingPhi0;
siblingPhi0->node()=root_;
}
@ -633,38 +646,40 @@ void Foam::binaryTree<CompType, ThermoType>::deleteLeaf(chP*& phi0)
deleteDemandDrivenData(phi0);
deleteDemandDrivenData(z);
}
size_--;
--size_;
}
template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::balance()
{
scalarField mean(chemistry_.nEqns(),0.0);
//1) walk through the entire tree by starting with the tree's most left
// chemPoint
chP* x = treeMin();
List<chP*> chemPoints(size_);
label chPi=0;
chemPoint* x = treeMin();
List<chemPoint*> chemPoints(size_);
label chemPointi = 0;
//2) compute the mean composition
label n = x->phi().size();
scalarField mean(n, 0.0);
while (x != nullptr)
{
const scalarField& phij = x->phi();
mean += phij;
chemPoints[chPi++] = x;
chemPoints[chemPointi++] = x;
x = treeSuccessor(x);
}
mean /= size_;
//3) compute the variance for each space direction
List<scalar> variance(chemistry_.nEqns(),0.0);
List<scalar> variance(n, 0.0);
forAll(chemPoints, j)
{
const scalarField& phij = chemPoints[j]->phi();
forAll(variance, vi)
{
variance[vi] += sqr(phij[vi]-mean[vi]);
variance[vi] += sqr(phij[vi] - mean[vi]);
}
}
@ -684,7 +699,7 @@ void Foam::binaryTree<CompType, ThermoType>::balance()
// in this direction if these extreme points were not deleted in the
// cleaning that come before the balance function they are still important
// and the tree should therefore take them into account
SortableList<scalar> phiMaxDir(chemPoints.size(),0.0);
SortableList<scalar> phiMaxDir(chemPoints.size(), 0.0);
forAll(chemPoints, j)
{
phiMaxDir[j] = chemPoints[j]->phi()[maxDir];
@ -696,7 +711,7 @@ void Foam::binaryTree<CompType, ThermoType>::balance()
root_ = nullptr;
// add the node for the two extremum
bn* newNode = new bn
node* newNode = new node
(
chemPoints[phiMaxDir.indices()[0]],
chemPoints[phiMaxDir.indices()[phiMaxDir.size()-1]],
@ -707,9 +722,9 @@ void Foam::binaryTree<CompType, ThermoType>::balance()
chemPoints[phiMaxDir.indices()[0]]->node() = newNode;
chemPoints[phiMaxDir.indices()[phiMaxDir.size()-1]]->node() = newNode;
for (label cpi=1; cpi<chemPoints.size()-1; cpi++)
for (label cpi=1; cpi<chemPoints.size()-1; ++cpi)
{
chP* phi0;
chemPoint* phi0;
binaryTreeSearch
(
chemPoints[phiMaxDir.indices()[cpi]]->phi(),
@ -717,8 +732,8 @@ void Foam::binaryTree<CompType, ThermoType>::balance()
phi0
);
// add the chemPoint
bn* nodeToAdd =
new bn(phi0,chemPoints[phiMaxDir.indices()[cpi]], phi0->node());
node* nodeToAdd =
new node(phi0, chemPoints[phiMaxDir.indices()[cpi]], phi0->node());
// make the parent of phi0 point to the newly created node
insertNode(phi0, nodeToAdd);
phi0->node() = nodeToAdd;
@ -729,7 +744,7 @@ void Foam::binaryTree<CompType, ThermoType>::balance()
template<class CompType, class ThermoType>
Foam::chemPointISAT<CompType, ThermoType>*
Foam::binaryTree<CompType, ThermoType>::treeMin(bn* subTreeRoot)
Foam::binaryTree<CompType, ThermoType>::treeMin(node* subTreeRoot)
{
if (subTreeRoot!=nullptr)
{
@ -748,9 +763,9 @@ Foam::binaryTree<CompType, ThermoType>::treeMin(bn* subTreeRoot)
template<class CompType, class ThermoType>
Foam::chemPointISAT<CompType, ThermoType>*
Foam::binaryTree<CompType, ThermoType>::treeSuccessor(chP* x)
Foam::binaryTree<CompType, ThermoType>::treeSuccessor(chemPoint* x)
{
if (size_>1)
if (size_ > 1)
{
if (x == x->node()->leafLeft())
{
@ -765,8 +780,8 @@ Foam::binaryTree<CompType, ThermoType>::treeSuccessor(chP* x)
}
else if (x == x->node()->leafRight())
{
bn* y = x->node();
while((y->parent() !=nullptr))
node* y = x->node();
while ((y->parent() != nullptr))
{
if (y == y->parent()->nodeLeft())
{
@ -779,7 +794,7 @@ Foam::binaryTree<CompType, ThermoType>::treeSuccessor(chP* x)
return treeMin(y->parent()->nodeRight());
}
}
y=y->parent();
y = y->parent();
}
// when we reach this point, y points to the root and
// never entered in the if loop (coming from the right)
@ -791,6 +806,7 @@ Foam::binaryTree<CompType, ThermoType>::treeSuccessor(chP* x)
FatalErrorInFunction
<< "inconsistent structure of the tree, no leaf and no node"
<< exit(FatalError);
return nullptr;
}
}
@ -823,19 +839,19 @@ bool Foam::binaryTree<CompType, ThermoType>::isFull()
template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::resetNumRetrieve()
{
// Has to go along each chP of the tree
// Has to go along each chemPoint of the tree
if (size_ > 0)
{
// First finds the first leaf
chP* chP0 = treeMin();
chP0->resetNumRetrieve();
chemPoint* chemPoint0 = treeMin();
chemPoint0->resetNumRetrieve();
// Then go to each successor
chP* nextChP = treeSuccessor(chP0);
while (nextChP != nullptr)
chemPoint* nextchemPoint = treeSuccessor(chemPoint0);
while (nextchemPoint != nullptr)
{
nextChP->resetNumRetrieve();
nextChP = treeSuccessor(nextChP);
nextchemPoint->resetNumRetrieve();
nextchemPoint = treeSuccessor(nextchemPoint);
}
}
}

56
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryTree/binaryTree.H
View File

@ -59,8 +59,10 @@ class binaryTree
{
public:
typedef binaryNode<CompType, ThermoType> bn;
typedef chemPointISAT<CompType, ThermoType> chP;
typedef binaryNode<CompType, ThermoType> node;
typedef chemPointISAT<CompType, ThermoType> chemPoint;
private:
@ -68,7 +70,7 @@ private:
TDACChemistryModel<CompType, ThermoType>& chemistry_;
//- Root node of the binary tree
bn *root_;
node *root_;
//- Maximum number of elements in the binary tree
label maxNLeafs_;
@ -85,8 +87,8 @@ private:
// will be lost
void insertNode
(
chP*& phi0,
bn*& newNode
chemPoint*& phi0,
node*& newNode
);
//- Perform a search in the subtree starting from the subtree node y
@ -95,8 +97,8 @@ private:
bool inSubTree
(
const scalarField& phiq,
bn* y,
chP* x
node* y,
chemPoint* x
);
void deleteSubTree(binaryNode<CompType, ThermoType>* subTreeRoot);
@ -107,20 +109,22 @@ private:
}
//- Replace the binaryNode u with v
void transplant(bn* u, bn* v);
void transplant(node* u, node* v);
chP* chemPSibling(bn* y);
chP* chemPSibling(chP* x);
chemPoint* chemPSibling(node* y);
chemPoint* chemPSibling(chemPoint* x);
bn* nodeSibling(bn* y);
bn* nodeSibling(chP* x);
node* nodeSibling(node* y);
node* nodeSibling(chemPoint* x);
void deleteAllNode(bn* subTreeRoot);
void deleteAllNode(node* subTreeRoot);
dictionary coeffsDict_;
public:
//- Constructors
// Constructors
//- Construct from dictionary and chemistryOnLineLibrary
binaryTree
@ -129,21 +133,23 @@ public:
dictionary coeffsDict
);
//- Member functions
// Member functions
inline label size()
{
return size_;
}
//- Computes iteratively the depth of the subTree
label depth(bn* subTreeRoot);
label depth(node* subTreeRoot);
inline label depth()
{
return depth(root_);
}
inline bn* root()
inline node* root()
{
return root_;
}
@ -177,7 +183,7 @@ public:
const scalarField& scaleFactor,
const scalar& epsTol,
const label nCols,
chP*& phi0
chemPoint*& phi0
);
@ -187,20 +193,20 @@ public:
void binaryTreeSearch
(
const scalarField& phiq,
bn* node,
chP*& nearest
node* node,
chemPoint*& nearest
);
// Perform a secondary binary tree search starting from a failed
// chemPoint x, with a depth-first search algorithm
// If another candidate is found return true and x points to the chemP
bool secondaryBTSearch(const scalarField& phiq, chP*& x);
bool secondaryBTSearch(const scalarField& phiq, chemPoint*& x);
//- Delete a leaf from the binary tree and reshape the binary tree for
// the following binary tree search
// Return the index in the nodeList of the removed node
// (-1 when no node)
void deleteLeaf(chP*& phi0);
void deleteLeaf(chemPoint*& phi0);
//- Cheap balance function
// This function just roughly separate the space in two parts
@ -216,14 +222,14 @@ public:
deleteAllNode(root_);
}
chP* treeMin(bn* subTreeRoot);
chemPoint* treeMin(node* subTreeRoot);
inline chP* treeMin()
inline chemPoint* treeMin()
{
return treeMin(root_);
}
chP* treeSuccessor(chP* x);
chemPoint* treeSuccessor(chemPoint* x);
//- Removes every entries of the tree and delete the associated objects
void clear();

134
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C
View File

@ -46,25 +46,26 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
scalarField d(nCols);
scalar scale, sigma, sum;
for (label k=0; k<nCols-1; k++)
for (label k=0; k<nCols-1; ++k)
{
scale = 0;
for (label i=k; i<nCols; i++)
for (label i=k; i<nCols; ++i)
{
scale=max(scale, mag(R(i, k)));
scale = max(scale, mag(R(i, k)));
}
if (scale == 0)
{
c[k] = d[k] = 0;
}
else
{
for (label i=k; i<nCols; i++)
for (label i=k; i<nCols; ++i)
{
R(i, k) /= scale;
}
sum = 0;
for (label i=k; i<nCols; i++)
for (label i=k; i<nCols; ++i)
{
sum += sqr(R(i, k));
}
@ -72,15 +73,15 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
R(k, k) += sigma;
c[k] = sigma*R(k, k);
d[k] = -scale*sigma;
for (label j=k+1; j<nCols; j++)
for (label j=k+1; j<nCols; ++j)
{
sum = 0;
for ( label i=k; i<nCols; i++)
for ( label i=k; i<nCols; ++i)
{
sum += R(i, k)*R(i, j);
}
scalar tau = sum/c[k];
for ( label i=k; i<nCols; i++)
for ( label i=k; i<nCols; ++i)
{
R(i, j) -= tau*R(i, k);
}
@ -91,12 +92,12 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
d[nCols-1] = R(nCols-1, nCols-1);
// form R
for (label i=0; i<nCols; i++)
for (label i=0; i<nCols; ++i)
{
R(i, i) = d[i];
for ( label j=0; j<i; j++)
for ( label j=0; j<i; ++j)
{
R(i, j)=0;
R(i, j) = 0;
}
}
}
@ -114,7 +115,7 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrUpdate
label k;
scalarField w(u);
for (k=n-1; k>=0; k--)
for (k=n-1; k>=0; --k)
{
if (w[k] != 0)
{
@ -127,7 +128,7 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrUpdate
k = 0;
}
for (label i=k-1; i>=0; i--)
for (label i=k-1; i>=0; --i)
{
rotate(R, i, w[i],-w[i+1], n);
if (w[i] == 0)
@ -144,12 +145,12 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrUpdate
}
}
for (label i=0; i<n; i++)
for (label i=0; i<n; ++i)
{
R(0, i) += w[0]*v[i];
}
for (label i=0; i<k; i++)
for (label i=0; i<k; ++i)
{
rotate(R, i, R(i, i), -R(i+1, i), n);
}
@ -184,7 +185,8 @@ void Foam::chemPointISAT<CompType, ThermoType>::rotate
s = sign(b)/sqrt(1.0 + sqr(fact));
c = fact*s;
}
for (label j=i;j<n;j++)
for (label j=i; j<n ;++j)
{
y = R(i, j);
w = R(i+1, j);
@ -223,8 +225,8 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
timeTag_(chemistry_.timeSteps()),
lastTimeUsed_(chemistry_.timeSteps()),
toRemove_(false),
maxNumNewDim_(coeffsDict.lookupOrDefault("maxNumNewDim",0)),
printProportion_(coeffsDict.lookupOrDefault("printProportion",false)),
maxNumNewDim_(coeffsDict.lookupOrDefault("maxNumNewDim", 0)),
printProportion_(coeffsDict.lookupOrDefault("printProportion", false)),
numRetrieve_(0),
nLifeTime_(0),
completeToSimplifiedIndex_
@ -251,12 +253,12 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
bool isMechRedActive = chemistry_.mechRed()->active();
if (isMechRedActive)
{
for (label i=0; i<completeSpaceSize-nAdditionalEqns_; i++)
for (label i=0; i<completeSpaceSize-nAdditionalEqns_; ++i)
{
completeToSimplifiedIndex_[i] =
chemistry.completeToSimplifiedIndex()[i];
}
for (label i=0; i<nActiveSpecies_; i++)
for (label i=0; i<nActiveSpecies_; ++i)
{
simplifiedToCompleteIndex_[i] =
chemistry.simplifiedToCompleteIndex()[i];
@ -276,7 +278,7 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
const scalarDiagonalMatrix& S = svdA.S();
// Replace the value of vector D by max(D, 1/2), first ISAT paper
for (label i=0; i<reduOrCompDim; i++)
for (label i=0; i<reduOrCompDim; ++i)
{
D[i] = max(S[i], 0.5);
}
@ -287,9 +289,9 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
// Result stored in Atilde
multiply(Atilde, svdA.U(), D, svdA.V().T());
for (label i=0; i<reduOrCompDim; i++)
for (label i=0; i<reduOrCompDim; ++i)
{
for (label j=0; j<reduOrCompDim; j++)
for (label j=0; j<reduOrCompDim; ++j)
{
label compi = i;
@ -379,7 +381,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
scalar epsTemp = 0;
List<scalar> propEps(completeSpaceSize(), scalar(0));
for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; ++i)
{
scalar temp = 0;
@ -392,11 +394,11 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
||(isMechRedActive && completeToSimplifiedIndex_[i] != -1)
)
{
label si=(isMechRedActive) ? completeToSimplifiedIndex_[i] : i;
label si = isMechRedActive ? completeToSimplifiedIndex_[i] : i;
for (label j=si; j<dim; j++)// LT is upper triangular
for (label j=si; j<dim; ++j)// LT is upper triangular
{
label sj=(isMechRedActive) ? simplifiedToCompleteIndex_[j] : j;
label sj = isMechRedActive ? simplifiedToCompleteIndex_[j] : j;
temp += LT_(si, j)*dphi[sj];
}
@ -469,13 +471,11 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
+ LT_(dim, dim+1)*dphi[idp_]
);
propEps[idp_] =
sqr(LT_(dim+1, dim+1)*dphi[idp_]);
propEps[idp_] = sqr(LT_(dim+1, dim+1)*dphi[idp_]);
if (variableTimeStep())
{
propEps[iddeltaT_] =
sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
propEps[iddeltaT_] = sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
}
}
@ -486,7 +486,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
{
scalar max = -1;
label maxIndex = -1;
for (label i=0; i<completeSpaceSize(); i++)
for (label i=0; i<completeSpaceSize(); ++i)
{
if (max < propEps[i])
{
@ -514,9 +514,10 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
{
propName = chemistry_.Y()[maxIndex].name();
}
Info<< "Direction maximum impact to error in ellipsoid: "
<< propName << endl;
Info<< "Proportion to the total error on the retrieve: "
<< propName << nl
<< "Proportion to the total error on the retrieve: "
<< max/(epsTemp+SMALL) << endl;
}
return false;
@ -542,7 +543,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
const scalarSquareMatrix& Avar(A());
bool isMechRedActive = chemistry_.mechRed()->active();
scalar dRl = 0;
label dim = completeSpaceSize()-2;
label dim = completeSpaceSize() - 2;
if (isMechRedActive)
{
dim = nActiveSpecies_;
@ -550,7 +551,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
// Since we build only the solution for the species, T and p are not
// included
for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; ++i)
{
dRl = 0;
if (isMechRedActive)
@ -560,7 +561,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
// If this species is active
if (si != -1)
{
for (label j=0; j<dim; j++)
for (label j=0; j<dim; ++j)
{
label sj=simplifiedToCompleteIndex_[j];
dRl += Avar(si, j)*dphi[sj];
@ -579,7 +580,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
}
else
{
for (label j=0; j<completeSpaceSize(); j++)
for (label j=0; j<completeSpaceSize(); ++j)
{
dRl += Avar(i, j)*dphi[j];
}
@ -615,7 +616,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
// check if the difference of active species is lower than the maximum
// number of new dimensions allowed
for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; ++i)
{
// first test if the current chemPoint has an inactive species
// corresponding to an active one in the query point
@ -625,7 +626,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
&& chemistry_.completeToSimplifiedIndex()[i]!=-1
)
{
activeAdded++;
++activeAdded;
dimToAdd.append(i);
}
// then test if an active species in the current chemPoint
@ -636,7 +637,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
&& chemistry_.completeToSimplifiedIndex()[i] == -1
)
{
activeAdded++;
++activeAdded;
// we don't need to add a new dimension but we count it to have
// control on the difference through maxNumNewDim
}
@ -649,7 +650,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
&& dphi[i] != 0
)
{
activeAdded++;
++activeAdded;
dimToAdd.append(i);
}
}
@ -661,11 +662,12 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
}
// the number of added dimension to the current chemPoint
nActiveSpecies_ += dimToAdd.size();
const label nDimAdd = dimToAdd.size();
nActiveSpecies_ += nDimAdd;
simplifiedToCompleteIndex_.setSize(nActiveSpecies_);
forAll(dimToAdd, i)
{
label si = nActiveSpecies_ - dimToAdd.size() + i;
label si = nActiveSpecies_ - nDimAdd + i;
// add the new active species
simplifiedToCompleteIndex_[si] = dimToAdd[i];
completeToSimplifiedIndex_[dimToAdd[i]] = si;
@ -685,9 +687,9 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
A_ = scalarSquareMatrix(nActiveSpecies_+nAdditionalEqns_, Zero);
// write the initial active species
for (label i=0; i<initNActiveSpecies; i++)
for (label i=0; i<initNActiveSpecies; ++i)
{
for (label j=0; j<initNActiveSpecies; j++)
for (label j=0; j<initNActiveSpecies; ++j)
{
LT_(i, j) = LTvar(i, j);
A_(i, j) = Avar(i, j);
@ -695,37 +697,37 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
}
// write the columns for temperature and pressure
for (label i=0; i<initNActiveSpecies; i++)
for (label i=0; i<initNActiveSpecies; ++i)
{
for (label j=1; j>=0; j--)
for (label j=1; j>=0; --j)
{
LT_(i, nActiveSpecies_+j)=LTvar(i, initNActiveSpecies+j);
A_(i, nActiveSpecies_+j)=Avar(i, initNActiveSpecies+j);
LT_(nActiveSpecies_+j, i)=LTvar(initNActiveSpecies+j, i);
A_(nActiveSpecies_+j, i)=Avar(initNActiveSpecies+j, i);
LT_(i, nActiveSpecies_+j) = LTvar(i, initNActiveSpecies+j);
A_(i, nActiveSpecies_+j) = Avar(i, initNActiveSpecies+j);
LT_(nActiveSpecies_+j, i) = LTvar(initNActiveSpecies+j, i);
A_(nActiveSpecies_+j, i) = Avar(initNActiveSpecies+j, i);
}
}
// end with the diagonal elements for temperature and pressure
LT_(nActiveSpecies_, nActiveSpecies_)=
LT_(nActiveSpecies_, nActiveSpecies_) =
LTvar(initNActiveSpecies, initNActiveSpecies);
A_(nActiveSpecies_, nActiveSpecies_)=
A_(nActiveSpecies_, nActiveSpecies_) =
Avar(initNActiveSpecies, initNActiveSpecies);
LT_(nActiveSpecies_+1, nActiveSpecies_+1)=
LT_(nActiveSpecies_+1, nActiveSpecies_+1) =
LTvar(initNActiveSpecies+1, initNActiveSpecies+1);
A_(nActiveSpecies_+1, nActiveSpecies_+1)=
A_(nActiveSpecies_+1, nActiveSpecies_+1) =
Avar(initNActiveSpecies+1, initNActiveSpecies+1);
if (variableTimeStep())
{
LT_(nActiveSpecies_+2, nActiveSpecies_+2)=
LT_(nActiveSpecies_+2, nActiveSpecies_+2) =
LTvar(initNActiveSpecies+2, initNActiveSpecies+2);
A_(nActiveSpecies_+2, nActiveSpecies_+2)=
A_(nActiveSpecies_+2, nActiveSpecies_+2) =
Avar(initNActiveSpecies+2, initNActiveSpecies+2);
}
for (label i=initNActiveSpecies; i<nActiveSpecies_;i++)
for (label i=initNActiveSpecies; i<nActiveSpecies_; ++i)
{
LT_(i, i)=
LT_(i, i) =
1.0
/(tolerance_*scaleFactor_[simplifiedToCompleteIndex_[i]]);
A_(i, i) = 1;
@ -740,9 +742,9 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
scalar normPhiTilde = 0;
// p' = L^T.(p-phi)
for (label i=0; i<dim; i++)
for (label i=0; i<dim; ++i)
{
for (label j=i; j<dim-nAdditionalEqns_; j++)// LT is upper triangular
for (label j=i; j<dim-nAdditionalEqns_; ++j)// LT is upper triangular
{
label sj = j;
if (isMechRedActive)
@ -770,16 +772,16 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
scalarField u(gamma*phiTilde);
scalarField v(dim, 0);
for (label i=0; i<dim; i++)
for (label i=0; i<dim; ++i)
{
for (label j=0; j<=i;j++)
for (label j=0; j<=i; ++j)
{
v[i] += phiTilde[j]*LT_(j, i);
}
}
qrUpdate(LT_,dim, u, v);
nGrowth_++;
qrUpdate(LT_, dim, u, v);
++nGrowth_;
return true;
}

10
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
View File

@ -142,7 +142,7 @@ class chemPointISAT
{
// Private data
//- Pointer to the chemistryModel object
//- Reference to the chemistryModel object
TDACChemistryModel<CompType, ThermoType>& chemistry_;
//- Vector storing the composition, temperature and pressure
@ -195,7 +195,7 @@ class chemPointISAT
label numRetrieve_;
//- Variable to store the number of time steps the chempoint is allowed
// to still live according to the maxChPLifeTime_ parameter
// to still live according to the maxChPLifeTime_ parameter
label nLifeTime_;
List<label> completeToSimplifiedIndex_;
@ -402,9 +402,9 @@ public:
// ISAT functions
//- To RETRIEVE the mapping from the stored chemPoint phi, the query
// point phiq has to be in the EOA of phi.
// To test if phiq is in the ellipsoid:
// ||L^T.dphi|| <= 1
// point phiq has to be in the EOA of phi.
// To test if phiq is in the ellipsoid:
// ||L^T.dphi|| <= 1
bool inEOA(const scalarField& phiq);
//- More details about the minimum-volume ellipsoid covering an