diff --git a/dumux/boxmodels/common/boxdarcyfluxvariables.hh b/dumux/boxmodels/common/boxdarcyfluxvariables.hh
index 953d7a3b62db147c1499a5cbfd613ff4ef0848a8..c7707293676d7f96660be763dcdc6bf063424d78 100644
--- a/dumux/boxmodels/common/boxdarcyfluxvariables.hh
+++ b/dumux/boxmodels/common/boxdarcyfluxvariables.hh
@@ -1,3 +1,78 @@
-#warning This file is deprecated. Include dumux/implicit/box/boxdarcyfluxvariables.hh instead.
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * 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 This file contains the data which is required to calculate
+ * all fluxes of fluid phases over a face of a finite volume.
+ *
+ * This means pressure and temperature gradients, phase densities at
+ * the integration point, etc.
+ */
+#ifndef DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
+#define DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup BoxModel
+ * \ingroup BoxFluxVariables
+ * \brief Evaluates the normal component of the Darcy velocity
+ * on a (sub)control volume face.
+ */
+template <class TypeTag>
+class BoxDarcyFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
+{
+ typedef ImplicitDarcyFluxVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+public:
+ /*
+ * \brief The constructor
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry in the box scheme
+ * \param faceIdx The local index of the SCV (sub-control-volume) face
+ * \param elemVolVars The volume variables of the current element
+ * \param onBoundary A boolean variable to specify whether the flux variables
+ * are calculated for interior SCV faces or boundary faces, default=false
+ */
+ DUNE_DEPRECATED_MSG("Use ImplicitDarcyFluxVariables from "
+ "dumux/implicit/common/implicitdarcyfluxvariables.hh.")
+ BoxDarcyFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int faceIdx,
+ const ElementVolumeVariables &elemVolVars,
+ const bool onBoundary = false)
+ : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+ {}
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/boxmodels/common/boxforchheimerfluxvariables.hh b/dumux/boxmodels/common/boxforchheimerfluxvariables.hh
index 186d17e2c90598e68cf1c48c091540469b4eefea..422f53b93320e92fd11174dd0b48366d8d5a4fcc 100644
--- a/dumux/boxmodels/common/boxforchheimerfluxvariables.hh
+++ b/dumux/boxmodels/common/boxforchheimerfluxvariables.hh
@@ -1,3 +1,101 @@
-#warning This file is deprecated. Include dumux/implicit/box/boxforchheimerfluxvariables.hh instead.
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * 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 This file contains the data which is required to calculate
+ * all fluxes of fluid phases over a face of a finite volume,
+ * according to the Forchheimer-relation between velocity and pressure.
+ *
+ */
+#ifndef DUMUX_BOX_FORCHHEIMER_FLUX_VARIABLES_HH
+#define DUMUX_BOX_FORCHHEIMER_FLUX_VARIABLES_HH
-#include <dumux/implicit/box/boxforchheimerfluxvariables.hh>
+#include <dumux/implicit/common/implicitforchheimerfluxvariables.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup BoxModel
+ * \ingroup BoxFluxVariables
+ * \brief Evaluates the normal component of the Forchheimer velocity
+ * on a (sub)control volume face.
+ *
+ * The commonly used Darcy relation looses it's validity for \f$ Re < 1\f$.
+ * If one encounters flow velocities in porous media above this Reynolds number,
+ * the Forchheimer relation can be used. Like the Darcy relation, it relates
+ * the gradient in potential to velocity.
+ * However, this relation is not linear (as in the Darcy case) any more.
+ *
+ * Therefore, a Newton scheme is implemented in this class in order to calculate a
+ * velocity from the current set of variables. This velocity can subsequently be used
+ * e.g. by the localresidual.
+ *
+ * For Reynolds numbers above \f$ 500\f$ the (Standard) forchheimer relation also
+ * looses it's validity.
+ *
+ * The Forchheimer equation looks as follows:
+ * \f[ \nabla \left({p_\alpha + \rho_\alpha g z} \right)= - \frac{\mu_\alpha}{k_{r \alpha} K} \underline{v_\alpha} - \frac{c_F}{\eta_{\alpha r} \sqrt{K}} \rho_\alpha |\underline{v_\alpha}| \underline{v_\alpha} \f]
+ *
+ * For the formulation that is actually used in this class, see the documentation of the function calculating the residual.
+ *
+ * - This algorithm does not find a solution if the fluid is incompressible and the
+ * initial pressure distribution is uniform.
+ * - This algorithm assumes the relative passabilities to have the same functional
+ * form as the relative permeabilities.
+ */
+template <class TypeTag>
+class BoxForchheimerFluxVariables
+ : public ImplicitForchheimerFluxVariables<TypeTag>
+{
+ typedef ImplicitForchheimerFluxVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry in the box scheme
+ * \param faceIdx The local index of the SCV (sub-control-volume) face
+ * \param elemVolVars The volume variables of the current element
+ * \param onBoundary A boolean variable to specify whether the flux variables
+ * are calculated for interior SCV faces or boundary faces, default=false
+ */
+ DUNE_DEPRECATED_MSG("Use ImplicitForchheimerFluxVariables from "
+ "dumux/implicit/common/implicitforchheimerfluxvariables.hh.")
+ BoxForchheimerFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const unsigned int faceIdx,
+ const ElementVolumeVariables &elemVolVars,
+ const bool onBoundary = false)
+ : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+ {}
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/1p/1ppropertydefaults.hh b/dumux/implicit/1p/1ppropertydefaults.hh
index 34b3c254d4a93f329396d400c6e032c45759325e..7bacc08a40ce7e52aee7e14664cc36038abc9aa6 100644
--- a/dumux/implicit/1p/1ppropertydefaults.hh
+++ b/dumux/implicit/1p/1ppropertydefaults.hh
@@ -32,7 +32,7 @@
#include "1pmodel.hh"
#include "1plocalresidual.hh"
#include "1pvolumevariables.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include "1pindices.hh"
#include <dumux/material/fluidsystems/gasphase.hh>
@@ -64,7 +64,7 @@ SET_TYPE_PROP(OneP, Model, OnePBoxModel<TypeTag>);
SET_TYPE_PROP(OneP, VolumeVariables, OnePVolumeVariables<TypeTag>);
//! the FluxVariables property
-SET_TYPE_PROP(OneP, FluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(OneP, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
//! The indices required by the isothermal single-phase model
SET_TYPE_PROP(OneP, Indices, OnePIndices);
diff --git a/dumux/implicit/1p/1pvolumevariables.hh b/dumux/implicit/1p/1pvolumevariables.hh
index c2e4588845a51278c23c9b897b2ffc1a8eb73285..01cbeab00b7809f6c3dee6c1fcdd9856c2ea75d3 100644
--- a/dumux/implicit/1p/1pvolumevariables.hh
+++ b/dumux/implicit/1p/1pvolumevariables.hh
@@ -145,7 +145,7 @@ public:
/*!
* \brief Returns the mobility.
*
- * This function enables the use of BoxDarcyFluxVariables
+ * This function enables the use of ImplicitDarcyFluxVariables
* with the 1p box model, ALTHOUGH the term mobility is
* usually not employed in the one phase context.
*
diff --git a/dumux/implicit/1p2c/1p2cfluxvariables.hh b/dumux/implicit/1p2c/1p2cfluxvariables.hh
index 9dd73c15635172005f11af43b8b22105e5cabd7d..eca03684d0d2abb5d61c84759bb105e23b8ad148 100644
--- a/dumux/implicit/1p2c/1p2cfluxvariables.hh
+++ b/dumux/implicit/1p2c/1p2cfluxvariables.hh
@@ -363,9 +363,17 @@ protected:
}
else
{
+ const Element& elementI = *fvGeometry_.neighbors[face().i];
+ FVElementGeometry fvGeometryI;
+ fvGeometryI.subContVol[0].global = elementI.geometry().center();
+
+ const Element& elementJ = *fvGeometry_.neighbors[face().j];
+ FVElementGeometry fvGeometryJ;
+ fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
+
sp.meanK(K_,
- sp.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
- sp.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
+ sp.intrinsicPermeability(elementI, fvGeometryI, 0),
+ sp.intrinsicPermeability(elementJ, fvGeometryJ, 0));
}
}
diff --git a/dumux/implicit/2p/2ppropertydefaults.hh b/dumux/implicit/2p/2ppropertydefaults.hh
index ff139c3bf71d02da3336c09fcfccd285091a14c3..5fea8679838004dcbc5f7ae478eb2404ffc9db0a 100644
--- a/dumux/implicit/2p/2ppropertydefaults.hh
+++ b/dumux/implicit/2p/2ppropertydefaults.hh
@@ -30,7 +30,7 @@
#include "2pmodel.hh"
#include "2pindices.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include "2pvolumevariables.hh"
#include "2pproperties.hh"
@@ -68,7 +68,7 @@ SET_TYPE_PROP(TwoP, Model, TwoPModel<TypeTag>);
SET_TYPE_PROP(TwoP, VolumeVariables, TwoPVolumeVariables<TypeTag>);
//! the FluxVariables property
-SET_TYPE_PROP(TwoP, FluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(TwoP, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
//! the upwind weight for the mass conservation equations.
SET_SCALAR_PROP(TwoP, ImplicitMassUpwindWeight, 1.0);
diff --git a/dumux/implicit/2p2c/2p2cpropertydefaults.hh b/dumux/implicit/2p2c/2p2cpropertydefaults.hh
index 8a41d489879a0511f30b4be9af6e4ece8e49ed2e..7ecbdb6c4391c4f278e8a79ec790036e10b8a395 100644
--- a/dumux/implicit/2p2c/2p2cpropertydefaults.hh
+++ b/dumux/implicit/2p2c/2p2cpropertydefaults.hh
@@ -35,7 +35,7 @@
#include "2p2cproperties.hh"
#include "2p2cnewtoncontroller.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include <dumux/material/spatialparams/boxspatialparams.hh>
namespace Dumux
@@ -122,7 +122,7 @@ SET_TYPE_PROP(TwoPTwoC, VolumeVariables, TwoPTwoCVolumeVariables<TypeTag>);
SET_TYPE_PROP(TwoPTwoC, FluxVariables, TwoPTwoCFluxVariables<TypeTag>);
//! define the base flux variables to realize Darcy flow
-SET_TYPE_PROP(TwoPTwoC, BaseFluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(TwoPTwoC, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
//! the upwind weight for the mass conservation equations.
SET_SCALAR_PROP(TwoPTwoC, ImplicitMassUpwindWeight, 1.0);
diff --git a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
index 0dba379ddb530eb924c549a1678f3b84391d77d1..958a06f317c20e6ddb06e846bf87dd9f18f866c7 100644
--- a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
+++ b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
@@ -27,7 +27,7 @@
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include "2pdfmproperties.hh"
@@ -45,7 +45,7 @@ namespace Dumux
* the intergration point, etc.
*/
template <class TypeTag>
-class TwoPDFMFluxVariables : public BoxDarcyFluxVariables<TypeTag>
+class TwoPDFMFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
@@ -86,7 +86,7 @@ public:
int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
- : BoxDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary),
+ : ImplicitDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary),
fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
{
faceSCV_ = &this->face();
diff --git a/dumux/implicit/2pdfm/2pdfmpropertydefaults.hh b/dumux/implicit/2pdfm/2pdfmpropertydefaults.hh
index 1af745335d0c1638e22aed0a25dc58096b2566c7..0ea4c8d9bdeaa7ecbbc493584bf1dea1f214ba62 100644
--- a/dumux/implicit/2pdfm/2pdfmpropertydefaults.hh
+++ b/dumux/implicit/2pdfm/2pdfmpropertydefaults.hh
@@ -27,7 +27,7 @@
#ifndef DUMUX_BOXMODELS_2PDFM_PROPERTY_DEFAULTS_HH
#define DUMUX_BOXMODELS_2PDFM_PROPERTY_DEFAULTS_HH
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include <dumux/material/components/nullcomponent.hh>
#include <dumux/material/fluidstates/immisciblefluidstate.hh>
#include <dumux/material/fluidsystems/2pimmisciblefluidsystem.hh>
diff --git a/dumux/implicit/2pni/2pnifluxvariables.hh b/dumux/implicit/2pni/2pnifluxvariables.hh
index 2d331d6a89e7d3d7b7b53d4ad0fe70d90e4c6ef5..ca7496b340d5d17744e2901b3ee097361ddae110 100644
--- a/dumux/implicit/2pni/2pnifluxvariables.hh
+++ b/dumux/implicit/2pni/2pnifluxvariables.hh
@@ -29,7 +29,7 @@
#define DUMUX_2PNI_FLUX_VARIABLES_HH
#include <dumux/common/math.hh>
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
namespace Dumux
{
@@ -45,9 +45,9 @@ namespace Dumux
* the integration point, etc.
*/
template <class TypeTag>
-class TwoPNIFluxVariables : public BoxDarcyFluxVariables<TypeTag>
+class TwoPNIFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
{
- typedef BoxDarcyFluxVariables<TypeTag> ParentType;
+ typedef ImplicitDarcyFluxVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
diff --git a/dumux/implicit/3p3c/3p3cpropertydefaults.hh b/dumux/implicit/3p3c/3p3cpropertydefaults.hh
index a644c778c93a1d9fef2dfe6df7b2b7fe899fb891..4309856760af22d5e5169cc483b773c184b6292d 100644
--- a/dumux/implicit/3p3c/3p3cpropertydefaults.hh
+++ b/dumux/implicit/3p3c/3p3cpropertydefaults.hh
@@ -38,7 +38,7 @@
#include "3p3cnewtoncontroller.hh"
// #include "3p3cboundaryvariables.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include <dumux/material/spatialparams/boxspatialparams.hh>
namespace Dumux
@@ -108,7 +108,7 @@ SET_TYPE_PROP(ThreePThreeC, VolumeVariables, ThreePThreeCVolumeVariables<TypeTag
SET_TYPE_PROP(ThreePThreeC, FluxVariables, ThreePThreeCFluxVariables<TypeTag>);
//! define the base flux variables to realize Darcy flow
-SET_TYPE_PROP(ThreePThreeC, BaseFluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(ThreePThreeC, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
//! the upwind factor for the mobility.
SET_SCALAR_PROP(ThreePThreeC, ImplicitMassUpwindWeight, 1.0);
diff --git a/dumux/implicit/box/boxdarcyfluxvariables.hh b/dumux/implicit/box/boxdarcyfluxvariables.hh
deleted file mode 100644
index e950b5323ec5eea0ed052e5473f033e73b590843..0000000000000000000000000000000000000000
--- a/dumux/implicit/box/boxdarcyfluxvariables.hh
+++ /dev/null
@@ -1,324 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * 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 This file contains the data which is required to calculate
- * all fluxes of fluid phases over a face of a finite volume.
- *
- * This means pressure and temperature gradients, phase densities at
- * the integration point, etc.
- */
-#ifndef DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
-#define DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
-
-#include "boxproperties.hh"
-
-#include <dumux/common/parameters.hh>
-#include <dumux/common/math.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-NEW_PROP_TAG(ImplicitMobilityUpwindWeight);
-NEW_PROP_TAG(SpatialParams);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(ProblemEnableGravity);
-}
-
-/*!
- * \ingroup BoxModel
- * \ingroup BoxFluxVariables
- * \brief Evaluates the normal component of the Darcy velocity
- * on a (sub)control volume face.
- */
-template <class TypeTag>
-class BoxDarcyFluxVariables
-{
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- enum { dim = GridView::dimension} ;
- enum { dimWorld = GridView::dimensionworld} ;
- enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
- typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-
-public:
- /*
- * \brief The constructor
- *
- * \param problem The problem
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry in the box scheme
- * \param faceIdx The local index of the SCV (sub-control-volume) face
- * \param elemVolVars The volume variables of the current element
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- */
- BoxDarcyFluxVariables(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int faceIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
- {
- mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
- calculateGradients_(problem, element, elemVolVars);
- calculateNormalVelocity_(problem, element, elemVolVars);
- }
-
-public:
- /*!
- * \brief Return the volumetric flux over a face of a given phase.
- *
- * This is the calculated velocity multiplied by the unit normal
- * and the area of the face.
- * face().normal
- * has already the magnitude of the area.
- *
- * \param phaseIdx index of the phase
- */
- Scalar volumeFlux(const unsigned int phaseIdx) const
- { return volumeFlux_[phaseIdx]; }
-
- /*!
- * \brief Return the velocity of a given phase.
- *
- * This is the full velocity vector on the
- * face (without being multiplied with normal).
- *
- * \param phaseIdx index of the phase
- */
- DimVector velocity(const unsigned int phaseIdx) const
- { return velocity_[phaseIdx] ; }
-
- /*!
- * \brief Return intrinsic permeability multiplied with potential
- * gradient multiplied with normal.
- * I.e. everything that does not need upwind decisions.
- *
- * \param phaseIdx index of the phase
- */
- Scalar kGradPNormal(const unsigned int phaseIdx) const
- { return kGradPNormal_[phaseIdx] ; }
-
- /*!
- * \brief Return the local index of the downstream control volume
- * for a given phase.
- *
- * \param phaseIdx index of the phase
- */
- const unsigned int downstreamIdx(const unsigned phaseIdx) const
- { return downstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the local index of the upstream control volume
- * for a given phase.
- *
- * \param phaseIdx index of the phase
- */
- const unsigned int upstreamIdx(const unsigned phaseIdx) const
- { return upstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the SCV (sub-control-volume) face. This may be either
- * a face within the element or a face on the element boundary,
- * depending on the value of onBoundary_.
- */
- const SCVFace &face() const
- {
- if (onBoundary_)
- return fvGeometry_.boundaryFace[faceIdx_];
- else
- return fvGeometry_.subContVolFace[faceIdx_];
- }
-
-protected:
-
- /*
- * \brief Calculation of the potential gradients
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- * are calculated for interior SCV faces or boundary faces, default=false
- */
- void calculateGradients_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- gradPotential_[phaseIdx]= 0.0 ;
- for (int idx = 0;
- idx < fvGeometry_.numFAP;
- idx++) // loop over adjacent vertices
- {
- // FE gradient at vertex idx
- const DimVector &feGrad = face().grad[idx];
-
- // index for the element volume variables
- int volVarsIdx = face().fapIndices[idx];
-
- // the pressure gradient
- DimVector tmp(feGrad);
- tmp *= elemVolVars[volVarsIdx].fluidState().pressure(phaseIdx);
- gradPotential_[phaseIdx] += tmp;
- }
-
- // correct the pressure gradient by the gravitational acceleration
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
- {
- // estimate the gravitational acceleration at a given SCV face
- // using the arithmetic mean
- DimVector g(problem.boxGravity(element, fvGeometry_, face().i));
- g += problem.boxGravity(element, fvGeometry_, face().j);
- g /= 2;
-
- // calculate the phase density at the integration point. we
- // only do this if the wetting phase is present in both cells
- Scalar SI = elemVolVars[face().i].fluidState().saturation(phaseIdx);
- Scalar SJ = elemVolVars[face().j].fluidState().saturation(phaseIdx);
- Scalar rhoI = elemVolVars[face().i].fluidState().density(phaseIdx);
- Scalar rhoJ = elemVolVars[face().j].fluidState().density(phaseIdx);
- Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
- Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
- if (fI + fJ == 0)
- // doesn't matter because no wetting phase is present in
- // both cells!
- fI = fJ = 0.5;
- const Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
-
- // make gravity acceleration a force
- DimVector f(g);
- f *= density;
-
- // calculate the final potential gradient
- gradPotential_[phaseIdx] -= f;
- } // gravity
- } // loop over all phases
- }
-
- /*
- * \brief Actual calculation of the normal Darcy velocities.
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- */
- void calculateNormalVelocity_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // calculate the mean intrinsic permeability
- const SpatialParams &spatialParams = problem.spatialParams();
- Tensor K;
- if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
- {
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- fvGeometry_,
- face().i),
- spatialParams.intrinsicPermeability(element,
- fvGeometry_,
- face().j));
- }
- else
- {
- spatialParams.meanK(K,
- spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
- spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
- }
-
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- // calculate the flux in the normal direction of the
- // current sub control volume face:
- //
- // v = - (K_f grad phi) * n
- // with K_f = rho g / mu K
- //
- // Mind, that the normal has the length of it's area.
- // This means that we are actually calculating
- // Q = - (K grad phi) dot n /|n| * A
-
-
- K.mv(gradPotential_[phaseIdx], kGradP_[phaseIdx]);
- kGradPNormal_[phaseIdx] = kGradP_[phaseIdx]*face().normal;
-
- // determine the upwind direction
- if (kGradPNormal_[phaseIdx] < 0)
- {
- upstreamIdx_[phaseIdx] = face().i;
- downstreamIdx_[phaseIdx] = face().j;
- }
- else
- {
- upstreamIdx_[phaseIdx] = face().j;
- downstreamIdx_[phaseIdx] = face().i;
- }
-
- // obtain the upwind volume variables
- const VolumeVariables& upVolVars = elemVolVars[ upstreamIdx(phaseIdx) ];
- const VolumeVariables& downVolVars = elemVolVars[ downstreamIdx(phaseIdx) ];
-
- // the minus comes from the Darcy relation which states that
- // the flux is from high to low potentials.
- // set the velocity
- velocity_[phaseIdx] = kGradP_[phaseIdx];
- velocity_[phaseIdx] *= - ( mobilityUpwindWeight_*upVolVars.mobility(phaseIdx)
- + (1.0 - mobilityUpwindWeight_)*downVolVars.mobility(phaseIdx)) ;
-
- // set the volume flux
- volumeFlux_[phaseIdx] = velocity_[phaseIdx] * face().normal ;
- }// loop all phases
- }
-
- const FVElementGeometry &fvGeometry_; //!< Information about the geometry of discretization
- const unsigned int faceIdx_; //!< The index of the sub control volume face
- const bool onBoundary_; //!< Specifying whether we are currently on the boundary of the simulation domain
- unsigned int upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
- Scalar volumeFlux_[numPhases] ; //!< Velocity multiplied with normal (magnitude=area)
- DimVector velocity_[numPhases] ; //!< The velocity as determined by Darcy's law or by the Forchheimer relation
- Scalar kGradPNormal_[numPhases] ; //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
- DimVector kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
- DimVector gradPotential_[numPhases] ; //!< Gradient of potential, which drives flow
- Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/box/boxforchheimerfluxvariables.hh b/dumux/implicit/box/boxforchheimerfluxvariables.hh
deleted file mode 100644
index 78d008f20316903f1c6235e14cbc7bea4d4c3fd5..0000000000000000000000000000000000000000
--- a/dumux/implicit/box/boxforchheimerfluxvariables.hh
+++ /dev/null
@@ -1,367 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * 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 This file contains the data which is required to calculate
- * all fluxes of fluid phases over a face of a finite volume,
- * according to the Forchheimer-relation between velocity and pressure.
- *
- */
-#ifndef DUMUX_BOX_FORCHHEIMER_FLUX_VARIABLES_HH
-#define DUMUX_BOX_FORCHHEIMER_FLUX_VARIABLES_HH
-
-#include "boxproperties.hh"
-
-#include <dumux/common/parameters.hh>
-#include <dumux/common/math.hh>
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-NEW_PROP_TAG(MobilityUpwindWeight);
-NEW_PROP_TAG(SpatialParams);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(ProblemEnableGravity);
-}
-
-/*!
- * \ingroup BoxModel
- * \ingroup BoxFluxVariables
- * \brief Evaluates the normal component of the Forchheimer velocity
- * on a (sub)control volume face.
- *
- * The commonly used Darcy relation looses it's validity for \f$ Re < 1\f$.
- * If one encounters flow velocities in porous media above this Reynolds number,
- * the Forchheimer relation can be used. Like the Darcy relation, it relates
- * the gradient in potential to velocity.
- * However, this relation is not linear (as in the Darcy case) any more.
- *
- * Therefore, a Newton scheme is implemented in this class in order to calculate a
- * velocity from the current set of variables. This velocity can subsequently be used
- * e.g. by the localresidual.
- *
- * For Reynolds numbers above \f$ 500\f$ the (Standard) forchheimer relation also
- * looses it's validity.
- *
- * The Forchheimer equation looks as follows:
- * \f[ \nabla \left({p_\alpha + \rho_\alpha g z} \right)= - \frac{\mu_\alpha}{k_{r \alpha} K} \underline{v_\alpha} - \frac{c_F}{\eta_{\alpha r} \sqrt{K}} \rho_\alpha |\underline{v_\alpha}| \underline{v_\alpha} \f]
- *
- * For the formulation that is actually used in this class, see the documentation of the function calculating the residual.
- *
- * - This algorithm does not find a solution if the fluid is incompressible and the
- * initial pressure distribution is uniform.
- * - This algorithm assumes the relative passabilities to have the same functional
- * form as the relative permeabilities.
- */
-template <class TypeTag>
-class BoxForchheimerFluxVariables
- : public BoxDarcyFluxVariables<TypeTag>
-{
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- enum { dim = GridView::dimension} ;
- enum { dimWorld = GridView::dimensionworld} ;
- enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
- typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-
-public:
- /*!
- * \brief The constructor
- *
- * \param problem The problem
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry in the box scheme
- * \param faceIdx The local index of the SCV (sub-control-volume) face
- * \param elemVolVars The volume variables of the current element
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- */
- BoxForchheimerFluxVariables(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const unsigned int faceIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- : BoxDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
- {
- calculateNormalVelocity_(problem, element, elemVolVars);
- }
-
-protected:
- /*!
- * \brief Function for calculation of velocities.
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- */
- void calculateNormalVelocity_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // calculate the mean intrinsic permeability
- const SpatialParams &spatialParams = problem.spatialParams();
- Tensor K;
- if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
- {
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- fvGeometry_,
- face().i),
- spatialParams.intrinsicPermeability(element,
- fvGeometry_,
- face().j));
- }
- else
- {
- spatialParams.meanK(K,
- spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
- spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
- }
-
- // obtain the Forchheimer coefficient from the spatial parameters
- const Scalar forchCoeff = spatialParams.forchCoeff(element,
- this->fvGeometry_,
- this->face().i);
-
- // Make sure the permeability matrix does not have off-diagonal entries
- assert( isDiagonal_(K) );
-
- Tensor sqrtK(0.0);
- for (int i = 0; i < dim; ++i)
- sqrtK[i][i] = std::sqrt(K[i][i]);
-
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- // initial guess of velocity: Darcy relation
- // first taken from base class, later overwritten in base class
- DimVector velocity = this-> velocity(phaseIdx);
-
- DimVector deltaV; // the change in velocity between Newton iterations
- DimVector residual(10e10); // the residual (function value that is to be minimized ) of the equation that is to be fulfilled
- DimVector tmp; // temporary variable for numerical differentiation
- Tensor gradF; // slope of equation that is to be solved
-
- // search by means of the Newton method for a root of Forchheimer equation
- for (int k = 0; residual.two_norm() > 1e-12 ; ++k) {
-
- if (k >= 30)
- DUNE_THROW(NumericalProblem, "could not determine forchheimer velocity within "<< k <<" iterations");
-
- // calculate current value of Forchheimer equation
- forchheimerResidual_(residual,
- forchCoeff,
- sqrtK,
- K,
- velocity,
- elemVolVars,
- this->gradPotential_[phaseIdx],
- phaseIdx);
-
- // newton's method: slope (gradF) and current value (residual) of function is known,
- forchheimerDerivative_(gradF,
- forchCoeff,
- sqrtK,
- velocity,
- elemVolVars,
- phaseIdx);
-
- // solve for change in velocity ("x-Axis intercept")
- gradF.solve(deltaV, residual);
- velocity -= deltaV;
- }
-
- this->velocity_[phaseIdx] = velocity ; // store the converged velocity solution
- this->volumeFlux_[phaseIdx] = velocity * this->face().normal; // velocity dot normal times cross sectional area is volume flux
- }// end loop over all phases
- }
-
- /*! \brief Function for calculation of Forchheimer relation.
- *
- * The relative passability \f$ \eta_r\f$ is the "Forchheimer-equivalent" to the relative permeability \f$ k_r\f$.
- * We use the same function as for \f$ k_r \f$ (VG, BC, linear) other authors use a simple power law e.g.: \f$\eta_{rw} = S_w^3\f$
- *
- * Some rearrangements have been made to the original Forchheimer relation:
- *
- * \f[ \underline{v_\alpha} + c_F \sqrt{K} \frac{\rho_\alpha}{\mu_\alpha } |\underline{v_\alpha}| \underline{v_\alpha} + \frac{k_{r \alpha}}{\mu_\alpha} K \nabla \left(p_\alpha + \rho_\alpha g z \right)= 0 \f]
- *
- * This already includes the assumption \f$ k_r(S_w) = \eta_r(S_w)\f$:
- * - \f$\eta_{rw} = S_w^x\f$ looks very similar to e.g. Van Genuchten relative permeabilities
- * - Fichot, et al. (2006), Nuclear Engineering and Design, state that several authors claim that \f$ k_r(S_w), \eta_r(S_w)\f$ can be chosen equal
- * - It leads to the equation not degenerating for the case of \f$S_w=1\f$, because I do not need to multiply with two different functions, and therefore there are terms not being zero.
- * - If this assumption is not to be made: Some regularization needs to be introduced ensuring that not all terms become zero for \f$S_w=1\f$.
- *
- * As long as the correct velocity is not found yet, the above equation is not fulfilled, but
- * the left hand side yields some residual. This function calculates the left hand side of the
- * equation and returns the residual.
- *
- * \param residual The current function value for the given velocity
- * \param forchCoeff The Forchheimer coefficient
- * \param sqrtK A diagonal matrix whose entries are square roots of the permeability matrix entries
- * \param K The permeability matrix
- * \param velocity The current velocity approximation.
- * \param elemVolVars The volume variables of the current element
- * \param gradPotential The gradient in potential
- * \param phaseIdx The index of the currently considered phase
- */
- void forchheimerResidual_(DimVector & residual,
- const Scalar forchCoeff,
- const Tensor & sqrtK,
- const Tensor & K,
- const DimVector & velocity,
- const ElementVolumeVariables & elemVolVars,
- const DimVector & gradPotential ,
- const unsigned int phaseIdx) const
- {
- const VolumeVariables upVolVars = elemVolVars[this->upstreamIdx(phaseIdx)];
- const VolumeVariables downVolVars = elemVolVars[this->downstreamIdx(phaseIdx)];
-
- // Obtaining the physical quantities
- const Scalar mobility = this->mobilityUpwindWeight_ * upVolVars.mobility(phaseIdx)
- + (1.-this->mobilityUpwindWeight_)* downVolVars.mobility(phaseIdx) ;
-
- const Scalar viscosity = this->mobilityUpwindWeight_ * upVolVars.fluidState().viscosity(phaseIdx)
- + (1.-this->mobilityUpwindWeight_)* downVolVars.fluidState().viscosity(phaseIdx) ;
-
- const Scalar density = this->mobilityUpwindWeight_ * upVolVars.fluidState().density(phaseIdx)
- + (1.-this->mobilityUpwindWeight_) * downVolVars.fluidState().density(phaseIdx) ;
-
- // residual = v
- residual = velocity;
-
- // residual += k_r/mu K grad p
- K.usmv(mobility, gradPotential, residual);
-
- // residual += c_F rho / mu abs(v) sqrt(K) v
- sqrtK.usmv(forchCoeff * density / viscosity * velocity.two_norm(), velocity, residual);
- }
-
-
- /*! \brief Function for calculation of the gradient of Forchheimer relation with respect to velocity.
- *
- * This function already exploits that \f$ \sqrt{K}\f$ is a diagonal matrix. Therefore, we only have
- * to deal with the main entries.
- * The gradient of the Forchheimer relations looks as follows (mind that \f$ \sqrt{K}\f$ is a tensor):
- *
- * \f[ f\left(\underline{v_\alpha}\right) =
- * \left(
- * \begin{array}{ccc}
- * 1 & 0 &0 \\
- * 0 & 1 &0 \\
- * 0 & 0 &1 \\
- * \end{array}
- * \right)
- * +
- * c_F \frac{\rho_\alpha}{\mu_\alpha} |\underline{v}_\alpha| \sqrt{K}
- * +
- * c_F \frac{\rho_\alpha}{\mu_\alpha}\frac{1}{|\underline{v}_\alpha|} \sqrt{K}
- * \left(
- * \begin{array}{ccc}
- * v_x^2 & v_xv_y & v_xv_z \\
- * v_yv_x & v_{y}^2 & v_yv_z \\
- * v_zv_x & v_zv_y &v_{z}^2 \\
- * \end{array}
- * \right)
- * \f]
- *
- * \param derivative The gradient of the forchheimer equation
- * \param forchCoeff Forchheimer coefficient
- * \param sqrtK A diagonal matrix whose entries are square roots of the permeability matrix entries.
- * \param velocity The current velocity approximation.
- * \param elemVolVars The volume variables of the current element
- * \param phaseIdx The index of the currently considered phase
- */
- void forchheimerDerivative_(Tensor & derivative,
- const Scalar forchCoeff,
- const Tensor & sqrtK,
- const DimVector & velocity,
- const ElementVolumeVariables & elemVolVars,
- const unsigned int phaseIdx) const
- {
- const VolumeVariables upVolVars = elemVolVars[this->upstreamIdx(phaseIdx)];
- const VolumeVariables downVolVars = elemVolVars[this->downstreamIdx(phaseIdx)];
-
- // Obtaining the physical quantities
- const Scalar viscosity = this->mobilityUpwindWeight_ * upVolVars.fluidState().viscosity(phaseIdx)
- + (1.-this->mobilityUpwindWeight_)* downVolVars.fluidState().viscosity(phaseIdx) ;
-
- const Scalar density = this->mobilityUpwindWeight_ * upVolVars.fluidState().density(phaseIdx)
- + (1.-this->mobilityUpwindWeight_) * downVolVars.fluidState().density(phaseIdx) ;
-
- // Initialize because for low velocities, we add and do not set the entries.
- derivative = 0.;
-
- const Scalar absV = velocity.two_norm() ;
- // This part of the derivative is only calculated if the velocity is sufficiently small: do not divide by zero!
- // This should not be very bad because the derivative is only intended to give an approximation of the gradient of the
- // function that goes into the Newton scheme.
- // This is important in the case of a one-phase region in two-phase flow. The non-existing phase has a velocity of zero (kr=0).
- // f = sqrtK * vTimesV (this is a matrix)
- // f *= forchCoeff density / |velocity| / viscosity (multiply all entries with scalars)
- if(absV > 1e-20)
- for (int i=0; i<dim; i++)
- for (int k=0; k<dim; k++)
- derivative[i][k]= sqrtK[i][i] * velocity[i] * velocity[k] * forchCoeff * density / absV / viscosity;
-
- // add on the main diagonal:
- // 1 + sqrtK_i forchCoeff density |v| / viscosity
- for (int i=0; i<dim; i++)
- derivative[i][i] += (1.0 + ( sqrtK[i][i]*forchCoeff * density * absV / viscosity ) ) ;
- }
-
- /*!
- * \brief Check whether all off-diagonal entries of a tensor are zero.
- *
- * \param K the tensor that is to be checked.
- *
- * \return True iff all off-diagonals are zero.
- *
- */
- const bool isDiagonal_(const Tensor & K) const
- {
- for (int i = 0; i < dim; i++) {
- for (int k = 0; k < dim; k++) {
- if ((i != k) && (std::abs(K[i][k]) >= 1e-25)) {
- return false;
- }
- }
- }
- return true;
- }
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/common/implicitdarcyfluxvariables.hh b/dumux/implicit/common/implicitdarcyfluxvariables.hh
index 6e5d3df449f65a22a12966497c6a5a2c0a6e800c..9507d03a895b2bdea7db71dbabcd84d9b7189d71 100644
--- a/dumux/implicit/common/implicitdarcyfluxvariables.hh
+++ b/dumux/implicit/common/implicitdarcyfluxvariables.hh
@@ -25,10 +25,10 @@
* This means pressure and temperature gradients, phase densities at
* the integration point, etc.
*/
-#ifndef DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
-#define DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
+#ifndef DUMUX_IMPLICIT_DARCY_FLUX_VARIABLES_HH
+#define DUMUX_IMPLICIT_DARCY_FLUX_VARIABLES_HH
-#include "boxproperties.hh"
+#include "implicitproperties.hh"
#include <dumux/common/parameters.hh>
#include <dumux/common/math.hh>
@@ -46,13 +46,13 @@ NEW_PROP_TAG(ProblemEnableGravity);
}
/*!
- * \ingroup BoxModel
- * \ingroup BoxFluxVariables
+ * \ingroup ImplicitModel
+ * \ingroup ImplicitFluxVariables
* \brief Evaluates the normal component of the Darcy velocity
* on a (sub)control volume face.
*/
template <class TypeTag>
-class BoxDarcyFluxVariables
+class ImplicitDarcyFluxVariables
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
@@ -85,13 +85,13 @@ public:
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
- BoxDarcyFluxVariables(const Problem &problem,
+ ImplicitDarcyFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
- : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
+ : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
{
mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
calculateGradients_(problem, element, elemVolVars);
@@ -165,17 +165,7 @@ public:
}
protected:
- const FVElementGeometry &fvGeometry_; //!< Information about the geometry of discretization
- const unsigned int faceIdx_; //!< The index of the sub control volume face
- const bool onBoundary_; //!< Specifying whether we are currently on the boundary of the simulation domain
- unsigned int upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
- Scalar volumeFlux_[numPhases] ; //!< Velocity multiplied with normal (magnitude=area)
- DimVector velocity_[numPhases] ; //!< The velocity as determined by Darcy's law or by the Forchheimer relation
- Scalar kGradPNormal_[numPhases] ; //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
- DimVector kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
- DimVector gradPotential_[numPhases] ; //!< Gradient of potential, which drives flow
- Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
-
+
/*
* \brief Calculation of the potential gradients
*
@@ -255,13 +245,31 @@ protected:
// calculate the mean intrinsic permeability
const SpatialParams &spatialParams = problem.spatialParams();
Tensor K;
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- fvGeometry_,
- face().i),
- spatialParams.intrinsicPermeability(element,
- fvGeometry_,
- face().j));
+ if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+ {
+ spatialParams.meanK(K,
+ spatialParams.intrinsicPermeability(element,
+ fvGeometry_,
+ face().i),
+ spatialParams.intrinsicPermeability(element,
+ fvGeometry_,
+ face().j));
+ }
+ else
+ {
+ const Element& elementI = *fvGeometry_.neighbors[face().i];
+ FVElementGeometry fvGeometryI;
+ fvGeometryI.subContVol[0].global = elementI.geometry().center();
+
+ const Element& elementJ = *fvGeometry_.neighbors[face().j];
+ FVElementGeometry fvGeometryJ;
+ fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
+
+ spatialParams.meanK(K,
+ spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0),
+ spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0));
+ }
+
// loop over all phases
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
{
@@ -306,6 +314,17 @@ protected:
volumeFlux_[phaseIdx] = velocity_[phaseIdx] * face().normal ;
}// loop all phases
}
+
+ const FVElementGeometry &fvGeometry_; //!< Information about the geometry of discretization
+ const unsigned int faceIdx_; //!< The index of the sub control volume face
+ const bool onBoundary_; //!< Specifying whether we are currently on the boundary of the simulation domain
+ unsigned int upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
+ Scalar volumeFlux_[numPhases] ; //!< Velocity multiplied with normal (magnitude=area)
+ DimVector velocity_[numPhases] ; //!< The velocity as determined by Darcy's law or by the Forchheimer relation
+ Scalar kGradPNormal_[numPhases] ; //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
+ DimVector kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
+ DimVector gradPotential_[numPhases] ; //!< Gradient of potential, which drives flow
+ Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
};
} // end namespace
diff --git a/dumux/implicit/common/implicitforchheimerfluxvariables.hh b/dumux/implicit/common/implicitforchheimerfluxvariables.hh
index ec3117cb08f9215b4d3048aaec4f665abbb11de9..9956bae633fb9d6ecb94dc1c5d28f875c73e45a1 100644
--- a/dumux/implicit/common/implicitforchheimerfluxvariables.hh
+++ b/dumux/implicit/common/implicitforchheimerfluxvariables.hh
@@ -24,14 +24,14 @@
* according to the Forchheimer-relation between velocity and pressure.
*
*/
-#ifndef DUMUX_BOX_FORCHHEIMER_FLUX_VARIABLES_HH
-#define DUMUX_BOX_FORCHHEIMER_FLUX_VARIABLES_HH
+#ifndef DUMUX_IMPLICIT_FORCHHEIMER_FLUX_VARIABLES_HH
+#define DUMUX_IMPLICIT_FORCHHEIMER_FLUX_VARIABLES_HH
-#include "boxproperties.hh"
+#include "implicitproperties.hh"
#include <dumux/common/parameters.hh>
#include <dumux/common/math.hh>
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
namespace Dumux
{
@@ -46,8 +46,8 @@ NEW_PROP_TAG(ProblemEnableGravity);
}
/*!
- * \ingroup BoxModel
- * \ingroup BoxFluxVariables
+ * \ingroup ImplicitModel
+ * \ingroup ImplicitFluxVariables
* \brief Evaluates the normal component of the Forchheimer velocity
* on a (sub)control volume face.
*
@@ -75,8 +75,8 @@ NEW_PROP_TAG(ProblemEnableGravity);
* form as the relative permeabilities.
*/
template <class TypeTag>
-class BoxForchheimerFluxVariables
- : public BoxDarcyFluxVariables<TypeTag>
+class ImplicitForchheimerFluxVariables
+ : public ImplicitDarcyFluxVariables<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
@@ -107,13 +107,13 @@ public:
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
- BoxForchheimerFluxVariables(const Problem &problem,
+ ImplicitForchheimerFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const unsigned int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
- : BoxDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+ : ImplicitDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
{
calculateNormalVelocity_(problem, element, elemVolVars);
}
@@ -133,14 +133,31 @@ protected:
// calculate the mean intrinsic permeability
const SpatialParams &spatialParams = problem.spatialParams();
Tensor K;
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- this->fvGeometry_,
- this->face().i),
- spatialParams.intrinsicPermeability(element,
- this->fvGeometry_,
- this->face().j));
-
+ if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+ {
+ spatialParams.meanK(K,
+ spatialParams.intrinsicPermeability(element,
+ this->fvGeometry_,
+ this->face().i),
+ spatialParams.intrinsicPermeability(element,
+ this->fvGeometry_,
+ this->face().j));
+ }
+ else
+ {
+ const Element& elementI = *this->fvGeometry_.neighbors[this->face().i];
+ FVElementGeometry fvGeometryI;
+ fvGeometryI.subContVol[0].global = elementI.geometry().center();
+
+ const Element& elementJ = *this->fvGeometry_.neighbors[this->face().j];
+ FVElementGeometry fvGeometryJ;
+ fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
+
+ spatialParams.meanK(K,
+ spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0),
+ spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0));
+ }
+
// obtain the Forchheimer coefficient from the spatial parameters
const Scalar forchCoeff = spatialParams.forchCoeff(element,
this->fvGeometry_,
diff --git a/dumux/implicit/mpnc/mpncfluxvariables.hh b/dumux/implicit/mpnc/mpncfluxvariables.hh
index b557dcd1437c733b7282f8aca85afeea21355685..ac913caea42ffb83ecdbc0a6865ecc7154317d40 100644
--- a/dumux/implicit/mpnc/mpncfluxvariables.hh
+++ b/dumux/implicit/mpnc/mpncfluxvariables.hh
@@ -32,8 +32,8 @@
#include "diffusion/fluxvariables.hh"
#include "energy/mpncfluxvariablesenergy.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
-#include <dumux/implicit/box/boxforchheimerfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitforchheimerfluxvariables.hh>
namespace Dumux
{
diff --git a/dumux/implicit/mpnc/mpncmodel.hh b/dumux/implicit/mpnc/mpncmodel.hh
index 18f9cd2e71d43790a9f02a326b3dc7ef2eeb42a6..5b31ce887164a5bc595a869303c3836f0d5f03d0 100644
--- a/dumux/implicit/mpnc/mpncmodel.hh
+++ b/dumux/implicit/mpnc/mpncmodel.hh
@@ -38,7 +38,7 @@ namespace Dumux
* denoted by the upper index \f$\kappa \in \{ 1, \dots, N \} \f$.
*
* The momentum approximation can be selected via "BaseFluxVariables":
- * Darcy (BoxDarcyFluxVariables) and Forchheimer (BoxForchheimerFluxVariables)
+ * Darcy (ImplicitDarcyFluxVariables) and Forchheimer (ImplicitForchheimerFluxVariables)
* relations are available for all Box models.
*
* By inserting this into the equations for the conservation of the
diff --git a/dumux/implicit/mpnc/mpncpropertydefaults.hh b/dumux/implicit/mpnc/mpncpropertydefaults.hh
index 00f8c33b5f409fc3b59c4f4708b807ba517e1245..cd1f94ae84b76ae16fe7ab8128f3d3b1443134a8 100644
--- a/dumux/implicit/mpnc/mpncpropertydefaults.hh
+++ b/dumux/implicit/mpnc/mpncpropertydefaults.hh
@@ -173,7 +173,7 @@ SET_TYPE_PROP(MPNC, VolumeVariables, MPNCVolumeVariables<TypeTag>);
SET_TYPE_PROP(MPNC, FluxVariables, MPNCFluxVariables<TypeTag>);
//! the Base of the Fluxvariables property (Darcy / Forchheimer)
-SET_TYPE_PROP(MPNC, BaseFluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(MPNC, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
//! truncate the newton update for the first few Newton iterations of a time step
SET_BOOL_PROP(MPNC, NewtonEnableChop, true);
diff --git a/dumux/implicit/richards/richardspropertydefaults.hh b/dumux/implicit/richards/richardspropertydefaults.hh
index 2e312d08051bbb19dec0631e18087b387804a5af..2766eb2ffff0247b15ffb960c73fd62364772930 100644
--- a/dumux/implicit/richards/richardspropertydefaults.hh
+++ b/dumux/implicit/richards/richardspropertydefaults.hh
@@ -31,7 +31,7 @@
#include "richardsmodel.hh"
#include "richardsproblem.hh"
#include "richardsindices.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include "richardsvolumevariables.hh"
#include "richardsproperties.hh"
#include "richardsnewtoncontroller.hh"
@@ -65,7 +65,7 @@ SET_TYPE_PROP(Richards, Model, RichardsModel<TypeTag>);
SET_TYPE_PROP(Richards, VolumeVariables, RichardsVolumeVariables<TypeTag>);
//! The class for the quantities required for the flux calculation
-SET_TYPE_PROP(Richards, FluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(Richards, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
//! The class of the newton controller
SET_TYPE_PROP(Richards, NewtonController, RichardsNewtonController<TypeTag>);
diff --git a/dumux/material/spatialparams/boxspatialparams1p.hh b/dumux/material/spatialparams/boxspatialparams1p.hh
index b1c6e0e7060ac44a42eb260f29110408694dd288..bb73bcfd7b32ac89a5606949faa3e097a8a058b3 100644
--- a/dumux/material/spatialparams/boxspatialparams1p.hh
+++ b/dumux/material/spatialparams/boxspatialparams1p.hh
@@ -107,13 +107,6 @@ public:
result[i][j] = harmonicMean(K1[i][j], K2[i][j]);
}
- Scalar elementIntrinsicPermeability (const Element &element) const
- {
- FVElementGeometry fvGeometry;
- fvGeometry.subContVol[0].global = element.geometry().center();
- return asImp_().intrinsicPermeability(element, fvGeometry, /*scvIdx=*/0);
- }
-
/*!
* \brief Function for defining the intrinsic (absolute) permeability.
*
diff --git a/test/implicit/mpnc/forchheimer1pproblem.hh b/test/implicit/mpnc/forchheimer1pproblem.hh
index 638ddae5c1517dd5f63053589bd80aaabeb6b36a..c0aa59f82409d9f819198da7ba67d1809ac36676 100644
--- a/test/implicit/mpnc/forchheimer1pproblem.hh
+++ b/test/implicit/mpnc/forchheimer1pproblem.hh
@@ -98,7 +98,7 @@ SET_INT_PROP(Forchheimer1pProblem, ImplicitNumericDifferenceMethod, +1);
SET_TYPE_PROP(Forchheimer1pProblem, Scalar, double);
// decide which to use for velocity calculation: Darcy / Forchheimer
-SET_TYPE_PROP(Forchheimer1pProblem, BaseFluxVariables, BoxForchheimerFluxVariables<TypeTag>);
+SET_TYPE_PROP(Forchheimer1pProblem, BaseFluxVariables, ImplicitForchheimerFluxVariables<TypeTag>);
SET_BOOL_PROP(Forchheimer1pProblem, VtkAddVelocities, true);
}
diff --git a/test/implicit/mpnc/forchheimer2pproblem.hh b/test/implicit/mpnc/forchheimer2pproblem.hh
index 49a926b069c1bf27dba67a104e19231c0ae924cf..5061431cafd79d2ee95a396ab8c89f5d94967f34 100644
--- a/test/implicit/mpnc/forchheimer2pproblem.hh
+++ b/test/implicit/mpnc/forchheimer2pproblem.hh
@@ -98,7 +98,7 @@ SET_INT_PROP(Forchheimer2pProblem, ImplicitNumericDifferenceMethod, +1);
SET_TYPE_PROP(Forchheimer2pProblem, Scalar, double);
// decide how to calculate velocity: Darcy / Forchheimer
-SET_TYPE_PROP(Forchheimer2pProblem, BaseFluxVariables, BoxForchheimerFluxVariables<TypeTag>);
+SET_TYPE_PROP(Forchheimer2pProblem, BaseFluxVariables, ImplicitForchheimerFluxVariables<TypeTag>);
SET_BOOL_PROP(Forchheimer2pProblem, VtkAddVelocities, true);
}