Commit e026b920 authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[stokes]

fixed bug when massUpwindWeight < 1 is used in stokes

reviewed by bernd


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14916 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 49e72151
......@@ -190,21 +190,13 @@ protected:
if(!useMoles)
{
// mass balance with upwinded density
if (Dune::FloatCmp::eq<Scalar>(massUpwindWeight_, 1.0)) // fully upwind
massBalanceResidual *= up.density();
else
massBalanceResidual *= massUpwindWeight_ * up.density() +
(1.-massUpwindWeight_) * dn.density();
massBalanceResidual *= (massUpwindWeight_ * up.density()
+ (1.-massUpwindWeight_) * dn.density());
}
else
{
// mass balance with upwinded molarDensity
if (Dune::FloatCmp::eq<Scalar>(massUpwindWeight_, 1.0)) // fully upwind
massBalanceResidual *= up.molarDensity();
else
massBalanceResidual *= massUpwindWeight_ * up.molarDensity() +
(1.-massUpwindWeight_) * dn.molarDensity();
massBalanceResidual *= (massUpwindWeight_ * up.molarDensity()
+ (1.-massUpwindWeight_) * dn.molarDensity());
}
if (!fluxVars.onBoundary())
......
......@@ -163,23 +163,15 @@ public:
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
Scalar tmp = fluxVars.normalVelocity();
Scalar tmp = fluxVars.normalVelocity();
if(!useMoles)
{
// for transport equations
if (this->massUpwindWeight_ > 0.0)
tmp *= this->massUpwindWeight_ // upwind data
* up.density()
* up.massFraction(transportCompIdx);
if (this->massUpwindWeight_ < 1.0)
tmp += (1.0 - this->massUpwindWeight_) // rest
* dn.density()
* dn.massFraction(transportCompIdx);
flux[transportEqIdx] += tmp;
Valgrind::CheckDefined(flux[transportEqIdx]);
tmp *= (this->massUpwindWeight_ * up.density() * up.massFraction(transportCompIdx)
+ (1.-this->massUpwindWeight_) * dn.density() * dn.massFraction(transportCompIdx));
flux[transportEqIdx] += tmp;
Valgrind::CheckDefined(flux[transportEqIdx]);
}
else
{
......@@ -188,14 +180,8 @@ public:
{
if (conti0EqIdx+compIdx != massBalanceIdx) //mass balance is calculated above
{
if (this->massUpwindWeight_ > 0.0)
tmp *= this->massUpwindWeight_ // upwind data
* up.molarDensity()
* up.moleFraction(compIdx);
if (this->massUpwindWeight_ < 1.0)
tmp += (1.0 - this->massUpwindWeight_) // rest
* dn.molarDensity()
* dn.moleFraction(compIdx);
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]);
......
......@@ -109,11 +109,8 @@ public:
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
Scalar tmp = fluxVars.normalVelocity();
tmp *= this->massUpwindWeight_ * // upwind data
up.density() * up.enthalpy() +
(1.0 - this->massUpwindWeight_) * // rest
dn.density() * dn.enthalpy();
tmp *= (this->massUpwindWeight_ * up.density() * up.enthalpy()
+ (1.0 - this->massUpwindWeight_) * dn.density() * dn.enthalpy());
flux[energyEqIdx] += tmp;
Valgrind::CheckDefined(flux[energyEqIdx]);
......
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment