From f78e26935620f25a0ee8fce2f6a8b9ea9d5da1bd Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 17 Apr 2018 18:11:14 +0200
Subject: [PATCH] [tpfa] Free darcy's law from TypeTag
---
dumux/common/properties.hh | 4 -
.../cellcentered/tpfa/darcyslaw.hh | 140 +++++++++---------
dumux/material/spatialparams/fv1p.hh | 8 +
dumux/porousmediumflow/properties.hh | 3 -
4 files changed, 80 insertions(+), 75 deletions(-)
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 6fb5a53ff6..b6a57311eb 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -120,10 +120,6 @@ NEW_PROP_TAG(Formulation); //!< The formulation of the m
NEW_PROP_TAG(UseConstraintSolver); //!< Whether to use a contraint solver for computing the secondary variables
NEW_PROP_TAG(UseKelvinEquation); //!< If we use Kelvin equation to lower the vapor pressure as a function of capillary pressure, temperature
-// specify if we evaluate the permeability in the volume (for discontinuous fields)
-// or at the scvf center for analytical permeability fields (e.g. convergence studies)
-NEW_PROP_TAG(EvaluatePermeabilityAtScvfIP);
-
// When using the box method in a multi-phase context, an interface solver might be necessary
NEW_PROP_TAG(EnableBoxInterfaceSolver);
diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
index cd9c5bcb62..6a2393293c 100644
--- a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
@@ -31,13 +31,21 @@
#include <dumux/discretization/methods.hh>
#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
-namespace Dumux
-{
+namespace Dumux {
+
// forward declarations
template<class TypeTag, DiscretizationMethod discMethod>
class DarcysLawImplementation;
-template<class TypeTag, bool isNetwork>
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation
+ * \note Darcy's law is speialized for network and surface grids (i.e. if grid dim < dimWorld)
+ * \tparam Scalar the scalar type for scalar physical quantities
+ * \tparam FVGridGeometry the grid geometry
+ * \tparam isNetwork whether we are computing on a network grid embedded in a higher world dimension
+ */
+template<class Scalar, class FVGridGeometry, bool isNetwork>
class CCTpfaDarcysLaw;
/*!
@@ -47,27 +55,27 @@ class CCTpfaDarcysLaw;
*/
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethod::cctpfa>
-: public CCTpfaDarcysLaw<TypeTag, (GET_PROP_TYPE(TypeTag, Grid)::dimension < GET_PROP_TYPE(TypeTag, Grid)::dimensionworld) >
+: public CCTpfaDarcysLaw<typename GET_PROP_TYPE(TypeTag, Scalar),
+ typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ (GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView::dimension < GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView::dimensionworld)>
{};
/*!
* \ingroup CCTpfaDiscretization
* \brief Class that fills the cache corresponding to tpfa Darcy's Law
*/
-template<class TypeTag>
+template<class FVGridGeometry>
class TpfaDarcysLawCacheFiller
{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
- using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+ using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
public:
//! Function to fill a TpfaDarcysLawCache of a given scvf
//! This interface has to be met by any advection-related cache filler class
- template<class FluxVariablesCacheFiller>
+ //! TODO: Probably get cache type out of the filler
+ template<class FluxVariablesCache, class Problem, class ElementVolumeVariables, class FluxVariablesCacheFiller>
static void fill(FluxVariablesCache& scvfFluxVarsCache,
const Problem& problem,
const Element& element,
@@ -84,20 +92,18 @@ public:
* \ingroup CCTpfaDiscretization
* \brief The cache corresponding to tpfa Darcy's Law
*/
-template<class TypeTag>
+template<class AdvectionType, class FVGridGeometry>
class TpfaDarcysLawCache
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using Scalar = typename AdvectionType::Scalar;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+ using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
public:
- using Filler = TpfaDarcysLawCacheFiller<TypeTag>;
+ using Filler = TpfaDarcysLawCacheFiller<FVGridGeometry>;
+ template<class Problem, class ElementVolumeVariables>
void updateAdvection(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -118,37 +124,33 @@ private:
* \ingroup CCTpfaDiscretization
* \brief Specialization of the CCTpfaDarcysLaw grids where dim=dimWorld
*/
-template<class TypeTag>
-class CCTpfaDarcysLaw<TypeTag, /*isNetwork*/ false>
+template<class ScalarType, class FVGridGeometry>
+class CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
{
- using Implementation = DarcysLawImplementation<TypeTag, DiscretizationMethod::cctpfa>;
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
- using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
- using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using ThisType = CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+ using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
+ //! state the scalar type of the law
+ using Scalar = ScalarType;
+
//! state the discretization method this implementation belongs to
static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
//! state the type for the corresponding cache
- using Cache = TpfaDarcysLawCache<TypeTag>;
+ using Cache = TpfaDarcysLawCache<ThisType, FVGridGeometry>;
//! Compute the advective flux
+ template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -157,7 +159,7 @@ class CCTpfaDarcysLaw<TypeTag, /*isNetwork*/ false>
int phaseIdx,
const ElementFluxVarsCache& elemFluxVarsCache)
{
- static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
+ static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
@@ -210,6 +212,7 @@ class CCTpfaDarcysLaw<TypeTag, /*isNetwork*/ false>
// The flux variables cache has to be bound to an element prior to flux calculations
// During the binding, the transmissibility will be computed and stored using the method below.
+ template<class Problem, class ElementVolumeVariables>
static Scalar calculateTransmissibility(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -223,15 +226,17 @@ class CCTpfaDarcysLaw<TypeTag, /*isNetwork*/ false>
const auto& insideVolVars = elemVolVars[insideScvIdx];
// check if we evaluate the permeability in the volume (for discontinuous fields, default)
// or at the scvf center for analytical permeability fields (e.g. convergence studies)
+ using SpatialParams = typename Problem::SpatialParams;
+ using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
auto getPermeability = [&problem](const VolumeVariables& volVars,
const GlobalPosition& scvfIpGlobal) -> typename SpatialParams::PermeabilityType
{
- if (GET_PROP_VALUE(TypeTag, EvaluatePermeabilityAtScvfIP))
+ if (SpatialParams::evaluatePermeabilityAtScvfIP())
{
- // TODO: Make PermeabilityType a property!!
+ // If the permeability at pos method is not overloaded it might have a different return type
// We do an implicit cast to permeability type in case EvaluatePermeabilityAtScvfIP is not true so that it compiles
// If it is true make sure that no cast is necessary and the correct return type is specified in the spatial params
- static_assert(!GET_PROP_VALUE(TypeTag, EvaluatePermeabilityAtScvfIP)
+ static_assert(!SpatialParams::evaluatePermeabilityAtScvfIP()
|| std::is_same<std::decay_t<typename SpatialParams::PermeabilityType>, std::decay_t<decltype(problem.spatialParams().permeabilityAtPos(scvfIpGlobal))>>::value,
"permeabilityAtPos doesn't return PermeabilityType stated in the spatial params!");
return problem.spatialParams().permeabilityAtPos(scvfIpGlobal);
@@ -274,37 +279,33 @@ class CCTpfaDarcysLaw<TypeTag, /*isNetwork*/ false>
* \ingroup CCTpfaDiscretization
* \brief Specialization of the CCTpfaDarcysLaw grids where dim < dimWorld (network/surface grids)
*/
-template<class TypeTag>
-class CCTpfaDarcysLaw<TypeTag, /*isNetwork*/ true>
+template<class ScalarType, class FVGridGeometry>
+class CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ true>
{
- using Implementation = DarcysLawImplementation<TypeTag, DiscretizationMethod::cctpfa>;
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
- using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
- using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using ThisType = CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ true>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+ using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
+ //! state the scalar type of the law
+ using Scalar = ScalarType;
+
//! state the discretization method this implementation belongs to
static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
//! state the type for the corresponding cache
- using Cache = TpfaDarcysLawCache<TypeTag>;
+ using Cache = TpfaDarcysLawCache<ThisType, FVGridGeometry>;
//! Compute the advective flux
+ template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -313,7 +314,7 @@ public:
int phaseIdx,
const ElementFluxVarsCache& elemFluxVarsCache)
{
- static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
+ static const bool gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
@@ -445,6 +446,7 @@ public:
// The flux variables cache has to be bound to an element prior to flux calculations
// During the binding, the transmissibility will be computed and stored using the method below.
+ template<class Problem, class ElementVolumeVariables>
static Scalar calculateTransmissibility(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -458,15 +460,17 @@ public:
const auto& insideVolVars = elemVolVars[insideScvIdx];
// check if we evaluate the permeability in the volume (for discontinuous fields, default)
// or at the scvf center for analytical permeability fields (e.g. convergence studies)
+ using SpatialParams = typename Problem::SpatialParams;
+ using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
auto getPermeability = [&problem](const VolumeVariables& volVars,
const GlobalPosition& scvfIpGlobal) -> typename SpatialParams::PermeabilityType
{
- if (GET_PROP_VALUE(TypeTag, EvaluatePermeabilityAtScvfIP))
+ if (SpatialParams::evaluatePermeabilityAtScvfIP())
{
- // TODO: Make PermeabilityType a property!!
- // We do an implicit cast to permeability type in case EvaluatePermeabilityAtScvfIP is not true so that it compiles
+ // If the permeability at pos method is not overloaded it might have a different return type
+ // We do an implicit cast to permeability type in case evaluatePermeabilityAtScvfIP is not true so that it compiles
// If it is true make sure that no cast is necessary and the correct return type is specified in the spatial params
- static_assert(!GET_PROP_VALUE(TypeTag, EvaluatePermeabilityAtScvfIP)
+ static_assert(!SpatialParams::evaluatePermeabilityAtScvfIP()
|| std::is_same<std::decay_t<typename SpatialParams::PermeabilityType>, std::decay_t<decltype(problem.spatialParams().permeabilityAtPos(scvfIpGlobal))>>::value,
"permeabilityAtPos doesn't return PermeabilityType stated in the spatial params!");
return problem.spatialParams().permeabilityAtPos(scvfIpGlobal);
diff --git a/dumux/material/spatialparams/fv1p.hh b/dumux/material/spatialparams/fv1p.hh
index c72a00466d..72fa841872 100644
--- a/dumux/material/spatialparams/fv1p.hh
+++ b/dumux/material/spatialparams/fv1p.hh
@@ -143,6 +143,14 @@ public:
"a permeability() or permeabilityAtPos() method.");
}
+ /*!
+ * \brief If the permeability should be evaluated directly at the scvf integration point
+ * (for convergence tests with analytical and continuous perm functions) or is evaluated
+ * at the scvs (for permeability fields with discontinuities) -> default
+ */
+ static constexpr bool evaluatePermeabilityAtScvfIP()
+ { return false; }
+
/*!
* \brief Function for defining the porosity.
* That is possibly solution dependent.
diff --git a/dumux/porousmediumflow/properties.hh b/dumux/porousmediumflow/properties.hh
index 10a5765aee..7cfee9b501 100644
--- a/dumux/porousmediumflow/properties.hh
+++ b/dumux/porousmediumflow/properties.hh
@@ -66,9 +66,6 @@ SET_BOOL_PROP(PorousMediumFlow, SolutionDependentAdvection, true);
SET_BOOL_PROP(PorousMediumFlow, SolutionDependentMolecularDiffusion, true);
SET_BOOL_PROP(PorousMediumFlow, SolutionDependentHeatConduction, true);
-//! By default, we evaluate the permeability in the volume
-SET_BOOL_PROP(PorousMediumFlow, EvaluatePermeabilityAtScvfIP, false);
-
//! The default implementation of the energy balance equation for flow problems in porous media.
SET_TYPE_PROP(PorousMediumFlow, EnergyLocalResidual, EnergyLocalResidual<TypeTag> );
--
GitLab