From b01c42faf9a0b22598672424de036bae8ae42781 Mon Sep 17 00:00:00 2001
From: Markus Wolff <markus.wolff@twt-gmbh.de>
Date: Mon, 20 Aug 2012 11:40:54 +0000
Subject: [PATCH] improved doxygen documentation
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8900 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../lmethod/fvmpfal2pfaboundpressure2p.hh | 199 +++++++++++-----
.../fvmpfal2pfaboundpressure2padaptive.hh | 218 +++++++++++++-----
.../lmethod/fvmpfal2pfaboundvelocity2p.hh | 64 ++++-
.../fvmpfal2pfaboundvelocity2padaptive.hh | 60 ++++-
.../omethod/fvmpfao2pfaboundpressure2p.hh | 172 +++++++++-----
.../omethod/fvmpfao2pfaboundvelocity2p.hh | 58 ++++-
.../fvmpfa/omethod/fvmpfaopressure2p.hh | 5 +-
.../fvmpfa/omethod/fvmpfaovelocity2p.hh | 2 +-
dumux/decoupled/common/fv/fvpressure.hh | 2 +
.../common/fv/mpfa/mpfalinteractionvolume.hh | 149 +++++++++++-
.../common/fv/mpfa/mpfaointeractionvolume.hh | 180 ++++++++++++++-
test/decoupled/2p/test_mpfa2pproblem.hh | 36 ++-
12 files changed, 914 insertions(+), 231 deletions(-)
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2p.hh
index 4d5021dafc..46a799b075 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2p.hh
@@ -28,25 +28,33 @@
/**
* @file
- * @brief Finite Volume-MPFAL implementation of a two-phase pressure equation
- * @brief Remark1: only for 2-D quadrilateral grid.
- * @brief Remark2: can use UGGrid or SGrid/YaspGrid
+ * @brief Finite volume MPFA L-method discretization of a two-phase pressure equation of the sequential IMPES model.
* @author Markus Wolff
*/
namespace Dumux
{
-//! \ingroup diffusion
-/*! Finite Volume-MPFAO-Implementation of the equation
- * \f$ - \text{div}\, \mathbf_{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_w + f_n \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap} ) = 0, \f$, or
- * \f$ - \text{div}\, \mathbf_{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_n - f_w \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap} ) = 0, \f$.
- * \f$\Phi = g\f$ on \f$\Gamma_1\f$, and
- * \f$-\text{div}\, \mathbf_{v}_t \cdot \mathbf{n} = J\f$
- * on \f$\Gamma_2\f$. Here,
- * \f$Phi_\alpha \f$ denotes the potential of phase \f$\alpha\f$, \f$K\f$ the intrinsic permeability,
- * \f$\lambda_t\f$ the total mobility, \f$this->f_\alpha\f$ the phase fractional flow function.
- * More details on the equations can be found in
- * H. Hoteit, A. Firoozabadi. Numerical modeling of thwo-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Res (31), 2008
+//! \ingroup FVPressure2p
+/*! \brief Finite volume MPFA L-method discretization of a two-phase flow pressure equation of the sequential IMPES model.
+ *
+ * Finite volume MPFA L-method discretization of the equations
+ * \f[ - \text{div}\, \boldsymbol v_t = - \text{div}\, (\lambda_t \boldsymbol K \text{grad}\, \Phi_w + f_n \lambda_t \boldsymbol K \text{grad}\, \Phi_{cap} ) = 0, \f]
+ * or
+ * \f[ - \text{div}\, \boldsymbol v_t = - \text{div}\, (\lambda_t \boldsymbol K \text{grad}\, \Phi_n - f_w \lambda_t \boldsymbol K \text{grad}\, \Phi_{cap} ) = 0. \f]
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * \f[ \Phi = g \; \text{on} \; \Gamma_1, \quad \text{and} \quad
+ * -\text{div}\, \boldsymbol v_t \cdot \mathbf{n} = J \; \text{on} \; \Gamma_2. \f]
+ * Here, \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * \f$ \lambda_t \f$ the total mobility, \f$ f_\alpha \f$ the phase fractional flow function.
+ *
+ * More details on the equations can be found in H. Hoteit, A. Firoozabadi.
+ * Numerical modeling of thwo-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Res (31), 2008
+ *
+ * Remark1: only for 2-D quadrilateral grid
+ *
+ * Remark2: implemented for UGGrid, ALUGrid, or SGrid/YaspGrid
+ *
+ *\tparam TypeTag The problem Type Tag
*/
template<class TypeTag>
class FVMPFAL2PFABoundPressure2P: public FVPressure<TypeTag>
@@ -74,9 +82,12 @@ class FVMPFAL2PFABoundPressure2P: public FVPressure<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolutionType;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename SolutionTypes::ScalarSolution ScalarSolutionType;
typedef typename GET_PROP_TYPE(TypeTag, GridTypeIndices) GridTypeIndices;
@@ -125,6 +136,12 @@ class FVMPFAL2PFABoundPressure2P: public FVPressure<TypeTag>
typedef Dune::FieldVector<Scalar, dim> DimVector;
public:
+ /*! \brief Type of the interaction volume objects
+ *
+ * Type of the interaction volume objects used to store the geometric information which is needed
+ * to calculated the transmissibility matrices of one MPFA interaction volume.
+ *
+ */
typedef Dumux::FVMPFALInteractionVolume<TypeTag> InteractionVolume;
private:
@@ -145,6 +162,11 @@ public:
//constitutive functions are initialized and stored in the variables object
void updateMaterialLaws();
+ /*! \brief Updates interaction volumes
+ *
+ * Globally rebuilds the MPFA interaction volumes.
+ *
+ */
void updateInteractionVolumeInfo()
{
interactionVolumes_.clear();
@@ -156,7 +178,11 @@ public:
storeInteractionVolumeInfo();
}
- void initialize(bool solveTwice = true)
+ /*! \brief Initializes the pressure model
+ *
+ * \copydetails ParentType::initialize()
+ */
+ void initialize()
{
ParentType::initialize();
@@ -175,6 +201,9 @@ public:
return;
}
+ /*! \brief Globally stores the pressure solution
+ *
+ */
void storePressureSolution()
{
// iterate through leaf grid an evaluate c0 at cell center
@@ -185,6 +214,10 @@ public:
}
}
+ /*! \brief Stores the pressure solution of a cell
+ *
+ * \param element Dune grid element
+ */
void storePressureSolution(const Element& element)
{
int globalIdx = problem_.variables().index(element);
@@ -219,7 +252,11 @@ public:
cellData.fluxData().resetVelocity();
}
- //! updates the pressure field (analog to update function in Dumux::IMPET)
+ /*! \brief Pressure update
+ *
+ * \copydetails ParentType::update()
+ *
+ */
void update()
{
//error bounds for error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
@@ -256,40 +293,15 @@ public:
return;
}
- Scalar evaluateErrorTerm(CellData& cellData)
- {
- //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
- // error reduction routine: volumetric error is damped and inserted to right hand side
- Scalar sat = 0;
- switch (saturationType_)
- {
- case Sw:
- sat = cellData.saturation(wPhaseIdx);
- break;
- case Sn:
- sat = cellData.saturation(nPhaseIdx);
- break;
- }
- Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
- if (sat < 0.0)
- {
- error = sat;
- }
- error /= timeStep_;
-
- Scalar errorAbs = std::abs(error);
-
- if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
- && (!problem_.timeManager().willBeFinished()))
- {
- return ErrorTermFactor_ * error;
- }
- return 0.0;
- }
-
- //! \brief Write data files
- /* \param name file name */
+ /*! \brief Adds pressure output to the output file
+ *
+ * Adds the pressure, the potential and the capillary pressure to the output.
+ *
+ * \tparam MultiWriter Class defining the output writer
+ * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+ *
+ */
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
@@ -335,7 +347,10 @@ public:
return;
}
-
+ //! Constructs a FVMPFAL2PFABoundPressure2P object
+ /**
+ * \param problem A problem class object
+ */
FVMPFAL2PFABoundPressure2P(Problem& problem) :
ParentType(problem), problem_(problem), R_(0),
gravity_(problem.gravity()),
@@ -387,14 +402,16 @@ private:
DimMatrix R_;
protected:
+
+ // Calculates tranmissibility matrix
bool calculateTransmissibility(
Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
InteractionVolume& interactionVolume,
std::vector<DimVector >& lambda,
int idx1, int idx2, int idx3, int idx4);
- GlobalInteractionVolumeVector interactionVolumes_;
- InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_;
+ GlobalInteractionVolumeVector interactionVolumes_;//!< Global Vector of interaction volumes
+ InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_; //!< Vector marking faces which intersect the boundary
private:
const GlobalPosition& gravity_; //!< vector including the gravity constant
@@ -413,6 +430,50 @@ private:
static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation); //!< gives kind of saturation used (\f$ 0 = S_w\f$, \f$ 1 = S_n\f$)
static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation); //!< gives kind of velocity used (\f$ 0 = v_w\f$, \f$ 1 = v_n\f$, \f$ 2 = v_t\f$)
+ /* Volume correction term to correct for unphysical saturation overshoots/undershoots.
+ * These can occur if the estimated time step for the explicit transport was too large. Correction by an artificial source term allows to correct
+ * this errors due to wrong time-stepping without losing mass conservation. The error term looks as follows:
+ * \f[
+ * q_{error} = \begin{cases}
+ * S < 0 & a_{error} \frac{S}{\Delta t} V \\
+ * S > 1 & a_{error} \frac{(S - 1)}{\Delta t} V \\
+ * 0 \le S \le 1 & 0
+ * \end{cases}
+ * \f]
+ * where \f$a_{error}\f$ is a weighting factor (default: \f$a_{error} = 0.5\f$)
+ */
+ Scalar evaluateErrorTerm_(CellData& cellData)
+ {
+ //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
+ // error reduction routine: volumetric error is damped and inserted to right hand side
+ Scalar sat = 0;
+ switch (saturationType_)
+ {
+ case Sw:
+ sat = cellData.saturation(wPhaseIdx);
+ break;
+ case Sn:
+ sat = cellData.saturation(nPhaseIdx);
+ break;
+ }
+
+ Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
+ if (sat < 0.0)
+ {
+ error = sat;
+ }
+ error /= timeStep_;
+
+ Scalar errorAbs = std::abs(error);
+
+ if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
+ && (!problem_.timeManager().willBeFinished()))
+ {
+ return ErrorTermFactor_ * error;
+ }
+ return 0.0;
+ }
+
};
template<class TypeTag>
@@ -1300,10 +1361,10 @@ void FVMPFAL2PFABoundPressure2P<TypeTag>::assemble()
problem_.source(source, *elementPointer4);
this->f_[globalIdx4] += volume4 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
- this->f_[globalIdx1] += evaluateErrorTerm(cellData1) * volume1 / (4.0);
- this->f_[globalIdx2] += evaluateErrorTerm(cellData2) * volume2 / (4.0);
- this->f_[globalIdx3] += evaluateErrorTerm(cellData3) * volume3 / (4.0);
- this->f_[globalIdx4] += evaluateErrorTerm(cellData4) * volume4 / (4.0);
+ this->f_[globalIdx1] += evaluateErrorTerm_(cellData1) * volume1 / (4.0);
+ this->f_[globalIdx2] += evaluateErrorTerm_(cellData2) * volume2 / (4.0);
+ this->f_[globalIdx3] += evaluateErrorTerm_(cellData3) * volume3 / (4.0);
+ this->f_[globalIdx4] += evaluateErrorTerm_(cellData4) * volume4 / (4.0);
//get mobilities of the phases
Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
@@ -1689,7 +1750,7 @@ void FVMPFAL2PFABoundPressure2P<TypeTag>::assemble()
problem_.source(source, *elementPointer);
this->f_[globalIdx] += volume / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
- this->f_[globalIdx] += evaluateErrorTerm(cellData) * volume / (4.0);
+ this->f_[globalIdx] += evaluateErrorTerm_(cellData) * volume / (4.0);
//get mobilities of the phases
Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
@@ -1889,6 +1950,19 @@ void FVMPFAL2PFABoundPressure2P<TypeTag>::assemble()
return;
}
+/*! \brief Calculates tranmissibility matrix
+ *
+ * Calculates tranmissibility matrix of an L-shape for a certain flux face.
+ * Automatically selects one of the two possible L-shape (left, or right).
+ *
+ * \param transmissibility Matrix for the resulting transmissibility
+ * \param interactionVolume The interaction volume object (includes geometric information)
+ * \param lambda Mobilities of cells 1-4
+ * \param idx1 Index of cell 1 of the L-stencil
+ * \param idx2 Index of cell 2 of the L-stencil
+ * \param idx3 Index of cell 3 of the L-stencil
+ * \param idx4 Index of cell 4 of the L-stencil
+ */
template<class TypeTag>
bool FVMPFAL2PFABoundPressure2P<TypeTag>::calculateTransmissibility(
Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
@@ -2186,9 +2260,10 @@ bool FVMPFAL2PFABoundPressure2P<TypeTag>::calculateTransmissibility(
}
}
-
-
-//constitutive functions are updated once if new saturations are calculated and stored in the variables object
+/*! \brief Updates constitutive relations and stores them in the variable class
+ *
+ * Stores mobility, fractional flow function and capillary pressure for all grid cells.
+ */
template<class TypeTag>
void FVMPFAL2PFABoundPressure2P<TypeTag>::updateMaterialLaws()
{
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2padaptive.hh
index dddd2463f6..c12186771c 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2padaptive.hh
@@ -28,25 +28,35 @@
/**
* @file
- * @brief Finite Volume-MPFAL implementation of a two-phase pressure equation on an adaptive grid
- * @brief Remark1: only for 2-D quadrilateral grid.
- * @brief Remark2: can use UGGrid or SGrid/YaspGrid
+ * @brief Grid adaptive finite volume MPFA L-method discretization of a two-phase pressure equation of the sequential IMPES model.
* @author Markus Wolff
*/
namespace Dumux
{
-//! \ingroup diffusion
-/*! Finite Volume-MPFAO-Implementation of the equation
- * \f$ - \text{div}\, \mathbf_{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_w + f_n \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap} ) = 0, \f$, or
- * \f$ - \text{div}\, \mathbf_{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_n - f_w \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap} ) = 0, \f$.
- * \f$\Phi = g\f$ on \f$\Gamma_1\f$, and
- * \f$-\text{div}\, \mathbf_{v}_t \cdot \mathbf{n} = J\f$
- * on \f$\Gamma_2\f$. Here,
- * \f$Phi_\alpha \f$ denotes the potential of phase \f$\alpha\f$, \f$K\f$ the intrinsic permeability,
- * \f$\lambda_t\f$ the total mobility, \f$this->f_\alpha\f$ the phase fractional flow function.
- * More details on the equations can be found in
- * H. Hoteit, A. Firoozabadi. Numerical modeling of thwo-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Res (31), 2008
+//! \ingroup FVPressure2p
+/*! \brief Grid adaptive finite volume MPFA L-method discretization of a two-phase flow pressure equation of the sequential IMPES model.
+ *
+ * Grid adaptive finite volume MPFA L-method discretization of the equations
+ * \f[ - \text{div}\, \boldsymbol v_t = - \text{div}\, (\lambda_t \boldsymbol K \text{grad}\, \Phi_w + f_n \lambda_t \boldsymbol K \text{grad}\, \Phi_{cap} ) = 0, \f]
+ * or
+ * \f[ - \text{div}\, \boldsymbol v_t = - \text{div}\, (\lambda_t \boldsymbol K \text{grad}\, \Phi_n - f_w \lambda_t \boldsymbol K \text{grad}\, \Phi_{cap} ) = 0. \f]
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * \f[ \Phi = g \; \text{on} \; \Gamma_1, \quad \text{and} \quad
+ * -\text{div}\, \boldsymbol v_t \cdot \mathbf{n} = J \; \text{on} \; \Gamma_2. \f]
+ * Here, \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * \f$ \lambda_t \f$ the total mobility, \f$ f_\alpha \f$ the phase fractional flow function.
+ *
+ * More details on the equations can be found in H. Hoteit, A. Firoozabadi.
+ * Numerical modeling of thwo-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Res (31), 2008
+ *
+ * Remark1: only for 2-D quadrilateral grid
+ *
+ * Remark2: implemented for UGGrid, ALUGrid, or SGrid/YaspGrid
+ *
+ * Remark3: Allowed difference in grid levels of two neighboring cells: 1
+ *
+ *\tparam TypeTag The problem Type Tag
*/
template<class TypeTag>
class FVMPFAL2PFABoundPressure2PAdaptive: public FVPressure<TypeTag>
@@ -74,9 +84,12 @@ class FVMPFAL2PFABoundPressure2PAdaptive: public FVPressure<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolutionType;
+
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename SolutionTypes::ScalarSolution ScalarSolutionType;
typedef typename GET_PROP_TYPE(TypeTag, GridTypeIndices) GridTypeIndices;
@@ -132,6 +145,13 @@ public:
noTransmissibility = 0,
rightTriangle = 1
};
+
+ /*! \brief Type of the interaction volume objects
+ *
+ * Type of the interaction volume objects used to store the geometric information which is needed
+ * to calculated the transmissibility matrices of one MPFA interaction volume.
+ *
+ */
typedef Dumux::FVMPFALInteractionVolume<TypeTag> InteractionVolume;
private:
@@ -154,6 +174,11 @@ public:
//constitutive functions are initialized and stored in the variables object
void updateMaterialLaws();
+ /*! \brief Updates interaction volumes
+ *
+ * Globally rebuilds the MPFA interaction volumes.
+ *
+ */
void updateInteractionVolumeInfo()
{
interactionVolumes_.clear();
@@ -166,7 +191,11 @@ public:
// printInteractionVolumes();
}
- void initialize(bool solveTwice = true)
+ /*! \brief Initializes the pressure model
+ *
+ * \copydetails ParentType::initialize()
+ */
+ void initialize()
{
ParentType::initialize();
@@ -182,6 +211,9 @@ public:
return;
}
+ /*! \brief Globally stores the pressure solution
+ *
+ */
void storePressureSolution()
{
// iterate through leaf grid an evaluate c0 at cell center
@@ -192,6 +224,10 @@ public:
}
}
+ /*! \brief Stores the pressure solution of a cell
+ *
+ * \param element Dune grid element
+ */
void storePressureSolution(const Element& element)
{
int globalIdx = problem_.variables().index(element);
@@ -226,7 +262,11 @@ public:
cellData.fluxData().resetVelocity();
}
- //! updates the pressure field (analog to update function in Dumux::IMPET)
+ /*! \brief Pressure update
+ *
+ * \copydetails ParentType::update()
+ *
+ */
void update()
{
int gridSize = problem_.gridView().size(0);
@@ -275,41 +315,14 @@ public:
return;
}
- Scalar evaluateErrorTerm(CellData& cellData)
- {
- //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
- // error reduction routine: volumetric error is damped and inserted to right hand side
- Scalar sat = 0;
- switch (saturationType_)
- {
- case Sw:
- sat = cellData.saturation(wPhaseIdx);
- break;
- case Sn:
- sat = cellData.saturation(nPhaseIdx);
- break;
- }
-
- Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
- if (sat < 0.0)
- {
- error = sat;
- }
- error /= timeStep_;
-
- Scalar errorAbs = std::abs(error);
-
- if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
- && (!problem_.timeManager().willBeFinished()))
- {
- return ErrorTermFactor_ * error;
- }
- return 0.0;
- }
-
-
- //! \brief Write data files
- /* \param name file name */
+ /*! \brief Adds pressure output to the output file
+ *
+ * Adds the pressure, the potential and the capillary pressure to the output.
+ *
+ * \tparam MultiWriter Class defining the output writer
+ * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+ *
+ */
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
@@ -361,6 +374,10 @@ public:
return;
}
+ //! Constructs a FVMPFAL2PFABoundPressure2PAdaptive object
+ /**
+ * \param problem A problem class object
+ */
FVMPFAL2PFABoundPressure2PAdaptive(Problem& problem) :
ParentType(problem), problem_(problem), R_(0),
gravity_(problem.gravity()),
@@ -418,24 +435,27 @@ private:
DimMatrix R_;
public:
+ // Calculates tranmissibility matrix
int calculateTransmissibility(Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1>& transmissibility,
InteractionVolume& interactionVolume,
std::vector<DimVector >& lambda,
int idx1, int idx2, int idx3, int idx4);
+ // Calculates tranmissibility matrix of left L-shape
int calculateLeftHNTransmissibility(Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1>& transmissibility,
InteractionVolume& interactionVolume,
std::vector<DimVector >& lambda,
int idx1, int idx2, int idx3);
+ // Calculates tranmissibility matrix of right L-shape
int calculateRightHNTransmissibility(Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1>& transmissibility,
InteractionVolume& interactionVolume,
std::vector<DimVector >& lambda,
int idx1, int idx2, int idx3);
protected:
- GlobalInteractionVolumeVector interactionVolumes_;
- InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_;
+ GlobalInteractionVolumeVector interactionVolumes_;//!< Global Vector of interaction volumes
+ InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_;//!< Vector marking faces which intersect the boundary
private:
const GlobalPosition& gravity_; //!< vector including the gravity constant
@@ -454,6 +474,38 @@ private:
static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation); //!< gives kind of saturation used (\f$ 0 = S_w\f$, \f$ 1 = S_n\f$)
static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation); //!< gives kind of velocity used (\f$ 0 = v_w\f$, \f$ 1 = v_n\f$, \f$ 2 = v_t\f$)
+ Scalar evaluateErrorTerm_(CellData& cellData)
+ {
+ //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
+ // error reduction routine: volumetric error is damped and inserted to right hand side
+ Scalar sat = 0;
+ switch (saturationType_)
+ {
+ case Sw:
+ sat = cellData.saturation(wPhaseIdx);
+ break;
+ case Sn:
+ sat = cellData.saturation(nPhaseIdx);
+ break;
+ }
+
+ Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
+ if (sat < 0.0)
+ {
+ error = sat;
+ }
+ error /= timeStep_;
+
+ Scalar errorAbs = std::abs(error);
+
+ if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
+ && (!problem_.timeManager().willBeFinished()))
+ {
+ return ErrorTermFactor_ * error;
+ }
+ return 0.0;
+ }
+
};
template<class TypeTag>
@@ -1842,10 +1894,10 @@ void FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::assemble()
problem_.source(source, *elementPointer4);
this->f_[globalIdx4] += volume4 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
- this->f_[globalIdx1] += evaluateErrorTerm(cellData1) * volume1 / (4.0);
- this->f_[globalIdx2] += evaluateErrorTerm(cellData2) * volume2 / (4.0);
- this->f_[globalIdx3] += evaluateErrorTerm(cellData3) * volume3 / (4.0);
- this->f_[globalIdx4] += evaluateErrorTerm(cellData4) * volume4 / (4.0);
+ this->f_[globalIdx1] += evaluateErrorTerm_(cellData1) * volume1 / (4.0);
+ this->f_[globalIdx2] += evaluateErrorTerm_(cellData2) * volume2 / (4.0);
+ this->f_[globalIdx3] += evaluateErrorTerm_(cellData3) * volume3 / (4.0);
+ this->f_[globalIdx4] += evaluateErrorTerm_(cellData4) * volume4 / (4.0);
//get mobilities of the phases
Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
@@ -2216,8 +2268,8 @@ void FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::assemble()
problem_.source(source, *elementPointer2);
this->f_[globalIdx2] += volume2 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
- this->f_[globalIdx1] += evaluateErrorTerm(cellData1) * volume1 / (4.0);
- this->f_[globalIdx2] += evaluateErrorTerm(cellData2) * volume2 / (4.0);
+ this->f_[globalIdx1] += evaluateErrorTerm_(cellData1) * volume1 / (4.0);
+ this->f_[globalIdx2] += evaluateErrorTerm_(cellData2) * volume2 / (4.0);
//get mobilities of the phases
Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
@@ -2499,7 +2551,7 @@ void FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::assemble()
problem_.source(source, *elementPointer);
this->f_[globalIdx] += volume / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
- this->f_[globalIdx] += evaluateErrorTerm(cellData) * volume / (4.0);
+ this->f_[globalIdx] += evaluateErrorTerm_(cellData) * volume / (4.0);
//get mobilities of the phases
Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
@@ -2699,6 +2751,19 @@ void FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::assemble()
return;
}
+/*! \brief Calculates tranmissibility matrix
+ *
+ * Calculates tranmissibility matrix of an L-shape for a certain flux face.
+ * Automatically selects one of the two possible L-shape (left, or right).
+ *
+ * \param transmissibility Matrix for the resulting transmissibility
+ * \param interactionVolume The interaction volume object (includes geometric information)
+ * \param lambda Mobilities of cells 1-4
+ * \param idx1 Index of cell 1 of the L-stencil
+ * \param idx2 Index of cell 2 of the L-stencil
+ * \param idx3 Index of cell 3 of the L-stencil
+ * \param idx4 Index of cell 4 of the L-stencil
+ */
template<class TypeTag>
int FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::calculateTransmissibility(
Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1>& transmissibility,
@@ -2941,6 +3006,18 @@ int FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::calculateTransmissibility(
}
}
+/*! \brief Calculates tranmissibility matrix
+ *
+ * Calculates tranmissibility matrix of an L-shape for a certain flux face.
+ * Calculates only the transmissibility of the left L-shape (needed at hanging nodes HN).
+ *
+ * \param transmissibilityLeft Matrix for the resulting transmissibility
+ * \param interactionVolume The interaction volume object (includes geometric information)
+ * \param lambda Mobilities of cells 1-3
+ * \param idx1 Index of cell 1 of the L-stencil
+ * \param idx2 Index of cell 2 of the L-stencil
+ * \param idx3 Index of cell 3 of the L-stencil
+ */
template<class TypeTag>
int FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::calculateLeftHNTransmissibility(
Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1>& transmissibilityLeft,
@@ -3078,6 +3155,18 @@ int FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::calculateLeftHNTransmissibility
return leftTriangle;
}
+/*! \brief Calculates tranmissibility matrix
+ *
+ * Calculates tranmissibility matrix of an L-shape for a certain flux face.
+ * Calculates only the transmissibility of the right L-shape (needed at hanging nodes HN).
+ *
+ * \param transmissibilityRight Matrix for the resulting transmissibility
+ * \param interactionVolume The interaction volume object (includes geometric information)
+ * \param lambda Mobilities of cells 1-3
+ * \param idx1 Index of cell 1 of the L-stencil
+ * \param idx2 Index of cell 2 of the L-stencil
+ * \param idx3 Index of cell 3 of the L-stencil
+ */
template<class TypeTag>
int FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::calculateRightHNTransmissibility(
Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1>& transmissibilityRight,
@@ -3206,7 +3295,10 @@ int FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::calculateRightHNTransmissibilit
return rightTriangle;
}
-//constitutive functions are updated once if new saturations are calculated and stored in the variables object
+/*! \brief Updates constitutive relations and stores them in the variable class
+ *
+ * Stores mobility, fractional flow function and capillary pressure for all grid cells.
+ */
template<class TypeTag>
void FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>::updateMaterialLaws()
{
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
index 0d31e740b0..03a3e284d9 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2p.hh
@@ -25,13 +25,33 @@
/**
* @file
- * @brief Base class for defining an instance of a numerical diffusion model
+ * @brief Velocity Field from a finite volume solution of a pressure equation using a MPFA L-method.
* @author Markus Wolff
*/
namespace Dumux
{
+//! \ingroup FVPressure2p
+/*! \brief Determines the velocity from a finite volume solution of the pressure equation of a sequential model (IMPES).
+ * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA L-method.
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * The pressure has to be given as piecewise constant cell values.
+ * The velocities are calculated as
+ *
+ *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \text{grad}\, \Phi_\alpha, \f]
+ * and,
+ * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
+ *
+ * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * and \f$ \lambda_\alpha \f$ a phase mobility.
+ *
+ * Remark1: only for 2-D quadrilateral grids!
+ *
+ * Remark2: can use UGGrid, ALUGrid or SGrid/YaspGrid!
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
template<class TypeTag> class FVMPFAL2PFABoundVelocity2P: public FVMPFAL2PFABoundPressure2P<TypeTag>
{
typedef FVMPFAL2PFABoundPressure2P<TypeTag> ParentType;
@@ -58,7 +78,8 @@ template<class TypeTag> class FVMPFAL2PFABoundVelocity2P: public FVMPFAL2PFABoun
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
@@ -111,6 +132,10 @@ template<class TypeTag> class FVMPFAL2PFABoundVelocity2P: public FVMPFAL2PFABoun
typedef Dune::FieldVector<Scalar, dim> DimVector;
public:
+ //! Constructs a FVMPFAL2PFABoundVelocity2P object
+ /*!
+ * \param problem A problem class object
+ */
FVMPFAL2PFABoundVelocity2P(Problem& problem) :
ParentType(problem), problem_(problem), gravity_(problem.gravity())
{
@@ -127,17 +152,27 @@ public:
viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
}
+ //Calculates the velocities at all cell-cell interfaces.
void calculateVelocity();
- void initialize(bool solveTwice = true)
+ /*! \brief Initializes pressure and velocity
+ *
+ * \copydetails ParentType::initialize()
+ */
+ void initialize()
{
- ParentType::initialize(solveTwice);
+ ParentType::initialize();
calculateVelocity();
return;
}
+ /*! \brief Pressure and velocity update
+ *
+ * \copydetails ParentType::update()
+ *
+ */
void update()
{
ParentType::update();
@@ -147,8 +182,14 @@ public:
return;
}
- //! \brief Write data files
- /* \param name file name */
+ /*! \brief Adds velocity output to the output file
+ *
+ * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
+ *
+ * \tparam MultiWriter Class defining the output writer
+ * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+ *
+ */
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
@@ -233,6 +274,11 @@ private:
};
// end of template
+/*! \brief Calculates the velocities at a cell-cell interfaces.
+ *
+ * Calculates the velocities at a cell-cell interfaces from a given pressure field.
+ *
+ */
template<class TypeTag>
void FVMPFAL2PFABoundVelocity2P<TypeTag>::calculateVelocity()
{
@@ -263,12 +309,6 @@ void FVMPFAL2PFABoundVelocity2P<TypeTag>::calculateVelocity()
CellData& cellData3 = problem_.variables().cellData(globalIdx3);
CellData& cellData4 = problem_.variables().cellData(globalIdx4);
- // get global coordinate of cell centers
- const GlobalPosition& globalPos1 = elementPointer1->geometry().center();
- const GlobalPosition& globalPos2 = elementPointer2->geometry().center();
- const GlobalPosition& globalPos3 = elementPointer3->geometry().center();
- const GlobalPosition& globalPos4 = elementPointer4->geometry().center();
-
// get pressure values
Dune::FieldVector < Scalar, 2 * dim > pW(0);
Dune::FieldVector < Scalar, 2 * dim > pN(0);
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
index 387ed571ec..b88b21d7bf 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundvelocity2padaptive.hh
@@ -25,13 +25,35 @@
/**
* @file
- * @brief Base class for defining an instance of a numerical diffusion model
+ * @brief Velocity Field from a finite volume solution of a pressure equation using a grid adaptive MPFA L-method.
* @author Markus Wolff
*/
namespace Dumux
{
+//! \ingroup FVPressure2p
+/*! \brief Determines the velocity from a grid adaptive finite volume solution of the pressure equation of a sequential model (IMPES).
+ * Calculates phase velocities or total velocity from a known pressure field applying a grid adaptive finite volume discretization and a MPFA L-method.
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * The pressure has to be given as piecewise constant cell values.
+ * The velocities are calculated as
+ *
+ *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \text{grad}\, \Phi_\alpha, \f]
+ * and,
+ * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
+ *
+ * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * and \f$ \lambda_\alpha \f$ a phase mobility.
+ *
+ * Remark1: only for 2-D quadrilateral grids!
+ *
+ * Remark2: can use UGGrid, ALUGrid or SGrid/YaspGrid!
+ *
+ * Remark3: Allowed difference in grid levels of two neighboring cells: 1
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
template<class TypeTag> class FVMPFAL2PFABoundVelocity2PAdaptive: public FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>
{
typedef FVMPFAL2PFABoundPressure2PAdaptive<TypeTag> ParentType;
@@ -58,7 +80,8 @@ template<class TypeTag> class FVMPFAL2PFABoundVelocity2PAdaptive: public FVMPFAL
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
@@ -117,6 +140,10 @@ template<class TypeTag> class FVMPFAL2PFABoundVelocity2PAdaptive: public FVMPFAL
typedef Dune::FieldVector<Scalar, dim> DimVector;
public:
+ //! Constructs a FVMPFAO2PFABoundVelocity2P object
+ /*!
+ * \param problem A problem class object
+ */
FVMPFAL2PFABoundVelocity2PAdaptive(Problem& problem) :
ParentType(problem), problem_(problem), gravity_(problem.gravity())
{
@@ -133,17 +160,27 @@ public:
viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
}
+ //Calculates the velocities at all cell-cell interfaces.
void calculateVelocity();
- void initialize(bool solveTwice = true)
+ /*! \brief Initializes pressure and velocity
+ *
+ * \copydetails ParentType::initialize()
+ */
+ void initialize()
{
- ParentType::initialize(solveTwice);
+ ParentType::initialize();
calculateVelocity();
return;
}
+ /*! \brief Pressure and velocity update
+ *
+ * \copydetails ParentType::update()
+ *
+ */
void update()
{
ParentType::update();
@@ -153,8 +190,14 @@ public:
return;
}
- //! \brief Write data files
- /* \param name file name */
+ /*! \brief Adds velocity output to the output file
+ *
+ * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
+ *
+ * \tparam MultiWriter Class defining the output writer
+ * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+ *
+ */
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
@@ -239,6 +282,11 @@ private:
};
// end of template
+/*! \brief Calculates the velocities at a cell-cell interfaces.
+ *
+ * Calculates the velocities at a cell-cell interfaces from a given pressure field.
+ *
+ */
template<class TypeTag>
void FVMPFAL2PFABoundVelocity2PAdaptive<TypeTag>::calculateVelocity()
{
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundpressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundpressure2p.hh
index ef61767c48..a2734ced08 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundpressure2p.hh
@@ -28,25 +28,33 @@
/**
* @file
- * @brief Finite Volume-MPFAO implementation of a two-phase pressure equation
- * @brief Remark1: only for 2-D quadrilateral grid.
- * @brief Remark2: can use UGGrid or SGrid/YaspGrid
+ * @brief Finite volume MPFA O-method discretization of a two-phase pressure equation of the sequential IMPES model.
* @author Markus Wolff
*/
namespace Dumux
{
-//! \ingroup diffusion
-/*! Finite Volume-MPFAO-Implementation of the equation
- * \f$ - \text{div}\, \mathbf_{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_w + f_n \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap} ) = 0, \f$, or
- * \f$ - \text{div}\, \mathbf_{v}_t = - \text{div}\, (\lambda_t \mathbf{K} \text{grad}\, \Phi_n - f_w \lambda_t \mathbf{K} \text{grad}\, \Phi_{cap} ) = 0, \f$.
- * \f$\Phi = g\f$ on \f$\Gamma_1\f$, and
- * \f$-\text{div}\, \mathbf_{v}_t \cdot \mathbf{n} = J\f$
- * on \f$\Gamma_2\f$. Here,
- * \f$Phi_\alpha \f$ denotes the potential of phase \f$\alpha\f$, \f$K\f$ the intrinsic permeability,
- * \f$\lambda_t\f$ the total mobility, \f$this->f_\alpha\f$ the phase fractional flow function.
- * More details on the equations can be found in
- * H. Hoteit, A. Firoozabadi. Numerical modeling of thwo-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Res (31), 2008
+//! \ingroup FVPressure2p
+/*! \brief Finite volume MPFA O-method discretization of a two-phase flow pressure equation of the sequential IMPES model.
+ *
+ * Finite volume MPFA O-method discretization of the equations
+ * \f[ - \text{div}\, \boldsymbol v_t = - \text{div}\, (\lambda_t \boldsymbol K \text{grad}\, \Phi_w + f_n \lambda_t \boldsymbol K \text{grad}\, \Phi_{cap} ) = 0, \f]
+ * or
+ * \f[ - \text{div}\, \boldsymbol v_t = - \text{div}\, (\lambda_t \boldsymbol K \text{grad}\, \Phi_n - f_w \lambda_t \boldsymbol K \text{grad}\, \Phi_{cap} ) = 0. \f]
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * \f[ \Phi = g \; \text{on} \; \Gamma_1, \quad \text{and} \quad
+ * -\text{div}\, \boldsymbol v_t \cdot \mathbf{n} = J \; \text{on} \; \Gamma_2. \f]
+ * Here, \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * \f$ \lambda_t \f$ the total mobility, \f$ f_\alpha \f$ the phase fractional flow function.
+ *
+ * More details on the equations can be found in H. Hoteit, A. Firoozabadi.
+ * Numerical modeling of thwo-phase flow in heterogeneous permeable media with different capillarity pressures. Adv Water Res (31), 2008
+ *
+ * Remark1: only for 2-D quadrilateral grid
+ *
+ * Remark2: implemented for UGGrid, ALUGrid, or SGrid/YaspGrid
+ *
+ *\tparam TypeTag The problem Type Tag
*/
template<class TypeTag>
class FVMPFAO2PFABoundPressure2P: public FVPressure<TypeTag>
@@ -74,9 +82,10 @@ class FVMPFAO2PFABoundPressure2P: public FVPressure<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolutionType;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename SolutionTypes::ScalarSolution ScalarSolutionType;
typedef typename GET_PROP_TYPE(TypeTag, GridTypeIndices) GridTypeIndices;
@@ -142,6 +151,11 @@ public:
//constitutive functions are initialized and stored in the variables object
void updateMaterialLaws();
+ /*! \brief Updates interaction volumes
+ *
+ * Globally rebuilds the MPFA interaction volumes.
+ *
+ */
void updateInteractionVolumeInfo()
{
interactionVolumes_.clear();
@@ -153,7 +167,12 @@ public:
storeInteractionVolumeInfo();
}
- void initialize(bool solveTwice = true)
+ /*! \brief Initializes the pressure model
+ *
+ * \copydetails ParentType::initialize()
+ *
+ */
+ void initialize()
{
ParentType::initialize();
@@ -172,7 +191,11 @@ public:
return;
}
- //! updates the pressure field (analog to update function in Dumux::IMPET)
+ /*! \brief Pressure update
+ *
+ * \copydetails ParentType::update()
+ *
+ */
void update()
{
//error bounds for error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
@@ -209,6 +232,9 @@ public:
return;
}
+ /*! \brief Globally stores the pressure solution
+ *
+ */
void storePressureSolution()
{
// iterate through leaf grid an evaluate c0 at cell center
@@ -219,6 +245,10 @@ public:
}
}
+ /*! \brief Stores the pressure solution of a cell
+ *
+ * \param element Dune grid element
+ */
void storePressureSolution(const Element& element)
{
int globalIdx = problem_.variables().index(element);
@@ -253,8 +283,14 @@ public:
cellData.fluxData().resetVelocity();
}
- //! \brief Write data files
- /* \param name file name */
+ /*! \brief Adds pressure output to the output file
+ *
+ * Adds the pressure, the potential and the capillary pressure to the output.
+ *
+ * \tparam MultiWriter Class defining the output writer
+ * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+ *
+ */
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
@@ -300,38 +336,10 @@ public:
return;
}
- Scalar evaluateErrorTerm(CellData& cellData)
- {
- //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
- // error reduction routine: volumetric error is damped and inserted to right hand side
- Scalar sat = 0;
- switch (saturationType_)
- {
- case Sw:
- sat = cellData.saturation(wPhaseIdx);
- break;
- case Sn:
- sat = cellData.saturation(nPhaseIdx);
- break;
- }
-
- Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
- if (sat < 0.0)
- {
- error = sat;
- }
- error /= timeStep_;
-
- Scalar errorAbs = std::abs(error);
-
- if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
- && (!problem_.timeManager().willBeFinished()))
- {
- return ErrorTermFactor_ * error;
- }
- return 0.0;
- }
-
+ //! Constructs a FVMPFAO2PFABoundPressure2P object
+ /**
+ * \param problem A problem class object
+ */
FVMPFAO2PFABoundPressure2P(Problem& problem) :
ParentType(problem), problem_(problem), gravity_(problem.gravity()), maxError_(
0.), timeStep_(1.)
@@ -373,6 +381,50 @@ public:
private:
Problem& problem_;
+ /* Volume correction term to correct for unphysical saturation overshoots/undershoots.
+ * These can occur if the estimated time step for the explicit transport was too large. Correction by an artificial source term allows to correct
+ * this errors due to wrong time-stepping without losing mass conservation. The error term looks as follows:
+ * \f[
+ * q_{error} = \begin{cases}
+ * S < 0 & a_{error} \frac{S}{\Delta t} V \\
+ * S > 1 & a_{error} \frac{(S - 1)}{\Delta t} V \\
+ * 0 \le S \le 1 & 0
+ * \end{cases}
+ * \f]
+ * where \f$a_{error}\f$ is a weighting factor (default: \f$a_{error} = 0.5\f$)
+ */
+ Scalar evaluateErrorTerm_(CellData& cellData)
+ {
+ //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
+ // error reduction routine: volumetric error is damped and inserted to right hand side
+ Scalar sat = 0;
+ switch (saturationType_)
+ {
+ case Sw:
+ sat = cellData.saturation(wPhaseIdx);
+ break;
+ case Sn:
+ sat = cellData.saturation(nPhaseIdx);
+ break;
+ }
+
+ Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
+ if (sat < 0.0)
+ {
+ error = sat;
+ }
+ error /= timeStep_;
+
+ Scalar errorAbs = std::abs(error);
+
+ if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
+ && (!problem_.timeManager().willBeFinished()))
+ {
+ return ErrorTermFactor_ * error;
+ }
+ return 0.0;
+ }
+
protected:
GlobalInteractionVolumeVector interactionVolumes_;
InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_;
@@ -1378,10 +1430,10 @@ void FVMPFAO2PFABoundPressure2P<TypeTag>::assemble()
this->f_[globalIdx4] += volume4 / (4.0)
* (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
- this->f_[globalIdx1] += evaluateErrorTerm(cellData1) * volume1 / (4.0);
- this->f_[globalIdx2] += evaluateErrorTerm(cellData2) * volume2 / (4.0);
- this->f_[globalIdx3] += evaluateErrorTerm(cellData3) * volume3 / (4.0);
- this->f_[globalIdx4] += evaluateErrorTerm(cellData4) * volume4 / (4.0);
+ this->f_[globalIdx1] += evaluateErrorTerm_(cellData1) * volume1 / (4.0);
+ this->f_[globalIdx2] += evaluateErrorTerm_(cellData2) * volume2 / (4.0);
+ this->f_[globalIdx3] += evaluateErrorTerm_(cellData3) * volume3 / (4.0);
+ this->f_[globalIdx4] += evaluateErrorTerm_(cellData4) * volume4 / (4.0);
//get mobilities of the phases
Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
@@ -1771,7 +1823,7 @@ void FVMPFAO2PFABoundPressure2P<TypeTag>::assemble()
this->f_[globalIdx] += volume / (4.0)
* (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
- this->f_[globalIdx] += evaluateErrorTerm(cellData) * volume / (4.0);
+ this->f_[globalIdx] += evaluateErrorTerm_(cellData) * volume / (4.0);
//get mobilities of the phases
Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
@@ -1972,7 +2024,11 @@ void FVMPFAO2PFABoundPressure2P<TypeTag>::assemble()
return;
}
-//constitutive functions are updated once if new saturations are calculated and stored in the variables object
+
+/*! \brief Updates constitutive relations and stores them in the variable class
+ *
+ * Stores mobility, fractional flow function and capillary pressure for all grid cells.
+ */
template<class TypeTag>
void FVMPFAO2PFABoundPressure2P<TypeTag>::updateMaterialLaws()
{
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
index 787fe0eb46..944f51a21a 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2pfaboundvelocity2p.hh
@@ -25,13 +25,33 @@
/**
* @file
- * @brief Base class for defining an instance of a numerical diffusion model
+ * @brief Velocity Field from a finite volume solution of a pressure equation using a MPFA O-method.
* @author Markus Wolff
*/
namespace Dumux
{
+//! \ingroup FVPressure2p
+/*! \brief Determines the velocity from a finite volume solution of the pressure equation of a sequential model (IMPES).
+ * Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization and a MPFA O-method.
+ * At Dirichlet boundaries a two-point flux approximation is used.
+ * The pressure has to be given as piecewise constant cell values.
+ * The velocities are calculated as
+ *
+ *\f[ \boldsymbol v_\alpha = - \lambda_\alpha \boldsymbol K \text{grad}\, \Phi_\alpha, \f]
+ * and,
+ * \f[ \boldsymbol v_t = \boldsymbol v_w + \boldsymbol v_n,\f]
+ *
+ * where \f$ \Phi_\alpha \f$ denotes the potential of phase \f$ \alpha \f$, \f$ \boldsymbol K \f$ the intrinsic permeability,
+ * and \f$ \lambda_\alpha \f$ a phase mobility.
+ *
+ * Remark1: only for 2-D quadrilateral grids!
+ *
+ * Remark2: can use UGGrid, ALUGrid or SGrid/YaspGrid!
+ *
+ * \tparam TypeTag The problem Type Tag
+ */
template<class TypeTag> class FVMPFAO2PFABoundVelocity2P: public FVMPFAO2PFABoundPressure2P<TypeTag>
{
typedef FVMPFAO2PFABoundPressure2P<TypeTag> ParentType;
@@ -58,7 +78,8 @@ template<class TypeTag> class FVMPFAO2PFABoundVelocity2P: public FVMPFAO2PFABoun
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
@@ -111,6 +132,10 @@ template<class TypeTag> class FVMPFAO2PFABoundVelocity2P: public FVMPFAO2PFABoun
typedef Dune::FieldVector<Scalar, dim> DimVector;
public:
+ //! Constructs a FVMPFAO2PFABoundVelocity2P object
+ /*!
+ * \param problem A problem class object
+ */
FVMPFAO2PFABoundVelocity2P(Problem& problem) :
ParentType(problem), problem_(problem), gravity_(problem.gravity())
{
@@ -127,17 +152,27 @@ public:
viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
}
+ //Calculates the velocities at all cell-cell interfaces.
void calculateVelocity();
- void initialize(bool solveTwice = true)
+ /*! \brief Initializes pressure and velocity
+ *
+ * \copydetails ParentType::initialize()
+ */
+ void initialize()
{
- ParentType::initialize(solveTwice);
+ ParentType::initialize();
calculateVelocity();
return;
}
+ /*! \brief Pressure and velocity update
+ *
+ * \copydetails ParentType::update()
+ *
+ */
void update()
{
ParentType::update();
@@ -147,8 +182,14 @@ public:
return;
}
- //! \brief Write data files
- /* \param name file name */
+ /*! \brief Adds velocity output to the output file
+ *
+ * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
+ *
+ * \tparam MultiWriter Class defining the output writer
+ * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
+ *
+ */
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
@@ -233,6 +274,11 @@ private:
};
// end of template
+/*! \brief Calculates the velocities at a cell-cell interfaces.
+ *
+ * Calculates the velocities at a cell-cell interfaces from a given pressure field.
+ *
+ */
template<class TypeTag>
void FVMPFAO2PFABoundVelocity2P<TypeTag>::calculateVelocity()
{
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaopressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaopressure2p.hh
index 670b563f18..c473828d60 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaopressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaopressure2p.hh
@@ -36,9 +36,8 @@
namespace Dumux
{
-/*! \ingroup FVPressure2p
- *
- * \brief Finite Volume MPFA O-method discretization of a two-phase pressure equation of the sequential IMPES model.
+//!\ingroup FVPressure2p
+/*! \brief Finite Volume MPFA O-method discretization of a two-phase pressure equation of the sequential IMPES model.
*
* This class provides a finite volume (FV) implementation using the multi-point flux approximation O-method (MPFA O-method) for solving equations of the form
* \f[
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaovelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaovelocity2p.hh
index 7d5bc2a0dc..8d4a7b1da2 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaovelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaovelocity2p.hh
@@ -50,7 +50,7 @@ class FVMPFAOPressure2P;
* Remark2: can use UGGrid or SGrid/YaspGrid!
* Remark3: gravity is neglected!
*
- * \tparam TypeTag The Type Tag
+ * \tparam TypeTag The problem Type Tag
*/
template<class TypeTag> class FVMPFAOVelocity2P:public FVMPFAOPressure2P<TypeTag>
{
diff --git a/dumux/decoupled/common/fv/fvpressure.hh b/dumux/decoupled/common/fv/fvpressure.hh
index f4812c7049..54af178f9b 100644
--- a/dumux/decoupled/common/fv/fvpressure.hh
+++ b/dumux/decoupled/common/fv/fvpressure.hh
@@ -173,11 +173,13 @@ public:
const Scalar pressure(const int globalIdx) const
{ return pressure_[globalIdx];}
+ //!Returns the global matrix of the last pressure solution step.
const Matrix& globalMatrix()
{
return A_;
}
+ //!Returns the right hand side vector of the last pressure solution step.
const RHSVector& rightHandSide()
{
return f_;
diff --git a/dumux/decoupled/common/fv/mpfa/mpfalinteractionvolume.hh b/dumux/decoupled/common/fv/mpfa/mpfalinteractionvolume.hh
index bd55d4a1d4..3a9a87dcc2 100644
--- a/dumux/decoupled/common/fv/mpfa/mpfalinteractionvolume.hh
+++ b/dumux/decoupled/common/fv/mpfa/mpfalinteractionvolume.hh
@@ -22,7 +22,7 @@
/**
* @file
- * @brief Class including the information of an interaction volume of a MPFA method that does not change with time
+ * @brief Class including the information of an interaction volume of a MPFA L-method that does not change with time.
* @author Markus Wolff
*/
@@ -31,6 +31,12 @@
namespace Dumux
{
+//! \ingroup IMPET
+/*! \brief Class including the information of an interaction volume of a MPFA L-method that does not change with time.
+ *
+ * Includes information needed to calculate the transmissibility matrix of an L-interaction-volume.
+ *
+ */
template<class TypeTag>
class FVMPFALInteractionVolume
{
@@ -47,7 +53,9 @@ private:
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ ///@cond 0
typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+ ///@endcond
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef Dune::FieldVector<DimVector, dim> FieldVectorVector;
@@ -63,9 +71,7 @@ public:
outside = -1
};
- //! Constructs a FVMPFALInteractionVolumeInfo object
- /**
- */
+ //! Constructs a FVMPFALInteractionVolume object
FVMPFALInteractionVolume() :
stored_(false), normal_(FieldVectorVector(DimVector(0.0))),
facePos_(FieldVectorVector(DimVector(0.0))),
@@ -86,6 +92,7 @@ public:
faceIndexOnSubVolume_[3][1] = 2;
}
+ //! Delete stored information
void reset()
{
stored_ = false;
@@ -101,42 +108,73 @@ public:
elementNum_ = 0;
}
+ //! Mark storage as completed
void setStored()
{
stored_ = true;
}
+ //! Returns true if information has already been stored
bool isStored() const
{
return stored_;
}
+ //! Store position of the center vertex of the interaction volume
void setCenterPosition(DimVector ¢erVertexPos)
{
centerVertexPos_ = centerVertexPos;
}
+ //! Store an element of the interaction volume
+ /*!
+ * \param pointer The Dune::EntityPointer to the element
+ * \param subVolumeIdx The local element index in the interaction volume
+ */
void setSubVolumeElement(ElementPointer pointer, int subVolumeIdx)
{
elements_[subVolumeIdx].push_back(pointer);
elementNum_++;
}
+ //! Store position of the center of a flux face
+ /*!
+ * \param pos Position of face center
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ */
void setFacePosition(const DimVector& pos, int subVolumeIdx, int subVolumeFaceIdxInInside)
{
facePos_[subVolumeIdx][subVolumeFaceIdxInInside] = pos;
}
+ //! Store a flux face area
+ /*!
+ * \param faceArea The flux face area
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ */
void setFaceArea(Scalar& faceArea, int subVolumeIdx, int subVolumeFaceIdxInInside)
{
faceArea_[subVolumeIdx][subVolumeFaceIdxInInside] = faceArea;
}
+ //! Store a flux face normal
+ /*!
+ * \param normal The normal vector
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ */
void setNormal(DimVector& normal, int subVolumeIdx, int subVolumeFaceIdxInInside)
{
normal_[subVolumeIdx][subVolumeFaceIdxInInside] = normal;
}
+ //! Store boundary condtion types for a flux face
+ /*!
+ * \param boundaryTypes BoundaryTypes object
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
void setBoundary(BoundaryTypes& boundaryTypes, int subVolumeFaceIdx)
{
if (boundaryTypes_.size() == 0)
@@ -147,11 +185,20 @@ public:
faceType_[subVolumeFaceIdx] = boundary;
}
+ //! Mark a flux face to be outside the model domain
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
void setOutsideFace(int subVolumeFaceIdx)
{
faceType_[subVolumeFaceIdx] = outside;
}
+ //! Store Dirichlet boundary condtions for a flux face
+ /*!
+ * \param condition Vector of primary variables
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
void setDirichletCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
{
if (dirichletValues_.size() == 0)
@@ -161,6 +208,11 @@ public:
dirichletValues_[subVolumeFaceIdx] = condition;
}
+ //! Store Neumann boundary condtions for a flux face
+ /*!
+ * \param condition Vector phase fluxes
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
void setNeumannCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
{
if (neumannValues_.size() == 0)
@@ -170,51 +222,94 @@ public:
neumannValues_[subVolumeFaceIdx] = condition;
}
- DimVector& getCenterPosition()
+ //! Store map from local interaction volume numbering to numbering of the Dune reference element.
+ /*!
+ * \param indexInInside Face index of the Dune reference element
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ */
+ void setIndexOnElement(int indexInInside, int subVolumeIdx, int subVolumeFaceIdxInInside)
{
- return centerVertexPos_;
+ indexOnElement_[subVolumeIdx][subVolumeFaceIdxInInside] = indexInInside;
}
- void setIndexOnElement(int indexInInside, int subVolumeIdx, int subVolumeFaceIdxInInside)
+ //! Get position vector of central vertex
+ DimVector& getCenterPosition()
{
- indexOnElement_[subVolumeIdx][subVolumeFaceIdxInInside] = indexInInside;
+ return centerVertexPos_;
}
+ //! Get number of stored elements
int getElementNumber()
{
return elementNum_;
}
+ //! Map from local interaction volume numbering to numbering of the Dune reference element.
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdx The local face index in the interaction volume element
+ *
+ * \return Face index of the Dune reference element
+ */
int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
{
return indexOnElement_[subVolumeIdx][subVolumeFaceIdx];
}
+ //! Map from local interaction volume numbering on element to numbering on interaction volume.
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdx The local face index in the interaction volume element
+ *
+ * \return Local face index int the interaction volume
+ */
int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
{
return faceIndexOnSubVolume_[subVolumeIdx][subVolumeFaceIdx];
}
+ //! Get an element of the interaction volume.
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ *
+ * \return Dune::EntityPointer to the interaction volume sub-element.
+ */
ElementPointer& getSubVolumeElement(int subVolumeIdx)
{
return elements_[subVolumeIdx][0];
}
+ //! Get boundary condtion types for a flux face
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ *
+ * \return Object containing information about the boundary types.
+ */
BoundaryTypes& getBoundaryType(int subVolumeFaceIdx)
{
return boundaryTypes_[subVolumeFaceIdx];
}
+ //! Returns true if an interaction volume flux face is outside the model domain.
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
bool isOutsideFace(int subVolumeFaceIdx)
{
return (faceType_[subVolumeFaceIdx] == outside);
}
+ //! Returns true if an interaction volume flux face is inside the model domain and no boundary face.
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
bool isInsideFace(int subVolumeFaceIdx)
{
return (faceType_[subVolumeFaceIdx] == inside);
}
+ //! Returns true if the interaction volume is completely inside the model domain.
bool isInnerVolume()
{
for (int i = 0; i < faceType_.size(); i++)
@@ -225,36 +320,74 @@ public:
return true;
}
+ //! Returns true if an interaction volume flux face is a boundary face.
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
bool isBoundaryFace(int subVolumeFaceIdx)
{
return (faceType_[subVolumeFaceIdx] == boundary);
}
+ //! Get the Dirichlet boundary condtions for a flux face
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ *
+ * \return Vector of primary variables.
+ */
PrimaryVariables& getDirichletValues(int subVolumeFaceIdx)
{
return dirichletValues_[subVolumeFaceIdx];
}
+ //! Get the Neumann boundary condtions for a flux face
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ *
+ * \return Vector of phase fluxes.
+ */
PrimaryVariables& getNeumannValues(int subVolumeFaceIdx)
{
return neumannValues_[subVolumeFaceIdx];
}
+ //! Get a flux face normal
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ *
+ * \return Normal vector
+ */
DimVector& getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
{
return normal_[subVolumeIdx][subVolumeFaceIdxInInside];
}
+ //! Get position of the center of a flux face
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ *
+ * \return Position of face center
+ */
DimVector& getFacePosition(int subVolumeIdx, int subVolumeFaceIdxInInside)
{
return facePos_[subVolumeIdx][subVolumeFaceIdxInInside];
}
+ //! Get a flux face area
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ *
+ * \return Area of the flux face
+ */
Scalar& getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
{
return faceArea_[subVolumeIdx][subVolumeFaceIdxInInside];
}
+ //! Outputs stored information
void printInteractionVolumeInfo()
{
std::cout<<"\nNumber of stored elements: "<<elementNum_<<"\n";
diff --git a/dumux/decoupled/common/fv/mpfa/mpfaointeractionvolume.hh b/dumux/decoupled/common/fv/mpfa/mpfaointeractionvolume.hh
index bc5e5fe66c..2e615ac058 100644
--- a/dumux/decoupled/common/fv/mpfa/mpfaointeractionvolume.hh
+++ b/dumux/decoupled/common/fv/mpfa/mpfaointeractionvolume.hh
@@ -22,7 +22,7 @@
/**
* @file
- * @brief Class including the information of an interaction volume of a MPFA method that does not change with time
+ * @brief Class including the information of an interaction volume of a MPFA O-method that does not change with time.
* @author Markus Wolff
*/
@@ -31,6 +31,12 @@
namespace Dumux
{
+//! \ingroup IMPET
+/*! \brief Class including the information of an interaction volume of a MPFA O-method that does not change with time.
+ *
+ * Includes information needed to calculate the transmissibility matrix of an O-interaction-volume.
+ *
+ */
template<class TypeTag>
class FVMPFAOInteractionVolume
{
@@ -47,7 +53,9 @@ private:
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ ///@cond 0
typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+ ///@endcond
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
@@ -64,9 +72,7 @@ public:
outside = -1
};
- //! Constructs a InteractionVolumeInfo object
- /**
- */
+ //! Constructs a FVMPFAOInteractionVolume object
FVMPFAOInteractionVolume() :
stored_(false), permTimesNu_(FieldVectorVector(DimVector(0.0))),
nu_(FieldVectorVector(DimVector(0.0))), normal_(FieldVectorVector(DimVector(0.0))),
@@ -85,31 +91,56 @@ public:
faceIndexOnSubVolume_[3][1] = 2;
}
+ //! Mark storage as completed
void setStored()
{
stored_ = true;
}
+ //! Returns true if information has already been stored
bool isStored() const
{
return stored_;
}
+ //! Store an element of the interaction volume
+ /*!
+ * \param pointer The Dune::EntityPointer to the element
+ * \param subVolumeIdx The local element index in the interaction volume
+ */
void setSubVolumeElement(ElementPointer pointer, int subVolumeIdx)
{
elements_[subVolumeIdx].push_back(pointer);
}
+ //! Store the \f$ dF \f$ for the transmissiblity calculation
+ /*!
+ * \param dF Value of \f$ dF \f$
+ * \param subVolumeIdx The local element index in the interaction volume
+ */
void setDF(Scalar dF, int subVolumeIdx)
{
dF_[subVolumeIdx] = dF;
}
+ //! Store a flux face area
+ /*!
+ * \param faceArea The flux face area
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ */
void setFaceArea(Scalar& faceArea, int subVolumeIdx, int subVolumeFaceIdxInInside)
{
faceArea_[subVolumeIdx][subVolumeFaceIdxInInside] = faceArea;
}
+ //! Store \f$ \boldsymbol K \boldsymbol \nu\f$ for the transmissiblity calculation
+ /*!
+ * \param nu \f$ \boldsymbol \nu \f$ vector
+ * \param perm Permeability matrix
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ */
void setPermTimesNu(DimVector& nu, DimMatrix& perm, int subVolumeIdx, int subVolumeFaceIdxInInside)
{
permTimesNu_[subVolumeIdx][subVolumeFaceIdxInInside] = 0;
@@ -119,11 +150,22 @@ public:
nu_[subVolumeIdx][subVolumeFaceIdxInInside] = nu;
}
+ //! Store a flux face normal
+ /*!
+ * \param normal The normal vector
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ */
void setNormal(DimVector& normal, int subVolumeIdx, int subVolumeFaceIdxInInside)
{
normal_[subVolumeIdx][subVolumeFaceIdxInInside] = normal;
}
+ //! Store boundary condtion types for a flux face
+ /*!
+ * \param boundaryTypes BoundaryTypes object
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
void setBoundary(BoundaryTypes& boundaryTypes, int subVolumeFaceIdx)
{
// std::cout<<"old BCType = "<<boundaryType_[subVolumeFaceIdx]<<"\n";
@@ -136,11 +178,20 @@ public:
faceType_[subVolumeFaceIdx] = boundary;
}
+ //! Mark a flux face to be outside the model domain
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
void setOutsideFace(int subVolumeFaceIdx)
{
faceType_[subVolumeFaceIdx] = outside;
}
+ //! Store Dirichlet boundary condtions for a flux face
+ /*!
+ * \param condition Vector of primary variables
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
void setDirichletCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
{
if (dirichletValues_.size() == 0)
@@ -150,6 +201,11 @@ public:
dirichletValues_[subVolumeFaceIdx] = condition;
}
+ //! Store Neumann boundary condtions for a flux face
+ /*!
+ * \param condition Vector phase fluxes
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
void setNeumannCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
{
if (neumannValues_.size() == 0)
@@ -159,41 +215,82 @@ public:
neumannValues_[subVolumeFaceIdx] = condition;
}
+ //! Store map from local interaction volume numbering to numbering of the Dune reference element.
+ /*!
+ * \param indexInInside Face index of the Dune reference element
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ */
void setIndexOnElement(int indexInInside, int subVolumeIdx, int subVolumeFaceIdxInInside)
{
indexOnElement_[subVolumeIdx][subVolumeFaceIdxInInside] = indexInInside;
}
+ //! Map from local interaction volume numbering to numbering of the Dune reference element.
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdx The local face index in the interaction volume element
+ *
+ * \return Face index of the Dune reference element
+ */
int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
{
return indexOnElement_[subVolumeIdx][subVolumeFaceIdx];
}
+ //! Map from local interaction volume numbering on element to numbering on interaction volume.
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdx The local face index in the interaction volume element
+ *
+ * \return Local face index int the interaction volume
+ */
int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
{
return faceIndexOnSubVolume_[subVolumeIdx][subVolumeFaceIdx];
}
+ //! Get an element of the interaction volume.
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ *
+ * \return Dune::EntityPointer to the interaction volume sub-element.
+ */
ElementPointer& getSubVolumeElement(int subVolumeIdx)
{
return elements_[subVolumeIdx][0];
}
+ //! Get boundary condtion types for a flux face
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ *
+ * \return Object containing information about the boundary types.
+ */
BoundaryTypes& getBoundaryType(int subVolumeFaceIdx)
{
return boundaryTypes_[subVolumeFaceIdx];
}
+ //! Returns true if an interaction volume flux face is outside the model domain.
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
bool isOutsideFace(int subVolumeFaceIdx)
{
return (faceType_[subVolumeFaceIdx] == outside);
}
+ //! Returns true if an interaction volume flux face is inside the model domain and no boundary face.
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
bool isInsideFace(int subVolumeFaceIdx)
{
return (faceType_[subVolumeFaceIdx] == inside);
}
+ //! Returns true if the interaction volume is completely inside the model domain.
bool isInnerVolume()
{
for (int i = 0; i < faceType_.size(); i++)
@@ -204,43 +301,95 @@ public:
return true;
}
+ //! Returns true if an interaction volume flux face is a boundary face.
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ */
bool isBoundaryFace(int subVolumeFaceIdx)
{
return (faceType_[subVolumeFaceIdx] == boundary);
}
+ //! Get the Dirichlet boundary condtions for a flux face
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ *
+ * \return Vector of primary variables.
+ */
PrimaryVariables& getDirichletValues(int subVolumeFaceIdx)
{
return dirichletValues_[subVolumeFaceIdx];
}
+ //! Get the Neumann boundary condtions for a flux face
+ /*!
+ * \param subVolumeFaceIdx The local face index in the interaction volume
+ *
+ * \return Vector of phase fluxes.
+ */
PrimaryVariables& getNeumannValues(int subVolumeFaceIdx)
{
return neumannValues_[subVolumeFaceIdx];
}
+ //! Get a flux face normal
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ *
+ * \return Normal vector
+ */
DimVector& getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
{
return normal_[subVolumeIdx][subVolumeFaceIdxInInside];
}
+ //! Get a flux face area
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ *
+ * \return Area of the flux face
+ */
Scalar& getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
{
return faceArea_[subVolumeIdx][subVolumeFaceIdxInInside];
}
+ //! Get \f$ \boldsymbol \nu \f$ vector used for the transmissiblity calculation
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
+ *
+ * \return \f$ \boldsymbol \nu \f$ vector
+ */
DimVector& getNu(int subVolumeIdx, int subVolumeFaceIdxInInside)
{
return nu_[subVolumeIdx][subVolumeFaceIdxInInside];
}
- //returns n^TKK_rnu
+ //! Get \f$ \boldsymbol n^\text{T} \boldsymbol K \boldsymbol \nu \f$ for the transmissiblity calculation
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInsideN The local face index in the interaction volume element for the normal \f$ \boldsymbol n \f$
+ * \param subVolumeFaceIdxInInsideNu The local face index in the interaction volume element for the vector \f$ \boldsymbol \nu \f$
+ *
+ * \return \f$ \boldsymbol n^\text{T} \boldsymbol K \boldsymbol \nu \f$
+ */
Scalar getNTKNu(int subVolumeIdx, int subVolumeFaceIdxInInsideN, int subVolumeFaceIdxInInsideNu) const
{
return normal_[subVolumeIdx][subVolumeFaceIdxInInsideN] * permTimesNu_[subVolumeIdx][subVolumeFaceIdxInInsideNu];
}
- //returns n^TKK_rnu
+ //! Get \f$ \boldsymbol n^\text{T} k_{r\alpha} \boldsymbol K \boldsymbol \nu\f$ for the transmissiblity calculation
+ /*!
+ * \param relPerm relative permeability value (\f$ \boldsymbol n^\text{T} k_{r\alpha} \f$)
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInsideN The local face index in the interaction volume element for the normal \f$ \boldsymbol n \f$
+ * \param subVolumeFaceIdxInInsideNu The local face index in the interaction volume element for the vector \f$ \boldsymbol \nu \f$
+ *
+ * \return \f$ \boldsymbol n^\text{T} k_{r\alpha} \boldsymbol K \boldsymbol \nu \f$
+ */
Scalar getNTKrKNu(Scalar& relPerm, int subVolumeIdx, int subVolumeFaceIdxInInsideN, int subVolumeFaceIdxInInsideNu) const
{
DimVector krKNu(permTimesNu_[subVolumeIdx][subVolumeFaceIdxInInsideNu]);
@@ -250,13 +399,28 @@ public:
return normal_[subVolumeIdx][subVolumeFaceIdxInInsideN] * krKNu;
}
- //returns n^TKK_rnu/dF
+ //! Get \f$ \frac{1}{dF} \left(\boldsymbol n^\text{T} \boldsymbol K \boldsymbol \nu \right) \f$ for the transmissiblity calculation
+ /*!
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInsideN The local face index in the interaction volume element for the normal \f$ \boldsymbol n \f$
+ * \param subVolumeFaceIdxInInsideNu The local face index in the interaction volume element for the vector \f$ \boldsymbol \nu \f$
+ *
+ * \return \f$ \frac{1}{dF} \left(\boldsymbol n^\text{T} \boldsymbol K \boldsymbol \nu \right) \f$
+ */
Scalar getNTKNu_by_dF(int subVolumeIdx, int subVolumeFaceIdxInInsideN, int subVolumeFaceIdxInInsideNu) const
{
return faceArea_[subVolumeIdx][subVolumeFaceIdxInInsideN]*getNTKNu(subVolumeIdx, subVolumeFaceIdxInInsideN, subVolumeFaceIdxInInsideNu) / dF_[subVolumeIdx];
}
- //returns n^TKK_rnu/dF
+ //! Get \f$ \frac{1}{dF} \left(\boldsymbol n^\text{T} k_{r\alpha} \boldsymbol K \boldsymbol \nu \right) \f$ for the transmissiblity calculation
+ /*!
+ * \param relPerm relative permeability value (\f$ \boldsymbol n^\text{T} k_{r\alpha} \f$)
+ * \param subVolumeIdx The local element index in the interaction volume
+ * \param subVolumeFaceIdxInInsideN The local face index in the interaction volume element for the normal \f$ \boldsymbol n \f$
+ * \param subVolumeFaceIdxInInsideNu The local face index in the interaction volume element for the vector \f$ \boldsymbol \nu \f$
+ *
+ * \return \f$ \frac{1}{dF} \left(\boldsymbol n^\text{T} k_{r\alpha} \boldsymbol K \boldsymbol \nu \right) \f$
+ */
Scalar getNTKrKNu_by_dF(Scalar& relPerm, int subVolumeIdx, int subVolumeFaceIdxInInsideN, int subVolumeFaceIdxInInsideNu) const
{
return faceArea_[subVolumeIdx][subVolumeFaceIdxInInsideN]*getNTKrKNu(relPerm, subVolumeIdx, subVolumeFaceIdxInInsideN, subVolumeFaceIdxInInsideNu) / dF_[subVolumeIdx];
diff --git a/test/decoupled/2p/test_mpfa2pproblem.hh b/test/decoupled/2p/test_mpfa2pproblem.hh
index 53fa2f3f8f..638636654a 100644
--- a/test/decoupled/2p/test_mpfa2pproblem.hh
+++ b/test/decoupled/2p/test_mpfa2pproblem.hh
@@ -17,6 +17,13 @@
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
+
+/*!
+ * \file
+ *
+ * \brief test problem for sequential 2p models
+ */
+
#ifndef DUMUX_TEST_MPFA2P_PROBLEM_HH
#define DUMUX_TEST_MPFA2P_PROBLEM_HH
@@ -107,7 +114,17 @@ NEW_TYPE_TAG(MPFALAdaptiveTwoPTestProblem, INHERITS_FROM(FVMPFALPressureTwoPAdap
}
/*!
- * \ingroup DecoupledProblems
+ * \ingroup IMPETtests
+ *
+ * \brief test problem for sequential 2p models
+ *
+ * A DNAPL is injected from the top into a rectangular 2D domain saturated by water.
+ * The remaining upper and the lower boundary is closed (Neumann = 0). At the sides a hydrostatic pressure condition
+ * and free outflow for saturation are set. The domain is heterogeneous with a backround material and three lenses.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_mpfa2p -parameterFile ./test_mpfa2p.input</tt>,
+ * where the arguments define the parameter file..
*/
template<class TypeTag>
class MPFATwoPTestProblem: public IMPESProblem2P<TypeTag>
@@ -124,8 +141,9 @@ typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
-typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
-typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolutionType;
+typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+typedef typename SolutionTypes::ScalarSolution ScalarSolutionType;
typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
@@ -234,7 +252,7 @@ Scalar temperatureAtPos(const GlobalPosition& globalPos) const
// \}
-
+//! Returns the reference pressure for evaluation of constitutive relations
Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
{
return 1e5; // -> 10°C
@@ -245,6 +263,13 @@ void source(PrimaryVariables &values,const Element& element) const
values = 0;
}
+/*!
+* \brief Returns the type of boundary condition.
+*
+* BC for pressure equation can be dirichlet (pressure) or neumann (flux).
+*
+* BC for saturation equation can be dirichlet (saturation), neumann (flux), or outflow.
+*/
void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos) const
{
if (isInlet(globalPos))
@@ -263,6 +288,7 @@ void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos)
}
}
+//! set dirichlet condition (pressure [Pa], saturation [-])
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
{
Scalar pRef = referencePressureAtPos(globalPos);
@@ -277,6 +303,7 @@ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) c
}
}
+//! set neumann condition for phases (flux, [kg/(m^2 s)])
void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
{
values = 0;
@@ -287,6 +314,7 @@ void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) con
}
+//! return initial solution -> only saturation values have to be given!
void initialAtPos(PrimaryVariables &values,
const GlobalPosition& globalPos) const
{
--
GitLab