Commit 884f2f86 authored by Andrew Heather's avatar Andrew Heather
Browse files

STYLE: Code clean-up

parent 94b62c28
......@@ -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];
......
......@@ -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;
}
}
......
......@@ -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;
}
}
}
......
......@@ -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
......
......@@ -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;