diff --git a/dumux/boxmodels/common/boxproblem.hh b/dumux/boxmodels/common/boxproblem.hh
index c455443807b5b990d59e43d0bad635b32fd87f98..3a52ac48383aa0955cf9f9ef101f55d31cc9738b 100644
--- a/dumux/boxmodels/common/boxproblem.hh
+++ b/dumux/boxmodels/common/boxproblem.hh
@@ -217,7 +217,7 @@ public:
* potentially solution dependent and requires some box method
* specific things.
*
- * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
+ * \param values The neumann values for the conservation equations in units of \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param is The intersection between element and boundary
@@ -249,7 +249,7 @@ public:
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
- * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
+ * \param values The neumann values for the conservation equations in units of \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param is The intersection between element and boundary
@@ -274,7 +274,7 @@ public:
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
- * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
+ * \param values The neumann values for the conservation equations in units of \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
* \param pos The position of the boundary face's integration point in global coordinates
*
* For this method, the \a values parameter stores the mass flux
@@ -299,7 +299,7 @@ public:
* potentially solution dependent and requires some box method
* specific things.
*
- * \param values The source and sink values for the conservation equations
+ * \param values The source and sink values for the conservation equations in units of \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
@@ -323,7 +323,7 @@ public:
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
- * \param values The source and sink values for the conservation equations
+ * \param values The source and sink values for the conservation equations in units of \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
@@ -345,7 +345,7 @@ public:
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
- * \param values The source and sink values for the conservation equations
+ * \param values The source and sink values for the conservation equations in units of \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$
* \param pos The position of the center of the finite volume
* for which the source term ought to be
* specified in global coordinates
diff --git a/dumux/boxmodels/mpnc/mpncindices.hh b/dumux/boxmodels/mpnc/mpncindices.hh
index eb1cc5b0736d3b741ed03c48ad5a955f0226df30..2e80219300592faf203037f813bae460ba5902b4 100644
--- a/dumux/boxmodels/mpnc/mpncindices.hh
+++ b/dumux/boxmodels/mpnc/mpncindices.hh
@@ -26,6 +26,20 @@
namespace Dumux
{
+
+/*!
+ * \ingroup MPNCModel
+ * \ingroup BoxIndices
+ * \brief Enumerates the formulations which the 2p2c model accepts.
+ */
+struct MpNcPressureFormulation
+{
+ enum {
+ mostWettingFirst,
+ leastWettingFirst
+ };
+};
+
/*!
* \ingroup MPNCModel
* \ingroup BoxIndices
diff --git a/dumux/boxmodels/mpnc/mpncproperties.hh b/dumux/boxmodels/mpnc/mpncproperties.hh
index a8833bb7cd2051e7e4556a5584f23847c6930d6e..3a4dc02af599115d5ab38dee6e963db8983c71c0 100644
--- a/dumux/boxmodels/mpnc/mpncproperties.hh
+++ b/dumux/boxmodels/mpnc/mpncproperties.hh
@@ -51,6 +51,8 @@ NEW_PROP_TAG(Indices); //!< Enumerations for the model
NEW_PROP_TAG(BaseFluxVariables); //!< The type of velocity calculation that is to be used
+NEW_PROP_TAG(PressureFormulation); //!< The formulation of the model
+
NEW_PROP_TAG(MPNCVtkCommonModule); //!< Vtk writer module for writing the common quantities into the VTK output file
NEW_PROP_TAG(MPNCVtkMassModule); //!< Vtk writer module for writing the mass related quantities into the VTK output file
NEW_PROP_TAG(MPNCVtkEnergyModule); //!< Vtk writer module for writing the energy related quantities into the VTK output file
@@ -60,6 +62,7 @@ NEW_PROP_TAG(VelocityAveragingInModel);//!< Should the averaging of velocities b
//! specify which quantities are written to the vtk output files
NEW_PROP_TAG(VtkAddPorosity);
+NEW_PROP_TAG(VtkAddPermeability);
NEW_PROP_TAG(VtkAddBoundaryTypes);
NEW_PROP_TAG(VtkAddSaturations);
NEW_PROP_TAG(VtkAddPressures);
diff --git a/dumux/boxmodels/mpnc/mpncpropertydefaults.hh b/dumux/boxmodels/mpnc/mpncpropertydefaults.hh
index 0ead27ef49ade35d95b3420fa2eb967a236a0744..578332fc70d9769d438e681d18e0a194cebd152d 100644
--- a/dumux/boxmodels/mpnc/mpncpropertydefaults.hh
+++ b/dumux/boxmodels/mpnc/mpncpropertydefaults.hh
@@ -52,6 +52,12 @@ namespace Properties
// default property values
//////////////////////////////////////////////////////////////////
+
+//! Set the default pressure formulation to the pressure of the (most) wetting phase
+SET_INT_PROP(BoxMPNC,
+ PressureFormulation,
+ MpNcPressureFormulation::mostWettingFirst);
+
/*!
* \brief Set the property for the number of components.
*
@@ -247,6 +253,7 @@ SET_BOOL_PROP(BoxMPNC, VelocityAveragingInModel, false);
//! Specify what to add to the VTK output by default
SET_BOOL_PROP(BoxMPNC, VtkAddPorosity, true);
+SET_BOOL_PROP(BoxMPNC, VtkAddPermeability, false);
SET_BOOL_PROP(BoxMPNC, VtkAddBoundaryTypes, false);
SET_BOOL_PROP(BoxMPNC, VtkAddSaturations, true);
SET_BOOL_PROP(BoxMPNC, VtkAddPressures, true);
diff --git a/dumux/boxmodels/mpnc/mpncvolumevariables.hh b/dumux/boxmodels/mpnc/mpncvolumevariables.hh
index 409218898f666c9855d2deb26067b5cec5b21b0b..bba2e1b0f7550cfbb5354d0dcbcd35352bce1269 100644
--- a/dumux/boxmodels/mpnc/mpncvolumevariables.hh
+++ b/dumux/boxmodels/mpnc/mpncvolumevariables.hh
@@ -64,6 +64,19 @@ class MPNCVolumeVariables
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum {dimWorld=GridView::dimensionworld};
+ typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition;
+ enum {wPhaseIdx = FluidSystem::wPhaseIdx};
+ enum {nPhaseIdx = FluidSystem::nPhaseIdx};
+
+ // formulations
+ enum {
+ pressureFormulation = GET_PROP_VALUE(TypeTag, PressureFormulation),
+ mostWettingFirst = MpNcPressureFormulation::mostWettingFirst,
+ leastWettingFirst = MpNcPressureFormulation::leastWettingFirst
+ };
+
enum {numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
enum {numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
enum {enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
@@ -138,7 +151,6 @@ public:
/////////////
// set the phase pressures
/////////////
-
// capillary pressure parameters
const MaterialLawParams &materialParams =
problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
@@ -146,9 +158,23 @@ public:
Scalar capPress[numPhases];
MaterialLaw::capillaryPressures(capPress, materialParams, fluidState_);
// add to the pressure of the first fluid phase
- Scalar p0 = priVars[p0Idx];
- for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
- fluidState_.setPressure(phaseIdx, p0 - capPress[0] + capPress[phaseIdx]);
+
+ // depending on which pressure is stored in the primary variables
+ if(pressureFormulation == mostWettingFirst){
+ // This means that the pressures are sorted from the most wetting to the least wetting-1 in the primary variables vector.
+ // For two phases this means that there is one pressure as primary variable: pw
+ const Scalar pw = priVars[p0Idx];
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ fluidState_.setPressure(phaseIdx, pw - capPress[0] + capPress[phaseIdx]);
+ }
+ else if(pressureFormulation == leastWettingFirst){
+ // This means that the pressures are sorted from the least wetting to the most wetting-1 in the primary variables vector.
+ // For two phases this means that there is one pressure as primary variable: pn
+ const Scalar pn = priVars[p0Idx];
+ for (int phaseIdx = numPhases-1; phaseIdx >= 0; --phaseIdx)
+ fluidState_.setPressure(phaseIdx, pn - capPress[numPhases-1] + capPress[phaseIdx]);
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << pressureFormulation << " is invalid.");
/////////////
// set the fluid compositions
@@ -240,7 +266,9 @@ public:
* the control volume.
*/
Scalar mobility(const unsigned int phaseIdx) const
- { return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx); }
+ {
+ return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx);
+ }
/*!
* \brief Returns the viscosity of a given phase within
diff --git a/dumux/boxmodels/mpnc/mpncvtkwritercommon.hh b/dumux/boxmodels/mpnc/mpncvtkwritercommon.hh
index c92536bf6dfce21f96e696ff48a13e115e038cad..6a83d7dae20e0c8e3caa9260acebfd1dda65c36c 100644
--- a/dumux/boxmodels/mpnc/mpncvtkwritercommon.hh
+++ b/dumux/boxmodels/mpnc/mpncvtkwritercommon.hh
@@ -68,6 +68,7 @@ public:
: ParentType(problem)
{
porosityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPorosity);
+ permeabilityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPermeability);
boundaryTypesOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddBoundaryTypes);
saturationOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddSaturations);
pressureOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddPressures);
@@ -89,6 +90,8 @@ public:
{
if (porosityOutput_)
this->resizeScalarBuffer_(porosity_);
+ if (permeabilityOutput_)
+ this->resizeScalarBuffer_(permeability_);
if (boundaryTypesOutput_)
this->resizeScalarBuffer_(boundaryTypes_);
@@ -127,6 +130,9 @@ public:
if (porosityOutput_) porosity_[globalIdx] = volVars.porosity();
+ // works for scalar permeability in spatialparameters
+ if (permeabilityOutput_) permeability_[globalIdx] = this->problem_.spatialParams().intrinsicPermeability(elem,fvGeometry,localVertexIdx);
+
// calculate a single value for the boundary type: use one
// bit for each equation and set it to 1 if the equation
// is used for a dirichlet condition
@@ -239,6 +245,7 @@ public:
private:
bool porosityOutput_;
+ bool permeabilityOutput_ ;
bool boundaryTypesOutput_;
bool saturationOutput_;
bool pressureOutput_;
@@ -260,6 +267,7 @@ private:
ScalarVector boxSurface_;
ScalarVector porosity_;
+ ScalarVector permeability_;
ScalarVector temperature_;
ScalarVector boundaryTypes_;
diff --git a/dumux/material/constraintsolvers/ncpflash.hh b/dumux/material/constraintsolvers/ncpflash.hh
index 3bff92de02ca84c0eea87020366d7021373e25f6..4534ac9628aa49a4475ea69a5a6ad763d3e86a4d 100644
--- a/dumux/material/constraintsolvers/ncpflash.hh
+++ b/dumux/material/constraintsolvers/ncpflash.hh
@@ -62,7 +62,7 @@ namespace Dumux {
* this also sums up to M*(N + 2).
*
* We use the following catches: Capillary pressures are taken
- * into accout expicitly, so that only the pressure of the first
+ * into account explicitly, so that only the pressure of the first
* phase is solved implicitly, also the closure condition for the
* saturations is taken into account explicitly, which means, that
* we don't need to implicitly solve for the last
diff --git a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
index 73f6dee33012a86b4239f7d9f2a2cb44fc31a9f6..e346dda92c6639bf692ac04f766120b2bcb66971 100644
--- a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
@@ -136,7 +136,7 @@ public:
static Scalar krw(const Params ¶ms, Scalar Swe)
{
return std::max(std::min(Swe,1.0),0.0);
- };
+ }
/*!
* \brief The relative permeability for the non-wetting phase.