diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index b81375dca7b05a028d76825f2a1937d6f8f6ceb9..10b3b9aceaa535a74f149cc426781ee0ea133896 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -80,6 +80,9 @@ class FVProblem using BoundaryTypes = Dumux::BoundaryTypes; public: + //! Export spatial parameter type + using SpatialParams = GetPropType; + //! export traits of this problem struct Traits { @@ -91,16 +94,33 @@ public: /*! * \brief Constructor * \param gridGeometry The finite volume grid geometry + * \param spatialParams Spatially varying parameters * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") */ - FVProblem(std::shared_ptr gridGeometry, const std::string& paramGroup = "") + FVProblem(std::shared_ptr gridGeometry, + std::shared_ptr spatialParams, + const std::string& paramGroup = "") : gridGeometry_(gridGeometry) + , spatialParams_(spatialParams) , paramGroup_(paramGroup) { // set a default name for the problem problemName_ = getParamFromGroup(paramGroup, "Problem.Name"); } + /*! + * \brief Constructor + * \param gridGeometry The finite volume grid geometry + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + * \note This constructor assumes the spatial parameters to be constructible from a grid geometry + */ + FVProblem(std::shared_ptr gridGeometry, + const std::string& paramGroup = "") + : FVProblem(gridGeometry, + std::make_shared(gridGeometry), + paramGroup) + {} + /*! * \brief The problem name. * @@ -545,6 +565,7 @@ public: * are assumed to extend 1 m to the back. */ template + [[deprecated("Use problem.spatialParams().extrusionFactor() instead!")]] Scalar extrusionFactor(const Element& element, const SubControlVolume& scv, const ElementSolution& elemSol) const @@ -562,6 +583,7 @@ public: * thought as pipes with a cross section of 1 m^2 and 2D problems * are assumed to extend 1 m to the back. */ + [[deprecated("Use problem.spatialParams().extrusionFactorAtPos() instead!")]] Scalar extrusionFactorAtPos(const GlobalPosition &globalPos) const { // As a default, i.e. if the user's problem does not overload @@ -579,6 +601,14 @@ public: const std::string& paramGroup() const { return paramGroup_; } + //! Return the spatial parameters + SpatialParams& spatialParams() + { return *spatialParams_; } + + //! Return the spatial parameters + const SpatialParams& spatialParams() const + { return *spatialParams_; } + protected: //! Returns the implementation of the problem (i.e. static polymorphism) Implementation &asImp_() @@ -592,6 +622,9 @@ private: //! The finite volume grid geometry std::shared_ptr gridGeometry_; + //! Spatially varying parameters + std::shared_ptr spatialParams_; + //! The parameter group in which to retrieve runtime parameters std::string paramGroup_; diff --git a/dumux/common/fvspatialparams.hh b/dumux/common/fvspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..156c8e3477289aca62db5269a645c39db1757007 --- /dev/null +++ b/dumux/common/fvspatialparams.hh @@ -0,0 +1,141 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Common + * \ingroup SpatialParameters + * \brief Basic spatial parameters + */ +#ifndef DUMUX_COMMON_FV_SPATIAL_PARAMS_BASE_HH +#define DUMUX_COMMON_FV_SPATIAL_PARAMS_BASE_HH + +#include + +#include +#include + +namespace Dumux { + +/*! + * \ingroup SpatialParameters + * \brief The base class for spatial parameters. + */ +template +class FVSpatialParamsBase +{ + using GridView = typename GridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using SubControlVolume = typename GridGeometry::SubControlVolume; + + static constexpr int dimWorld = GridView::dimensionworld; + + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + using GravityVector = Dune::FieldVector; + +public: + FVSpatialParamsBase(std::shared_ptr gridGeometry) + : gridGeometry_(gridGeometry) + , gravity_(0.0) + { + if (getParam("Problem.EnableGravity")) + gravity_[dimWorld-1] = -9.81; + } + + /*! + * \brief Return how much the domain is extruded at a given sub-control volume. + * + * This means the factor by which a lower-dimensional (1D or 2D) + * entity needs to be expanded to get a full dimensional cell. The + * default is 1.0 which means that 1D problems are actually + * thought as pipes with a cross section of 1 m^2 and 2D problems + * are assumed to extend 1 m to the back. + */ + template + Scalar extrusionFactor(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + // forward to generic interface + return asImp_().extrusionFactorAtPos(scv.center()); + } + + /*! + * \brief Return how much the domain is extruded at a given position. + */ + Scalar extrusionFactorAtPos(const GlobalPosition& globalPos) const + { + // As a default, i.e. if the user's problem does not overload + // any extrusion factor method, return 1.0 + return 1.0; + } + + /*! + * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$. + * + * The default behaviour is a constant gravity vector; + * if the Problem.EnableGravity parameter is true, + * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, + * else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$. + * + * \param pos the spatial position at which to evaulate the gravity vector + */ + const GravityVector& gravity(const GlobalPosition& pos) const + { return gravity_; } + + //! The finite volume grid geometry + const GridGeometry& gridGeometry() const + { return *gridGeometry_; } + +protected: + //! Returns the implementation of the spatial parameters (static polymorphism) + Implementation &asImp_() + { return *static_cast(this); } + + //! \copydoc asImp_() + const Implementation &asImp_() const + { return *static_cast(this); } + +private: + std::shared_ptr gridGeometry_; + GravityVector gravity_; +}; + +/*! + * \ingroup SpatialParameters + * \brief Default spatial parameters class to be reused in models + * that solely need to define gravity and extrusion. + */ +template +class DefaultFVSpatialParams +: public FVSpatialParamsBase> +{ + using ThisType = DefaultFVSpatialParams; + using ParentType = FVSpatialParamsBase; + +public: + using ParentType::ParentType; +}; + +} // namespace Dumux + +#endif diff --git a/dumux/freeflow/navierstokes/mass/1p/model.hh b/dumux/freeflow/navierstokes/mass/1p/model.hh index a518e299c7b66e35c832d3ce0aa61cf151c27c3b..3cf152bb57f369775ba4c48d1bea86ed13c7ea3e 100644 --- a/dumux/freeflow/navierstokes/mass/1p/model.hh +++ b/dumux/freeflow/navierstokes/mass/1p/model.hh @@ -146,6 +146,14 @@ template struct ModelTraits { using type = NavierStokesMassOnePModelTraits; }; +//! Per default, we use the base spatial parameters +template +struct SpatialParams +{ + using type = DefaultFVSpatialParams, + GetPropType>; +}; + /*! * \brief The fluid state which is used by the volume variables to * store the thermodynamic state. This should be chosen diff --git a/dumux/freeflow/navierstokes/momentum/localresidual.hh b/dumux/freeflow/navierstokes/momentum/localresidual.hh index 5db4a54c6229d27f9ae6aa240f6967d9aadf7e22..a023f0f1c4a5ebc0ad6d614bb417e3518545b101 100644 --- a/dumux/freeflow/navierstokes/momentum/localresidual.hh +++ b/dumux/freeflow/navierstokes/momentum/localresidual.hh @@ -123,7 +123,7 @@ public: const SubControlVolume& scv) const { NumEqVector source = ParentType::computeSource(problem, element, fvGeometry, elemVolVars, scv); - source += problem.gravity()[scv.dofAxis()] * problem.density(element, scv); + source += problem.spatialParams().gravity(scv.dofPosition())[scv.dofAxis()] * problem.density(element, scv); // Axisymmetric problems in 2D feature an extra source terms arising from the transformation to cylindrical coordinates. // See Ferziger/Peric: Computational methods for fluid dynamics chapter 8. diff --git a/dumux/freeflow/navierstokes/momentum/model.hh b/dumux/freeflow/navierstokes/momentum/model.hh index a3110591375fa454a730ff288c43c0147a256586..e37fad35afd3340a8fe8a79d3da0abaa2bd874cf 100644 --- a/dumux/freeflow/navierstokes/momentum/model.hh +++ b/dumux/freeflow/navierstokes/momentum/model.hh @@ -221,6 +221,14 @@ public: using type = EmptyCouplingManager; }; +//! Per default, we use the base spatial parameters +template +struct SpatialParams +{ + using type = DefaultFVSpatialParams, + GetPropType>; +}; + //! The specific I/O fields // template // struct IOFields { using type = NavierStokesIOFields; }; diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 092603f0324a6b1c326f92ede9dd6bb86291e7b8..d20f9a88b2aafabf503437a76dcd67350be7ee24 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -82,13 +82,13 @@ class NavierStokesProblemImpl using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; using VelocityVector = Dune::FieldVector; - using GravityVector = Dune::FieldVector; using CouplingManager = GetPropType; using ModelTraits = GetPropType; static constexpr bool isCoupled_ = !std::is_empty_v; public: + using typename ParentType::SpatialParams; //! Export traits of this problem struct Traits @@ -121,12 +121,25 @@ public: std::shared_ptr couplingManager, const std::string& paramGroup = "") : ParentType(gridGeometry, paramGroup) - , gravity_(0.0) , couplingManager_(couplingManager) { - if (getParamFromGroup(paramGroup, "Problem.EnableGravity")) - gravity_[dim-1] = -9.81; + enableInertiaTerms_ = getParamFromGroup(paramGroup, "Problem.EnableInertiaTerms"); + } + /*! + * \brief The constructor + * \param gridGeometry The finite volume grid geometry + * \param spatialParams Spatially varying parameters + * \param couplingManager The coupling manager (couples mass and momentum equations) + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ + NavierStokesProblemImpl(std::shared_ptr gridGeometry, + std::shared_ptr spatialParams, + std::shared_ptr couplingManager, + const std::string& paramGroup = "") + : ParentType(gridGeometry, spatialParams, paramGroup) + , couplingManager_(couplingManager) + { enableInertiaTerms_ = getParamFromGroup(paramGroup, "Problem.EnableInertiaTerms"); } @@ -140,6 +153,18 @@ public: : NavierStokesProblemImpl(gridGeometry, {}, paramGroup) {} + /*! + * \brief The constructor for usage without a coupling manager + * \param gridGeometry The finite volume grid geometry + * \param spatialParams Spatially varying parameters + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ + NavierStokesProblemImpl(std::shared_ptr gridGeometry, + std::shared_ptr spatialParams, + const std::string& paramGroup = "") + : NavierStokesProblemImpl(gridGeometry, spatialParams, {}, paramGroup) + {} + /*! * \brief Evaluate the source term for all phases within a given * sub-control-volume. @@ -223,8 +248,9 @@ public: * If the Problem.EnableGravity parameter is true, this means * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$ */ - const GravityVector& gravity() const - { return gravity_; } + [[deprecated("Use problem.spatialParams().gravity(globalPos) instead!")]] + decltype(auto) gravity() const + { return this->spatialParams().gravity(GlobalPosition{}); } /*! * \brief Returns whether intertia terms should be considered. @@ -527,7 +553,6 @@ private: const Implementation& asImp_() const { return *static_cast(this); } - GravityVector gravity_; bool enableInertiaTerms_; std::shared_ptr couplingManager_; }; @@ -553,6 +578,8 @@ class NavierStokesProblemImpl static constexpr bool isCoupled_ = !std::is_empty_v; public: + using typename ParentType::SpatialParams; + //! Export traits of this problem struct Traits { @@ -586,12 +613,39 @@ public: , couplingManager_(couplingManager) {} + /*! + * \brief The constructor + * \param gridGeometry The finite volume grid geometry + * \param spatialParams Spatially varying parameters + * \param couplingManager The coupling manager (couples mass and momentum equations) + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ + NavierStokesProblemImpl(std::shared_ptr gridGeometry, + std::shared_ptr spatialParams, + std::shared_ptr couplingManager, + const std::string& paramGroup = "") + : ParentType(gridGeometry, spatialParams, paramGroup) + , couplingManager_(couplingManager) + {} + + /*! + * \brief The constructor for usage without a coupling manager + * \param gridGeometry The finite volume grid geometry + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ + NavierStokesProblemImpl(std::shared_ptr gridGeometry, + const std::string& paramGroup = "") + : NavierStokesProblemImpl(gridGeometry, {}, paramGroup) + {} + /*! * \brief The constructor for usage without a coupling manager * \param gridGeometry The finite volume grid geometry + * \param spatialParams Spatially varying parameters * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") */ NavierStokesProblemImpl(std::shared_ptr gridGeometry, + std::shared_ptr spatialParams, const std::string& paramGroup = "") : NavierStokesProblemImpl(gridGeometry, {}, paramGroup) {} @@ -664,7 +718,7 @@ private: * \ingroup NavierStokesModel * \brief Navier-Stokes problem base class. * - * This implements gravity (if desired) and a function returning the temperature. + * This implements a function returning the temperature. * Includes a specialized method used only by the staggered grid discretization. */ template @@ -699,9 +753,10 @@ class NavierStokesProblemImpl using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; using VelocityVector = Dune::FieldVector; - using GravityVector = Dune::FieldVector; public: + using typename ParentType::SpatialParams; + /*! * \brief The constructor * \param gridGeometry The finite volume grid geometry @@ -709,11 +764,21 @@ public: */ NavierStokesProblemImpl(std::shared_ptr gridGeometry, const std::string& paramGroup = "") : ParentType(gridGeometry, paramGroup) - , gravity_(0.0) { - if (getParamFromGroup(paramGroup, "Problem.EnableGravity")) - gravity_[dim-1] = -9.81; + enableInertiaTerms_ = getParamFromGroup(paramGroup, "Problem.EnableInertiaTerms"); + } + /*! + * \brief The constructor + * \param gridGeometry The finite volume grid geometry + * \param spatialParams Spatially varying parameters + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ + NavierStokesProblemImpl(std::shared_ptr gridGeometry, + std::shared_ptr spatialParams, + const std::string& paramGroup = "") + : ParentType(gridGeometry, spatialParams, paramGroup) + { enableInertiaTerms_ = getParamFromGroup(paramGroup, "Problem.EnableInertiaTerms"); } @@ -742,8 +807,9 @@ public: * If the Problem.EnableGravity parameter is true, this means * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$ */ - const GravityVector& gravity() const - { return gravity_; } + [[deprecated("Use problem.spatialParams().gravity(globalPos) instead")]] + decltype(auto) gravity() const + { return this->spatialParams().gravity(GlobalPosition{}); } /*! * \brief Returns whether interia terms should be considered. @@ -886,7 +952,6 @@ private: const Implementation &asImp_() const { return *static_cast(this); } - GravityVector gravity_; bool enableInertiaTerms_; }; diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh index 2e68a79bf58792ad69a42385b2119a9a1910a883..aa2d555345efb5d611c5f33ce182361587914b99 100644 --- a/dumux/freeflow/navierstokes/staggered/localresidual.hh +++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh @@ -177,7 +177,7 @@ public: { FacePrimaryVariables source(0.0); const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; - source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density(); + source += problem.spatialParams().gravity(scvf.ipGlobal())[scvf.directionIndex()] * insideVolVars.density(); source += problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]; // Axisymmetric problems in 2D feature an extra source terms arising from the transformation to cylindrical coordinates. diff --git a/dumux/freeflow/properties.hh b/dumux/freeflow/properties.hh index d918192b27a378bedae6980a6e430d4d641c58a7..c6df54fe47a80cd5f36625cf19a6a99c20d0b4f0 100644 --- a/dumux/freeflow/properties.hh +++ b/dumux/freeflow/properties.hh @@ -27,6 +27,8 @@ #include #include +#include + #include #include #include "turbulencemodel.hh" @@ -54,6 +56,14 @@ struct FluxVariablesCacheFiller { using type = FluxVari template struct HeatConductionType { using type = FouriersLaw; }; +//! Per default, we use the base spatial parameters +template +struct SpatialParams +{ + using type = DefaultFVSpatialParams, + GetPropType>; +}; + } // namespace Properties } // namespace Dumux diff --git a/dumux/geomechanics/elastic/fvspatialparams.hh b/dumux/geomechanics/elastic/fvspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..a1c98ccb59fe9cb2d184d1392abcdc9186c72091 --- /dev/null +++ b/dumux/geomechanics/elastic/fvspatialparams.hh @@ -0,0 +1,142 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup SpatialParameters + * \brief The base class for spatial parameters of linear elastic geomechanical problems + */ +#ifndef DUMUX_GEOMECHANICS_ELASTIC_FV_SPATIAL_PARAMS_HH +#define DUMUX_GEOMECHANICS_ELASTIC_FV_SPATIAL_PARAMS_HH + +#include + +#include + +#include +#include +#include + +namespace Dumux { + +#ifndef DOXYGEN +namespace Detail { +// helper struct detecting if the user-defined spatial params class has a lameParamsAtPos function +// for g++ > 5.3, this can be replaced by a lambda +template +struct hasLameParamsAtPos +{ + template + auto operator()(const SpatialParams& a) + -> decltype(a.lameParamsAtPos(std::declval())) + {} +}; + +} // end namespace Detail +#endif + +/*! + * \ingroup SpatialParameters + * \brief The base class for spatial parameters of linear elastic geomechanical problems + */ +template +class FVSpatialParamsElastic +: public FVSpatialParamsBase +{ + using ParentType = FVSpatialParamsBase + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + using ParentType::ParentType; + + /*! + * \brief Function for defining the solid volume fraction. + * That is possibly solution dependent. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \param compIdx The solid component index + * \return the volume fraction of the inert solid component with index compIdx + */ + template + Scalar inertVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { + static_assert(SolidSystem::isInert(), "Elastic model can only be used with inert solid systems"); + + // when there is only one component, the volume fraction must be one + if (SolidSystem::numInertComponents == 1) + return 1.0; + + // otherwise we require the user to define the solid composition + return this->asImp_().template inertVolumeFractionAtPos(scv.center(), compIdx); + } + + /*! + * \brief Function for defining the solid volume fraction. + * That is possibly solution dependent. + * + * \param globalPos The global position + * \param compIdx The solid component index + * \return the volume fraction of the inert solid component with index compIdx + */ + template + Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const + { DUNE_THROW(Dune::InvalidStateException, "The spatial parameters do not provide inertVolumeFractionAtPos() method."); } + + /*! + * \brief Define the Lame parameters + * \note These are possibly solution dependent and are evaluated + * for an integration point inside the element. Therefore, + * a flux variables cache object is passed to this function + * containing data on shape functions at the integration point. + * + * \param element The current element + * \param fvGeometry The local finite volume geometry + * \param elemVolVars Primary/Secondary variables inside the element + * \param fluxVarsCache Contains data on shape functions at the integration point + * \return lame parameters + */ + template + decltype(auto) lameParams(const Element& element, + const FVElementGeometry& fvGeometry, + const ElemVolVars& elemVolVars, + const FluxVarsCache& fluxVarsCache) const + { + static_assert(decltype(isValid(Detail::hasLameParamsAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const\n\n" + " or overload this function\n\n" + " template\n" + " const LameParams& lameParams(const Element& element,\n" + " const FVElementGeometry& fvGeometry,\n" + " const ElemVolVars& elemVolVars,\n" + " const FluxVarsCache& fluxVarsCache) const\n\n"); + + return this->asImp_().lameParamsAtPos(fluxVarsCache.ipGlobal()); + } +}; + +} // end namespace Dumux +#endif diff --git a/dumux/geomechanics/poroelastic/fvspatialparams.hh b/dumux/geomechanics/poroelastic/fvspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..92ca7c25b845d37290c6c97598e0b6c76d8e2015 --- /dev/null +++ b/dumux/geomechanics/poroelastic/fvspatialparams.hh @@ -0,0 +1,298 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup SpatialParameters + * \brief The base class for spatial parameters of poro-elastic geomechanical problems + */ +#ifndef DUMUX_GEOMECHANICS_POROELASTIC_FV_SPATIAL_PARAMS_HH +#define DUMUX_GEOMECHANICS_POROELASTIC_FV_SPATIAL_PARAMS_HH + +#include + +#include +#include +#include +#include + +namespace Dumux { + +#ifndef DOXYGEN +namespace Detail { +// helper struct detecting if the user-defined spatial params class has a lameParamsAtPos function +// for g++ > 5.3, this can be replaced by a lambda +template +struct hasLameParamsAtPos +{ + template + auto operator()(const SpatialParams& a) + -> decltype(a.lameParamsAtPos(std::declval())) + {} +}; + +// helper struct detecting if the user-defined spatial params class has a reactiveVolumeFractionAtPos function +// for g++ > 5.3, this can be replaced by a lambda +template +struct hasReactiveVolumeFractionAtPos +{ + template + auto operator()(const SpatialParams& a) + -> decltype(a.template reactiveVolumeFractionAtPos(std::declval(), 0)) + {} +}; + +// helper struct detecting if the user-defined spatial params class has a biotCoefficientAtPos function +// for g++ > 5.3, this can be replaced by a lambda +template +struct hasBiotCoeffAtPos +{ + template + auto operator()(const SpatialParams& a) + -> decltype(a.biotCoefficientAtPos(std::declval())) + {} +}; + +} // end namespace Detail +#endif + +/*! + * \ingroup SpatialParameters + * \brief The base class for spatial parameters of poro-elastic geomechanical problems + */ +template +class FVSpatialParamsPoroElastic +: public FVSpatialParamsBase +{ + using ParentType = FVSpatialParamsBase; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + using ParentType::ParentType; + + /*! + * \brief Function for defining the porosity. + * That is possibly solution dependent. + * \note this can only be used for solids with one inert component + * (see inertVolumeFraction for the more general interface) + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return the porosity + */ + template + Scalar porosity(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + static_assert(decltype(isValid(Detail::hasPorosityAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " Scalar porosityAtPos(const GlobalPosition& globalPos) const\n\n" + " or overload this function\n\n" + " template\n" + " Scalar porosity(const Element& element,\n" + " const SubControlVolume& scv,\n" + " const ElementSolution& elemSol) const\n\n"); + + return this->asImp_().porosityAtPos(scv.center()); + } + + /*! + * \brief Function for defining the solid volume fraction. + * That is possibly solution dependent. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \param compIdx The solid component index + * \return the volume fraction of the inert solid component with index compIdx + * + * \note this overload is enabled if there is only one inert solid component and the + * user didn't choose to implement an inertVolumeFractionAtPos overload. + * It then forwards to the simpler porosity interface. + * With more than one solid components or active solid components (i.e. dissolution) + * please overload the more general inertVolumeFraction/inertVolumeFractionAtPos interface. + */ + template()) + (std::declval()))::value, int> = 0> + Scalar inertVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { return 1.0 - this->asImp_().porosity(element, scv, elemSol); } + + // specialization if there are no inert components at all + template = 0> + Scalar inertVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { return 0.0; } + + // the more general interface forwarding to inertVolumeFractionAtPos + template 1) || + ( + (SolidSystem::numInertComponents > 0) && + ( + !SolidSystem::isInert() + || decltype(isValid(Detail::hasInertVolumeFractionAtPos()) + (std::declval()))::value + ) + ), int> = 0> + Scalar inertVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { + static_assert(decltype(isValid(Detail::hasInertVolumeFractionAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " template\n" + " Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const\n\n" + " or overload this function\n\n" + " template\n" + " Scalar inertVolumeFraction(const Element& element,\n" + " const SubControlVolume& scv,\n" + " const ElementSolution& elemSol,\n" + " int compIdx) const\n\n"); + + return this->asImp_().template inertVolumeFractionAtPos(scv.center(), compIdx); + } + + /*! + * \brief Function for defining the solid volume fraction of a solid + * component that takes part in some sort of reaction. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \param compIdx The solid component index + * \return the volume fraction of the inert solid component with index compIdx + * + * \note This overload is enabled if there are only inert solid components + * and the user did not choose to implement a reactiveVolumeFractionAtPos + * function. The reactive volume fraction is zero in this case. + */ + template()) + (std::declval()))::value, int > = 0 > + Scalar reactiveVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { return 0.0; } + + //! overload for the case of reactive solids or user-provided overload + template()) + (std::declval()))::value, int > = 0 > + Scalar reactiveVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { + static_assert(decltype(isValid(Detail::hasReactiveVolumeFractionAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " template\n" + " Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const\n\n" + " or overload this function\n\n" + " template\n" + " Scalar inertVolumeFraction(const Element& element,\n" + " const SubControlVolume& scv,\n" + " const ElementSolution& elemSol,\n" + " int compIdx) const\n\n"); + + return this->asImp_().template reactiveVolumeFractionAtPos(scv.center(), compIdx); + } + + /*! + * \brief Define the Lame parameters + * \note These are possibly solution dependent and are evaluated + * for an integration point inside the element. Therefore, + * a flux variables cache object is passed to this function + * containing data on shape functions at the integration point. + * + * \param element The current element + * \param fvGeometry The local finite volume geometry + * \param elemVolVars Primary/Secondary variables inside the element + * \param fluxVarsCache Contains data on shape functions at the integration point + * \return lame parameters + */ + template + decltype(auto) lameParams(const Element& element, + const FVElementGeometry& fvGeometry, + const ElemVolVars& elemVolVars, + const FluxVarsCache& fluxVarsCache) const + { + static_assert(decltype(isValid(Detail::hasLameParamsAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const\n\n" + " or overload this function\n\n" + " template\n" + " const LameParams& lameParams(const Element& element,\n" + " const FVElementGeometry& fvGeometry,\n" + " const ElemVolVars& elemVolVars,\n" + " const FluxVarsCache& fluxVarsCache) const\n\n"); + + return this->asImp_().lameParamsAtPos(fluxVarsCache.ipGlobal()); + } + + /*! + * \brief Returns the Biot coefficient in an element + * \note This is possibly solution dependent and is evaluated + * for an integration point inside the element. Therefore, + * a flux variables cache object is passed to this function + * containing data on shape functions at the integration point. + * + * \param element The current element + * \param fvGeometry The local finite volume geometry + * \param elemVolVars Primary/Secondary variables inside the element + * \param fluxVarsCache Contains data on shape functions at the integration point + * \return Biot coefficient + */ + template + Scalar biotCoefficient(const Element& element, + const FVElementGeometry& fvGeometry, + const ElemVolVars& elemVolVars, + const FluxVarsCache& fluxVarsCache) const + { + static_assert(decltype(isValid(Detail::hasBiotCoeffAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " const LameParams& biotCoefficientAtPos(const GlobalPosition& globalPos) const\n\n" + " or overload this function\n\n" + " template\n" + " const LameParams& biotCoefficient(const Element& element,\n" + " const FVElementGeometry& fvGeometry,\n" + " const ElemVolVars& elemVolVars,\n" + " const FluxVarsCache& fluxVarsCache) const\n\n"); + + return this->asImp_().biotCoefficientAtPos(fluxVarsCache.ipGlobal()); + } +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/material/spatialparams/fv.hh b/dumux/material/spatialparams/fv.hh index 520098da684204945ff0fd08e8b481976905de5f..ce593b9b7cc026ee3fc9136a53b9d1ee6af983e9 100644 --- a/dumux/material/spatialparams/fv.hh +++ b/dumux/material/spatialparams/fv.hh @@ -25,102 +25,15 @@ #ifndef DUMUX_FV_SPATIAL_PARAMS_HH #define DUMUX_FV_SPATIAL_PARAMS_HH -#include -#include -#include "fv1p.hh" +#warning "This header will be removed after 3.5 in favor of dumux/porousmediumflow/fvspatialparams.hh" +#include namespace Dumux { -#ifndef DOXYGEN -namespace Detail { -// helper struct detecting if the user-defined spatial params class -// has a fluidMatrixInteractionAtPos function -template -struct hasFluidMatrixInteractionAtPos -{ - template - auto operator()(const SpatialParams& a) - -> decltype(a.fluidMatrixInteractionAtPos(std::declval())) - {} -}; -} // end namespace Detail -#endif - -/*! - * \ingroup SpatialParameters - * \brief The base class for spatial parameters of multi-phase problems - * using a fully implicit discretization method. - */ template -class FVSpatialParams : public FVSpatialParamsOneP -{ - using ParentType = FVSpatialParamsOneP; - using GridView = typename GridGeometry::GridView; - using FVElementGeometry = typename GridGeometry::LocalView; - using SubControlVolume = typename GridGeometry::SubControlVolume; - using Element = typename GridView::template Codim<0>::Entity; - - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - -public: - FVSpatialParams(std::shared_ptr gridGeometry) - : ParentType(gridGeometry) - {} - - /*! - * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.). - * - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - */ - template - decltype(auto) fluidMatrixInteraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const - { - static_assert(decltype(isValid(Detail::hasFluidMatrixInteractionAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const\n\n" - " or overload this function\n\n" - " template\n" - " auto fluidMatrixInteraction(const Element& element,\n" - " const SubControlVolume& scv,\n" - " const ElementSolution& elemSol) const\n\n"); - - return this->asImp_().fluidMatrixInteractionAtPos(scv.center()); - } - - /*! - * \brief Function for defining which phase is to be considered as the wetting phase. - * - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - * \return the wetting phase index - */ - template - int wettingPhase(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const - { - return this->asImp_().template wettingPhaseAtPos(scv.center()); - } - - /*! - * \brief Function for defining which phase is to be considered as the wetting phase. - * - * \return the wetting phase index - * \param globalPos The global position - */ - template - int wettingPhaseAtPos(const GlobalPosition& globalPos) const - { - DUNE_THROW(Dune::InvalidStateException, - "The spatial parameters do not provide " - "a wettingPhaseAtPos() method."); - } -}; +using FVSpatialParams +[[deprecated("Use FVPorousMediumSpatialParams in dumux/porousmediumflow/fvspatialparams.hh instead!")]] += FVPorousMediumSpatialParams; } // namespace Dumux diff --git a/dumux/material/spatialparams/fv1p.hh b/dumux/material/spatialparams/fv1p.hh index d638f82598642d7d15af2e7b39f8540dc5f79d79..a6a1a34e02c592171831a0a7bbb4da43c505af40 100644 --- a/dumux/material/spatialparams/fv1p.hh +++ b/dumux/material/spatialparams/fv1p.hh @@ -25,323 +25,15 @@ #ifndef DUMUX_FV_SPATIAL_PARAMS_ONE_P_HH #define DUMUX_FV_SPATIAL_PARAMS_ONE_P_HH -#include -#include - -#include -#include -#include +#warning "This header will be removed after 3.5 in favor of dumux/porousmediumflow/fvspatialparams.hh" +#include namespace Dumux { -#ifndef DOXYGEN -namespace Detail { -// helper struct detecting if the user-defined spatial params class has a permeabilityAtPos function -// for g++ > 5.3, this can be replaced by a lambda -template -struct hasPermeabilityAtPos -{ - template - auto operator()(const SpatialParams& a) - -> decltype(a.permeabilityAtPos(std::declval())) - {} -}; - -template -struct hasInertVolumeFractionAtPos -{ - template - auto operator()(const SpatialParams& a) - -> decltype(a.template inertVolumeFractionAtPos(std::declval(), 0)) - {} -}; - -template -struct hasPorosityAtPos -{ - template - auto operator()(const SpatialParams& a) - -> decltype(a.porosityAtPos(std::declval())) - {} -}; -} // end namespace Detail -#endif - -/*! - * \ingroup SpatialParameters - * \brief The base class for spatial parameters of one-phase problems - * using a fully implicit discretization method. - */ template -class FVSpatialParamsOneP -{ - using GridView = typename GridGeometry::GridView; - using FVElementGeometry = typename GridGeometry::LocalView; - using SubControlVolume = typename GridGeometry::SubControlVolume; - using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; - using Element = typename GridView::template Codim<0>::Entity; - - enum { dim = GridView::dimension }; - enum { dimWorld = GridView::dimensionworld }; - using DimWorldMatrix = Dune::FieldMatrix; - - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - -public: - FVSpatialParamsOneP(std::shared_ptr gridGeometry) - : gridGeometry_(gridGeometry) - , gravity_(0.0) - { - const bool enableGravity = getParam("Problem.EnableGravity"); - if (enableGravity) - gravity_[dimWorld-1] = -9.81; - - /* \brief default forchheimer coefficient - * Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90 \cite ward1964 . - * Actually the Forchheimer coefficient is also a function of the dimensions of the - * porous medium. Taking it as a constant is only a first approximation - * (Nield, Bejan, Convection in porous media, 2006, p. 10 \cite nield2006 ) - */ - forchCoeffDefault_ = getParam("SpatialParams.ForchCoeff", 0.55); - } - - /*! - * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$. - * - * The default behaviour is a constant gravity vector; - * if the Problem.EnableGravity parameter is true, - * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, - * else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$. - * - * \param pos the spatial position at which to evaulate the gravity vector - */ - const GlobalPosition& gravity(const GlobalPosition &pos) const - { return gravity_; } - - /*! - * \brief Harmonic average of a discontinuous scalar field at discontinuity interface - * (for compatibility reasons with the function below) - * \return the averaged scalar - * \param T1 first scalar parameter - * \param T2 second scalar parameter - * \param normal The unit normal vector of the interface - */ - Scalar harmonicMean(const Scalar T1, - const Scalar T2, - const GlobalPosition& normal) const - { return Dumux::harmonicMean(T1, T2); } - - /*! - * \brief Harmonic average of a discontinuous tensorial field at discontinuity interface - * \note We do a harmonic average of the part normal to the interface (alpha*I) and - * an arithmetic average of the tangential part (T - alpha*I). - * \return the averaged tensor - * \param T1 first tensor - * \param T2 second tensor - * \param normal The unit normal vector of the interface - */ - DimWorldMatrix harmonicMean(const DimWorldMatrix& T1, - const DimWorldMatrix& T2, - const GlobalPosition& normal) const - { - // determine nT*k*n - GlobalPosition tmp; - GlobalPosition tmp2; - T1.mv(normal, tmp); - T2.mv(normal, tmp2); - const Scalar alpha1 = tmp*normal; - const Scalar alpha2 = tmp2*normal; - - const Scalar alphaHarmonic = Dumux::harmonicMean(alpha1, alpha2); - const Scalar alphaAverage = 0.5*(alpha1 + alpha2); - - DimWorldMatrix T(0.0); - for (int i = 0; i < dimWorld; ++i) - { - for (int j = 0; j < dimWorld; ++j) - { - T[i][j] += 0.5*(T1[i][j] + T2[i][j]); - if (i == j) - T[i][j] += alphaHarmonic - alphaAverage; - } - } - - return T; - } - - /*! - * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$ - * \note It is possibly solution dependent. - * - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - * \return permeability - */ - template - decltype(auto) permeability(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const - { - static_assert(decltype(isValid(Detail::hasPermeabilityAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " const PermeabilityType& permeabilityAtPos(const GlobalPosition& globalPos) const\n\n" - " or overload this function\n\n" - " template\n" - " const PermeabilityType& permeability(const Element& element,\n" - " const SubControlVolume& scv,\n" - " const ElementSolution& elemSol) const\n\n"); - - return asImp_().permeabilityAtPos(scv.center()); - } - - /*! - * \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. - * \note this can only be used for solids with one inert component - * (see inertVolumeFraction for the more general interface) - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - * \return the porosity - */ - template - Scalar porosity(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const - { - static_assert(decltype(isValid(Detail::hasPorosityAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " Scalar porosityAtPos(const GlobalPosition& globalPos) const\n\n" - " or overload this function\n\n" - " template\n" - " Scalar porosity(const Element& element,\n" - " const SubControlVolume& scv,\n" - " const ElementSolution& elemSol) const\n\n"); - - return asImp_().porosityAtPos(scv.center()); - } - - /*! - * \brief Function for defining the solid volume fraction. - * That is possibly solution dependent. - * - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - * \param compIdx The solid component index - * \return the volume fraction of the inert solid component with index compIdx - * - * \note this overload is enable if there is only one inert solid component and the - * user didn't choose to implement a inertVolumeFractionAtPos overload. - * It then forwards to the simpler porosity interface. - * With more than one solid components or active solid components (i.e. dissolution) - * please overload the more general inertVolumeFraction/inertVolumeFractionAtPos interface. - */ - template())(std::declval()))::value, - int> = 0> - Scalar inertVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { - return 1.0 - asImp_().porosity(element, scv, elemSol); - } - - // specialization if there are no inert components at all - template = 0> - Scalar inertVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { - return 0.0; - } - - // the more general interface forwarding to inertVolumeFractionAtPos - template 1) || - ( - (SolidSystem::numInertComponents > 0) && - ( - !SolidSystem::isInert() - || decltype(isValid(Detail::hasInertVolumeFractionAtPos()) - (std::declval()))::value - ) - ), - int> = 0> - Scalar inertVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { - static_assert(decltype(isValid(Detail::hasInertVolumeFractionAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " template\n" - " Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const\n\n" - " or overload this function\n\n" - " template\n" - " Scalar inertVolumeFraction(const Element& element,\n" - " const SubControlVolume& scv,\n" - " const ElementSolution& elemSol,\n" - " int compIdx) const\n\n"); - - return asImp_().template inertVolumeFractionAtPos(scv.center(), compIdx); - } - - /*! - * \brief Function for defining the Beavers-Joseph coefficient for multidomain - * problems\f$\mathrm{[-]}\f$. - * - * \return Beavers-Joseph coefficient \f$\mathrm{[-]}\f$ - * \param globalPos The global position - */ - Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const - { - DUNE_THROW(Dune::InvalidStateException, - "The spatial parameters do not provide a beaversJosephCoeffAtPos() method."); - } - - /*! - * \brief Apply the Forchheimer coefficient for inertial forces - * calculation. - * \param scvf The sub-control volume face where the - * intrinsic velocity ought to be calculated. - */ - Scalar forchCoeff(const SubControlVolumeFace &scvf) const - { - return forchCoeffDefault_; - } - - //! The finite volume grid geometry - const GridGeometry& gridGeometry() const - { return *gridGeometry_; } - - -protected: - Implementation &asImp_() - { return *static_cast(this); } - - const Implementation &asImp_() const - { return *static_cast(this); } - -private: - std::shared_ptr gridGeometry_; - GlobalPosition gravity_; //!< The gravity vector - Scalar forchCoeffDefault_; -}; +using FVSpatialParamsOneP +[[deprecated("Use FVPorousMediumSpatialParamsOneP in dumux/porousmediumflow/fvspatialparams.hh instead!")]] + = FVPorousMediumSpatialParamsOneP; } // namespace Dumux diff --git a/dumux/material/spatialparams/fv1pconstant.hh b/dumux/material/spatialparams/fv1pconstant.hh index aed91024053f6a9028956e1f9ea113d6cc8eeb0e..0cd4ae366ccccf1303d9f70bfdfb5e1ef869190e 100644 --- a/dumux/material/spatialparams/fv1pconstant.hh +++ b/dumux/material/spatialparams/fv1pconstant.hh @@ -24,48 +24,15 @@ #ifndef DUMUX_FV_CONSTANT_SPATIAL_PARAMS_ONE_P_HH #define DUMUX_FV_CONSTANT_SPATIAL_PARAMS_ONE_P_HH -#include -#include +#warning "This header will be removed after 3.5 in favor of dumux/porousmediumflow/fv1pconstantspatialparams.hh" +#include namespace Dumux { -/*! - * \ingroup SpatialParameters - * \brief A spatial params implementation for 1p problem with constant properties - */ template -class FVSpatialParamsOnePConstant -: public FVSpatialParamsOneP> -{ - using ThisType = FVSpatialParamsOnePConstant; - using ParentType = FVSpatialParamsOneP; - using GlobalPosition = typename GridGeometry::GridView::template Codim<0>::Geometry::GlobalCoordinate; - -public: - using PermeabilityType = Scalar; - - FVSpatialParamsOnePConstant(std::shared_ptr gridGeometry) - : ParentType(gridGeometry) - , porosity_(getParam("SpatialParams.Porosity")) - , permeability_(getParam("SpatialParams.Permeability")) - {} - - /*! - * \brief The (intrinsic) permeability \f$[m^2]\f$ - */ - PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const - { return permeability_; } - - /*! - * \brief The porosity \f$[-]\f$ - */ - Scalar porosityAtPos(const GlobalPosition& globalPos) const - { return porosity_; } - -private: - const Scalar porosity_; - const Scalar permeability_; -}; +using FVSpatialParamsOnePConstant +[[deprecated("Use FVSpatialParamsOnePConstant in dumux/porousmediumflow/fv1pconstantspatialparams.hh instead!")]] += FVPorousMediumOnePConstantSpatialParams; } // end namespace Dumux diff --git a/dumux/material/spatialparams/fvelastic.hh b/dumux/material/spatialparams/fvelastic.hh index ed89b81e165aabdeb7eeabe00d79d05a7043d78f..a0aead1de91b569bb74abfffa5b1c8e1d57cdf7d 100644 --- a/dumux/material/spatialparams/fvelastic.hh +++ b/dumux/material/spatialparams/fvelastic.hh @@ -21,157 +21,10 @@ * \ingroup SpatialParameters * \brief The base class for spatial parameters of linear elastic geomechanical problems */ -#ifndef DUMUX_GEOMECHANICS_ELASTIC_FV_SPATIAL_PARAMS_HH -#define DUMUX_GEOMECHANICS_ELASTIC_FV_SPATIAL_PARAMS_HH +#ifndef DUMUX_MATERIAL_ELASTIC_SPATIAL_PARAMS_HH +#define DUMUX_MATERIAL_ELASTIC_SPATIAL_PARAMS_HH -#include +#warning "This header will be removed after 3.5 in favor of dumux/geomechanics/elastic/fvspatialparams.hh" +#include -#include - -#include -#include - -namespace Dumux { - -#ifndef DOXYGEN -namespace Detail { -// helper struct detecting if the user-defined spatial params class has a lameParamsAtPos function -// for g++ > 5.3, this can be replaced by a lambda -template -struct hasLameParamsAtPos -{ - template - auto operator()(const SpatialParams& a) - -> decltype(a.lameParamsAtPos(std::declval())) - {} -}; - -} // end namespace Detail -#endif - -/*! - * \ingroup SpatialParameters - * \brief The base class for spatial parameters of linear elastic geomechanical problems - */ -template -class FVSpatialParamsElastic -{ - using FVElementGeometry = typename GridGeometry::LocalView; - using SubControlVolume = typename GridGeometry::SubControlVolume; - using GridView = typename GridGeometry::GridView; - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - - enum { dimWorld = GridView::dimensionworld }; - -public: - //! The constructor - FVSpatialParamsElastic(std::shared_ptr gridGeometry) - : gridGeometry_(gridGeometry) - , gravity_(0.0) - { - const bool enableGravity = getParam("Problem.EnableGravity"); - if (enableGravity) - gravity_[dimWorld-1] = -9.81; - } - - /*! - * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$. - * - * The default behaviour is a constant gravity vector; - * if the Problem.EnableGravity parameter is true, - * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, - * else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$. - * - * \param pos the spatial position at which to evaulate the gravity vector - */ - const GlobalPosition& gravity(const GlobalPosition &pos) const - { return gravity_; } - - /*! - * \brief Function for defining the solid volume fraction. - * That is possibly solution dependent. - * - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - * \param compIdx The solid component index - * \return the volume fraction of the inert solid component with index compIdx - */ - template - Scalar inertVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { - static_assert(SolidSystem::isInert(), "Elastic model can only be used with inert solid systems"); - - // when there is only one component, the volume fraction must be one - if (SolidSystem::numInertComponents == 1) - return 1.0; - - // otherwise we require the user to define the solid composition - return asImp_().template inertVolumeFractionAtPos(scv.center(), compIdx); - } - - /*! - * \brief Function for defining the solid volume fraction. - * That is possibly solution dependent. - * - * \param globalPos The global position - * \param compIdx The solid component index - * \return the volume fraction of the inert solid component with index compIdx - */ - template - Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const - { DUNE_THROW(Dune::InvalidStateException, "The spatial parameters do not provide inertVolumeFractionAtPos() method."); } - - /*! - * \brief Define the Lame parameters - * \note These are possibly solution dependent and are evaluated - * for an integration point inside the element. Therefore, - * a flux variables cache object is passed to this function - * containing data on shape functions at the integration point. - * - * \param element The current element - * \param fvGeometry The local finite volume geometry - * \param elemVolVars Primary/Secondary variables inside the element - * \param fluxVarsCache Contains data on shape functions at the integration point - * \return lame parameters - */ - template - decltype(auto) lameParams(const Element& element, - const FVElementGeometry& fvGeometry, - const ElemVolVars& elemVolVars, - const FluxVarsCache& fluxVarsCache) const - { - static_assert(decltype(isValid(Detail::hasLameParamsAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const\n\n" - " or overload this function\n\n" - " template\n" - " const LameParams& lameParams(const Element& element,\n" - " const FVElementGeometry& fvGeometry,\n" - " const ElemVolVars& elemVolVars,\n" - " const FluxVarsCache& fluxVarsCache) const\n\n"); - - return asImp_().lameParamsAtPos(fluxVarsCache.ipGlobal()); - } - - //! The finite volume grid geometry - const GridGeometry& gridGeometry() const - { return *gridGeometry_; } - -protected: - Implementation &asImp_() - { return *static_cast(this); } - - const Implementation &asImp_() const - { return *static_cast(this); } - -private: - std::shared_ptr gridGeometry_; - GlobalPosition gravity_; //!< The gravity vector -}; -} // end namespace Dumuxs #endif diff --git a/dumux/material/spatialparams/fvporoelastic.hh b/dumux/material/spatialparams/fvporoelastic.hh index 22e7532cef78eb7c02174cc7ac7dba6bbd72806b..b1d6df31161513511f04cc4ad4fc70366f3a242d 100644 --- a/dumux/material/spatialparams/fvporoelastic.hh +++ b/dumux/material/spatialparams/fvporoelastic.hh @@ -21,312 +21,10 @@ * \ingroup SpatialParameters * \brief The base class for spatial parameters of poro-elastic geomechanical problems */ -#ifndef DUMUX_GEOMECHANICS_POROELASTIC_FV_SPATIAL_PARAMS_HH -#define DUMUX_GEOMECHANICS_POROELASTIC_FV_SPATIAL_PARAMS_HH +#ifndef DUMUX_MATERIAL_POROELASTIC_SPATIAL_PARAMS_HH +#define DUMUX_MATERIAL_POROELASTIC_SPATIAL_PARAMS_HH -#include +#warning "This header will be removed after 3.5 in favor of dumux/geomechanics/poroelastic/fvspatialparams.hh" +#include -#include -#include -#include - -namespace Dumux { - -#ifndef DOXYGEN -namespace Detail { -// helper struct detecting if the user-defined spatial params class has a lameParamsAtPos function -// for g++ > 5.3, this can be replaced by a lambda -template -struct hasLameParamsAtPos -{ - template - auto operator()(const SpatialParams& a) - -> decltype(a.lameParamsAtPos(std::declval())) - {} -}; - -// helper struct detecting if the user-defined spatial params class has a reactiveVolumeFractionAtPos function -// for g++ > 5.3, this can be replaced by a lambda -template -struct hasReactiveVolumeFractionAtPos -{ - template - auto operator()(const SpatialParams& a) - -> decltype(a.template reactiveVolumeFractionAtPos(std::declval(), 0)) - {} -}; - -// helper struct detecting if the user-defined spatial params class has a biotCoefficientAtPos function -// for g++ > 5.3, this can be replaced by a lambda -template -struct hasBiotCoeffAtPos -{ - template - auto operator()(const SpatialParams& a) - -> decltype(a.biotCoefficientAtPos(std::declval())) - {} -}; - -} // end namespace Detail -#endif - -/*! - * \ingroup SpatialParameters - * \brief The base class for spatial parameters of poro-elastic geomechanical problems - */ -template -class FVSpatialParamsPoroElastic -{ - using FVElementGeometry = typename GridGeometry::LocalView; - using SubControlVolume = typename GridGeometry::SubControlVolume; - using GridView = typename GridGeometry::GridView; - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - - enum { dimWorld = GridView::dimensionworld }; - -public: - //! The constructor - FVSpatialParamsPoroElastic(std::shared_ptr gridGeometry) - : gridGeometry_(gridGeometry) - , gravity_(0.0) - { - const bool enableGravity = getParam("Problem.EnableGravity"); - if (enableGravity) - gravity_[dimWorld-1] = -9.81; - } - - /*! - * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$. - * - * The default behaviour is a constant gravity vector; - * if the Problem.EnableGravity parameter is true, - * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, - * else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$. - * - * \param pos the spatial position at which to evaulate the gravity vector - */ - const GlobalPosition& gravity(const GlobalPosition &pos) const - { return gravity_; } - - /*! - * \brief Function for defining the porosity. - * That is possibly solution dependent. - * \note this can only be used for solids with one inert component - * (see inertVolumeFraction for the more general interface) - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - * \return the porosity - */ - template - Scalar porosity(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol) const - { - static_assert(decltype(isValid(Detail::hasPorosityAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " Scalar porosityAtPos(const GlobalPosition& globalPos) const\n\n" - " or overload this function\n\n" - " template\n" - " Scalar porosity(const Element& element,\n" - " const SubControlVolume& scv,\n" - " const ElementSolution& elemSol) const\n\n"); - - return asImp_().porosityAtPos(scv.center()); - } - - /*! - * \brief Function for defining the solid volume fraction. - * That is possibly solution dependent. - * - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - * \param compIdx The solid component index - * \return the volume fraction of the inert solid component with index compIdx - * - * \note this overload is enabled if there is only one inert solid component and the - * user didn't choose to implement an inertVolumeFractionAtPos overload. - * It then forwards to the simpler porosity interface. - * With more than one solid components or active solid components (i.e. dissolution) - * please overload the more general inertVolumeFraction/inertVolumeFractionAtPos interface. - */ - template()) - (std::declval()))::value, int> = 0> - Scalar inertVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { return 1.0 - asImp_().porosity(element, scv, elemSol); } - - // specialization if there are no inert components at all - template = 0> - Scalar inertVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { return 0.0; } - - // the more general interface forwarding to inertVolumeFractionAtPos - template 1) || - ( - (SolidSystem::numInertComponents > 0) && - ( - !SolidSystem::isInert() - || decltype(isValid(Detail::hasInertVolumeFractionAtPos()) - (std::declval()))::value - ) - ), int> = 0> - Scalar inertVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { - static_assert(decltype(isValid(Detail::hasInertVolumeFractionAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " template\n" - " Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const\n\n" - " or overload this function\n\n" - " template\n" - " Scalar inertVolumeFraction(const Element& element,\n" - " const SubControlVolume& scv,\n" - " const ElementSolution& elemSol,\n" - " int compIdx) const\n\n"); - - return asImp_().template inertVolumeFractionAtPos(scv.center(), compIdx); - } - - /*! - * \brief Function for defining the solid volume fraction of a solid - * component that takes part in some sort of reaction. - * - * \param element The current element - * \param scv The sub-control volume inside the element. - * \param elemSol The solution at the dofs connected to the element. - * \param compIdx The solid component index - * \return the volume fraction of the inert solid component with index compIdx - * - * \note This overload is enabled if there are only inert solid components - * and the user did not choose to implement a reactiveVolumeFractionAtPos - * function. The reactive volume fraction is zero in this case. - */ - template()) - (std::declval()))::value, int > = 0 > - Scalar reactiveVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { return 0.0; } - - //! overload for the case of reactive solids or user-provided overload - template()) - (std::declval()))::value, int > = 0 > - Scalar reactiveVolumeFraction(const Element& element, - const SubControlVolume& scv, - const ElementSolution& elemSol, - int compIdx) const - { - static_assert(decltype(isValid(Detail::hasReactiveVolumeFractionAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " template\n" - " Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const\n\n" - " or overload this function\n\n" - " template\n" - " Scalar inertVolumeFraction(const Element& element,\n" - " const SubControlVolume& scv,\n" - " const ElementSolution& elemSol,\n" - " int compIdx) const\n\n"); - - return asImp_().template reactiveVolumeFractionAtPos(scv.center(), compIdx); - } - - /*! - * \brief Define the Lame parameters - * \note These are possibly solution dependent and are evaluated - * for an integration point inside the element. Therefore, - * a flux variables cache object is passed to this function - * containing data on shape functions at the integration point. - * - * \param element The current element - * \param fvGeometry The local finite volume geometry - * \param elemVolVars Primary/Secondary variables inside the element - * \param fluxVarsCache Contains data on shape functions at the integration point - * \return lame parameters - */ - template - decltype(auto) lameParams(const Element& element, - const FVElementGeometry& fvGeometry, - const ElemVolVars& elemVolVars, - const FluxVarsCache& fluxVarsCache) const - { - static_assert(decltype(isValid(Detail::hasLameParamsAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const\n\n" - " or overload this function\n\n" - " template\n" - " const LameParams& lameParams(const Element& element,\n" - " const FVElementGeometry& fvGeometry,\n" - " const ElemVolVars& elemVolVars,\n" - " const FluxVarsCache& fluxVarsCache) const\n\n"); - - return asImp_().lameParamsAtPos(fluxVarsCache.ipGlobal()); - } - - /*! - * \brief Returns the Biot coefficient in an element - * \note This is possibly solution dependent and is evaluated - * for an integration point inside the element. Therefore, - * a flux variables cache object is passed to this function - * containing data on shape functions at the integration point. - * - * \param element The current element - * \param fvGeometry The local finite volume geometry - * \param elemVolVars Primary/Secondary variables inside the element - * \param fluxVarsCache Contains data on shape functions at the integration point - * \return Biot coefficient - */ - template - Scalar biotCoefficient(const Element& element, - const FVElementGeometry& fvGeometry, - const ElemVolVars& elemVolVars, - const FluxVarsCache& fluxVarsCache) const - { - static_assert(decltype(isValid(Detail::hasBiotCoeffAtPos())(this->asImp_()))::value," \n\n" - " Your spatial params class has to either implement\n\n" - " const LameParams& biotCoefficientAtPos(const GlobalPosition& globalPos) const\n\n" - " or overload this function\n\n" - " template\n" - " const LameParams& biotCoefficient(const Element& element,\n" - " const FVElementGeometry& fvGeometry,\n" - " const ElemVolVars& elemVolVars,\n" - " const FluxVarsCache& fluxVarsCache) const\n\n"); - - return asImp_().biotCoefficientAtPos(fluxVarsCache.ipGlobal()); - } - - //! The finite volume grid geometry - const GridGeometry& gridGeometry() const - { return *gridGeometry_; } - -protected: - Implementation &asImp_() - { return *static_cast(this); } - - const Implementation &asImp_() const - { return *static_cast(this); } - -private: - std::shared_ptr gridGeometry_; - GlobalPosition gravity_; //!< The gravity vector -}; -} // end namespace Dumux #endif diff --git a/dumux/porousmediumflow/fv1pconstantspatialparams.hh b/dumux/porousmediumflow/fv1pconstantspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..1631219f129d0dd746d6a9daab97caf627203d84 --- /dev/null +++ b/dumux/porousmediumflow/fv1pconstantspatialparams.hh @@ -0,0 +1,73 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup PorousMediumFlow + * \ingroup SpatialParameters + * \brief A spatial params implementation for 1p problem with constant properties + */ +#ifndef DUMUX_POROUS_MEDIUM_FLOW_FV_CONSTANT_SPATIAL_PARAMS_ONE_P_HH +#define DUMUX_POROUS_MEDIUM_FLOW_FV_CONSTANT_SPATIAL_PARAMS_ONE_P_HH + +#include +#include "fvspatialparams.hh" + +namespace Dumux { + +/*! + * \ingroup SpatialParameters + * \brief A spatial params implementation for 1p problem with constant properties + */ +template +class FVPorousMediumOnePConstantSpatialParams +: public FVPorousMediumSpatialParamsOneP> +{ + using ThisType = FVPorousMediumOnePConstantSpatialParams; + using ParentType = FVPorousMediumSpatialParamsOneP; + using GlobalPosition = typename GridGeometry::GridView::template Codim<0>::Geometry::GlobalCoordinate; + +public: + using PermeabilityType = Scalar; + + FVPorousMediumOnePConstantSpatialParams(std::shared_ptr gridGeometry) + : ParentType(gridGeometry) + , porosity_(getParam("SpatialParams.Porosity")) + , permeability_(getParam("SpatialParams.Permeability")) + {} + + /*! + * \brief The (intrinsic) permeability \f$[m^2]\f$ + */ + PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const + { return permeability_; } + + /*! + * \brief The porosity \f$[-]\f$ + */ + Scalar porosityAtPos(const GlobalPosition& globalPos) const + { return porosity_; } + +private: + const Scalar porosity_; + const Scalar permeability_; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/porousmediumflow/fvspatialparams.hh b/dumux/porousmediumflow/fvspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..d1fd3293dd2ca405292c306f753c70a34ad8fe14 --- /dev/null +++ b/dumux/porousmediumflow/fvspatialparams.hh @@ -0,0 +1,411 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup SpatialParameters + * \brief The base class for spatial parameters in porous-medium-flow problems. + */ +#ifndef DUMUX_POROUS_MEDIUM_FLOW_FV_SPATIAL_PARAMS_HH +#define DUMUX_POROUS_MEDIUM_FLOW_FV_SPATIAL_PARAMS_HH + +#include +#include + +#include +#include +#include +#include + +namespace Dumux { + +#ifndef DOXYGEN +namespace Detail { +// helper struct detecting if the user-defined spatial params class has a permeabilityAtPos function +// for g++ > 5.3, this can be replaced by a lambda +template +struct hasPermeabilityAtPos +{ + template + auto operator()(const SpatialParams& a) + -> decltype(a.permeabilityAtPos(std::declval())) + {} +}; + +template +struct hasInertVolumeFractionAtPos +{ + template + auto operator()(const SpatialParams& a) + -> decltype(a.template inertVolumeFractionAtPos(std::declval(), 0)) + {} +}; + +template +struct hasPorosityAtPos +{ + template + auto operator()(const SpatialParams& a) + -> decltype(a.porosityAtPos(std::declval())) + {} +}; + +// helper struct detecting if the user-defined spatial params class +// has a fluidMatrixInteractionAtPos function +template +struct hasFluidMatrixInteractionAtPos +{ + template + auto operator()(const SpatialParams& a) + -> decltype(a.fluidMatrixInteractionAtPos(std::declval())) + {} +}; +} // end namespace Detail +#endif + +/*! + * \ingroup SpatialParameters + * \brief The base class for spatial parameters of single-phase problems + */ +template +class FVPorousMediumSpatialParamsOneP +: public FVSpatialParamsBase +{ + using ParentType = FVSpatialParamsBase; + using GridView = typename GridGeometry::GridView; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Element = typename GridView::template Codim<0>::Entity; + + enum { dim = GridView::dimension }; + enum { dimWorld = GridView::dimensionworld }; + using DimWorldMatrix = Dune::FieldMatrix; + + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + FVPorousMediumSpatialParamsOneP(std::shared_ptr gridGeometry) + : ParentType(gridGeometry) + { + /* \brief default forchheimer coefficient + * Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90 \cite ward1964 . + * Actually the Forchheimer coefficient is also a function of the dimensions of the + * porous medium. Taking it as a constant is only a first approximation + * (Nield, Bejan, Convection in porous media, 2006, p. 10 \cite nield2006 ) + */ + forchCoeffDefault_ = getParam("SpatialParams.ForchCoeff", 0.55); + } + + /*! + * \brief Harmonic average of a discontinuous scalar field at discontinuity interface + * (for compatibility reasons with the function below) + * \return the averaged scalar + * \param T1 first scalar parameter + * \param T2 second scalar parameter + * \param normal The unit normal vector of the interface + */ + Scalar harmonicMean(const Scalar T1, + const Scalar T2, + const GlobalPosition& normal) const + { return Dumux::harmonicMean(T1, T2); } + + /*! + * \brief Harmonic average of a discontinuous tensorial field at discontinuity interface + * \note We do a harmonic average of the part normal to the interface (alpha*I) and + * an arithmetic average of the tangential part (T - alpha*I). + * \return the averaged tensor + * \param T1 first tensor + * \param T2 second tensor + * \param normal The unit normal vector of the interface + */ + DimWorldMatrix harmonicMean(const DimWorldMatrix& T1, + const DimWorldMatrix& T2, + const GlobalPosition& normal) const + { + // determine nT*k*n + GlobalPosition tmp; + GlobalPosition tmp2; + T1.mv(normal, tmp); + T2.mv(normal, tmp2); + const Scalar alpha1 = tmp*normal; + const Scalar alpha2 = tmp2*normal; + + const Scalar alphaHarmonic = Dumux::harmonicMean(alpha1, alpha2); + const Scalar alphaAverage = 0.5*(alpha1 + alpha2); + + DimWorldMatrix T(0.0); + for (int i = 0; i < dimWorld; ++i) + { + for (int j = 0; j < dimWorld; ++j) + { + T[i][j] += 0.5*(T1[i][j] + T2[i][j]); + if (i == j) + T[i][j] += alphaHarmonic - alphaAverage; + } + } + + return T; + } + + /*! + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$ + * \note It is possibly solution dependent. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return permeability + */ + template + decltype(auto) permeability(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + static_assert(decltype(isValid(Detail::hasPermeabilityAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " const PermeabilityType& permeabilityAtPos(const GlobalPosition& globalPos) const\n\n" + " or overload this function\n\n" + " template\n" + " const PermeabilityType& permeability(const Element& element,\n" + " const SubControlVolume& scv,\n" + " const ElementSolution& elemSol) const\n\n"); + + return asImp_().permeabilityAtPos(scv.center()); + } + + /*! + * \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. + * \note this can only be used for solids with one inert component + * (see inertVolumeFraction for the more general interface) + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return the porosity + */ + template + Scalar porosity(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + static_assert(decltype(isValid(Detail::hasPorosityAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " Scalar porosityAtPos(const GlobalPosition& globalPos) const\n\n" + " or overload this function\n\n" + " template\n" + " Scalar porosity(const Element& element,\n" + " const SubControlVolume& scv,\n" + " const ElementSolution& elemSol) const\n\n"); + + return asImp_().porosityAtPos(scv.center()); + } + + /*! + * \brief Function for defining the solid volume fraction. + * That is possibly solution dependent. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \param compIdx The solid component index + * \return the volume fraction of the inert solid component with index compIdx + * + * \note this overload is enable if there is only one inert solid component and the + * user didn't choose to implement a inertVolumeFractionAtPos overload. + * It then forwards to the simpler porosity interface. + * With more than one solid components or active solid components (i.e. dissolution) + * please overload the more general inertVolumeFraction/inertVolumeFractionAtPos interface. + */ + template())(std::declval()))::value, + int> = 0> + Scalar inertVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { + return 1.0 - asImp_().porosity(element, scv, elemSol); + } + + // specialization if there are no inert components at all + template = 0> + Scalar inertVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { + return 0.0; + } + + // the more general interface forwarding to inertVolumeFractionAtPos + template 1) || + ( + (SolidSystem::numInertComponents > 0) && + ( + !SolidSystem::isInert() + || decltype(isValid(Detail::hasInertVolumeFractionAtPos()) + (std::declval()))::value + ) + ), + int> = 0> + Scalar inertVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { + static_assert(decltype(isValid(Detail::hasInertVolumeFractionAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " template\n" + " Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const\n\n" + " or overload this function\n\n" + " template\n" + " Scalar inertVolumeFraction(const Element& element,\n" + " const SubControlVolume& scv,\n" + " const ElementSolution& elemSol,\n" + " int compIdx) const\n\n"); + + return asImp_().template inertVolumeFractionAtPos(scv.center(), compIdx); + } + + /*! + * \brief Function for defining the Beavers-Joseph coefficient for multidomain + * problems\f$\mathrm{[-]}\f$. + * + * \return Beavers-Joseph coefficient \f$\mathrm{[-]}\f$ + * \param globalPos The global position + */ + Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const + { + DUNE_THROW(Dune::InvalidStateException, + "The spatial parameters do not provide a beaversJosephCoeffAtPos() method."); + } + + /*! + * \brief Apply the Forchheimer coefficient for inertial forces + * calculation. + * \param scvf The sub-control volume face where the + * intrinsic velocity ought to be calculated. + */ + Scalar forchCoeff(const SubControlVolumeFace &scvf) const + { + return forchCoeffDefault_; + } + +protected: + Implementation &asImp_() + { return *static_cast(this); } + + const Implementation &asImp_() const + { return *static_cast(this); } + +private: + Scalar forchCoeffDefault_; +}; + +/*! + * \ingroup SpatialParameters + * \brief The base class for spatial parameters of multi-phase problems + */ +template +class FVPorousMediumSpatialParams +: public FVPorousMediumSpatialParamsOneP +{ + using ParentType = FVPorousMediumSpatialParamsOneP; + using GridView = typename GridGeometry::GridView; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using Element = typename GridView::template Codim<0>::Entity; + + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + FVPorousMediumSpatialParams(std::shared_ptr gridGeometry) + : ParentType(gridGeometry) + {} + + /*! + * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.). + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + */ + template + decltype(auto) fluidMatrixInteraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + static_assert(decltype(isValid(Detail::hasFluidMatrixInteractionAtPos())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const\n\n" + " or overload this function\n\n" + " template\n" + " auto fluidMatrixInteraction(const Element& element,\n" + " const SubControlVolume& scv,\n" + " const ElementSolution& elemSol) const\n\n"); + + return this->asImp_().fluidMatrixInteractionAtPos(scv.center()); + } + + /*! + * \brief Function for defining which phase is to be considered as the wetting phase. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return the wetting phase index + */ + template + int wettingPhase(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + return this->asImp_().template wettingPhaseAtPos(scv.center()); + } + + /*! + * \brief Function for defining which phase is to be considered as the wetting phase. + * + * \return the wetting phase index + * \param globalPos The global position + */ + template + int wettingPhaseAtPos(const GlobalPosition& globalPos) const + { + DUNE_THROW(Dune::InvalidStateException, + "The spatial parameters do not provide " + "a wettingPhaseAtPos() method."); + } +}; + +} // namespace Dumux + +#endif diff --git a/dumux/porousmediumflow/problem.hh b/dumux/porousmediumflow/problem.hh index 801f8463403e3953ab4d7749220c1e3ee91ad7c3..baecbe6eb8155cd6d660c8942d01f9ab92d38e34 100644 --- a/dumux/porousmediumflow/problem.hh +++ b/dumux/porousmediumflow/problem.hh @@ -50,44 +50,11 @@ class PorousMediumFlowProblem : public FVProblem using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - using GravityVector = Dune::FieldVector; - public: - //! Export spatial parameter type - using SpatialParams = GetPropType; - /*! - * \brief Constructor, passing the spatial parameters. - * - * \param gridGeometry The finite volume grid geometry - * \param spatialParams The spatial parameter class - * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") - */ - PorousMediumFlowProblem(std::shared_ptr gridGeometry, - std::shared_ptr spatialParams, - const std::string& paramGroup = "") - : ParentType(gridGeometry, paramGroup) - , gravity_(0.0) - , spatialParams_(spatialParams) - { - const bool enableGravity = getParamFromGroup(paramGroup, "Problem.EnableGravity"); - if (enableGravity) - gravity_[dimWorld-1] = -9.81; - } - - /*! - * \brief Constructor, constructing the spatial parameters. - * - * \param gridGeometry The finite volume grid geometry - * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") - */ - PorousMediumFlowProblem(std::shared_ptr gridGeometry, - const std::string& paramGroup = "") - : PorousMediumFlowProblem(gridGeometry, - std::make_shared(gridGeometry), - paramGroup) - {} + //! Use constructors of the base class + using ParentType::ParentType; /*! * \name Physical parameters for porous media problems @@ -117,26 +84,7 @@ public: DUNE_THROW(Dune::NotImplemented, "temperature() method not implemented by the user problem"); } - /*! - * \brief Returns the spatial parameters object. - */ - SpatialParams &spatialParams() - { return *spatialParams_; } - - /*! - * \brief Returns the spatial parameters object. - */ - const SpatialParams &spatialParams() const - { return *spatialParams_; } - // \} - -protected: - //! The gravity acceleration vector - GravityVector gravity_; - - // material properties of the porous medium - std::shared_ptr spatialParams_; }; } // end namespace Dumux diff --git a/test/discretization/test_fvgridvariables.cc b/test/discretization/test_fvgridvariables.cc index 5190170e01bc2e56a6556663dedbf5cafc07c7fb..91c0f8ce4fa5e2dc3e516829d71e63ba24da934a 100644 --- a/test/discretization/test_fvgridvariables.cc +++ b/test/discretization/test_fvgridvariables.cc @@ -27,6 +27,7 @@ #include #include #include +#include // we use the 1p type tag here in order not to be obliged // to define grid flux vars cache & vol vars cache... @@ -98,7 +99,10 @@ int main (int argc, char *argv[]) Dune::MPIHelper::instance(argc, argv); using namespace Dumux; - Dumux::Parameters::init(argc, argv); + Dumux::Parameters::init(argc, argv, [] (Dune::ParameterTree& tree) { + tree["SpatialParams.Porosity"] = "0.1"; + tree["SpatialParams.Permeability"] = "1e-12"; + }); using TypeTag = Properties::TTag::GridVariablesTestBox; using Grid = GetPropType;