Commit f78e2693 authored by Timo Koch's avatar Timo Koch
Browse files

[tpfa] Free darcy's law from TypeTag

parent 9a467dc5
......@@ -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);
......
......@@ -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);
......
......@@ -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.
......
......@@ -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> );
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment