STYLE: code clean-up/restructure

This commit is contained in:
andy
2010-08-06 12:09:06 +01:00
parent f243736c78
commit 02cd5c05da
5 changed files with 58 additions and 52 deletions

View File

@ -151,7 +151,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
) const ) const
{ {
scalarField c2(nSpecie_, 0.0); scalarField c2(nSpecie_, 0.0);
for (label i=0; i<nSpecie_; i++) for (label i = 0; i < nSpecie_; i++)
{ {
c2[i] = max(0.0, c[i]); c2[i] = max(0.0, c[i]);
} }
@ -169,7 +169,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
lRef = R.lhs()[slRef].index; lRef = R.lhs()[slRef].index;
pf = kf; pf = kf;
for (label s=1; s<Nl; s++) for (label s = 1; s < Nl; s++)
{ {
const label si = R.lhs()[s].index; const label si = R.lhs()[s].index;
@ -212,7 +212,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
// find the matrix element and element position for the rhs // find the matrix element and element position for the rhs
pr = kr; pr = kr;
for (label s=1; s<Nr; s++) for (label s = 1; s < Nr; s++)
{ {
const label si = R.rhs()[s].index; const label si = R.rhs()[s].index;
if (c[si] < c[rRef]) if (c[si] < c[rRef])
@ -270,7 +270,7 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
// dT/dt = ... // dT/dt = ...
scalar rho = 0.0; scalar rho = 0.0;
scalar cSum = 0.0; scalar cSum = 0.0;
for (label i=0; i<nSpecie_; i++) for (label i = 0; i < nSpecie_; i++)
{ {
const scalar W = specieThermo_[i].W(); const scalar W = specieThermo_[i].W();
cSum += c[i]; cSum += c[i];
@ -287,7 +287,7 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
cp /= mw; cp /= mw;
scalar dT = 0.0; scalar dT = 0.0;
for (label i=0; i<nSpecie_; i++) for (label i = 0; i < nSpecie_; i++)
{ {
const scalar hi = specieThermo_[i].h(T); const scalar hi = specieThermo_[i].h(T);
dT += hi*dcdt[i]; dT += hi*dcdt[i];
@ -552,6 +552,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::Sh() const
) )
); );
if (this->chemistry_) if (this->chemistry_)
{ {
scalarField& Sh = tSh(); scalarField& Sh = tSh();
@ -560,7 +561,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::Sh() const
{ {
forAll(Sh, cellI) forAll(Sh, cellI)
{ {
scalar hi = specieThermo_[i].Hc(); const scalar hi = specieThermo_[i].Hc();
Sh[cellI] -= hi*RR_[i][cellI]; Sh[cellI] -= hi*RR_[i][cellI];
} }
} }
@ -734,7 +735,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
t += dt; t += dt;
// update the temperature // update the temperature
scalar cTot = sum(c); const scalar cTot = sum(c);
ThermoType mixture(0.0*specieThermo_[0]); ThermoType mixture(0.0*specieThermo_[0]);
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {

View File

@ -39,7 +39,7 @@ Foam::EulerImplicit<CompType, ThermoType>::EulerImplicit
chemistrySolver<CompType, ThermoType>(model, modelName), chemistrySolver<CompType, ThermoType>(model, modelName),
coeffsDict_(model.subDict(modelName + "Coeffs")), coeffsDict_(model.subDict(modelName + "Coeffs")),
cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))), cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
equil_(coeffsDict_.lookup("equilibriumRateLimiter")) eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
{} {}
@ -65,12 +65,12 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
const label nSpecie = this->model_.nSpecie(); const label nSpecie = this->model_.nSpecie();
simpleMatrix<scalar> RR(nSpecie, 0, 0); simpleMatrix<scalar> RR(nSpecie, 0, 0);
for (label i=0; i<nSpecie; i++) for (label i = 0; i < nSpecie; i++)
{ {
c[i] = max(0.0, c[i]); c[i] = max(0.0, c[i]);
} }
for (label i=0; i<nSpecie; i++) for (label i = 0; i < nSpecie; i++)
{ {
RR.source()[i] = c[i]/dt; RR.source()[i] = c[i]/dt;
} }
@ -88,9 +88,9 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
); );
scalar corr = 1.0; scalar corr = 1.0;
if (equil_) if (eqRateLimiter_)
{ {
if (omegai<0.0) if (omegai < 0.0)
{ {
corr = 1.0/(1.0 + pr*dt); corr = 1.0/(1.0 + pr*dt);
} }
@ -100,31 +100,31 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
} }
} }
forAll(R.lhs(), s) forAll(R.lhs(), specieI)
{ {
label si = R.lhs()[s].index; const label id = R.lhs()[specieI].index;
scalar sl = R.lhs()[s].stoichCoeff; const scalar sc = R.lhs()[specieI].stoichCoeff;
RR[si][rRef] -= sl*pr*corr; RR[id][rRef] -= sc*pr*corr;
RR[si][lRef] += sl*pf*corr; RR[id][lRef] += sc*pf*corr;
} }
forAll(R.rhs(), s) forAll(R.rhs(), specieI)
{ {
label si = R.rhs()[s].index; const label id = R.rhs()[specieI].index;
scalar sr = R.rhs()[s].stoichCoeff; const scalar sc = R.rhs()[specieI].stoichCoeff;
RR[si][lRef] -= sr*pf*corr; RR[id][lRef] -= sc*pf*corr;
RR[si][rRef] += sr*pr*corr; RR[id][rRef] += sc*pr*corr;
} }
} }
for (label i=0; i<nSpecie; i++) for (label i = 0; i < nSpecie; i++)
{ {
RR[i][i] += 1.0/dt; RR[i][i] += 1.0/dt;
} }
c = RR.LUsolve(); c = RR.LUsolve();
for (label i=0; i<nSpecie; i++) for (label i = 0; i < nSpecie; i++)
{ {
c[i] = max(0.0, c[i]); c[i] = max(0.0, c[i]);
} }
@ -134,7 +134,7 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
const label nEqns = this->model_.nEqns(); const label nEqns = this->model_.nEqns();
scalarField c1(nEqns, 0.0); scalarField c1(nEqns, 0.0);
for (label i=0; i<nSpecie; i++) for (label i = 0; i < nSpecie; i++)
{ {
c1[i] = c[i]; c1[i] = c[i];
} }
@ -144,9 +144,9 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
scalarField dcdt(nEqns, 0.0); scalarField dcdt(nEqns, 0.0);
this->model_.derivatives(0.0, c1, dcdt); this->model_.derivatives(0.0, c1, dcdt);
scalar sumC = sum(c); const scalar sumC = sum(c);
for (label i=0; i<nSpecie; i++) for (label i = 0; i < nSpecie; i++)
{ {
scalar d = dcdt[i]; scalar d = dcdt[i];
if (d < -SMALL) if (d < -SMALL)

View File

@ -57,12 +57,17 @@ class EulerImplicit
{ {
// Private data // Private data
//- Coefficients dictionary
dictionary coeffsDict_; dictionary coeffsDict_;
// Model constants // Model constants
//- Chemistry timescale
scalar cTauChem_; scalar cTauChem_;
Switch equil_;
//- Equilibrium rate limiter flag (on/off)
Switch eqRateLimiter_;
public: public:

View File

@ -38,7 +38,7 @@ Foam::sequential<CompType, ThermoType>::sequential
chemistrySolver<CompType, ThermoType>(model, modelName), chemistrySolver<CompType, ThermoType>(model, modelName),
coeffsDict_(model.subDict(modelName + "Coeffs")), coeffsDict_(model.subDict(modelName + "Coeffs")),
cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))), cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
equil_(coeffsDict_.lookup("equilibriumRateLimiter")) eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
{} {}
@ -70,45 +70,41 @@ Foam::scalar Foam::sequential<CompType, ThermoType>::solve
{ {
const Reaction<ThermoType>& R = this->model_.reactions()[i]; const Reaction<ThermoType>& R = this->model_.reactions()[i];
scalar om0 = this->model_.omega scalar omega = this->model_.omega
( (
R, c, T, p, pf, cf, lRef, pb, cb, rRef R, c, T, p, pf, cf, lRef, pb, cb, rRef
); );
scalar omeg = 0.0; if (eqRateLimiter_)
if (!equil_)
{ {
omeg = om0; if (omega < 0.0)
}
else
{
if (om0<0.0)
{ {
omeg = om0/(1.0 + pb*dt); omega /= 1.0 + pb*dt;
} }
else else
{ {
omeg = om0/(1.0 + pf*dt); omega /= 1.0 + pf*dt;
} }
} }
tChemInv = max(tChemInv, mag(omeg));
tChemInv = max(tChemInv, mag(omega));
// update species // update species
forAll(R.lhs(), s) forAll(R.lhs(), specieI)
{ {
label si = R.lhs()[s].index; const label id = R.lhs()[specieI].index;
scalar sl = R.lhs()[s].stoichCoeff; const scalar sc = R.lhs()[specieI].stoichCoeff;
c[si] -= dt*sl*omeg; c[id] -= dt*sc*omega;
c[si] = max(0.0, c[si]); c[id] = max(0.0, c[id]);
} }
forAll(R.rhs(), s) forAll(R.rhs(), specieI)
{ {
label si = R.rhs()[s].index; const label id = R.rhs()[specieI].index;
scalar sr = R.rhs()[s].stoichCoeff; const scalar sc = R.rhs()[specieI].stoichCoeff;
c[si] += dt*sr*omeg; c[id] += dt*sc*omega;
c[si] = max(0.0, c[si]); c[id] = max(0.0, c[id]);
} }
} }

View File

@ -59,12 +59,17 @@ class sequential
{ {
// Private data // Private data
//- Coefficients dictionary
dictionary coeffsDict_; dictionary coeffsDict_;
// Model constants // Model constants
//- Chemistry time scale
scalar cTauChem_; scalar cTauChem_;
Switch equil_;
//- Equilibrium rate limiter flag (on/off)
Switch eqRateLimiter_;
public: public:
@ -75,7 +80,6 @@ public:
// Constructors // Constructors
//- Construct from components //- Construct from components
sequential sequential
( (