diff --git a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
index 0ad7c66f1b2c423c7f2fbce9e1cfb60eb61a2279..10cfd0ec825d2589a816d854ae1d255ff74eaee1 100644
--- a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
@@ -50,8 +50,9 @@ namespace Dumux
* the integration point, etc.
*/
template <class TypeTag>
-class TwoPTwoCFluxVariables
+class TwoPTwoCFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariables)
{
+ typedef typename GET_PROP_TYPE(TypeTag, BaseFluxVariables) BaseFluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
@@ -84,7 +85,7 @@ public:
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the box scheme
- * \param faceIdx The local index of the SCV (sub-control-volume) face
+ * \param faceIdx The local index of the SCV (sub-control-volume) face
* \param elemVolVars The volume variables of the current element
* \param onBoundary Distinguishes if we are on a SCV face or on a boundary face
*/
@@ -94,12 +95,11 @@ public:
const int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
- : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
+ : BaseFluxVariables(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
{
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
density_[phaseIdx] = Scalar(0);
molarDensity_[phaseIdx] = Scalar(0);
- potentialGrad_[phaseIdx] = Scalar(0);
massFractionGrad_[phaseIdx] = Scalar(0);
moleFractionGrad_[phaseIdx] = Scalar(0);
}
@@ -115,23 +115,22 @@ protected:
// calculate densities at the integration points of the face
DimVector tmp(0.0);
for (int idx = 0;
- idx < fvGeometry_.numFAP;
+ idx < this->fvGeometry_.numFAP;
idx++) // loop over adjacent vertices
{
// index for the element volume variables
- int volVarsIdx = face().fapIndices[idx];
+ int volVarsIdx = this->face().fapIndices[idx];
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
{
density_[phaseIdx] += elemVolVars[volVarsIdx].density(phaseIdx)*
- face().shapeValue[idx];
+ this->face().shapeValue[idx];
molarDensity_[phaseIdx] += elemVolVars[volVarsIdx].molarDensity(phaseIdx)*
- face().shapeValue[idx];
+ this->face().shapeValue[idx];
}
}
calculateGradients_(problem, element, elemVolVars);
- calculateVelocities_(problem, element, elemVolVars);
calculatePorousDiffCoeff_(problem, element, elemVolVars);
}
@@ -142,23 +141,14 @@ protected:
// calculate gradients
DimVector tmp(0.0);
for (int idx = 0;
- idx < fvGeometry_.numFAP;
+ idx < this->fvGeometry_.numFAP;
idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
- const DimVector &feGrad = face().grad[idx];
+ const DimVector &feGrad = this->face().grad[idx];
// index for the element volume variables
- int volVarsIdx = face().fapIndices[idx];
-
- // compute sum of pressure gradients for each phase
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- // the pressure gradient
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].pressure(phaseIdx);
- potentialGrad_[phaseIdx] += tmp;
- }
+ int volVarsIdx = this->face().fapIndices[idx];
// the concentration gradient of the non-wetting
// component in the wetting phase
@@ -180,41 +170,6 @@ protected:
tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
moleFractionGrad_[nPhaseIdx] += tmp;
}
-
- ///////////////
- // correct the pressure gradients by the gravitational acceleration
- ///////////////
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity)) {
- // estimate the gravitational acceleration at a given SCV face
- // using the arithmetic mean
- DimVector g(problem.boxGravity(element, fvGeometry_, face().i));
- g += problem.boxGravity(element, fvGeometry_, face().j);
- g /= 2;
-
- for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
- {
- // calculate the phase density at the integration point. we
- // only do this if the wetting phase is present in both cells
- Scalar SI = elemVolVars[face().i].saturation(phaseIdx);
- Scalar SJ = elemVolVars[face().j].saturation(phaseIdx);
- Scalar rhoI = elemVolVars[face().i].density(phaseIdx);
- Scalar rhoJ = elemVolVars[face().j].density(phaseIdx);
- Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
- Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
- if (fI + fJ == 0)
- // doesn't matter because no wetting phase is present in
- // both cells!
- fI = fJ = 0.5;
- Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
-
- // make gravity acceleration a force
- DimVector f(g);
- f *= density;
-
- // calculate the final potential gradient
- potentialGrad_[phaseIdx] -= f;
- }
- }
}
Scalar rhoFactor_(int phaseIdx, int scvIdx, const ElementVolumeVariables &vDat)
@@ -232,45 +187,12 @@ protected:
return sp.eval(sat);
}
- void calculateVelocities_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- const SpatialParams &spatialParams = problem.spatialParams();
- // multiply the pressure potential by the intrinsic permeability
- DimMatrix K;
- for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
- {
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- fvGeometry_,
- face().i),
- spatialParams.intrinsicPermeability(element,
- fvGeometry_,
- face().j));
- K.mv(potentialGrad_[phaseIdx], Kmvp_[phaseIdx]);
- KmvpNormal_[phaseIdx] = -(Kmvp_[phaseIdx]*face().normal);
- }
-
- // set the upstream and downstream vertices
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- upstreamIdx_[phaseIdx] = face().i;
- downstreamIdx_[phaseIdx] = face().j;
-
- if (KmvpNormal_[phaseIdx] < 0) {
- std::swap(upstreamIdx_[phaseIdx],
- downstreamIdx_[phaseIdx]);
- }
- }
- }
-
void calculatePorousDiffCoeff_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
- const VolumeVariables &volVarsI = elemVolVars[face().i];
- const VolumeVariables &volVarsJ = elemVolVars[face().j];
+ const VolumeVariables &volVarsI = elemVolVars[this->face().i];
+ const VolumeVariables &volVarsJ = elemVolVars[this->face().j];
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
@@ -312,29 +234,17 @@ public:
* decision on the upwinding approach (which is done in the
* actual model).
*/
+ DUNE_DEPRECATED_MSG("use the methods of the base flux variables instead")
Scalar KmvpNormal(int phaseIdx) const
- { return KmvpNormal_[phaseIdx]; }
+ { return -this->kGradPNormal(phaseIdx); }
/*!
* \brief Return the pressure potential multiplied with the
* intrinsic permeability as vector (for velocity output)
*/
+ DUNE_DEPRECATED_MSG("use the methods of the base flux variables instead")
DimVector Kmvp(int phaseIdx) const
- { return Kmvp_[phaseIdx]; }
-
- /*!
- * \brief Return the local index of the upstream control volume
- * for a given phase.
- */
- int upstreamIdx(int phaseIdx) const
- { return upstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the local index of the downstream control volume
- * for a given phase.
- */
- int downstreamIdx(int phaseIdx) const
- { return downstreamIdx_[phaseIdx]; }
+ { return this->kGradP_[phaseIdx]; }
/*!
* \brief The binary diffusion coefficient for each fluid phase.
@@ -397,25 +307,8 @@ public:
const DimVector &moleFractionGrad(int phaseIdx) const
{ return moleFractionGrad_[phaseIdx]; };
- /*!
- * \brief The face of the current sub-control volume. This may be either
- * an inner sub-control-volume face or a face on the boundary.
- */
- const SCVFace &face() const
- {
- if (onBoundary_)
- return fvGeometry_.boundaryFace[faceIdx_];
- else
- return fvGeometry_.subContVolFace[faceIdx_];
- }
-
protected:
- const FVElementGeometry &fvGeometry_;
- const int faceIdx_;
- const bool onBoundary_;
-
// gradients
- DimVector potentialGrad_[numPhases];
DimVector massFractionGrad_[numPhases];
DimVector moleFractionGrad_[numPhases];
@@ -424,13 +317,6 @@ protected:
// intrinsic permeability times pressure potential gradient
DimVector Kmvp_[numPhases];
- // projected on the face normal
- Scalar KmvpNormal_[numPhases];
-
- // local index of the upwind vertex for each phase
- int upstreamIdx_[numPhases];
- // local index of the downwind vertex for each phase
- int downstreamIdx_[numPhases];
// the diffusion coefficient for the porous medium
Scalar porousDiffCoeff_[numPhases];
diff --git a/dumux/boxmodels/2p2c/2p2clocalresidual.hh b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
index d0a9ab42ef556f698c91b291f32df25f9e2725d1..8d52b3390236afe511a28cd1057a2e3f9064624f 100644
--- a/dumux/boxmodels/2p2c/2p2clocalresidual.hh
+++ b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
@@ -229,26 +229,22 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
if (massUpwindWeight_ > 0.0)
// upstream vertex
flux[eqIdx] +=
- fluxVars.KmvpNormal(phaseIdx)
+ fluxVars.volumeFlux(phaseIdx)
* massUpwindWeight_
* up.density(phaseIdx)
- * up.mobility(phaseIdx)
* up.fluidState().massFraction(phaseIdx, compIdx);
if (massUpwindWeight_ < 1.0)
// downstream vertex
flux[eqIdx] +=
- fluxVars.KmvpNormal(phaseIdx)
+ fluxVars.volumeFlux(phaseIdx)
* (1 - massUpwindWeight_)
* dn.density(phaseIdx)
- * dn.mobility(phaseIdx)
* dn.fluidState().massFraction(phaseIdx, compIdx);
- Valgrind::CheckDefined(fluxVars.KmvpNormal(phaseIdx));
+ Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx));
Valgrind::CheckDefined(up.density(phaseIdx));
- Valgrind::CheckDefined(up.mobility(phaseIdx));
Valgrind::CheckDefined(up.fluidState().massFraction(phaseIdx, compIdx));
Valgrind::CheckDefined(dn.density(phaseIdx));
- Valgrind::CheckDefined(dn.mobility(phaseIdx));
Valgrind::CheckDefined(dn.fluidState().massFraction(phaseIdx, compIdx));
}
// flux of the total mass balance;
@@ -259,22 +255,18 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
// upstream vertex
if (massUpwindWeight_ > 0.0)
flux[replaceCompEqIdx] +=
- fluxVars.KmvpNormal(phaseIdx)
+ fluxVars.volumeFlux(phaseIdx)
* massUpwindWeight_
- * up.density(phaseIdx)
- * up.mobility(phaseIdx);
+ * up.density(phaseIdx);
// downstream vertex
if (massUpwindWeight_ < 1.0)
flux[replaceCompEqIdx] +=
- fluxVars.KmvpNormal(phaseIdx)
+ fluxVars.volumeFlux(phaseIdx)
* (1 - massUpwindWeight_)
- * dn.density(phaseIdx)
- * dn.mobility(phaseIdx);
- Valgrind::CheckDefined(fluxVars.KmvpNormal(phaseIdx));
+ * dn.density(phaseIdx);
+ Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx));
Valgrind::CheckDefined(up.density(phaseIdx));
- Valgrind::CheckDefined(up.mobility(phaseIdx));
Valgrind::CheckDefined(dn.density(phaseIdx));
- Valgrind::CheckDefined(dn.mobility(phaseIdx));
}
diff --git a/dumux/boxmodels/2p2c/2p2cmodel.hh b/dumux/boxmodels/2p2c/2p2cmodel.hh
index 58138be58c2f9cc712916131b366504f23afe150..ca230d8694177cbd134d1dd60ff931a26a079d24 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -181,8 +181,6 @@ public:
staticVertexDat_[globalIdx].oldPhasePresence
= staticVertexDat_[globalIdx].phasePresence;
}
-
- massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
}
/*!
@@ -432,13 +430,8 @@ public:
// Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolume
// face in the reference element.
- massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
PhasesVector flux;
- flux[phaseIdx] = fluxVars.KmvpNormal(phaseIdx)
- * (massUpwindWeight_
- * up.mobility(phaseIdx)
- + (1- massUpwindWeight_)
- * dn.mobility(phaseIdx)) / localArea;
+ flux[phaseIdx] = fluxVars.volumeFlux(phaseIdx) / localArea;
// transform the normal Darcy velocity into a vector
tmpVelocity[phaseIdx] = localNormal;
@@ -781,7 +774,6 @@ protected:
std::vector<StaticVars> staticVertexDat_;
bool switchFlag_;
bool velocityOutput_;
- Scalar massUpwindWeight_;
};
}
diff --git a/dumux/boxmodels/2p2c/2p2cproperties.hh b/dumux/boxmodels/2p2c/2p2cproperties.hh
index 599ac05a4986d6e98eb7ba85fab9c44107fdcc4b..61e0c8bef18692ac216ade943484e9b04f1f82e6 100644
--- a/dumux/boxmodels/2p2c/2p2cproperties.hh
+++ b/dumux/boxmodels/2p2c/2p2cproperties.hh
@@ -69,9 +69,11 @@ NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered i
NEW_PROP_TAG(EnableGravity); //!< DEPRECATED Returns whether gravity is considered in the problem
NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind weight for the mass conservation equations
NEW_PROP_TAG(MassUpwindWeight); //!< DEPRECATED The value of the upwind weight for the mass conservation equations
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
NEW_PROP_TAG(ReplaceCompEqIdx); //!< The index of the total mass balance equation, if one component balance is replaced (ReplaceCompEqIdx < NumComponents)
NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
NEW_PROP_TAG(EnableVelocityOutput); //!< DEPRECATED Returns whether vertex velocity vectors are written into the vtk output
+NEW_PROP_TAG(BaseFluxVariables); //! The base flux variables
}
}
diff --git a/dumux/boxmodels/2p2c/2p2cpropertydefaults.hh b/dumux/boxmodels/2p2c/2p2cpropertydefaults.hh
index 7067782dbbe6b2b141fda84c68003dcc80f8ce25..692387a12e4bbdb5669c35a764edfe2024fd512d 100644
--- a/dumux/boxmodels/2p2c/2p2cpropertydefaults.hh
+++ b/dumux/boxmodels/2p2c/2p2cpropertydefaults.hh
@@ -39,6 +39,8 @@
#include "2p2cproperties.hh"
#include "2p2cnewtoncontroller.hh"
+#include <dumux/boxmodels/common/boxdarcyfluxvariables.hh>
+
namespace Dumux
{
@@ -122,10 +124,16 @@ SET_TYPE_PROP(BoxTwoPTwoC, VolumeVariables, TwoPTwoCVolumeVariables<TypeTag>);
//! the FluxVariables property
SET_TYPE_PROP(BoxTwoPTwoC, FluxVariables, TwoPTwoCFluxVariables<TypeTag>);
+//! define the base flux variables to realize Darcy flow
+SET_TYPE_PROP(BoxTwoPTwoC, BaseFluxVariables, BoxDarcyFluxVariables<TypeTag>);
+
//! the upwind weight for the mass conservation equations.
SET_SCALAR_PROP(BoxTwoPTwoC, ImplicitMassUpwindWeight, GET_PROP_VALUE(TypeTag, MassUpwindWeight));
SET_SCALAR_PROP(BoxTwoPTwoC, MassUpwindWeight, 1.0);//DEPRECATED
+//! set default mobility upwind weight to 1.0, i.e. fully upwind
+SET_SCALAR_PROP(BoxTwoPTwoC, ImplicitMobilityUpwindWeight, 1.0);
+
//! The indices required by the isothermal 2p2c model
SET_PROP(BoxTwoPTwoC,
Indices)
diff --git a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
index 9686c77e10872b7b49e21afd2377831d8e72bd8f..f241d5b48fa80f99bfd789409dcf6168dbdf0780 100644
--- a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
@@ -141,14 +141,12 @@ public:
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
flux[energyEqIdx] +=
- fluxVars.KmvpNormal(phaseIdx) * (massUpwindWeight_ * // upstream vertex
+ fluxVars.volumeFlux(phaseIdx) * (massUpwindWeight_ * // upstream vertex
( up.density(phaseIdx) *
- up.mobility(phaseIdx) *
up.enthalpy(phaseIdx))
+
(1-massUpwindWeight_) * // downstream vertex
( dn.density(phaseIdx) *
- dn.mobility(phaseIdx) *
dn.enthalpy(phaseIdx)) );
}
}