From c6f89f953afd20938742010451f4a7186df0dbd8 Mon Sep 17 00:00:00 2001 From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de> Date: Fri, 22 Dec 2017 15:31:23 +0100 Subject: [PATCH] [mpfa] template classes by interaction volume instead of MpfaMethod We now do not specialize the interaction volume depending on the method chosen but extract the chosen method from the interaction volume type that has been set as property. This enables us to hand in the traits class to the interaction volume from the outside, which is actually what you want. Now, for certain grids one can use static matrix/vector types to speed up simulations. --- .../cellcentered/mpfa/CMakeLists.txt | 1 - .../cellcentered/mpfa/helper.hh | 2 +- .../cellcentered/mpfa/interactionvolume.hh | 50 ----------- .../mpfa/interactionvolumebase.hh | 56 +++++++----- .../mpfa/interactionvolumedatahandle.hh | 68 ++++++++------- .../cellcentered/mpfa/localassembler.hh | 10 +-- .../mpfa/omethod/interactionvolume.hh | 87 +++++++++---------- .../mpfa/omethod/localassembler.hh | 74 ++++++++-------- .../mpfa/omethod/localsubcontrolentities.hh | 8 +- .../cellcentered/mpfa/properties.hh | 26 +++--- .../cellcentered/mpfa/subcontrolvolumeface.hh | 42 ++++----- 11 files changed, 193 insertions(+), 231 deletions(-) delete mode 100644 dumux/discretization/cellcentered/mpfa/interactionvolume.hh diff --git a/dumux/discretization/cellcentered/mpfa/CMakeLists.txt b/dumux/discretization/cellcentered/mpfa/CMakeLists.txt index 4390abd49b..86256a3acb 100644 --- a/dumux/discretization/cellcentered/mpfa/CMakeLists.txt +++ b/dumux/discretization/cellcentered/mpfa/CMakeLists.txt @@ -18,7 +18,6 @@ gridinteractionvolumeindexsets.hh helper.hh interactionvolumebase.hh interactionvolumedatahandle.hh -interactionvolume.hh localassembler.hh localfacedata.hh methods.hh diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh index 3de020e6b8..e4ce6f1a1c 100644 --- a/dumux/discretization/cellcentered/mpfa/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/helper.hh @@ -567,7 +567,7 @@ class CCMpfaHelperImplementation : public MpfaDimensionHelper<TypeTag, dim, dimW using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType; - using GlobalPosition = typename InteractionVolume::Traits::GlobalPosition; + using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; using CoordScalar = typename GridView::ctype; using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>; diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/interactionvolume.hh deleted file mode 100644 index 0923b59f2d..0000000000 --- a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh +++ /dev/null @@ -1,50 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- -// vi: set et ts=4 sw=4 sts=4: -/***************************************************************************** - * See the file COPYING for full copying permissions. * - * * - * This program is free software: you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation, either version 2 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program. If not, see <http://www.gnu.org/licenses/>. * - *****************************************************************************/ -/*! - * \file - * \ingroup CCMpfaDiscretization - * \brief Class used for interaction volumes in mpfa schemes. - */ -#ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUME_HH -#define DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUME_HH - -#include <dumux/common/properties.hh> -#include "methods.hh" - -namespace Dumux -{ - -// forward declaration of the method-specific implementation -template<class TypeTag, MpfaMethods MpfaMethod> -class CCMpfaInteractionVolumeImplementation; - -/*! - * \ingroup CCMpfaDiscretization - * \brief Alias to select the correct implementation of the interactionvolume - * volume. The implementations for the schemes have to be included below. - */ -template<class TypeTag> -using CCMpfaInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>; - -} // end namespace Dumux - -// the specializations of this class for the available methods have to be included here -#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh> - -#endif diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh index 38bf007cca..381d501f30 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh @@ -31,27 +31,35 @@ namespace Dumux /* * \ingroup CCMpfaDiscretization - * \brief Type Traits to retrieve types associated with an implementation of a Dumux::CCMpfaInteractionVolumeBase. - * You have to specialize this class for every implementation of Dumux::CCMpfaInteractionVolumeBase. + * \brief Type Traits to retrieve types associated with an implementation of a Dumux::CCMpfaInteractionVolume. + * You have to provide a traits class for every implementation of Dumux::CCMpfaInteractionVolume. Also, + * make sure to publicly export the traits class type in your interaction volume implementation. + * The traits should contain the following type definitions: * * \code - * //! export the interaction volume type to be used on boundaries etc. - * using SecondaryInteractionVolume = ...; + * //! export the problem type (needed for iv-local assembly) + * using Problem = ...; + * //! export the type of the local view on the finite volume grid geometry + * using FVElementGeometry = ...; + * //! export the type of the local view on the grid volume variables + * using ElementVolumeVariables = ...; + * //! export the type of grid view + * using GridView = ...; + * //! export the type used for scalar values + * using ScalarType = ...; + * //! export the type of the mpfa helper class + * using MpfaHelper = ...; * //! export the type used for local indices * using LocalIndexType = ...; - * //! export the type used for indices on the grid - * using GridIndexType = ...; * //! export the type for the interaction volume index set * using IndexSet = ...; - * //! export the type used for global coordinates - * using GlobalPosition = ...; * //! export the type of interaction-volume local scvs * using LocalScvType = ...; * //! export the type of interaction-volume local scvfs * using LocalScvfType = ...; * //! export the type of used for the iv-local face data * using LocalFaceData = ...; - * //! export the type of face data container (use dynamic container here) + * //! export the type of face data container * using LocalFaceDataContainer = ...; * //! export the type used for iv-local matrices * using Matrix = ...; @@ -60,35 +68,35 @@ namespace Dumux * //! export the type used for the iv-stencils * using Stencil = ...; * //! export the data handle type for this iv - * using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >; + * using DataHandle = ...; * \endcode */ -template< class InteractionVolume > -struct CCMpfaInteractionVolumeTraits {}; /*! * \ingroup CCMpfaDiscretization * \brief Base class for the interaction volumes of mpfa methods. It defines * the interface and actual implementations should derive from this class. + * + * \tparam Impl The actual implementation of the interaction volume + * \tparam T The traits class to be used */ -template<class TypeTag, class Implementation> +template< class Impl, class T> class CCMpfaInteractionVolumeBase { // Curiously recurring template pattern - Implementation & asImp() { return static_cast<Implementation&>(*this); } - const Implementation & asImp() const { return static_cast<const Implementation&>(*this); } + Impl & asImp() { return static_cast<Impl&>(*this); } + const Impl & asImp() const { return static_cast<const Impl&>(*this); } + + using Problem = typename T::Problem; + using FVElementGeometry = typename T::FVElementGeometry; + using ElementVolumeVariables = typename T::ElementVolumeVariables; - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using GridView = typename T::GridView; using Element = typename GridView::template Codim<0>::Entity; public: - // Types required in the assembly of the local eq system - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); - //! state the traits type publicly - using Traits = CCMpfaInteractionVolumeTraits< Implementation >; + using Traits = T; //! Prepares everything for the assembly void setUpLocalScope(const typename Traits::IndexSet& indexSet, @@ -132,7 +140,7 @@ public: //! returns the number of interaction volumes living around a vertex template< class NodalIndexSet > static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) - { return Implementation::numInteractionVolumesAtVertex(nodalIndexSet); } + { return Impl::numInteractionVolumesAtVertex(nodalIndexSet); } //! adds the iv index sets living around a vertex to a given container //! and stores the the corresponding index in a map for each scvf @@ -140,7 +148,7 @@ public: static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer, ScvfIndexMap& scvfIndexMap, const NodalIndexSet& nodalIndexSet) - { Implementation::addInteractionVolumeIndexSets(ivIndexSetContainer, scvfIndexMap, nodalIndexSet); } + { Impl::addInteractionVolumeIndexSets(ivIndexSetContainer, scvfIndexMap, nodalIndexSet); } }; } // end namespace Dumux diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh index d1d032e088..f91f5e9e97 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh @@ -33,22 +33,22 @@ namespace Dumux { //! Empty data handle class -template<class InteractionVolume> class EmptyDataHandle { public: + template< class InteractionVolume > void resize(const InteractionVolume& iv) {} }; //! Data handle for quantities related to advection -template<class TypeTag, class InteractionVolume, bool EnableAdvection> +template<class TypeTag, class M, class V, class LI, bool EnableAdvection> class AdvectionDataHandle { // export matrix & vector types from interaction volume - using Matrix = typename InteractionVolume::Traits::Matrix; - using Vector = typename InteractionVolume::Traits::Vector; + using Matrix = M; + using Vector = V; + using LocalIndexType = LI; using Scalar = typename Vector::value_type; - using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType; using OutsideDataContainer = std::vector< std::vector<Vector> >; @@ -59,7 +59,7 @@ class AdvectionDataHandle public: //! prepare data handle for subsequent fill (normal grids) - template< int d = dim, int dw = dimWorld, std::enable_if_t<(d==dw), int> = 0> + template< class InteractionVolume, int d = dim, int dw = dimWorld, std::enable_if_t<(d==dw), int> = 0> void resize(const InteractionVolume& iv) { // resize transmissibility matrix & pressure vectors @@ -79,7 +79,7 @@ public: //! prepare data handle for subsequent fill (surface grids) - template< int d = dim, int dw = dimWorld, std::enable_if_t<(d<dw), int> = 0> + template< class InteractionVolume, int d = dim, int dw = dimWorld, std::enable_if_t<(d<dw), int> = 0> void resize(const InteractionVolume& iv) { static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"); @@ -177,12 +177,12 @@ private: }; //! Data handle for quantities related to diffusion -template<class TypeTag, class InteractionVolume, bool EnableDiffusion> +template<class TypeTag, class M, class V, class LI, bool EnableDiffusion> class DiffusionDataHandle { // export matrix & vector types from interaction volume - using Matrix = typename InteractionVolume::Traits::Matrix; - using Vector = typename InteractionVolume::Traits::Vector; + using Matrix = M; + using Vector = V; using OutsideTContainer = std::vector< std::vector<Vector> >; static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension; @@ -197,6 +197,7 @@ public: void setComponentIndex(unsigned int compIdx) { compIdx_ = compIdx; } //! prepare data handle for subsequent fill + template< class InteractionVolume > void resize(const InteractionVolume& iv) { for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) @@ -245,18 +246,18 @@ private: }; //! Data handle for quantities related to heat conduction -template<class TypeTag, class InteractionVolume, bool EnableHeatConduction> +template<class TypeTag, class M, class V, class LI, bool EnableHeatConduction> class HeatConductionDataHandle { - //! export matrix & vector types from interaction volume - using Matrix = typename InteractionVolume::Traits::Matrix; - using Vector = typename InteractionVolume::Traits::Vector; + using Matrix = M; + using Vector = V; static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension; static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld; public: //! prepare data handle for subsequent fill + template< class InteractionVolume > void resize(const InteractionVolume& iv) { //! resize transmissibility matrix & temperature vector @@ -294,25 +295,34 @@ private: }; //! Process-dependent data handles when related process is disabled -template<class TypeTag, class InteractionVolume> -class AdvectionDataHandle<TypeTag, InteractionVolume, false> : public EmptyDataHandle<InteractionVolume> {}; -template<class TypeTag, class InteractionVolume> -class DiffusionDataHandle<TypeTag, InteractionVolume, false> : public EmptyDataHandle<InteractionVolume> {}; -template<class TypeTag, class InteractionVolume> -class HeatConductionDataHandle<TypeTag, InteractionVolume, false> : public EmptyDataHandle<InteractionVolume> {}; - -//! Interaction volume data handle class -template<class TypeTag, class InteractionVolume> -class InteractionVolumeDataHandle : public AdvectionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableAdvection)>, - public DiffusionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>, - public HeatConductionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)> +template<class TypeTag, class M, class V, class LI> +class AdvectionDataHandle<TypeTag, M, V, LI, false> : public EmptyDataHandle {}; +template<class TypeTag, class M, class V, class LI> +class DiffusionDataHandle<TypeTag, M, V, LI, false> : public EmptyDataHandle {}; +template<class TypeTag, class M, class V, class LI> +class HeatConductionDataHandle<TypeTag, M, V, LI, false> : public EmptyDataHandle {}; + +/*! + * \ingroup CCMpfaDiscretization + * \brief Class for the interaction volume data handle. + * + * \tparam TypeTag The problem TypeTag + * \tparam M The type used for iv-local matrices + * \tparam V The type used for iv-local vectors + * \tparam LI The type used for iv-local indexing + */ +template<class TypeTag, class M, class V, class LI> +class InteractionVolumeDataHandle : public AdvectionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableAdvection)>, + public DiffusionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>, + public HeatConductionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)> { - using AdvectionHandle = AdvectionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableAdvection)>; - using DiffusionHandle = DiffusionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>; - using HeatConductionHandle = HeatConductionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>; + using AdvectionHandle = AdvectionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableAdvection)>; + using DiffusionHandle = DiffusionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>; + using HeatConductionHandle = HeatConductionDataHandle<TypeTag, M, V, LI, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>; public: //! prepare data handles for subsequent fills + template< class InteractionVolume > void resize(const InteractionVolume& iv) { AdvectionHandle::resize(iv); diff --git a/dumux/discretization/cellcentered/mpfa/localassembler.hh b/dumux/discretization/cellcentered/mpfa/localassembler.hh index 3f690d7a0d..2c4aec9157 100644 --- a/dumux/discretization/cellcentered/mpfa/localassembler.hh +++ b/dumux/discretization/cellcentered/mpfa/localassembler.hh @@ -48,17 +48,17 @@ using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< InteractionVo * \tparam Interaction Volume The interaction volume implementation * used by the mpfa scheme */ -template<class InteractionVolume> +template< class InteractionVolume > class InteractionVolumeAssemblerBase { - using Problem = typename InteractionVolume::Problem; - using FVElementGeometry = typename InteractionVolume::FVElementGeometry; - using ElementVolumeVariables = typename InteractionVolume::ElementVolumeVariables; - using Traits = typename InteractionVolume::Traits; using Matrix = typename Traits::Matrix; using Vector = typename Traits::Vector; + using Problem = typename Traits::Problem; + using FVElementGeometry = typename Traits::FVElementGeometry; + using ElementVolumeVariables = typename Traits::ElementVolumeVariables; + public: /*! * \brief The constructor. diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh index 0730d11aca..0e8cd03584 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh @@ -33,7 +33,6 @@ #include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh> #include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh> -#include <dumux/discretization/cellcentered/mpfa/interactionvolume.hh> #include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh> #include <dumux/discretization/cellcentered/mpfa/localfacedata.hh> #include <dumux/discretization/cellcentered/mpfa/methods.hh> @@ -43,37 +42,37 @@ namespace Dumux { -//! Forward declaration of the interaction volume specialization for the mpfa-o scheme -template< class TypeTag > -class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>; +//! Forward declaration of the o-method's interaction volume +template< class Traits > class CCMpfaOInteractionVolume; -//! Specialization of the interaction volume traits class for the mpfa-o method +//! Specialization of the default interaction volume traits class for the mpfa-o method template< class TypeTag > -struct CCMpfaInteractionVolumeTraits< CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> > +struct CCMpfaODefaultInteractionVolumeTraits { private: - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - - using InteractionVolumeType = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>; + using GridIndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType; + static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension; + static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld; public: - //! Per default, we use the same interaction volume everywhere (also on boundaries etc...) - using SecondaryInteractionVolume = InteractionVolumeType; - + //! export the problem type (needed for iv-local assembly) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + //! export the type of the local view on the finite volume grid geometry + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + //! export the type of the local view on the grid volume variables + using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); + //! export the type of grid view + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + //! export the type used for scalar values + using ScalarType = typename GET_PROP_TYPE(TypeTag, Scalar); + //! export the type of the mpfa helper class + using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); //! export the type used for local indices using LocalIndexType = std::uint8_t; - //! export the type used for indices on the grid - using GridIndexType = typename GridView::IndexSet::IndexType; //! export the type for the interaction volume index set using IndexSet = CCMpfaOInteractionVolumeIndexSet< DualGridNodalIndexSet<GridIndexType, LocalIndexType, dim> >; - //! export the type used for global coordinates - using GlobalPosition = Dune::FieldVector< typename GridView::ctype, dimWorld >; //! export the type of interaction-volume local scvs - using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, GlobalPosition, dim >; + using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, ScalarType, dim, dimWorld >; //! export the type of interaction-volume local scvfs using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf< IndexSet >; //! export the type of used for the iv-local face data @@ -81,45 +80,45 @@ public: //! export the type of face data container (use dynamic container here) using LocalFaceDataContainer = std::vector< LocalFaceData >; //! export the type used for iv-local matrices - using Matrix = Dune::DynamicMatrix< Scalar >; + using Matrix = Dune::DynamicMatrix< ScalarType >; //! export the type used for iv-local vectors - using Vector = Dune::DynamicVector< Scalar >; + using Vector = Dune::DynamicVector< ScalarType >; //! export the type used for the iv-stencils using Stencil = std::vector< GridIndexType >; //! export the data handle type for this iv - using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >; + using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >; }; /*! * \ingroup CCMpfaDiscretization * \brief Class for the interaction volume of the mpfa-o method. */ -template<class TypeTag> -class CCMpfaInteractionVolumeImplementation< TypeTag, MpfaMethods::oMethod > - : public CCMpfaInteractionVolumeBase< TypeTag, - CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> > +template< class Traits > +class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits > { - using ThisType = CCMpfaInteractionVolumeImplementation< TypeTag, MpfaMethods::oMethod >; - using ParentType = CCMpfaInteractionVolumeBase< TypeTag, ThisType >; + using ThisType = CCMpfaOInteractionVolume< Traits >; + using ParentType = CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >; - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename Traits::ScalarType; + using Helper = typename Traits::MpfaHelper; + using Problem = typename Traits::Problem; + using FVElementGeometry = typename Traits::FVElementGeometry; + + using GridView = typename Traits::GridView; using Element = typename GridView::template Codim<0>::Entity; static constexpr int dim = GridView::dimension; using DimVector = Dune::FieldVector<Scalar, dim>; - using TraitsType = CCMpfaInteractionVolumeTraits< ThisType >; - using Matrix = typename TraitsType::Matrix; - using LocalScvType = typename TraitsType::LocalScvType; - using LocalScvfType = typename TraitsType::LocalScvfType; - using LocalFaceData = typename TraitsType::LocalFaceData; + using Matrix = typename Traits::Matrix; + using LocalScvType = typename Traits::LocalScvType; + using LocalScvfType = typename Traits::LocalScvfType; + using LocalFaceData = typename Traits::LocalFaceData; - using IndexSet = typename TraitsType::IndexSet; - using GridIndexType = typename TraitsType::GridIndexType; - using LocalIndexType = typename TraitsType::LocalIndexType; - using Stencil = typename TraitsType::Stencil; + using IndexSet = typename Traits::IndexSet; + using GridIndexType = typename GridView::IndexSet::IndexType; + using LocalIndexType = typename Traits::LocalIndexType; + using Stencil = typename Traits::Stencil; //! Data attached to scvf touching Dirichlet boundaries. //! For the default o-scheme, we only store the corresponding vol vars index. @@ -139,10 +138,6 @@ public: //! publicly state the mpfa-scheme this interaction volume is associated with static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod; - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); - //! Sets up the local scope for a given iv index set void setUpLocalScope(const IndexSet& indexSet, const Problem& problem, diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh index 26f922da34..aabba86066 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh @@ -28,10 +28,9 @@ #include <dumux/common/math.hh> #include <dumux/common/properties.hh> -#include <dumux/discretization/cellcentered/mpfa/methods.hh> #include <dumux/discretization/cellcentered/mpfa/localassembler.hh> #include <dumux/discretization/cellcentered/mpfa/computetransmissibility.hh> -#include <dumux/discretization/cellcentered/mpfa/interactionvolume.hh> +#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh> namespace Dumux { @@ -40,24 +39,22 @@ namespace Dumux * \ingroup CCMpfaDiscretization * \brief Specialization of the interaction volume-local * assembler class for the mpfa-o scheme. + * + * \tparam IVTraits The traits class of the interaction volume */ -template< class TypeTag > -class InteractionVolumeAssemblerImpl< CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> > - : public InteractionVolumeAssemblerBase< CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> > +template< class IVTraits > +class InteractionVolumeAssemblerImpl< CCMpfaOInteractionVolume<IVTraits> > + : public InteractionVolumeAssemblerBase< CCMpfaOInteractionVolume<IVTraits> > { - using InteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>; + using InteractionVolume = CCMpfaOInteractionVolume< IVTraits >; using ParentType = InteractionVolumeAssemblerBase< InteractionVolume >; - using Traits = typename InteractionVolume::Traits; - using LocalIndexType = typename Traits::LocalIndexType; - using Matrix = typename Traits::Matrix; - using Vector = typename Traits::Vector; + using LocalIndexType = typename IVTraits::LocalIndexType; + using Matrix = typename IVTraits::Matrix; + using Vector = typename IVTraits::Vector; - static constexpr int dim = Traits::LocalScvType::myDimension; - static constexpr int dimWorld = Traits::LocalScvType::worldDimension; - - // obtain the number of phases via the property system - static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); + static constexpr int dim = IVTraits::LocalScvType::myDimension; + static constexpr int dimWorld = IVTraits::LocalScvType::worldDimension; public: //! Use the constructor of the base class @@ -315,9 +312,13 @@ public: const Matrix& CA, const GetTensorFunction& getTensor) { - // we require the CA matrix and the g vector to have the correct size already - assert(g.size() == numPhases && "Provided gravity container does not have NumPhases entries"); - assert(g[0].size() == iv.numFaces() && "Gravitation vector g does not have the correct size"); + //! We require the gravity container to be a two-dimensional vector/array type, structured as follows: + //! - first index adresses the respective phases + //! - second index adresses the face within the interaction volume + + // make sure CA matrix and the g vector to have the correct size already + assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }) + && "size of gravity vector does not match the number of iv-local faces!"); assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns() && "Matrix CA does not have the correct size"); //! For each face, we... @@ -326,12 +327,10 @@ public: //! - compute \f$ \alpha^* = \alpha_{outside} - \alpha_{inside} \f$ using Scalar = typename Vector::value_type; - std::array< Vector, numPhases > sum_alphas; - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - { - g[pIdx] = 0.0; - sum_alphas[pIdx].resize(iv.numUnknowns(), 0.0); - } + // reset gravity containers to zero + const auto numPhases = g.size(); + std::vector< Vector > sum_alphas(numPhases, Vector(iv.numUnknowns(), 0.0)); + std::for_each(g.begin(), g.end(), [&iv](auto& v) { v = 0.0; }); for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) { @@ -352,7 +351,7 @@ public: // for which there is the specialization of this function below assert(neighborScvIndices.size() <= 2 && "Scvf seems to have more than one outside scv!"); - std::array< Scalar, numPhases > rho; + std::vector< Scalar > rho(numPhases); const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), tensor, gravity); if (!curLocalScvf.isDirichlet()) { @@ -431,14 +430,20 @@ public: const Matrix& A, const GetTensorFunction& getTensor) { + //! We require the gravity container to be a two-dimensional vector/array type, structured as follows: + //! - first index adresses the respective phases + //! - second index adresses the face within the interaction volume + //! We require the outside gravity container to be a three-dimensional vector/array type, structured as follows: + //! - first index adresses the respective phases + //! - second index adresses the face within the interaction volume + //! - third index adresses the i-th "outside" face of the current face + // we require the CA matrix and the gravity containers to have the correct size already assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns() && "Matrix CA does not have the correct size"); - assert(g.size() == numPhases && "Provided gravity container does not have NumPhases entries"); - assert(outsideG.size() == numPhases && "Provided outside gravity container does not have NumPhases entries"); - assert(std::all_of(g.cbegin(), g.cend(), [&iv](const auto& d) { return d.size() == iv.numFaces(); }) - && "Gravitation vector g does not have the correct size"); - assert(std::all_of(outsideG.cbegin(), outsideG.cend(), [&iv](const auto& d) { return d.size() == iv.numFaces(); }) - && "Outside gravity container does not have the correct size"); + assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }) + && "size of gravity container does not match the number of iv-local faces!"); + assert(std::all_of(outsideG.begin(), outsideG.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }) + && "Outside gravity container does not match the number of iv-local faces!"); //! For each face, we... //! - arithmetically average the phase densities @@ -447,12 +452,13 @@ public: using Scalar = typename Vector::value_type; // reset everything to zero - std::array< Vector, numPhases > sum_alphas; + const auto numPhases = g.size(); + std::vector< Vector > sum_alphas(numPhases, Vector(iv.numUnknowns(), 0.0)); + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) { g[pIdx] = 0.0; std::for_each(outsideG[pIdx].begin(), outsideG[pIdx].end(), [] (auto& v) { v = 0.0; }); - sum_alphas[pIdx].resize(iv.numUnknowns(), 0.0); } for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) @@ -473,7 +479,7 @@ public: const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), tensor, gravity); const auto numOutsideFaces = curGlobalScvf.boundary() ? 0 : curGlobalScvf.numOutsideScvs(); std::vector< Scalar > alpha_outside(numOutsideFaces); - std::array< Scalar, numPhases > rho; + std::vector< Scalar > rho(numPhases); if (!curLocalScvf.isDirichlet()) { diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh index a0553b3ef7..51ab083d0f 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh @@ -36,10 +36,10 @@ namespace Dumux * in the mpfa-o scheme. * * \tparam IvIndexSet The type used for index sets within interaction volumes - * \tparam GC The type used for global coordinates * \tparam dim The dimensionality of the grid + * \tparam dimWorld The dimension of the world the grid is embedded in */ -template< class IvIndexSet, class GC, int dim > +template< class IvIndexSet, class Scalar, int dim, int dimWorld> class CCMpfaOInteractionVolumeLocalScv { @@ -47,12 +47,12 @@ public: // export some types using GridIndexType = typename IvIndexSet::GridIndexType; using LocalIndexType = typename IvIndexSet::LocalIndexType; - using GlobalCoordinate = GC; + using GlobalCoordinate = Dune::FieldVector<Scalar, dimWorld>; using ctype = typename GlobalCoordinate::value_type; using LocalBasis = std::array< GlobalCoordinate, dim >; static constexpr int myDimension = dim; - static constexpr int worldDimension = GlobalCoordinate::dimension; + static constexpr int worldDimension = dimWorld; /*! * \brief The constructor diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh index f306b74653..a2f6f6cba8 100644 --- a/dumux/discretization/cellcentered/mpfa/properties.hh +++ b/dumux/discretization/cellcentered/mpfa/properties.hh @@ -40,7 +40,6 @@ #include <dumux/discretization/cellcentered/elementsolution.hh> #include <dumux/discretization/cellcentered/elementboundarytypes.hh> -#include <dumux/discretization/cellcentered/mpfa/methods.hh> #include <dumux/discretization/cellcentered/mpfa/connectivitymap.hh> #include <dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh> #include <dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh> @@ -49,7 +48,8 @@ #include <dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh> #include <dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh> #include <dumux/discretization/cellcentered/mpfa/helper.hh> -#include <dumux/discretization/cellcentered/mpfa/interactionvolume.hh> + +#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh> namespace Dumux { @@ -64,22 +64,29 @@ SET_PROP(CCMpfaModel, DiscretizationMethod) static const DiscretizationMethods value = DiscretizationMethods::CCMpfa; }; -//! By default we set the o-method as the Mpfa method of choice +//! Extract the used mpfa method from the primary interaction volume SET_PROP(CCMpfaModel, MpfaMethod) { - static const MpfaMethods value = MpfaMethods::oMethod; + static const MpfaMethods value = GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume)::MpfaMethod; }; //! The mpfa helper class SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>); -//! The interaction volume class -SET_TYPE_PROP(CCMpfaModel, PrimaryInteractionVolume, CCMpfaInteractionVolume<TypeTag>); +//! Per default, we use the mpfa-o interaction volume +SET_PROP(CCMpfaModel, PrimaryInteractionVolume) +{ +private: + //! use the default traits + using Traits = CCMpfaODefaultInteractionVolumeTraits< TypeTag >; +public: + using type = CCMpfaOInteractionVolume< Traits >; +}; -//! The secondary interaction volume class (exported from the traits) +//! Per default, we use the mpfa-o interaction volume everywhere (pure mpfa-o scheme) SET_TYPE_PROP(CCMpfaModel, SecondaryInteractionVolume, - typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume)::Traits::SecondaryInteractionVolume); + typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume)); //! Set the default for the global finite volume geometry SET_TYPE_PROP(CCMpfaModel, @@ -134,7 +141,6 @@ private: using Grid = typename GET_PROP_TYPE(TypeTag, Grid); static const int dim = Grid::dimension; static const int dimWorld = Grid::dimensionworld; - static const MpfaMethods method = GET_PROP_VALUE(TypeTag, MpfaMethod); // we use geometry traits that use static corner vectors to and a fixed geometry type template <class ct> @@ -168,7 +174,7 @@ private: }; public: - using type = Dumux::CCMpfaSubControlVolumeFace<method, ScvfGeometryTraits>; + using type = Dumux::CCMpfaSubControlVolumeFace< ScvfGeometryTraits>; }; //! Set the solution vector type for an element diff --git a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh index 7698ff82cd..45bd3cc55a 100644 --- a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh +++ b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh @@ -28,17 +28,18 @@ #include <dune/common/version.hh> #include <dune/geometry/type.hh> -#include "methods.hh" - namespace Dumux { /*! * \ingroup CCMpfaDiscretization - * \brief Default implementation of the class for a sub-control volume face in mpfa methods. + * \brief Class for a sub control volume face in mpfa methods, i.e a part of the boundary + * of a control volume we compute fluxes on. + * + * \param ScvfGeometryTraits the traits class for the geometry type */ template<class ScvfGeometryTraits> -class CCMpfaDefaultSubControlVolumeFace +class CCMpfaSubControlVolumeFace { using GridIndexType = typename ScvfGeometryTraits::GridIndexType; using Scalar = typename ScvfGeometryTraits::Scalar; @@ -65,16 +66,16 @@ public: * \param boundary Boolean to specify whether or not the scvf is on a boundary */ template<class MpfaHelper> - CCMpfaDefaultSubControlVolumeFace(const MpfaHelper& helper, - CornerStorage&& corners, - GlobalPosition&& unitOuterNormal, - GridIndexType vIdxGlobal, - unsigned int vIdxLocal, - GridIndexType scvfIndex, - GridIndexType insideScvIdx, - const std::vector<GridIndexType>& outsideScvIndices, - Scalar q, - bool boundary) + CCMpfaSubControlVolumeFace(const MpfaHelper& helper, + CornerStorage&& corners, + GlobalPosition&& unitOuterNormal, + GridIndexType vIdxGlobal, + unsigned int vIdxLocal, + GridIndexType scvfIndex, + GridIndexType insideScvIdx, + const std::vector<GridIndexType>& outsideScvIndices, + Scalar q, + bool boundary) : boundary_(boundary) , vertexIndex_(vIdxGlobal) , scvfIndex_(scvfIndex) @@ -170,19 +171,6 @@ private: Scalar area_; }; -/*! - * \ingroup CCMpfaDiscretization - * \brief Class for a sub control volume face in mpfa methods, i.e a part of the boundary - * of a control volume we compute fluxes on. Per default, we use the default - * implementation of the mpfa scvf class. If a scheme requires a different implementation, - * provide a specialization for it. - * - * \param M the mpfa method used - * \param GT the traits class for the geometry type - */ -template<MpfaMethods M, class GT> -using CCMpfaSubControlVolumeFace = CCMpfaDefaultSubControlVolumeFace< GT >; - } // end namespace Dumux #endif -- GitLab