Commit c3261ef3 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[implicit] enhance readability for useMoles if-statements

Replace

if (!useMoles)
  // mass stuff
else
  // mole stuff

by

if (useMoles)
  // mole stuff
else
  // mass stuff 

Reviewed by fetzer.



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@15200 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent d9044035
......@@ -123,12 +123,12 @@ protected:
storage = 0.0;
if(!useMoles)
// mass balance mass fraction based
storage[massBalanceIdx] = volVars.density();
else
if(useMoles)
// mass balance mole fraction based
storage[massBalanceIdx] = volVars.molarDensity();
else
// mass balance mass fraction based
storage[massBalanceIdx] = volVars.density();
// momentum balance
for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; ++momentumIdx)
......@@ -188,15 +188,15 @@ protected:
DimVector massBalanceResidual = fluxVars.velocity();
if(!useMoles)
if(useMoles)
{
massBalanceResidual *= (massUpwindWeight_ * up.density()
+ (1.-massUpwindWeight_) * dn.density());
massBalanceResidual *= (massUpwindWeight_ * up.molarDensity()
+ (1.-massUpwindWeight_) * dn.molarDensity());
}
else
{
massBalanceResidual *= (massUpwindWeight_ * up.molarDensity()
+ (1.-massUpwindWeight_) * dn.molarDensity());
massBalanceResidual *= (massUpwindWeight_ * up.density()
+ (1.-massUpwindWeight_) * dn.density());
}
if (!fluxVars.onBoundary())
......
......@@ -112,35 +112,35 @@ public:
this->prevVolVars_() : this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
if (!useMoles)
{
/* works for a maximum of two components, for more components
mole fractions must be used (property useMoles => true) */
if (useMoles)
{
// mass balance and transport equations
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx)
//else // transport equations
{
storage[conti0EqIdx+compIdx] = volVars.molarDensity()
* volVars.moleFraction(compIdx);
//storage of transported component
storage[transportEqIdx] = volVars.density()
* volVars.massFraction(transportCompIdx);
Valgrind::CheckDefined(volVars.molarDensity());
Valgrind::CheckDefined(volVars.moleFraction(compIdx));
}
}
}
else
{
/* works for a maximum of two components, for more components
mole fractions must be used (set property UseMoles to true) */
Valgrind::CheckDefined(volVars.density());
Valgrind::CheckDefined(volVars.massFraction(transportCompIdx));
//storage of transported component
storage[transportEqIdx] = volVars.density()
* volVars.massFraction(transportCompIdx);
}
else
{
// mass balance and transport equations
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx)
//else // transport equations
{
storage[conti0EqIdx+compIdx] = volVars.molarDensity()
* volVars.moleFraction(compIdx);
Valgrind::CheckDefined(volVars.molarDensity());
Valgrind::CheckDefined(volVars.moleFraction(compIdx));
}
}
}
Valgrind::CheckDefined(volVars.density());
Valgrind::CheckDefined(volVars.massFraction(transportCompIdx));
}
}
/*!
......@@ -165,31 +165,30 @@ public:
Scalar tmp = fluxVars.normalVelocity();
if(!useMoles)
{
if(useMoles)
{
//transport equations
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx) // mass balance is calculated above
{
tmp *= (this->massUpwindWeight_ * up.molarDensity() * up.moleFraction(compIdx)
+ (1.-this->massUpwindWeight_) * dn.molarDensity() * dn.moleFraction(compIdx));
flux[conti0EqIdx+compIdx] += tmp;
Valgrind::CheckDefined(flux[conti0EqIdx+compIdx]);
}
}
}
else
{
tmp *= (this->massUpwindWeight_ * up.density() * up.massFraction(transportCompIdx)
+ (1.-this->massUpwindWeight_) * dn.density() * dn.massFraction(transportCompIdx));
flux[transportEqIdx] += tmp;
Valgrind::CheckDefined(flux[transportEqIdx]);
}
else
{
//transport equations
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx) //mass balance is calculated above
{
tmp *= (this->massUpwindWeight_ * up.molarDensity() * up.moleFraction(compIdx)
+ (1.-this->massUpwindWeight_) * dn.molarDensity() * dn.moleFraction(compIdx));
flux[conti0EqIdx+compIdx] += tmp;
Valgrind::CheckDefined(flux[conti0EqIdx+compIdx]);
}
}
}
}
}
/*!
* \brief Adds the diffusive component flux to the flux vector over
......@@ -202,16 +201,7 @@ public:
const FluxVariables &fluxVars) const
{
// diffusive component flux
if(!useMoles)
{
flux[transportEqIdx] -= fluxVars.moleFractionGrad(transportCompIdx)
* fluxVars.face().normal
* (fluxVars.diffusionCoeff(transportCompIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity()
* FluidSystem::molarMass(transportCompIdx);// Multiplied by molarMass [kg/mol] to convert form [mol/m^3 s] to [kg/m^3 s]
Valgrind::CheckDefined(flux[transportEqIdx]);
}
else
if(useMoles)
{
//loop over secondary components
for (int compIdx=0; compIdx<numComponents; compIdx++)
......@@ -226,6 +216,15 @@ public:
}
}
}
else
{
flux[transportEqIdx] -= fluxVars.moleFractionGrad(transportCompIdx)
* fluxVars.face().normal
* (fluxVars.diffusionCoeff(transportCompIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity()
* FluidSystem::molarMass(transportCompIdx);// Multiplied by molarMass [kg/mol] to convert form [mol/m^3 s] to [kg/m^3 s]
Valgrind::CheckDefined(flux[transportEqIdx]);
}
}
};
......
......@@ -136,39 +136,39 @@ public:
const bool isOldSol = false)
{
if(!useMoles) //mass-fraction formulation
{
Scalar massOrMoleFrac[numComponents];
massOrMoleFrac[transportCompIdx] = priVars[massOrMoleFracIdx];
massOrMoleFrac[phaseCompIdx] = 1.0 - massOrMoleFrac[transportCompIdx];
// calculate average molar mass of the gas phase
Scalar M1 = FluidSystem::molarMass(transportCompIdx);
Scalar M2 = FluidSystem::molarMass(phaseCompIdx);
Scalar X2 = massOrMoleFrac[phaseCompIdx];
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
// convert mass to mole fractions and set the fluid state
fluidState.setMoleFraction(phaseIdx, transportCompIdx, massOrMoleFrac[transportCompIdx]*avgMolarMass/M1);
fluidState.setMoleFraction(phaseIdx, phaseCompIdx, massOrMoleFrac[phaseCompIdx]*avgMolarMass/M2);
}
else
{
Scalar moleFracPhase, sumMoleFrac(0.0);
//for components
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx)
{
sumMoleFrac += priVars[conti0EqIdx+compIdx];
fluidState.setMoleFraction(phaseIdx, compIdx, priVars[conti0EqIdx+compIdx]);
}
}
//molefraction for the main component (no primary variable)
moleFracPhase = 1 - sumMoleFrac;
fluidState.setMoleFraction(phaseIdx, phaseCompIdx, moleFracPhase);
}
if(useMoles) // mole-fraction formulation
{
Scalar moleFracPhase, sumMoleFrac(0.0);
// loop over components
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx)
{
sumMoleFrac += priVars[conti0EqIdx+compIdx];
fluidState.setMoleFraction(phaseIdx, compIdx, priVars[conti0EqIdx+compIdx]);
}
}
// mole fraction for the main component (no primary variable)
moleFracPhase = 1 - sumMoleFrac;
fluidState.setMoleFraction(phaseIdx, phaseCompIdx, moleFracPhase);
}
else // mass-fraction formulation
{
Scalar massOrMoleFrac[numComponents];
massOrMoleFrac[transportCompIdx] = priVars[massOrMoleFracIdx];
massOrMoleFrac[phaseCompIdx] = 1.0 - massOrMoleFrac[transportCompIdx];
// calculate average molar mass of the gas phase
Scalar M1 = FluidSystem::molarMass(transportCompIdx);
Scalar M2 = FluidSystem::molarMass(phaseCompIdx);
Scalar X2 = massOrMoleFrac[phaseCompIdx];
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
// convert mass to mole fractions and set the fluid state
fluidState.setMoleFraction(phaseIdx, transportCompIdx, massOrMoleFrac[transportCompIdx]*avgMolarMass/M1);
fluidState.setMoleFraction(phaseIdx, phaseCompIdx, massOrMoleFrac[phaseCompIdx]*avgMolarMass/M2);
}
}
/*!
......
......@@ -113,17 +113,17 @@ namespace Dumux
// storage term of continuity equation
storage[conti0EqIdx] += volVars.divU;
if(!useMoles)
if(useMoles)
{
//storage term of the transport equation - massfractions
storage[transportEqIdx] += volVars.fluidState().massFraction(phaseIdx, transportCompIdx)
* volVars.effPorosity;
// storage term of the transport equation - mole fractions
storage[transportEqIdx] += volVars.fluidState().moleFraction(phaseIdx, transportCompIdx)
* volVars.effPorosity;
}
else
{
// storage term of the transport equation - molefractions
storage[transportEqIdx] += volVars.fluidState().moleFraction(phaseIdx, transportCompIdx)
* volVars.effPorosity;
//storage term of the transport equation - mass fractions
storage[transportEqIdx] += volVars.fluidState().massFraction(phaseIdx, transportCompIdx)
* volVars.effPorosity;
}
}
......@@ -203,61 +203,60 @@ namespace Dumux
if(withStabilization_)
flux[conti0EqIdx] -= stabilizationTerm;
if(!useMoles)
if(useMoles)
{
// mass flux of the dissolved second component - massfraction
// advective flux of the component
flux[transportEqIdx] +=
fluxVars.KmvpNormal() *
(( upwindWeight_)* up.fluidState().massFraction(phaseIdx, transportCompIdx)/up.viscosity()
+
(1 - upwindWeight_)* dn.fluidState().massFraction(phaseIdx, transportCompIdx)/dn.viscosity());
// flux of the dissolved second component due to solid displacement
flux[transportEqIdx] +=
fluxVars.timeDerivUNormal() *
(( upwindWeight_)* up.fluidState().massFraction(phaseIdx, transportCompIdx)
* up.effPorosity
+
(1 - upwindWeight_)*dn.fluidState().massFraction(phaseIdx, transportCompIdx)
* up.effPorosity);
// stabilization term
if(withStabilization_)
flux[transportEqIdx] -=
stabilizationTerm *
(( upwindWeight_)* up.fluidState().massFraction(phaseIdx, transportCompIdx)
+
(1 - upwindWeight_)*dn.fluidState().massFraction(phaseIdx, transportCompIdx));
// mass flux of the dissolved second component - massfraction
// advective flux of the component
flux[transportEqIdx] +=
fluxVars.KmvpNormal() *
(( upwindWeight_)* up.fluidState().moleFraction(phaseIdx, transportCompIdx)/up.viscosity()
+
(1 - upwindWeight_)* dn.fluidState().moleFraction(phaseIdx, transportCompIdx)/dn.viscosity());
// flux of the dissolved second component due to solid displacement
flux[transportEqIdx] +=
fluxVars.timeDerivUNormal() *
(( upwindWeight_)* up.fluidState().moleFraction(phaseIdx, transportCompIdx)
* up.effPorosity
+
(1 - upwindWeight_)*dn.fluidState().moleFraction(phaseIdx, transportCompIdx)
* up.effPorosity);
// stabilization term
if(withStabilization_)
flux[transportEqIdx] -=
stabilizationTerm *
(( upwindWeight_)* up.fluidState().moleFraction(phaseIdx, transportCompIdx)
+
(1 - upwindWeight_)*dn.fluidState().moleFraction(phaseIdx, transportCompIdx));
}
else
{
// mass flux of the dissolved second component - massfraction
// advective flux of the component
flux[transportEqIdx] +=
fluxVars.KmvpNormal() *
(( upwindWeight_)* up.fluidState().moleFraction(phaseIdx, transportCompIdx)/up.viscosity()
+
(1 - upwindWeight_)* dn.fluidState().moleFraction(phaseIdx, transportCompIdx)/dn.viscosity());
// flux of the dissolved second component due to solid displacement
flux[transportEqIdx] +=
fluxVars.timeDerivUNormal() *
(( upwindWeight_)* up.fluidState().moleFraction(phaseIdx, transportCompIdx)
* up.effPorosity
+
(1 - upwindWeight_)*dn.fluidState().moleFraction(phaseIdx, transportCompIdx)
* up.effPorosity);
// stabilization term
if(withStabilization_)
flux[transportEqIdx] -=
stabilizationTerm *
(( upwindWeight_)* up.fluidState().moleFraction(phaseIdx, transportCompIdx)
+
(1 - upwindWeight_)*dn.fluidState().moleFraction(phaseIdx, transportCompIdx));
// mass flux of the dissolved second component - massfraction
// advective flux of the component
flux[transportEqIdx] +=
fluxVars.KmvpNormal() *
(( upwindWeight_)* up.fluidState().massFraction(phaseIdx, transportCompIdx)/up.viscosity()
+
(1 - upwindWeight_)* dn.fluidState().massFraction(phaseIdx, transportCompIdx)/dn.viscosity());
// flux of the dissolved second component due to solid displacement
flux[transportEqIdx] +=
fluxVars.timeDerivUNormal() *
(( upwindWeight_)* up.fluidState().massFraction(phaseIdx, transportCompIdx)
* up.effPorosity
+
(1 - upwindWeight_)*dn.fluidState().massFraction(phaseIdx, transportCompIdx)
* up.effPorosity);
// stabilization term
if(withStabilization_)
flux[transportEqIdx] -=
stabilizationTerm *
(( upwindWeight_)* up.fluidState().massFraction(phaseIdx, transportCompIdx)
+
(1 - upwindWeight_)*dn.fluidState().massFraction(phaseIdx, transportCompIdx));
}
}
/*!
......@@ -273,32 +272,32 @@ namespace Dumux
Scalar tmp(0);
// diffusive flux of second component
if(!useMoles)
if(useMoles)
{
// diffusive flux of the second component - massfraction
// diffusive flux of the second component - mole fraction
tmp = -(fluxVars.moleFractionGrad(transportCompIdx)*fluxVars.face().normal);
tmp *= fluxVars.diffCoeffPM();
// dispersive flux of second component - massfraction
DimVector normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
tmp -= (normalDisp * fluxVars.moleFractionGrad(transportCompIdx));
// dispersive flux of second component - mole fraction
DimVector normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
tmp -= (normalDisp * fluxVars.moleFractionGrad(transportCompIdx));
// convert it to a mass flux and add it
flux[transportEqIdx] += tmp * FluidSystem::molarMass(transportCompIdx);
flux[transportEqIdx] += tmp;
}
else
{
// diffusive flux of the second component - molefraction
// diffusive flux of the second component - mass fraction
tmp = -(fluxVars.moleFractionGrad(transportCompIdx)*fluxVars.face().normal);
tmp *= fluxVars.diffCoeffPM();
// dispersive flux of second component - molefraction
DimVector normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
tmp -= (normalDisp * fluxVars.moleFractionGrad(transportCompIdx));
// dispersive flux of second component - mass fraction
DimVector normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
tmp -= (normalDisp * fluxVars.moleFractionGrad(transportCompIdx));
flux[transportEqIdx] += tmp;
// convert it to a mass flux and add it
flux[transportEqIdx] += tmp * FluidSystem::molarMass(transportCompIdx);
}
}
......
......@@ -107,16 +107,7 @@ public:
const VolumeVariables &volVars = elemVolVars[scvIdx];
storage = 0;
if(!useMoles) //mass-fraction formulation
{
// storage term of continuity equation - massfractions
storage[conti0EqIdx] +=
volVars.density()*volVars.porosity();
//storage term of the transport equation - massfractions
storage[transportEqIdx] +=
volVars.density() * volVars.massFraction(transportCompIdx) * volVars.porosity();
}
else //mole-fraction formulation
if(useMoles) // mole-fraction formulation
{
// storage term of continuity equation- molefractions
//careful: molarDensity changes with moleFrac!
......@@ -126,7 +117,15 @@ public:
volVars.molarDensity()*volVars.moleFraction(transportCompIdx) *
volVars.porosity();
}
else // mass-fraction formulation
{
// storage term of continuity equation - massfractions
storage[conti0EqIdx] +=
volVars.density()*volVars.porosity();
//storage term of the transport equation - massfractions
storage[transportEqIdx] +=
volVars.density() * volVars.massFraction(transportCompIdx) * volVars.porosity();
}
}
/*!
......@@ -221,7 +220,21 @@ public:
Scalar tmp(0);
// diffusive flux of second component
if(!useMoles) //mass-fraction formulation
if(useMoles) // mole-fraction formulation
{
// diffusive flux of the second component - molefraction
tmp = -(fluxVars.moleFractionGrad(transportCompIdx)*fluxVars.face().normal);
tmp *= fluxVars.porousDiffCoeff() * fluxVars.molarDensity();
// dispersive flux of second component - molefraction
GlobalPosition normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
tmp -= fluxVars.molarDensity()*
(normalDisp * fluxVars.moleFractionGrad(transportCompIdx));
flux[transportEqIdx] += tmp;
}
else // mass-fraction formulation
{
// diffusive flux of the second component - massfraction
tmp = -(fluxVars.moleFractionGrad(transportCompIdx)*fluxVars.face().normal);
......@@ -236,20 +249,6 @@ public:
// convert it to a mass flux and add it
flux[transportEqIdx] += tmp * FluidSystem::molarMass(transportCompIdx);
}
else //mole-fraction formulation
{
// diffusive flux of the second component - molefraction
tmp = -(fluxVars.moleFractionGrad(transportCompIdx)*fluxVars.face().normal);
tmp *= fluxVars.porousDiffCoeff() * fluxVars.molarDensity();
// dispersive flux of second component - molefraction
GlobalPosition normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
tmp -= fluxVars.molarDensity()*
(normalDisp * fluxVars.moleFractionGrad(transportCompIdx));
flux[transportEqIdx] += tmp;
}
}
/*!
......
......@@ -132,7 +132,7 @@ public:
fluidState.setPressure(phaseIdx, priVars[pressureIdx]);
Scalar x1 = priVars[massOrMoleFracIdx]; //mole or mass fraction of component 1
if(!useMoles) //mass-fraction formulation
if(!useMoles) // mass-fraction formulation
{
// convert mass to mole fractions
Scalar M0 = FluidSystem::molarMass(phaseCompIdx);
......
......@@ -143,45 +143,45 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
// compute storage term of all components within all phases
storage = 0;
if(!useMoles) //mass fraction formulation
if(useMoles) // mole-fraction formulation
{
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
storage[eqIdx] += volVars.density(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.massFraction(phaseIdx, compIdx);
}
// this is only processed if one component mass balance equation
// is replaced by the total mass balance equation
if (replaceCompEqIdx < numComponents)
storage[replaceCompEqIdx] +=
volVars.density(phaseIdx)
* volVars.saturation(phaseIdx);
}
storage *= volVars.porosity();
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
storage[eqIdx] += volVars.molarDensity(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.moleFraction(phaseIdx, compIdx);
}
// this is only processed if one component mass balance equation
// is replaced by the total mass balance equation
if (replaceCompEqIdx < numComponents)
storage[replaceCompEqIdx] +=
volVars.molarDensity(phaseIdx)
* volVars.saturation(phaseIdx);
}
storage *= volVars.porosity();
}
else //mole-fraction formulation
else // mass-fraction formulation