diff --git a/dumux/discretization/cellcentered/mpfa/CMakeLists.txt b/dumux/discretization/cellcentered/mpfa/CMakeLists.txt index 74af590a45f2b5e6d037d176d406526cfac016e3..4390abd49b98bdfd050568567b3d9ba78ffedd0f 100644 --- a/dumux/discretization/cellcentered/mpfa/CMakeLists.txt +++ b/dumux/discretization/cellcentered/mpfa/CMakeLists.txt @@ -1,6 +1,7 @@ add_subdirectory("omethod") install(FILES +computetransmissibility connectivitymap.hh darcyslaw.hh dualgridindexset.hh @@ -18,6 +19,8 @@ helper.hh interactionvolumebase.hh interactionvolumedatahandle.hh interactionvolume.hh +localassembler.hh +localfacedata.hh methods.hh properties.hh subcontrolvolumeface.hh diff --git a/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh new file mode 100644 index 0000000000000000000000000000000000000000..673a9eb5590c7502ea9d737b24db5f1baf32664d --- /dev/null +++ b/dumux/discretization/cellcentered/mpfa/computetransmissibility.hh @@ -0,0 +1,101 @@ +// -*- 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 This file contains free functions to evaluate the transmissibilities + * associated with flux evaluations across sub-control volume faces + * in the context of cell-centered Mpfa schemes. + */ +#ifndef DUMUX_DISCRETIZATION_CC_MPFA_COMPUTE_TRANSMISSIBILITY_HH +#define DUMUX_DISCRETIZATION_CC_MPFA_COMPUTE_TRANSMISSIBILITY_HH + +#include <dune/common/typetraits.hh> +#include <dune/common/fmatrix.hh> +#include <dune/common/fvector.hh> + +#include <dumux/common/math.hh> + +namespace Dumux +{ + +/*! + * \ingroup CCMpfaDiscretization + * \brief Free function to evaluate the Mpfa transmissibility associated + * with the flux (in the form of flux = T*gradU) across a + * sub-control volume face stemming from a given sub-control + * volume with corresponding tensor T. + * + * \param scv The iv-local sub-control volume + * \param scvf The grid sub-control volume face + * \param T The tensor living in the scv + * \param extrusionFactor The extrusion factor of the scv + */ +template<class Scv, class Scvf, class Tensor> +Dune::FieldVector< typename Tensor::field_type, Scv::myDimension > +computeMpfaTransmissibility(const Scv& scv, + const Scvf& scvf, + const Tensor& T, + typename Scv::ctype extrusionFactor) +{ + Dune::FieldVector< typename Tensor::field_type, Scv::myDimension > wijk; + for (unsigned int dir = 0; dir < Scv::myDimension; ++dir) + wijk[dir] = vtmv(scvf.unitOuterNormal(), T, scv.nu(dir)); + + wijk *= scvf.area()*extrusionFactor; + wijk /= scv.detX(); + + return wijk; +} + +/*! + * \ingroup CCMpfaDiscretization + * \brief Free function to evaluate the Mpfa transmissibility associated + * with the flux (in the form of flux = T*gradU) across a + * sub-control volume face stemming from a given sub-control + * volume with corresponding tensor T, where T is a scalar. + * + * \param scv The iv-local sub-control volume + * \param scvf The grid sub-control volume face + * \param t the scalar quantity living in the scv + * \param extrusionFactor The extrusion factor of the scv + */ +template< class Scv, + class Scvf, + class Tensor, + std::enable_if_t< Dune::IsNumber<Tensor>::value, int > = 1 > +Dune::FieldVector<Tensor, Scv::myDimension> +computeMpfaTransmissibility(const Scv& scv, + const Scvf& scvf, + Tensor t, + typename Scv::ctype extrusionFactor) +{ + Dune::FieldVector<Tensor, Scv::myDimension> wijk; + for (unsigned int dir = 0; dir < Scv::myDimension; ++dir) + wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir)); + + wijk *= scvf.area()*extrusionFactor; + wijk /= scv.detX(); + + return wijk; +} + +} // end namespace Dumux + +#endif diff --git a/dumux/discretization/cellcentered/mpfa/connectivitymap.hh b/dumux/discretization/cellcentered/mpfa/connectivitymap.hh index 3c7152e75b9f48069ee935a476535f5718dffda8..95073ed462527cc3659b57b2d5fa9d4694790fc3 100644 --- a/dumux/discretization/cellcentered/mpfa/connectivitymap.hh +++ b/dumux/discretization/cellcentered/mpfa/connectivitymap.hh @@ -18,8 +18,11 @@ *****************************************************************************/ /*! * \file + * \ingroup CCMpfaDiscretization * \brief Stores the face indices corresponding to the neighbors of an element - * that contribute to the derivative calculation + * that contribute to the derivative calculation. Depending on if an + * mpfa scheme leads to a symmetric/unsymmetric sparsity pattern, the + * adequate implementation of the connectiviy map is chosen. */ #ifndef DUMUX_CC_MPFA_CONNECTIVITY_MAP_HH #define DUMUX_CC_MPFA_CONNECTIVITY_MAP_HH @@ -35,15 +38,15 @@ namespace Dumux template<class TypeTag, MpfaMethods method> class CCMpfaConnectivityMapImplementation; -// //! The Assembly map for models using mpfa methods +//! The Assembly map for models using mpfa methods template<class TypeTag> using CCMpfaConnectivityMap = CCMpfaConnectivityMapImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>; -//! The default is the general assembly map for mpfa schemes +//! The default is the general assembly map for mpfa schemes (non-symmetric schemes) template<class TypeTag, MpfaMethods method> class CCMpfaConnectivityMapImplementation : public CCMpfaGeneralConnectivityMap<TypeTag> {}; -//! The o-method can use the simple assembly map +//! The o-method can use the simple (symmetric) assembly map template<class TypeTag> class CCMpfaConnectivityMapImplementation<TypeTag, MpfaMethods::oMethod> : public CCSimpleConnectivityMap<TypeTag> {}; } diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh index 7aa203b3367c276ae4b0569c307d7454e746ae9d..1cbe5ee34f569b0da2d9bd7a7b52f83c13f04816 100644 --- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh @@ -18,27 +18,29 @@ *****************************************************************************/ /*! * \file - * \brief This file contains the class which is required to calculate - * volume and mass fluxes of fluid phases over a face of a finite volume by means - * of the Darcy approximation. This specializations is for cell-centered schemes - * using multi-point flux approximation. + * \ingroup CCMpfaDiscretization + * \brief Darcy's Law for cell-centered finite volume schemes + * with multi-point flux approximation. */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH #define DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH +#include <dune/common/densevector.hh> +#include <dune/common/densematrix.hh> + #include <dumux/common/properties.hh> #include <dumux/common/parameters.hh> #include <dumux/discretization/methods.hh> namespace Dumux { -// forward declarations +//! forward declaration of the method-specific implementation template<class TypeTag, DiscretizationMethods discMethod> class DarcysLawImplementation; /*! - * \ingroup Mpfa - * \brief Specialization of Darcy's Law for the CCMpfa method. + * \ingroup CCMpfaDiscretization + * \brief Darcy's law for cell-centered finite volume schemes with multi-point flux approximation. */ template <class TypeTag> class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> @@ -47,21 +49,13 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache); using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); - // Always use the dynamic type for vectors (compatibility with the boundary) - using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using CoefficientVector = typename PrimaryInteractionVolume::Traits::DynamicVector; - using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; - - static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - //! Class that fills the cache corresponding to mpfa Darcy's Law class MpfaDarcysLawCacheFiller { @@ -77,14 +71,16 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> const SubControlVolumeFace& scvf, const FluxVariablesCacheFiller& fluxVarsCacheFiller) { - // get interaction volume from the flux vars cache filler & upate the cache + // get interaction volume related data from the filler class & upate the cache if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.secondaryInteractionVolume(), - fluxVarsCacheFiller.dataHandle(), + fluxVarsCacheFiller.secondaryIvLocalFaceData(), + fluxVarsCacheFiller.secondaryIvDataHandle(), scvf); else scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.primaryInteractionVolume(), - fluxVarsCacheFiller.dataHandle(), + fluxVarsCacheFiller.primaryIvLocalFaceData(), + fluxVarsCacheFiller.primaryIvDataHandle(), scvf); } }; @@ -92,54 +88,124 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> //! The cache used in conjunction with the mpfa Darcy's Law class MpfaDarcysLawCache { - // We always use the dynamic types here to be compatible on the boundary - using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer; - using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); + + // In the current implementation of the flux variables cache we cannot + // make a disctinction between dynamic (mpfa-o) and static (mpfa-l) + // matrix and vector types, as currently the cache class can only be templated + // by a type tag (and there can only be one). We use a dynamic vector here to + // make sure it works in case one of the two used interaction volume types uses + // dynamic types performance is thus lowered for schemes using static types. + // TODO: this has to be checked thoroughly as soon as a scheme using static types + // is implemented. One idea to overcome the performance drop could be only + // storing the iv-local index here and obtain tij always from the datahandle + // of the fluxVarsCacheContainer + using Vector = Dune::DynamicVector< Scalar >; + using Matrix = Dune::DynamicMatrix< Scalar >; + + using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); + using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector; + using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix; + + static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value, + "The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" ); + static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value, + "The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" ); + + using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume); + using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector; + using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix; + + static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value, + "The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" ); + static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value, + "The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" ); public: // export the filler type using Filler = MpfaDarcysLawCacheFiller; - //! update cached objects - template<class InteractionVolume> - void updateAdvection(const InteractionVolume& iv, const DataHandle& dataHandle, const SubControlVolumeFace &scvf) + /*! + * \brief Update cached objects (transmissibilities and gravity) + * + * \tparam InteractionVolume The (mpfa scheme-specific) interaction volume + * \tparam LocalFaceData The object used to store iv-local info on an scvf + * \tparam DataHandle The object used to store transmissibility matrices etc. + * + * \param iv The interaction volume this scvf is embedded in + * \param localFaceData iv-local info on this scvf + * \param dataHandle Transmissibility matrix & gravity data of this iv + * \param scvf The sub-control volume face + */ + template< class InteractionVolume, class LocalFaceData, class DataHandle > + void updateAdvection(const InteractionVolume& iv, + const LocalFaceData& localFaceData, + const DataHandle& dataHandle, + const SubControlVolumeFace &scvf) { - const auto& localFaceData = iv.getLocalFaceData(scvf); + switchFluxSign_ = localFaceData.isOutside(); + + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + pj_[pIdx] = &dataHandle.pressures(pIdx); - // update the quantities that are equal for all phases - advectionSwitchFluxSign_ = localFaceData.isOutside(); - advectionVolVarsStencil_ = &dataHandle.volVarsStencil(); - advectionDirichletData_ = &dataHandle.dirichletData(); + static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"); + const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); - // the transmissibilities on surface grids have to be obtained from the outside + // standard grids if (dim == dimWorld) - advectionTij_ = &dataHandle.T()[localFaceData.ivLocalScvfIndex()]; + { + Tij_ = &dataHandle.advectionT()[ivLocalIdx]; + + if (enableGravity) + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx]; + } + + // surface grids else - advectionTij_ = localFaceData.isOutside() ? - &dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()] : - &dataHandle.T()[localFaceData.ivLocalScvfIndex()]; + { + if (!localFaceData.isOutside()) + { + Tij_ = &dataHandle.advectionT()[ivLocalIdx]; + + if (enableGravity) + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx]; + } + else + { + const auto idxInOutsideFaces = localFaceData.scvfLocalOutsideScvfIndex(); + Tij_ = &dataHandle.advectionTout()[ivLocalIdx][idxInOutsideFaces]; + + if (enableGravity) + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + g_[phaseIdx] = dataHandle.gravityOutside(phaseIdx)[ivLocalIdx][idxInOutsideFaces]; + } + } } - //! Returns the stencil for advective scvf flux computation - const Stencil& advectionVolVarsStencil() const { return *advectionVolVarsStencil_; } + //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions + const Vector& advectionTij() const { return *Tij_; } - //! Returns the transmissibilities associated with the volume variables - //! All phases flow through the same rock, thus, tij are equal for all phases - const CoefficientVector& advectionTij() const { return *advectionTij_; } + //! The cell (& Dirichlet) pressures within this interaction volume + const Vector& pressures(unsigned int phaseIdx) const { return *pj_[phaseIdx]; } - //! On faces that are "outside" w.r.t. a face in the interaction volume, - //! we have to take the negative value of the fluxes, i.e. multiply by -1.0 - bool advectionSwitchFluxSign() const { return advectionSwitchFluxSign_; } + //! The gravitational acceleration for a phase on this scvf + Scalar gravity(unsigned int phaseIdx) const { return g_[phaseIdx]; } - //! Returns the data on dirichlet boundary conditions affecting - //! the flux computation on this face - const DirichletDataContainer& advectionDirichletData() const { return *advectionDirichletData_; } + //! In the interaction volume-local system of eq we have one unknown per face. + //! On scvfs on this face, but in "outside" (neighbor) elements of it, we have + //! to take the negative value of the fluxes due to the flipped normal vector. + //! This function returns whether or not this scvf is an "outside" face in the iv. + bool advectionSwitchFluxSign() const { return switchFluxSign_; } private: - bool advectionSwitchFluxSign_; - const Stencil* advectionVolVarsStencil_; - const CoefficientVector* advectionTij_; - const DirichletDataContainer* advectionDirichletData_; + bool switchFluxSign_; + const Vector* Tij_; //!< The transmissibilities such that f = Tij*pj + std::array< Scalar, numPhases > g_; //!< Gravitational flux contribution on this face + std::array<const Vector*, numPhases> pj_; //!< The interaction-volume wide phase pressures pj }; public: @@ -149,6 +215,7 @@ public: // export the type for the corresponding cache using Cache = MpfaDarcysLawCache; + //! Compute the advective flux across an scvf static Scalar flux(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, @@ -157,78 +224,23 @@ public: const unsigned int phaseIdx, const ElementFluxVariablesCache& elemFluxVarsCache) { - static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"); - - // Calculate the interface density for gravity evaluation - const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx); - - // prepare computations - unsigned int i = 0; - Scalar scvfFlux(0.0); const auto& fluxVarsCache = elemFluxVarsCache[scvf]; const auto& tij = fluxVarsCache.advectionTij(); + const auto& pj = fluxVarsCache.pressures(phaseIdx); - // add contributions from cell-centered unknowns - for (const auto volVarIdx : fluxVarsCache.advectionVolVarsStencil()) - { - const auto& volVars = elemVolVars[volVarIdx]; - Scalar h = volVars.pressure(phaseIdx); + //! compute t_ij*p_j + Scalar scvfFlux = tij*pj; - // if gravity is enabled, add gravitational acceleration - if (gravity) - { - // gravitational acceleration in the center of the actual element - const auto x = fvGeometry.scv(volVarIdx).center(); - const auto g = problem.gravityAtPos(x); - h -= rho*(g*x); - } + //! maybe add gravitational acceleration + static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"); + if (enableGravity) + scvfFlux += fluxVarsCache.gravity(phaseIdx); - scvfFlux += tij[i++]*h; - } + //! switch the sign if necessary + if (fluxVarsCache.advectionSwitchFluxSign()) + scvfFlux *= -1.0; - // add contributions from possible dirichlet boundary conditions - for (const auto& d : fluxVarsCache.advectionDirichletData()) - { - const auto& volVars = elemVolVars[d.volVarIndex()]; - Scalar h = volVars.pressure(phaseIdx); - - // maybe add gravitational acceleration - if (gravity) - { - const auto x = d.ipGlobal(); - const auto g = problem.gravityAtPos(x); - h -= rho*(g*x); - } - - scvfFlux += tij[i++]*h; - } - - // return the flux (maybe adjust the sign) - return fluxVarsCache.advectionSwitchFluxSign() ? -scvfFlux : scvfFlux; - } - -private: - static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf, - const unsigned int phaseIdx) - { - static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"); - - if (!gravity) - return Scalar(0.0); - else - { - // use arithmetic mean of the densities around the scvf - if (!scvf.boundary()) - { - Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx); - for (const auto outsideIdx : scvf.outsideScvIndices()) - rho += elemVolVars[outsideIdx].density(phaseIdx); - return rho/(scvf.outsideScvIndices().size()+1); - } - else - return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx); - } + return scvfFlux; } }; diff --git a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh index d6e54463fd3a3fe7a51f615be072e9e33dd95c3f..47471c6415d13ac5600ae0eb32fd1ce62f4e5f3d 100644 --- a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh +++ b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh @@ -18,39 +18,40 @@ *****************************************************************************/ /*! * \file - * \brief Base class for the index sets of the dual grid in mpfa schemes. + * \ingroup CCMpfaDiscretization + * \brief Class for the index sets of the dual grid in mpfa schemes. */ #ifndef DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEXSET_BASE_HH #define DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEXSET_BASE_HH -#include <dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh> - namespace Dumux { /*! - * \ingroup Mpfa - * \brief Nodal index set for the dual grid of mpfa schemes. + * \ingroup CCMpfaDiscretization + * \brief Nodal index set for mpfa schemes, constructed + * around grid vertices. + * + * \tparam GI The type used for indices on the grid + * \tparam LI The type used for indexing in interaction volumes + * \tparam dim The dimension of the grid */ -template<class TypeTag> +template< class GI, class LI, int dim > class DualGridNodalIndexSet { - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; - using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using LocalIndexType = typename PrimaryInteractionVolume::Traits::LocalIndexType; - using LocalIndexContainer = typename PrimaryInteractionVolume::Traits::DynamicLocalIndexContainer; - using GlobalIndexContainer = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer; +public: + using GridIndexType = GI; + using LocalIndexType = LI; - static const int dim = GridView::dimension; + // we use dynamic containers to store indices here + using GridIndexContainer = std::vector<GridIndexType>; + using LocalIndexContainer = std::vector<LocalIndexType>; -public: - //! constructor + //! Constructor DualGridNodalIndexSet() : numBoundaryScvfs_(0) {} - // Inserts data for a given scvf. This should only called once per scvf! + //! Inserts data for a given scvf + template<typename SubControlVolumeFace> void insert(const SubControlVolumeFace& scvf) { insert(scvf.boundary(), @@ -59,20 +60,21 @@ public: scvf.outsideScvIndices()); } - // Inserts data for a given scvf. This should only called once per scvf! + //! Inserts scvf data void insert(const bool boundary, - const IndexType scvfIdx, - const IndexType insideScvIdx, - const std::vector<IndexType>& outsideScvIndices) + const GridIndexType scvfIdx, + const GridIndexType insideScvIdx, + const std::vector<GridIndexType>& outsideScvIndices) { - // this should always be called only once - assert( std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end() && "scvf has already been inserted!"); + // this should always be called only once per scvf + assert(std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end() + && "scvf has already been inserted!"); // the local index of the scvf data about to be inserted const LocalIndexType curScvfLocalIdx = scvfIndices_.size(); - // add global scvf data - GlobalIndexContainer scvIndices; + // add grid scvf data + GridIndexContainer scvIndices; scvIndices.reserve(outsideScvIndices.size()+1); scvIndices.push_back(insideScvIdx); scvIndices.insert(scvIndices.end(), outsideScvIndices.begin(), outsideScvIndices.end()); @@ -86,7 +88,7 @@ public: scvfIsOnBoundary_.push_back(boundary); scvfNeighborScvIndices_.emplace_back( std::move(scvIndices) ); - // if entry for the inside scv exists append scvf local index, create otherwise + // if entry for the inside scv exists, append scvf local index, create otherwise auto it = std::find( scvIndices_.begin(), scvIndices_.end(), insideScvIdx ); if (it != scvIndices_.end()) localScvfIndicesInScv_[ std::distance(scvIndices_.begin(), it) ].push_back(curScvfLocalIdx); @@ -100,105 +102,96 @@ public: } } - //! returns the number of scvs + //! returns the number of scvs around the node std::size_t numScvs() const { return scvIndices_.size(); } - //! returns the number of scvfs + //! returns the number of scvfs around the node std::size_t numScvfs() const { return scvfIndices_.size(); } - //! returns the number of boundary scvfs + //! returns the number of boundary scvfs around the node std::size_t numBoundaryScvfs() const { return numBoundaryScvfs_; } //! returns the global scv indices connected to this dual grid node - const GlobalIndexContainer& globalScvIndices() const { return scvIndices_; } + const GridIndexContainer& globalScvIndices() const { return scvIndices_; } //! returns the global scvf indices connected to this dual grid node - const GlobalIndexContainer& globalScvfIndices() const { return scvfIndices_; } + const GridIndexContainer& globalScvfIndices() const { return scvfIndices_; } + + //! returns whether or not the i-th scvf is on a domain boundary + bool scvfIsOnBoundary(unsigned int i) const { return scvfIsOnBoundary_[i]; } //! returns the global scv idx of the i-th scv - IndexType scvIdxGlobal(unsigned int i) const { return scvIndices_[i]; } + GridIndexType scvIdxGlobal(unsigned int i) const { return scvIndices_[i]; } //! returns the index of the i-th scvf - IndexType scvfIdxGlobal(unsigned int i) const { return scvfIndices_[i]; } + GridIndexType scvfIdxGlobal(unsigned int i) const { return scvfIndices_[i]; } //! returns the global index of the j-th scvf embedded in the i-th scv - IndexType scvfIdxGlobal(unsigned int i, unsigned int j) const + GridIndexType scvfIdxGlobal(unsigned int i, unsigned int j) const { assert(j < dim); // only dim faces can be embedded in an scv return scvfIndices_[ localScvfIndicesInScv_[i][j] ]; } - //! returns the nodel-local index of the j-th scvf embedded in the i-th scv - IndexType scvfIdxLocal(unsigned int i, unsigned int j) const + //! returns the node-local index of the j-th scvf embedded in the i-th scv + LocalIndexType scvfIdxLocal(unsigned int i, unsigned int j) const { assert(j < dim); // only dim faces can be embedded in an scv return localScvfIndicesInScv_[i][j]; } - //! returns whether or not the i-th scvf touches the boundary - bool scvfIsOnBoundary(unsigned int i) const { return scvfIsOnBoundary_[i]; } - //! returns the indices of the neighboring scvs of the i-th scvf - const GlobalIndexContainer& neighboringScvIndices(unsigned int i) const + const GridIndexContainer& neighboringScvIndices(unsigned int i) const { return scvfNeighborScvIndices_[i]; } private: - //! the indices of the scvs around a dual grid node - GlobalIndexContainer scvIndices_; - //! maps to each scv a list of scvf indices embedded in it - std::vector<LocalIndexContainer> localScvfIndicesInScv_; - - //! the indices of the scvfs around a dual grid node - GlobalIndexContainer scvfIndices_; - //! maps to each scvf a boolean to indicate if it is on the boundary - std::vector<bool> scvfIsOnBoundary_; - //! maps to each scvf a list of neighbouring scv indices - //! ordering: 0 - inside scv idx; 1..n - outside scv indices - std::vector<GlobalIndexContainer> scvfNeighborScvIndices_; - //! stores how many boundary scvfs are embedded in this dual grid node - std::size_t numBoundaryScvfs_; + GridIndexContainer scvIndices_; //!< The indices of the scvs around a dual grid node + std::vector<LocalIndexContainer> localScvfIndicesInScv_; //!< Maps to each scv a list of scvf indices embedded in it + + GridIndexContainer scvfIndices_; //!< the indices of the scvfs around a dual grid node + std::vector<bool> scvfIsOnBoundary_; //!< Maps to each scvf a boolean to indicate if it is on the boundary + std::vector<GridIndexContainer> scvfNeighborScvIndices_; //!< The neighboring scvs for the scvfs (order: 0 - inside, 1..n - outside) + std::size_t numBoundaryScvfs_; //!< stores how many boundary scvfs are embedded in this dual grid node }; /*! - * \ingroup Mpfa + * \ingroup CCMpfaDiscretization * \brief Class for the index sets of the dual grid in mpfa schemes. + * + * \tparam GridIndexType The type used for indices on the grid + * \tparam LocalIndexType The type used for indexing in interaction volumes + * \tparam dim The dimension of the grid */ -template<class TypeTag> +template< class GridIndexType, class LocalIndexType, int dim > class CCMpfaDualGridIndexSet { - using ParentType = std::vector<DualGridNodalIndexSet<TypeTag>>; - - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - - static const int dim = GridView::dimension; - public: - using NodalIndexSet = DualGridNodalIndexSet<TypeTag>; + using NodalIndexSet = DualGridNodalIndexSet< GridIndexType, LocalIndexType, dim >; - //! default constructor + //! Default constructor should not be used CCMpfaDualGridIndexSet() = delete; - //! constructor - explicit CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(dim)) {} + //! Constructor taking a grid view + template< class GridView > + CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(dim)) {} - //! access with an scvf + //! Access with an scvf + template< class SubControlVolumeFace > const NodalIndexSet& operator[] (const SubControlVolumeFace& scvf) const { return nodalIndexSets_[scvf.vertexIndex()]; } + + template< class SubControlVolumeFace > NodalIndexSet& operator[] (const SubControlVolumeFace& scvf) { return nodalIndexSets_[scvf.vertexIndex()]; } - //! access with an index - const NodalIndexSet& operator[] (IndexType i) const - { return nodalIndexSets_[i]; } - NodalIndexSet& operator[] (IndexType i) - { return nodalIndexSets_[i]; } + //! Access with an index + const NodalIndexSet& operator[] (GridIndexType i) const { return nodalIndexSets_[i]; } + NodalIndexSet& operator[] (GridIndexType i) { return nodalIndexSets_[i]; } private: std::vector<NodalIndexSet> nodalIndexSets_; }; -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh index 7e9281db78d17988fd21656319d4f718e491b07e..313cf34659b7ecd614d60f39fd8e517a13672483 100644 --- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh +++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh @@ -18,14 +18,17 @@ *****************************************************************************/ /*! * \file - * \brief The local object of flux var caches + * \ingroup CCMpfaDiscretization + * \brief The element-local object of flux variables caches */ #ifndef DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_FLUXVARSCACHE_HH #define DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_FLUXVARSCACHE_HH #include <dune/common/exceptions.hh> + #include <dumux/common/properties.hh> #include <dumux/common/parameters.hh> + #include "fluxvariablescachefiller.hh" #include "methods.hh" @@ -33,56 +36,61 @@ namespace Dumux { /*! - * \ingroup ImplicitModel - * \brief Base class for the local flux variables cache. - * Prepares the cache on all the faces in the stencil. + * \ingroup CCMpfaDiscretization + * \brief The flux variables caches for an element + * \note The class is specialized for a version with and without caching. + * If grid caching is enabled the flux caches are stored for the whole + * gridview in the corresponding GridFluxVariablesCache. This is memory + * intensive but faster. For caching disabled, the flux caches are locally + * computed for each element whenever needed. */ template<class TypeTag, bool EnableGridFluxVariablesCache> class CCMpfaElementFluxVariablesCache; /*! - * \ingroup ImplicitModel - * \brief Spezialization when caching globally + * \ingroup CCMpfaDiscretization + * \brief The flux variables caches for an element with caching enabled */ template<class TypeTag> class CCMpfaElementFluxVariablesCache<TypeTag, true> { - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); using GridFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache); - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + public: + //! The constructor CCMpfaElementFluxVariablesCache(const GridFluxVariablesCache& global) : gridFluxVarsCachePtr_(&global) {} - // Specialization for the global caching being enabled - do nothing here + //! Specialization for the global caching being enabled - do nothing here void bindElement(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) {} - // Specialization for the global caching being enabled - do nothing here + //! Specialization for the global caching being enabled - do nothing here void bind(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) {} - // Specialization for the global caching being enabled - do nothing here + //! Specialization for the global caching being enabled - do nothing here void bindScvf(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf) {} - // Specialization for the global caching being enabled - do nothing here + //! Specialization for the global caching being enabled - do nothing here void update(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) {} - // access operators in the case of caching + //! access operators in the case of caching const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const { return (*gridFluxVarsCachePtr_)[scvf]; } @@ -95,54 +103,55 @@ private: }; /*! - * \ingroup ImplicitModel - * \brief Spezialization when not using global caching + * \ingroup CCMpfaDiscretization + * \brief The flux variables caches for an element with caching disabled */ template<class TypeTag> class CCMpfaElementFluxVariablesCache<TypeTag, false> { - // the flux variables cache filler needs to be friend to fill - // the interaction volumes and data handles - friend CCMpfaFluxVariablesCacheFiller<TypeTag>; - - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; + using GridIndexType = typename GridView::IndexSet::IndexType; using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); using GridFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using FluxVariablesCacheFiller = CCMpfaFluxVariablesCacheFiller<TypeTag>; - using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); + using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume); - using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; - - static constexpr int dim = GridView::dimension; + using SecondaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; public: CCMpfaElementFluxVariablesCache(const GridFluxVariablesCache& global) : gridFluxVarsCachePtr_(&global) {} - // This function has to be called prior to flux calculations on the element. - // Prepares the transmissibilities of the scv faces in an element. The FvGeometry is assumed to be bound. + /*! + * \brief Prepares the transmissibilities of the scv faces in an element + * \note the fvGeometry is assumed to be bound to the same element + * \note this function has to be called prior to flux calculations on the element. + */ void bindElement(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) { - // TODO + // For mpfa schemes we will have to prepare the caches of all scvfs that are + // embedded in the interaction volumes in which the element-local scvfs are embedded DUNE_THROW(Dune::NotImplemented, "Local element binding of the flux variables cache in mpfa schemes"); } - // This function is called by the CCLocalResidual before flux calculations during assembly. - // Prepares the transmissibilities of the scv faces in the stencil. The FvGeometries are assumed to be bound. + /*! + * \brief Prepares the transmissibilities of the scv faces in the stencil of an element + * \note the fvGeometry is assumed to be bound to the same element + * \note this function has to be called prior to flux calculations on the element. + */ void bind(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) { - // clear data clear_(); // some references for convenience @@ -202,22 +211,31 @@ public: } } + /*! + * \brief Prepares the transmissibilities of a single scv face + * \note the fvGeometry is assumed to be bound to the same element + * \note this function has to be called prior to flux calculations on this scvf. + */ void bindScvf(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf) { - // TODO - DUNE_THROW(Dune::NotImplemented, "Local element binding of the flux variables cache in mpfa schemes"); + // For mpfa schemes we will have to prepare the caches of all + // scvfs that are embedded in the interaction volumes this scvf is embedded + DUNE_THROW(Dune::NotImplemented, "Scvf-local binding of the flux variables cache in mpfa schemes"); } - // This function is used to update the transmissibilities if the volume variables have changed - // Results in undefined behaviour if called before bind() or with a different element + /*! + * \brief Update the transmissibilities if the volume variables have changed + * \note Results in undefined behaviour if called before bind() or with a different element + */ void update(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) { - // update only if transmissibilities are solution-dependent + // Update only if the filler puts + // solution-dependent stuff into the caches if (FluxVariablesCacheFiller::isSolDependent) { const auto& problem = gridFluxVarsCache().problem(); @@ -253,19 +271,54 @@ public: } } - // access operators in the case of no caching + //! access operators in the case of no caching const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const { return fluxVarsCache_[getLocalScvfIdx_(scvf.index())]; } - const FluxVariablesCache& operator [](const IndexType scvfIdx) const + //! access operators in the case of no caching + const FluxVariablesCache& operator [](const GridIndexType scvfIdx) const { return fluxVarsCache_[getLocalScvfIdx_(scvfIdx)]; } + //! access operators in the case of no caching FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) { return fluxVarsCache_[getLocalScvfIdx_(scvf.index())]; } - FluxVariablesCache& operator [](const IndexType scvfIdx) + //! access operators in the case of no caching + FluxVariablesCache& operator [](const GridIndexType scvfIdx) { return fluxVarsCache_[getLocalScvfIdx_(scvfIdx)]; } + //! access to the stored interaction volumes + const std::vector<PrimaryInteractionVolume>& primaryInteractionVolumes() const + { return primaryInteractionVolumes_; } + + //! access to the stored interaction volumes + std::vector<PrimaryInteractionVolume>& primaryInteractionVolumes() + { return primaryInteractionVolumes_; } + + //! access to the stored data handles + const std::vector<PrimaryIvDataHandle>& primaryDataHandles() const + { return primaryIvDataHandles_; } + + //! access to the stored data handles + std::vector<PrimaryIvDataHandle>& primaryDataHandles() + { return primaryIvDataHandles_; } + + //! access to the stored interaction volumes + const std::vector<SecondaryInteractionVolume>& secondaryInteractionVolumes() const + { return secondaryInteractionVolumes_; } + + //! access to the stored interaction volumes + std::vector<SecondaryInteractionVolume>& secondaryInteractionVolumes() + { return secondaryInteractionVolumes_; } + + //! access to the stored data handles + const std::vector<SecondaryIvDataHandle>& secondaryDataHandles() const + { return secondaryIvDataHandles_; } + + //! access to the stored data handles + std::vector<SecondaryIvDataHandle>& secondaryDataHandles() + { return secondaryIvDataHandles_; } + //! The global object we are a restriction of const GridFluxVariablesCache& gridFluxVarsCache() const { return *gridFluxVarsCachePtr_; } @@ -273,6 +326,7 @@ public: private: const GridFluxVariablesCache* gridFluxVarsCachePtr_; + //! clears all containers void clear_() { fluxVarsCache_.clear(); @@ -283,35 +337,27 @@ private: secondaryIvDataHandles_.clear(); } - // get number of interaction volumes that are going to be required + //! get estimate of interaction volumes that are going to be required template<class AssemblyMap> std::size_t getNoInteractionVolumesEstimate_(const Element& element, const AssemblyMap& assemblyMap) { - // TODO: uncomment as soon as methods are re-implemented // if statements are optimized away by the compiler - if (GET_PROP_VALUE(TypeTag, MpfaMethod) == MpfaMethods::oMethod /*|| method == MpfaMethods::oMethodFps*/) + if (GET_PROP_VALUE(TypeTag, MpfaMethod) == MpfaMethods::oMethod) { - //! Reserve memory for the case of each facet having neighbors being 4 levels higher. Memory limitations - //! do not play an important role here as global caching is disabled. In the unlikely case you want - //! to use higher local differences in element levels set a higher value for the parameter below - //! in your input file (note that some grids might only support levels differences of one anyway) + // Reserve memory for the case of each facet having neighbors being 4 levels higher. Memory limitations + // do not play an important role here as global caching is disabled. In the unlikely case you want + // to use higher local differences in element levels set a higher value for the parameter below + // in your input file (note that some grids might only support levels differences of one anyway) static const unsigned int maxDiff = getParamFromGroup<unsigned int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Grid.MaxLocalElementLevelDifference", 4); - return element.subEntities(dim)*maxDiff; + return element.subEntities(GridView::dimension)*maxDiff; } - /* else if (method == MpfaMethods::lMethod) - { - std::size_t numInsideScvfs = MpfaHelper::getNumLocalScvfs(element.geometry().type()); - std::size_t numOutsideScvf = 0; - for (const auto& dataJ : assemblyMap) numOutsideScvf += dataJ.scvfsJ.size(); - return numOutsideScvf - numInsideScvfs; - }*/ else DUNE_THROW(Dune::NotImplemented, "number of interaction volumes estimate for chosen mpfa scheme"); } - // get index of an scvf in the local container + //! get index of an scvf in the local container unsigned int getLocalScvfIdx_(const int scvfIdx) const { auto it = std::find(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx); @@ -320,15 +366,15 @@ private: } - // the local flux vars caches and the index set + // the local flux vars caches and corresponding indices std::vector<FluxVariablesCache> fluxVarsCache_; - std::vector<IndexType> globalScvfIndices_; + std::vector<GridIndexType> globalScvfIndices_; - // store the interaction volumes and handles + // stored interaction volumes and handles std::vector<PrimaryInteractionVolume> primaryInteractionVolumes_; std::vector<SecondaryInteractionVolume> secondaryInteractionVolumes_; - std::vector<DataHandle> primaryIvDataHandles_; - std::vector<DataHandle> secondaryIvDataHandles_; + std::vector<PrimaryIvDataHandle> primaryIvDataHandles_; + std::vector<SecondaryIvDataHandle> secondaryIvDataHandles_; }; } // end namespace diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index c1098567c841bec79bad303b4d5ce3965f2c1e9d..b3c905c974f5970a19259b5ebf1cac807f824154 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -18,6 +18,7 @@ *****************************************************************************/ /*! * \file + * \ingroup CCMpfaDiscretization * \brief The local (stencil) volume variables class for cell centered mpfa models */ #ifndef DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH @@ -29,52 +30,56 @@ namespace Dumux { /*! - * \ingroup ImplicitModel - * \brief Base class for the local volume variables vector + * \ingroup CCMpfaDiscretization + * \brief The local (stencil) volume variables class for cell centered mpfa models + * \note The class is specilized for versions with and without caching */ template<class TypeTag, bool enableGridVolVarsCache> -class CCMpfaElementVolumeVariables -{}; +class CCMpfaElementVolumeVariables; -// specialization in case of storing the volume variables globally +/*! + * \ingroup CCMpfaDiscretization + * \brief The local (stencil) volume variables class for cell centered mpfa models with caching + * \note the volume variables are stored for the whole grid view in the corresponding GridVolumeVariables class + */ template<class TypeTag> class CCMpfaElementVolumeVariables<TypeTag, /*enableGridVolVarsCache*/true> { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables); using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); - using IndexType = typename GridView::IndexSet::IndexType; - - static const int dim = GridView::dimension; - using Element = typename GridView::template Codim<0>::Entity; + using GridIndexType = typename GridView::IndexSet::IndexType; public: //! Constructor CCMpfaElementVolumeVariables(const GridVolumeVariables& gridVolVars) : gridVolVarsPtr_(&gridVolVars) {} + //! operator for the access with an scv const VolumeVariables& operator [](const SubControlVolume& scv) const { return gridVolVars().volVars(scv.dofIndex()); } - // operator for the access with an index - // needed for cc methods for the access to the boundary volume variables - const VolumeVariables& operator [](const IndexType scvIdx) const + //! operator for the access with an index + const VolumeVariables& operator [](const GridIndexType scvIdx) const { return gridVolVars().volVars(scvIdx); } - // For compatibility reasons with the case of not storing the vol vars. - // function to be called before assembling an element, preparing the vol vars within the stencil + //! precompute all volume variables in a stencil of an element - do nothing volVars: are cached void bind(const Element& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) {} + const SolutionVector& sol) + {} - // function to prepare the vol vars within the element + //! precompute the volume variables of an element - do nothing: volVars are cached void bindElement(const Element& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) {} + const SolutionVector& sol) + {} //! The global volume variables object we are a restriction of const GridVolumeVariables& gridVolVars() const @@ -85,36 +90,38 @@ private: }; -// Specialization when the current volume variables are not stored +/*! + * \ingroup CCMpfaDiscretization + * \brief The local (stencil) volume variables class for cell centered tpfa models with caching + */ template<class TypeTag> class CCMpfaElementVolumeVariables<TypeTag, /*enableGridVolVarsCache*/false> { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); + using Element = typename GridView::template Codim<0>::Entity; + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables); using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using IndexType = typename GridView::IndexSet::IndexType; + using GridIndexType = typename GridView::IndexSet::IndexType; - static const int dim = GridView::dimension; - using Element = typename GridView::template Codim<0>::Entity; public: - //! Constructor CCMpfaElementVolumeVariables(const GridVolumeVariables& gridVolVars) : gridVolVarsPtr_(&gridVolVars) {} - // Binding of an element, prepares the volume variables within the element stencil - // called by the local jacobian to prepare element assembly + //! Prepares the volume variables within the element stencil void bind(const Element& element, const FVElementGeometry& fvGeometry, const SolutionVector& sol) { + clear(); + const auto& problem = gridVolVars().problem(); const auto& fvGridGeometry = fvGeometry.fvGridGeometry(); const auto& gridIvIndexSets = fvGridGeometry.gridInteractionVolumeIndexSets(); @@ -151,7 +158,7 @@ public: ++localIdx; } - // eventually prepare boundary volume variables + // maybe prepare boundary volume variables const auto maxNumBoundaryVolVars = maxNumBoundaryVolVars_(fvGeometry); if (maxNumBoundaryVolVars > 0) { @@ -167,7 +174,8 @@ public: const auto bcTypes = problem.boundaryTypes(element, scvf); - // on dirichlet boundaries use dirichlet values + // Only proceed on dirichlet boundaries. Fluxes across Neumann + // boundaries are never computed - the user-defined flux is taken. if (bcTypes.hasOnlyDirichlet()) { // boundary volume variables @@ -180,12 +188,6 @@ public: volumeVariables_.emplace_back(std::move(dirichletVolVars)); volVarIndices_.push_back(scvf.outsideScvIdx()); } - // use the inside volume variables for neumann boundaries - else - { - volumeVariables_.emplace_back(volumeVariables_[0]); - volVarIndices_.push_back(scvf.outsideScvIdx()); - } } // Update boundary volume variables in the neighbors @@ -194,16 +196,13 @@ public: if (fvGridGeometry.vertexUsesPrimaryInteractionVolume(scvf.vertexIndex())) { const auto& nodalIndexSet = gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet(); - // if present, insert boundary vol vars if (nodalIndexSet.numBoundaryScvfs() > 0) addBoundaryVolVars_(problem, fvGeometry, nodalIndexSet); - } else { const auto& nodalIndexSet = gridIvIndexSets.secondaryIndexSet(scvf).nodalIndexSet(); - // if present, insert boundary vol vars if (nodalIndexSet.numBoundaryScvfs() > 0) addBoundaryVolVars_(problem, fvGeometry, nodalIndexSet); @@ -235,11 +234,13 @@ public: // } } - // Binding of an element, prepares only the volume variables of the element + //! Prepares the volume variables of an element void bindElement(const Element& element, const FVElementGeometry& fvGeometry, const SolutionVector& sol) { + clear(); + auto eIdx = fvGeometry.fvGridGeometry().elementMapper().index(element); volumeVariables_.resize(1); volVarIndices_.resize(1); @@ -253,29 +254,40 @@ public: volVarIndices_[0] = scv.dofIndex(); } + //! access operator with scv const VolumeVariables& operator [](const SubControlVolume& scv) const { return volumeVariables_[getLocalIdx_(scv.dofIndex())]; } + //! access operator with scv VolumeVariables& operator [](const SubControlVolume& scv) { return volumeVariables_[getLocalIdx_(scv.dofIndex())]; } - const VolumeVariables& operator [](IndexType scvIdx) const + //! access operator with scv index + const VolumeVariables& operator [](GridIndexType scvIdx) const { return volumeVariables_[getLocalIdx_(scvIdx)]; } - VolumeVariables& operator [](IndexType scvIdx) + //! access operator with scv index + VolumeVariables& operator [](GridIndexType scvIdx) { return volumeVariables_[getLocalIdx_(scvIdx)]; } //! The global volume variables object we are a restriction of const GridVolumeVariables& gridVolVars() const { return *gridVolVarsPtr_; } + //! Clear all local storage + void clear() + { + volVarIndices_.clear(); + volumeVariables_.clear(); + } + private: const GridVolumeVariables* gridVolVarsPtr_; - //! Computes how many boundary vol vars come into play for flux calculations - //! on this element. This number is probably almost always higher than the actually - //! needed number of volume variables. However, memory is not an issue for the global - //! caching being deactivated and we want to make sure we reserve enough memory here. + // Computes how many boundary vol vars come into play for flux calculations + // on this element. This number here is probably always higher than the actually + // needed number of volume variables. However, memory is not an issue for the global + // caching being deactivated and we want to make sure we reserve enough memory here. // TODO: What about non-symmetric schemes? Is there a better way for estimating this? std::size_t maxNumBoundaryVolVars_(const FVElementGeometry& fvGeometry) { @@ -311,10 +323,10 @@ private: const auto insideElement = fvGeometry.fvGridGeometry().element(insideScvIdx); const auto bcTypes = problem.boundaryTypes(insideElement, ivScvf); - // on dirichlet boundaries use dirichlet values + // Only proceed on dirichlet boundaries. Fluxes across Neumann + // boundaries are never computed - the user-defined flux is taken. if (bcTypes.hasOnlyDirichlet()) { - // boundary volume variables VolumeVariables dirichletVolVars; const auto& ivScv = fvGeometry.scv(insideScvIdx); dirichletVolVars.update(ElementSolution(problem.dirichlet(insideElement, ivScvf)), @@ -325,15 +337,10 @@ private: volumeVariables_.emplace_back(std::move(dirichletVolVars)); volVarIndices_.push_back(ivScvf.outsideScvIdx()); } - // use the inside volume variables for neumann boundaries - else - { - volumeVariables_.emplace_back((*this)[insideScvIdx]); - volVarIndices_.push_back(ivScvf.outsideScvIdx()); - } } } + //! map a global scv index to the local storage index int getLocalIdx_(const int volVarIdx) const { auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx); @@ -341,7 +348,7 @@ private: return std::distance(volVarIndices_.begin(), it); } - std::vector<IndexType> volVarIndices_; + std::vector<GridIndexType> volVarIndices_; std::vector<VolumeVariables> volumeVariables_; }; diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh index e94332d2ec47e2eddbaed466444740628bf519a4..956b174e1637e540a98562996acff5e254e63c50 100644 --- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh @@ -18,26 +18,26 @@ *****************************************************************************/ /*! * \file - * \brief This file contains the class which is required to calculate - * molar and mass fluxes of a component in a fluid phase over a face of a finite volume by means - * of Fick's Law for cell-centered MPFA models. + * \ingroup CCMpfaDiscretization + * \brief Fick's law for cell-centered finite volume schemes with multi-point flux approximation */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_FICKS_LAW_HH #define DUMUX_DISCRETIZATION_CC_MPFA_FICKS_LAW_HH #include <dumux/common/math.hh> +#include <dumux/common/parameters.hh> #include <dumux/common/properties.hh> #include <dumux/discretization/methods.hh> -namespace Dumux { - -// forward declaration +namespace Dumux +{ +//! forward declaration of the method-specific implemetation template<class TypeTag, DiscretizationMethods discMethod> class FicksLawImplementation; /*! - * \ingroup Mpfa - * \brief Specialization of Fick's Law for the CCMpfa method. + * \ingroup CCMpfaDiscretization + * \brief Fick's law for cell-centered finite volume schemes with multi-point flux approximation */ template <class TypeTag> class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> @@ -46,6 +46,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); @@ -54,14 +55,6 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); using BalanceEqOpts = typename GET_PROP_TYPE(TypeTag, BalanceEqOpts); - // Always use the dynamic type for vectors (compatibility with the boundary) - using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using CoefficientVector = typename PrimaryInteractionVolume::Traits::DynamicVector; - using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; - - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); static constexpr int numComponents = GET_PROP_VALUE(TypeTag,NumComponents); using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>; @@ -81,78 +74,113 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> const SubControlVolumeFace& scvf, const FluxVariablesCacheFiller& fluxVarsCacheFiller) { - // get interaction volume from the flux vars cache filler & upate the cache - if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) - scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.secondaryInteractionVolume(), - fluxVarsCacheFiller.dataHandle(), - scvf, phaseIdx, compIdx); - else - scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.primaryInteractionVolume(), - fluxVarsCacheFiller.dataHandle(), - scvf, phaseIdx, compIdx); + // get interaction volume related data from the filler class & upate the cache + if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) + scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.secondaryInteractionVolume(), + fluxVarsCacheFiller.secondaryIvLocalFaceData(), + fluxVarsCacheFiller.secondaryIvDataHandle(), + scvf, phaseIdx, compIdx); + else + scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.primaryInteractionVolume(), + fluxVarsCacheFiller.primaryIvLocalFaceData(), + fluxVarsCacheFiller.primaryIvDataHandle(), + scvf, phaseIdx, compIdx); } }; //! The cache used in conjunction with the mpfa Fick's Law class MpfaFicksLawCache { - // We always use the dynamic types here to be compatible on the boundary - using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer; - using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer; + // In the current implementation of the flux variables cache we cannot + // make a disctinction between dynamic (mpfa-o) and static (mpfa-l) + // matrix and vector types, as currently the cache class can only be templated + // by a type tag (and there can only be one). We use a dynamic vector here to + // make sure it works in case one of the two used interaction volume types uses + // dynamic types performance is thus lowered for schemes using static types. + // TODO: this has to be checked thoroughly as soon as a scheme using static types + // is implemented. One idea to overcome the performance drop could be only + // storing the iv-local index here and obtain tij always from the datahandle + // of the fluxVarsCacheContainer + using Vector = Dune::DynamicVector< Scalar >; + using Matrix = Dune::DynamicMatrix< Scalar >; + + using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); + using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector; + using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix; + + static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value, + "The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" ); + static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value, + "The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" ); + + using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume); + using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector; + using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix; + + static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value, + "The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" ); + static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value, + "The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" ); + + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); public: // export filler type using Filler = MpfaFicksLawCacheFiller; - // update cached objects for the diffusive fluxes - template<class InteractionVolume> + /*! + * \brief Update cached objects (transmissibilities) + * + * \tparam InteractionVolume The (mpfa scheme-specific) interaction volume + * \tparam LocalFaceData The object used to store iv-local info on an scvf + * \tparam DataHandle The object used to store transmissibility matrices etc. + * + * \param iv The interaction volume this scvf is embedded in + * \param localFaceData iv-local info on this scvf + * \param dataHandle Transmissibility matrix & gravity data of this iv + * \param scvf The sub-control volume face + */ + template<class InteractionVolume, class LocalFaceData, class DataHandle> void updateDiffusion(const InteractionVolume& iv, + const LocalFaceData& localFaceData, const DataHandle& dataHandle, const SubControlVolumeFace &scvf, unsigned int phaseIdx, unsigned int compIdx) { - const auto& localFaceData = iv.getLocalFaceData(scvf); + switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside(); - // update the quantities that are equal for all phases - diffusionSwitchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside(); - diffusionVolVarsStencil_[phaseIdx][compIdx] = &dataHandle.volVarsStencil(); - diffusionDirichletData_[phaseIdx][compIdx] = &dataHandle.dirichletData(); + //! store pointer to the mole fraction vector of this iv + xj_[phaseIdx][compIdx] = &dataHandle.moleFractions(phaseIdx, compIdx); - // the transmissibilities on surface grids have to be obtained from the outside + const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); if (dim == dimWorld) - diffusionTij_[phaseIdx][compIdx] = &dataHandle.T()[localFaceData.ivLocalScvfIndex()]; + Tij_[phaseIdx][compIdx] = &dataHandle.diffusionT()[ivLocalIdx]; else - diffusionTij_[phaseIdx][compIdx] = localFaceData.isOutside() ? - &dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()] : - &dataHandle.T()[localFaceData.ivLocalScvfIndex()]; + Tij_[phaseIdx][compIdx] = localFaceData.isOutside() ? &dataHandle.diffusionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()] + : &dataHandle.diffusionT()[ivLocalIdx]; } - //! Returns the volume variables indices necessary for diffusive flux - //! computation. This includes all participating boundary volume variables - //! and it can be different for the phases & components. - const Stencil& diffusionVolVarsStencil(unsigned int phaseIdx, unsigned int compIdx) const - { return *diffusionVolVarsStencil_[phaseIdx][compIdx]; } - - //! On faces that are "outside" w.r.t. a face in the interaction volume, - //! we have to take the negative value of the fluxes, i.e. multiply by -1.0 + //! In the interaction volume-local system of eq we have one unknown per face. + //! On scvfs on this face, but in "outside" (neighbor) elements of it, we have + //! to take the negative value of the fluxes due to the flipped normal vector. + //! This function returns whether or not this scvf is an "outside" face in the iv. bool diffusionSwitchFluxSign(unsigned int phaseIdx, unsigned int compIdx) const - { return diffusionSwitchFluxSign_[phaseIdx][compIdx]; } + { return switchFluxSign_[phaseIdx][compIdx]; } - //! Returns the transmissibilities associated with the volume variables - //! This can be different for the phases & components. - const CoefficientVector& diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const - { return *diffusionTij_[phaseIdx][compIdx]; } + //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions + const Vector& diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const + { return *Tij_[phaseIdx][compIdx]; } - //! Returns the data on dirichlet boundary conditions affecting - //! the flux computation on this face - const DirichletDataContainer& diffusionDirichletData(unsigned int phaseIdx, unsigned int compIdx) const - { return *diffusionDirichletData_[phaseIdx][compIdx]; } + //! The cell (& Dirichlet) mole fractions within this interaction volume + const Vector& moleFractions(unsigned int phaseIdx, unsigned int compIdx) const + { return *xj_[phaseIdx][compIdx]; } private: - std::array< std::array<bool, numComponents>, numPhases> diffusionSwitchFluxSign_; - std::array< std::array<const Stencil*, numComponents>, numPhases> diffusionVolVarsStencil_; - std::array< std::array<const DirichletDataContainer*, numComponents>, numPhases> diffusionDirichletData_; - std::array< std::array<const CoefficientVector*, numComponents>, numPhases> diffusionTij_; + std::array< std::array<bool, numComponents>, numPhases > switchFluxSign_; + std::array< std::array<const Vector*, numComponents>, numPhases > Tij_; //!< The transmissibilities such that f = Tij*xj + std::array< std::array<const Vector*, numComponents>, numPhases > xj_; //!< The interaction-volume wide mole fractions xj }; public: @@ -162,13 +190,14 @@ public: // state the type for the corresponding cache and its filler using Cache = MpfaFicksLawCache; - static ComponentFluxVector flux (const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf, - const int phaseIdx, - const ElementFluxVariablesCache& elemFluxVarsCache) + //! Compute the diffusive flux across an scvf + static ComponentFluxVector flux(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf, + const int phaseIdx, + const ElementFluxVariablesCache& elemFluxVarsCache) { ComponentFluxVector componentFlux(0.0); for (int compIdx = 0; compIdx < numComponents; compIdx++) @@ -187,20 +216,16 @@ public: const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx); // prepare computations - Scalar flux(0.0); - unsigned int i = 0; const auto& fluxVarsCache = elemFluxVarsCache[scvf]; const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx); + const auto& xj = fluxVarsCache.moleFractions(phaseIdx, compIdx); // calculate Tij*xj - for (const auto volVarIdx : fluxVarsCache.diffusionVolVarsStencil(phaseIdx, compIdx)) - flux += tij[i++]*elemVolVars[volVarIdx].moleFraction(phaseIdx, compIdx); - - // add contributions from dirichlet BCs - for (const auto& d : fluxVarsCache.diffusionDirichletData(phaseIdx, compIdx)) - flux += tij[i++]*elemVolVars[d.volVarIndex()].moleFraction(phaseIdx, compIdx); + Scalar flux = tij*xj; + if (fluxVarsCache.diffusionSwitchFluxSign(phaseIdx, compIdx)) + flux *= -1.0; - componentFlux[compIdx] = fluxVarsCache.diffusionSwitchFluxSign(phaseIdx, compIdx) ? -1.0*flux*rho*effFactor : flux*rho*effFactor; + componentFlux[compIdx] = flux*rho*effFactor; } // accumulate the phase component flux @@ -212,6 +237,7 @@ public: } private: + //! compute the density at branching facets for network grids as arithmetic mean static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf, const unsigned int phaseIdx) diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh index 29935706928e45b69c71899ad44039ec5846f9ce..e62d4d03d4def3b00b36585041d5cbde2e670318 100644 --- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh +++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh @@ -16,37 +16,49 @@ * You should have received a copy of the GNU General Public License * * along with this program. If not, see <http://www.gnu.org/licenses/>. * *****************************************************************************/ -/*! + /*! * \file - * \brief The flux variables cache filler class + * \ingroup CCMpfaDiscretization + * \brief A helper class to fill the flux variable caches used in the flux constitutive laws */ #ifndef DUMUX_DISCRETIZATION_CCMPFA_FLUXVARSCACHE_FILLER_HH #define DUMUX_DISCRETIZATION_CCMPFA_FLUXVARSCACHE_FILLER_HH +#include <dumux/common/properties.hh> + #include <dumux/discretization/methods.hh> #include <dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh> +#include <dumux/discretization/cellcentered/mpfa/localassembler.hh> namespace Dumux { /*! - * \ingroup ImplicitModel - * \brief Helper class to fill the flux var caches + * \ingroup CCMpfaDiscretization + * \brief Helper class to fill the flux variables caches within + * the interaction volume around a given sub-control volume face. */ template<class TypeTag> class CCMpfaFluxVariablesCacheFiller { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); + using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); + using PrimaryDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; + using PrimaryLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData; using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume); - using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; + using SecondaryDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle; + using SecondaryLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData; - using Element = typename GridView::template Codim<0>::Entity; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; static constexpr bool doAdvection = GET_PROP_VALUE(TypeTag, EnableAdvection); static constexpr bool doDiffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion); @@ -57,11 +69,12 @@ class CCMpfaFluxVariablesCacheFiller static constexpr bool soldependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction); public: - static constexpr bool isSolDependent = (doAdvection && soldependentAdvection) || - (doDiffusion && soldependentDiffusion) || - (doHeatConduction && soldependentHeatConduction); + //! This cache filler is always solution-dependent, as it updates the + //! vectors of cell unknowns with which the transmissibilities have to be + //! multiplied in order to obtain the fluxes. + static constexpr bool isSolDependent = true; - //! The constructor. Sets the problem pointer + //! The constructor. Sets problem pointer. CCMpfaFluxVariablesCacheFiller(const Problem& problem) : problemPtr_(&problem) {} /*! @@ -71,7 +84,7 @@ public: * \param scvfFluxVarsCache The flux var cache to be updated corresponding to the given scvf * \param element The finite element * \param fvGeometry The finite volume geometry - * \param elemVolVars The element volume variables + * \param elemVolVars The element volume variables (primary/secondary variables) * \param scvf The corresponding sub-control volume face * \param forceUpdateAll if true, forces all caches to be updated (even the solution-independent ones) */ @@ -96,30 +109,30 @@ public: if (forceUpdateAll) { // the local index of the interaction volume to be created in its container - const auto ivIndexInContainer = fluxVarsCacheContainer.secondaryInteractionVolumes_.size(); + const auto ivIndexInContainer = fluxVarsCacheContainer.secondaryInteractionVolumes().size(); // prepare the locally cached boundary interaction volume const auto& indexSet = fvGridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf); - fluxVarsCacheContainer.secondaryInteractionVolumes_.emplace_back(); - secondaryIv_ = &fluxVarsCacheContainer.secondaryInteractionVolumes_.back(); + fluxVarsCacheContainer.secondaryInteractionVolumes().emplace_back(); + secondaryIv_ = &fluxVarsCacheContainer.secondaryInteractionVolumes().back(); secondaryIv_->setUpLocalScope(indexSet, problem(), fvGeometry); // prepare the corresponding data handle - fluxVarsCacheContainer.secondaryIvDataHandles_.emplace_back(); - ivDataHandle_ = &fluxVarsCacheContainer.secondaryIvDataHandles_.back(); - secondaryIv_->prepareDataHandle(*ivDataHandle_); + fluxVarsCacheContainer.secondaryDataHandles().emplace_back(); + secondaryIvDataHandle_ = &fluxVarsCacheContainer.secondaryDataHandles().back(); + secondaryIvDataHandle_->resize(*secondaryIv_); // fill the caches for all the scvfs in the interaction volume - fillCachesInInteractionVolume_(fluxVarsCacheContainer, *secondaryIv_, *ivDataHandle_, ivIndexInContainer, true); + fillCachesInInteractionVolume_(fluxVarsCacheContainer, *secondaryIv_, *secondaryIvDataHandle_, ivIndexInContainer, true); } else { const auto ivIndexInContainer = scvfFluxVarsCache.ivIndexInContainer(); - secondaryIv_ = &fluxVarsCacheContainer.secondaryInteractionVolumes_[ivIndexInContainer]; - ivDataHandle_ = &fluxVarsCacheContainer.secondaryIvDataHandles_[ivIndexInContainer]; + secondaryIv_ = &fluxVarsCacheContainer.secondaryInteractionVolumes()[ivIndexInContainer]; + secondaryIvDataHandle_ = &fluxVarsCacheContainer.secondaryDataHandles()[ivIndexInContainer]; // fill the caches for all the scvfs in the interaction volume - fillCachesInInteractionVolume_(fluxVarsCacheContainer, *secondaryIv_, *ivDataHandle_, ivIndexInContainer); + fillCachesInInteractionVolume_(fluxVarsCacheContainer, *secondaryIv_, *secondaryIvDataHandle_, ivIndexInContainer); } } else @@ -127,61 +140,69 @@ public: if (forceUpdateAll) { // the local index of the interaction volume to be created in its container - const auto ivIndexInContainer = fluxVarsCacheContainer.primaryInteractionVolumes_.size(); + const auto ivIndexInContainer = fluxVarsCacheContainer.primaryInteractionVolumes().size(); // prepare the locally cached boundary interaction volume const auto& indexSet = fvGridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf); - fluxVarsCacheContainer.primaryInteractionVolumes_.emplace_back(); - primaryIv_ = &fluxVarsCacheContainer.primaryInteractionVolumes_.back(); + fluxVarsCacheContainer.primaryInteractionVolumes().emplace_back(); + primaryIv_ = &fluxVarsCacheContainer.primaryInteractionVolumes().back(); primaryIv_->setUpLocalScope(indexSet, problem(), fvGeometry); // prepare the corresponding data handle - fluxVarsCacheContainer.primaryIvDataHandles_.emplace_back(); - ivDataHandle_ = &fluxVarsCacheContainer.primaryIvDataHandles_.back(); - primaryIv_->prepareDataHandle(*ivDataHandle_); + fluxVarsCacheContainer.primaryDataHandles().emplace_back(); + primaryIvDataHandle_ = &fluxVarsCacheContainer.primaryDataHandles().back(); + primaryIvDataHandle_->resize(*primaryIv_); // fill the caches for all the scvfs in the interaction volume - fillCachesInInteractionVolume_(fluxVarsCacheContainer, *primaryIv_, *ivDataHandle_, ivIndexInContainer, true); + fillCachesInInteractionVolume_(fluxVarsCacheContainer, *primaryIv_, *primaryIvDataHandle_, ivIndexInContainer, true); } else { const auto ivIndexInContainer = scvfFluxVarsCache.ivIndexInContainer(); - primaryIv_ = &fluxVarsCacheContainer.primaryInteractionVolumes_[ivIndexInContainer]; - ivDataHandle_ = &fluxVarsCacheContainer.primaryIvDataHandles_[ivIndexInContainer]; + primaryIv_ = &fluxVarsCacheContainer.primaryInteractionVolumes()[ivIndexInContainer]; + primaryIvDataHandle_ = &fluxVarsCacheContainer.primaryDataHandles()[ivIndexInContainer]; // fill the caches for all the scvfs in the interaction volume - fillCachesInInteractionVolume_(fluxVarsCacheContainer, *primaryIv_, *ivDataHandle_, ivIndexInContainer); + fillCachesInInteractionVolume_(fluxVarsCacheContainer, *primaryIv_, *primaryIvDataHandle_, ivIndexInContainer); } } } + //! returns the stored interaction volume pointer const PrimaryInteractionVolume& primaryInteractionVolume() const { return *primaryIv_; } + //! returns the stored interaction volume pointer const SecondaryInteractionVolume& secondaryInteractionVolume() const { return *secondaryIv_; } - const DataHandle& dataHandle() const - { return *ivDataHandle_; } + //! returns the stored data handle pointer + const PrimaryDataHandle& primaryIvDataHandle() const + { return *primaryIvDataHandle_; } -private: + //! returns the stored data handle pointer + const SecondaryDataHandle& secondaryIvDataHandle() const + { return *secondaryIvDataHandle_; } - const Problem& problem() const - { return *problemPtr_; } + //! returns the currently stored iv-local face data object + const PrimaryLocalFaceData& primaryIvLocalFaceData() const + { return *primaryLocalFaceData_; } - const Element& element() const - { return *elementPtr_; } + //! returns the currently stored iv-local face data object + const SecondaryLocalFaceData& secondaryIvLocalFaceData() const + { return *secondaryLocalFaceData_; } - const FVElementGeometry& fvGeometry() const - { return *fvGeometryPtr_; } +private: - const ElementVolumeVariables& elemVolVars() const - { return *elemVolVarsPtr_; } + const Problem& problem() const { return *problemPtr_; } + const Element& element() const { return *elementPtr_; } + const FVElementGeometry& fvGeometry() const { return *fvGeometryPtr_; } + const ElementVolumeVariables& elemVolVars() const { return *elemVolVarsPtr_; } //! Method to fill the flux var caches within an interaction volume - template<class FluxVariablesCacheContainer, class InteractionVolumeType> - void fillCachesInInteractionVolume_(FluxVariablesCacheContainer& fluxVarsCacheContainer, - InteractionVolumeType& iv, + template<class FluxVarsCacheContainer, class InteractionVolume, class DataHandle> + void fillCachesInInteractionVolume_(FluxVarsCacheContainer& fluxVarsCacheContainer, + InteractionVolume& iv, DataHandle& handle, unsigned int ivIndexInContainer, bool forceUpdateAll = false) @@ -205,31 +226,21 @@ private: i++; } - if (forceUpdateAll) - { - fillAdvection(fluxVarsCacheContainer, iv, handle, ivScvfs, ivFluxVarCaches); - fillDiffusion(fluxVarsCacheContainer, iv, handle, ivScvfs, ivFluxVarCaches); - fillHeatConduction(fluxVarsCacheContainer, iv, handle, ivScvfs, ivFluxVarCaches); - } - else - { - if (doAdvection && soldependentAdvection) - fillAdvection(fluxVarsCacheContainer, iv, handle, ivScvfs, ivFluxVarCaches); - if (doDiffusion && soldependentDiffusion) - fillDiffusion(fluxVarsCacheContainer, iv, handle, ivScvfs, ivFluxVarCaches); - if (doHeatConduction && soldependentHeatConduction) - fillHeatConduction(fluxVarsCacheContainer, iv, handle, ivScvfs, ivFluxVarCaches); - } + fillAdvection(iv, handle, ivScvfs, ivFluxVarCaches, forceUpdateAll); + fillDiffusion(iv, handle, ivScvfs, ivFluxVarCaches, forceUpdateAll); + fillHeatConduction(iv, handle, ivScvfs, ivFluxVarCaches, forceUpdateAll); } - //! method to fill the advective quantities - template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool advectionEnabled = doAdvection> - typename std::enable_if<advectionEnabled>::type - fillAdvection(FluxVariablesCacheContainer& fluxVarsCacheContainer, - InteractionVolumeType& iv, - DataHandle& handle, - const std::vector<const SubControlVolumeFace*>& ivScvfs, - const std::vector<FluxVariablesCache*>& ivFluxVarCaches) + //! fills the advective quantities (enabled advection) + template< class InteractionVolume, + class DataHandle, + bool enableAdvection = doAdvection, + typename std::enable_if_t<enableAdvection, int> = 0 > + void fillAdvection(InteractionVolume& iv, + DataHandle& handle, + const std::vector<const SubControlVolumeFace*>& ivScvfs, + const std::vector<FluxVariablesCache*>& ivFluxVarCaches, + bool forceUpdateAll = false) { using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); using AdvectionFiller = typename AdvectionType::Cache::Filler; @@ -237,15 +248,98 @@ private: static constexpr auto AdvectionMethod = AdvectionType::myDiscretizationMethod; using LambdaFactory = TensorLambdaFactory<TypeTag, AdvectionMethod>; - // set the advection context in the data handle - handle.setAdvectionContext(); - - // maybe solve the local system subject to K (if AdvectionType uses mpfa) + // skip the following if advection doesn't use mpfa if (AdvectionMethod == DiscretizationMethods::CCMpfa) - iv.solveLocalSystem(LambdaFactory::getAdvectionLambda(), problem(), fvGeometry(), elemVolVars(), handle); + { + // get instance of the interaction volume-local assembler + using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >; + IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars()); + + // Use different assembly if gravity is enabled + static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"); + + // Assemble T only if permeability is sol-dependent or if update is forced + if (forceUpdateAll || soldependentAdvection) + { + // distinguish between normal/surface grids (optimized away by compiler) + if (dim < dimWorld) + { + if (enableGravity) + localAssembler.assembleWithGravity( handle.advectionTout(), + handle.advectionT(), + handle.gravityOutside(), + handle.gravity(), + handle.advectionCA(), + handle.advectionA(), + iv, + LambdaFactory::getAdvectionLambda() ); + + else + localAssembler.assemble( handle.advectionTout(), + handle.advectionT(), + iv, + LambdaFactory::getAdvectionLambda() ); + } + + // normal grids + else + { + if (enableGravity) + localAssembler.assembleWithGravity( handle.advectionT(), + handle.gravity(), + handle.advectionCA(), + iv, + LambdaFactory::getAdvectionLambda() ); + else + localAssembler.assemble( handle.advectionT(), + iv, + LambdaFactory::getAdvectionLambda() ); + } + } + + // (maybe) only reassemble gravity vector + else if (enableGravity) + { + if (dim == dimWorld) + localAssembler.assembleGravity( handle.gravity(), + iv, + handle.advectionCA(), + LambdaFactory::getAdvectionLambda() ); + else + localAssembler.assembleGravity( handle.gravity(), + handle.gravityOutside(), + iv, + handle.advectionCA(), + handle.advectionA(), + LambdaFactory::getAdvectionLambda() ); + } + + // assemble pressure vectors + for (unsigned int pIdx = 0; pIdx < GET_PROP_VALUE(TypeTag, NumPhases); ++pIdx) + { + const auto& evv = &elemVolVars(); + auto getPressure = [&evv, pIdx] (auto volVarIdx) { return (evv->operator[](volVarIdx)).pressure(pIdx); }; + localAssembler.assemble(handle.pressures(pIdx), iv, getPressure); + } + } // fill advection caches for (unsigned int i = 0; i < iv.localFaceData().size(); ++i) + { + // set pointer to current local face data object + // ifs are evaluated at compile time and are optimized away + if (std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value) + { + // we cannot make a disctinction, thus we set both pointers + primaryLocalFaceData_ = &(iv.localFaceData()[i]); + secondaryLocalFaceData_ = &(iv.localFaceData()[i]); + } + else if (std::is_same<InteractionVolume, PrimaryInteractionVolume>::value) + primaryLocalFaceData_ = &(iv.localFaceData()[i]); + else + secondaryLocalFaceData_ = &(iv.localFaceData()[i]); + + // fill this scvfs cache AdvectionFiller::fill(*ivFluxVarCaches[i], problem(), iv.element(iv.localFaceData()[i].ivLocalInsideScvIndex()), @@ -253,25 +347,31 @@ private: elemVolVars(), *ivScvfs[i], *this); + } } //! do nothing if advection is not enabled - template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool advectionEnabled = doAdvection> - typename std::enable_if<!advectionEnabled>::type - fillAdvection(FluxVariablesCacheContainer& fluxVarsCacheContainer, - InteractionVolumeType& iv, - DataHandle& handle, - const std::vector<const SubControlVolumeFace*>& ivScvfs, - const std::vector<FluxVariablesCache*>& ivFluxVarCaches) {} - - //! method to fill the diffusive quantities - template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool diffusionEnabled = doDiffusion> - typename std::enable_if<diffusionEnabled>::type - fillDiffusion(FluxVariablesCacheContainer& fluxVarsCacheContainer, - InteractionVolumeType& iv, - DataHandle& handle, - const std::vector<const SubControlVolumeFace*>& ivScvfs, - const std::vector<FluxVariablesCache*>& ivFluxVarCaches) + template< class InteractionVolume, + class DataHandle, + bool enableAdvection = doAdvection, + typename std::enable_if_t<!enableAdvection, int> = 0 > + void fillAdvection(InteractionVolume& iv, + DataHandle& handle, + const std::vector<const SubControlVolumeFace*>& ivScvfs, + const std::vector<FluxVariablesCache*>& ivFluxVarCaches, + bool forceUpdateAll = false) + {} + + //! fills the diffusive quantities (diffusion enabled) + template< class InteractionVolume, + class DataHandle, + bool enableDiffusion = doDiffusion, + typename std::enable_if_t<enableDiffusion, int> = 0 > + void fillDiffusion(InteractionVolume& iv, + DataHandle& handle, + const std::vector<const SubControlVolumeFace*>& ivScvfs, + const std::vector<FluxVariablesCache*>& ivFluxVarCaches, + bool forceUpdateAll = false) { using DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType); using DiffusionFiller = typename DiffusionType::Cache::Filler; @@ -282,6 +382,10 @@ private: static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents); + // get instance of the interaction volume-local assembler + using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >; + IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars()); + for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx) @@ -289,15 +393,51 @@ private: if (phaseIdx == compIdx) continue; - // set the diffusion context in the data handle - handle.setDiffusionContext(phaseIdx, compIdx); - // solve the local system subject to the diffusion tensor (if uses mpfa) if (DiffusionMethod == DiscretizationMethods::CCMpfa) - iv.solveLocalSystem(LambdaFactory::getDiffusionLambda(phaseIdx, compIdx), problem(), fvGeometry(), elemVolVars(), handle); + { + // update the context in handle + handle.setPhaseIndex(phaseIdx); + handle.setComponentIndex(compIdx); + + // assemble T + if (forceUpdateAll || soldependentDiffusion) + { + if (dim < dimWorld) + localAssembler.assemble( handle.diffusionTout(), + handle.diffusionT(), + iv, + LambdaFactory::getDiffusionLambda(phaseIdx, compIdx) ); + else + localAssembler. assemble( handle.diffusionT(), + iv, + LambdaFactory::getDiffusionLambda(phaseIdx, compIdx) ); + } + + // assemble vector of mole fractions + const auto& evv = &elemVolVars(); + auto getMoleFraction = [&evv, phaseIdx, compIdx] (auto volVarIdx) + { return (evv->operator[](volVarIdx)).moleFraction(phaseIdx, compIdx); }; + localAssembler.assemble(handle.moleFractions(phaseIdx, compIdx), iv, getMoleFraction); + } // fill diffusion caches for (unsigned int i = 0; i < iv.localFaceData().size(); ++i) + { + // set pointer to current local face data object + // ifs are evaluated at compile time and are optimized away + if (std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value) + { + // we cannot make a disctinction, thus we set both pointers + primaryLocalFaceData_ = &(iv.localFaceData()[i]); + secondaryLocalFaceData_ = &(iv.localFaceData()[i]); + } + else if (std::is_same<InteractionVolume, PrimaryInteractionVolume>::value) + primaryLocalFaceData_ = &(iv.localFaceData()[i]); + else + secondaryLocalFaceData_ = &(iv.localFaceData()[i]); + + // fill this scvfs cache DiffusionFiller::fill(*ivFluxVarCaches[i], phaseIdx, compIdx, @@ -307,27 +447,33 @@ private: elemVolVars(), *ivScvfs[i], *this); + } } } } //! do nothing if diffusion is not enabled - template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool diffusionEnabled = doDiffusion> - typename std::enable_if<!diffusionEnabled>::type - fillDiffusion(FluxVariablesCacheContainer& fluxVarsCacheContainer, - InteractionVolumeType& iv, - DataHandle& handle, - const std::vector<const SubControlVolumeFace*>& ivScvfs, - const std::vector<FluxVariablesCache*>& ivFluxVarCaches) {} - - //! method to fill the quantities related to heat conduction - template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool heatConductionEnabled = doHeatConduction> - typename std::enable_if<heatConductionEnabled>::type - fillHeatConduction(FluxVariablesCacheContainer& fluxVarsCacheContainer, - InteractionVolumeType& iv, + template< class InteractionVolume, + class DataHandle, + bool enableDiffusion = doDiffusion, + typename std::enable_if_t<!enableDiffusion, int> = 0 > + void fillDiffusion(InteractionVolume& iv, DataHandle& handle, const std::vector<const SubControlVolumeFace*>& ivScvfs, - const std::vector<FluxVariablesCache*>& ivFluxVarCaches) + const std::vector<FluxVariablesCache*>& ivFluxVarCaches, + bool forceUpdateAll = false) + {} + + //! fills the quantities related to heat conduction (heat conduction enabled) + template< class InteractionVolume, + class DataHandle, + bool enableHeatConduction = doHeatConduction, + typename std::enable_if_t<enableHeatConduction, int> = 0 > + void fillHeatConduction(InteractionVolume& iv, + DataHandle& handle, + const std::vector<const SubControlVolumeFace*>& ivScvfs, + const std::vector<FluxVariablesCache*>& ivFluxVarCaches, + bool forceUpdateAll = false) { using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType); using HeatConductionFiller = typename HeatConductionType::Cache::Filler; @@ -335,15 +481,49 @@ private: static constexpr auto HeatConductionMethod = HeatConductionType::myDiscretizationMethod; using LambdaFactory = TensorLambdaFactory<TypeTag, HeatConductionMethod>; - // set the advection context in the data handle - handle.setHeatConductionContext(); - // maybe solve the local system subject to fourier coefficient if (HeatConductionMethod == DiscretizationMethods::CCMpfa) - iv.solveLocalSystem(LambdaFactory::getHeatConductionLambda(), problem(), fvGeometry(), elemVolVars(), handle); + { + // get instance of the interaction volume-local assembler + using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >; + IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars()); + + if (forceUpdateAll || soldependentAdvection) + { + if (dim < dimWorld) + localAssembler.assemble( handle.heatConductionTout(), + handle.heatConductionT(), + iv, + LambdaFactory::getHeatConductionLambda() ); + else + localAssembler.assemble( handle.heatConductionT(), + iv, + LambdaFactory::getHeatConductionLambda() ); + } + + // assemble vector of temperatures + const auto& evv = &elemVolVars(); + auto getMoleFraction = [&evv] (auto volVarIdx) { return (evv->operator[](volVarIdx)).temperature(); }; + localAssembler.assemble(handle.temperatures(), iv, getMoleFraction); + } // fill heat conduction caches for (unsigned int i = 0; i < iv.localFaceData().size(); ++i) + { + // set pointer to current local face data object + // ifs are evaluated at compile time and are optimized away + if (std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value) + { + // we cannot make a disctinction, thus we set both pointers + primaryLocalFaceData_ = &(iv.localFaceData()[i]); + secondaryLocalFaceData_ = &(iv.localFaceData()[i]); + } + else if (std::is_same<InteractionVolume, PrimaryInteractionVolume>::value) + primaryLocalFaceData_ = &(iv.localFaceData()[i]); + else + secondaryLocalFaceData_ = &(iv.localFaceData()[i]); + + // fill this scvfs cache HeatConductionFiller::fill(*ivFluxVarCaches[i], problem(), iv.element(iv.localFaceData()[i].ivLocalInsideScvIndex()), @@ -351,16 +531,20 @@ private: elemVolVars(), *ivScvfs[i], *this); + } } //! do nothing if heat conduction is disabled - template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool heatConductionEnabled = doHeatConduction> - typename std::enable_if<!heatConductionEnabled>::type - fillHeatConduction(FluxVariablesCacheContainer& fluxVarsCacheContainer, - InteractionVolumeType& iv, - DataHandle& handle, - const std::vector<const SubControlVolumeFace*>& ivScvfs, - const std::vector<FluxVariablesCache*>& ivFluxVarCaches) {} + template< class InteractionVolume, + class DataHandle, + bool enableHeatConduction = doHeatConduction, + typename std::enable_if_t<!enableHeatConduction, int> = 0 > + void fillHeatConduction(InteractionVolume& iv, + DataHandle& handle, + const std::vector<const SubControlVolumeFace*>& ivScvfs, + const std::vector<FluxVariablesCache*>& ivFluxVarCaches, + bool forceUpdateAll = false) + {} const Problem* problemPtr_; const Element* elementPtr_; @@ -374,9 +558,17 @@ private: SecondaryInteractionVolume* secondaryIv_; // pointer to the current interaction volume data handle - DataHandle* ivDataHandle_; + PrimaryDataHandle* primaryIvDataHandle_; + SecondaryDataHandle* secondaryIvDataHandle_; + + // We do an interaction volume-wise filling of the caches + // While filling, we store a pointer to the current localScvf + // face data object of the IV so that the individual caches + // can access it and don't have to retrieve it again + const PrimaryLocalFaceData* primaryLocalFaceData_; + const SecondaryLocalFaceData* secondaryLocalFaceData_; }; -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh index 14bb0fee6547a7abd637613d429799791899a469..b7a56d04e36c8535d87af5b94d1157241422e480 100644 --- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh @@ -18,8 +18,8 @@ *****************************************************************************/ /*! * \file - * \brief This file contains the class which is required to calculate - * heat conduction fluxes with Fourier's law for cell-centered MPFA models. + * \ingroup CCMpfaDiscretization + * \brief Fourier's law for cell-centered finite volume schemes with multi-point flux approximation */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH #define DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH @@ -29,14 +29,14 @@ namespace Dumux { -// forward declaration +//! forward declaration of the method-specific implementation template<class TypeTag, DiscretizationMethods discMethod> class FouriersLawImplementation; /*! - * \ingroup Mpfa - * \brief Specialization of Fourier's Law for the CCMpfa method. - */ +* \ingroup CCMpfaDiscretization +* \brief Fourier's law for cell-centered finite volume schemes with two-point flux approximation +*/ template <class TypeTag> class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> { @@ -51,15 +51,6 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel); - // Always use the dynamic type for vectors (compatibility with the boundary) - using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using CoefficientVector = typename PrimaryInteractionVolume::Traits::DynamicVector; - using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; - - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - static constexpr int energyEqIdx = GET_PROP_TYPE(TypeTag, Indices)::energyEqIdx; - //! Class that fills the cache corresponding to mpfa Darcy's Law class MpfaFouriersLawCacheFiller { @@ -78,11 +69,13 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> // get interaction volume from the flux vars cache filler & upate the cache if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.secondaryInteractionVolume(), - fluxVarsCacheFiller.dataHandle(), + fluxVarsCacheFiller.secondaryIvLocalFaceData(), + fluxVarsCacheFiller.secondaryIvDataHandle(), scvf); else scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.primaryInteractionVolume(), - fluxVarsCacheFiller.dataHandle(), + fluxVarsCacheFiller.primaryIvLocalFaceData(), + fluxVarsCacheFiller.primaryIvDataHandle(), scvf); } }; @@ -90,52 +83,91 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa> //! The cache used in conjunction with the mpfa Fourier's Law class MpfaFouriersLawCache { - // We always use the dynamic types here to be compatible on the boundary - using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer; - using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer; + // In the current implementation of the flux variables cache we cannot + // make a disctinction between dynamic (mpfa-o) and static (mpfa-l) + // matrix and vector types, as currently the cache class can only be templated + // by a type tag (and there can only be one). We use a dynamic vector here to + // make sure it works in case one of the two used interaction volume types uses + // dynamic types performance is thus lowered for schemes using static types. + // TODO: this has to be checked thoroughly as soon as a scheme using static types + // is implemented. One idea to overcome the performance drop could be only + // storing the iv-local index here and obtain tij always from the datahandle + // of the fluxVarsCacheContainer + using Vector = Dune::DynamicVector< Scalar >; + using Matrix = Dune::DynamicMatrix< Scalar >; + + using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); + using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector; + using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix; + + static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value, + "The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" ); + static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value, + "The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" ); + + using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume); + using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector; + using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix; + + static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value, + "The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" ); + static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value, + "The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" ); + + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + public: // export filler type using Filler = MpfaFouriersLawCacheFiller; - // update cached objects for heat conduction - template<class InteractionVolume> - void updateHeatConduction(const InteractionVolume& iv, const DataHandle& dataHandle, const SubControlVolumeFace &scvf) + /*! + * \brief Update cached objects (transmissibilities) + * + * \tparam InteractionVolume The (mpfa scheme-specific) interaction volume + * \tparam LocalFaceData The object used to store iv-local info on an scvf + * \tparam DataHandle The object used to store transmissibility matrices etc. + * + * \param iv The interaction volume this scvf is embedded in + * \param localFaceData iv-local info on this scvf + * \param dataHandle Transmissibility matrix & gravity data of this iv + * \param scvf The sub-control volume face + */ + template<class InteractionVolume, class LocalFaceData, class DataHandle> + void updateHeatConduction(const InteractionVolume& iv, + const LocalFaceData& localFaceData, + const DataHandle& dataHandle, + const SubControlVolumeFace &scvf) { - const auto& localFaceData = iv.getLocalFaceData(scvf); + switchFluxSign_ = localFaceData.isOutside(); - // update the quantities that are equal for all phases - heatConductionSwitchFluxSign_ = localFaceData.isOutside(); - heatConductionVolVarsStencil_ = &dataHandle.volVarsStencil(); - heatConductionDirichletData_ = &dataHandle.dirichletData(); + //! store pointer to the temperature vector of this iv + Tj_ = &dataHandle.temperatures(); - // the transmissibilities on surface grids have to be obtained from the outside + const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); if (dim == dimWorld) - heatConductionTij_ = &dataHandle.T()[localFaceData.ivLocalScvfIndex()]; + Tij_ = &dataHandle.heatConductionT()[ivLocalIdx]; else - heatConductionTij_ = localFaceData.isOutside() ? - &dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()] : - &dataHandle.T()[localFaceData.ivLocalScvfIndex()]; + Tij_ = localFaceData.isOutside() ? &dataHandle.heatConductionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()] + : &dataHandle.heatConductionT()[ivLocalIdx]; } - //! Returns the stencil for heat conduction flux computation on an scvf - const Stencil& heatConductionVolVarsStencil() const { return *heatConductionVolVarsStencil_; } - - //! Returns the transmissibilities associated with the volume variables - const CoefficientVector& heatConductionTij() const { return *heatConductionTij_; } + //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions + const Vector& heatConductionTij() const { return *Tij_; } - //! On faces that are "outside" w.r.t. a face in the interaction volume, - //! we have to take the negative value of the fluxes, i.e. multiply by -1.0 - bool heatConductionSwitchFluxSign() const { return heatConductionSwitchFluxSign_; } + //! In the interaction volume-local system of eq we have one unknown per face. + //! On scvfs on this face, but in "outside" (neighbor) elements of it, we have + //! to take the negative value of the fluxes due to the flipped normal vector. + //! This function returns whether or not this scvf is an "outside" face in the iv. + bool heatConductionSwitchFluxSign() const { return switchFluxSign_; } - //! Returns the data on dirichlet boundary conditions affecting - //! the flux computation on this face - const DirichletDataContainer& heatConductionDirichletData() const { return *heatConductionDirichletData_; } + //! The cell (& Dirichlet) temperatures within this interaction volume + const Vector& temperatures() const { return *Tj_; } private: - bool heatConductionSwitchFluxSign_; - const Stencil* heatConductionVolVarsStencil_; - const CoefficientVector* heatConductionTij_; - const DirichletDataContainer* heatConductionDirichletData_; + bool switchFluxSign_; + const Vector* Tij_; //!< The transmissibilities such that f = Tij*Tj + const Vector* Tj_; //!< The interaction-volume wide temperature Tj }; public: @@ -145,6 +177,7 @@ public: // state the type for the corresponding cache and its filler using Cache = MpfaFouriersLawCache; + //! Compute the conductive flux across an scvf static Scalar flux(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, @@ -152,22 +185,16 @@ public: const SubControlVolumeFace& scvf, const ElementFluxVarsCache& elemFluxVarsCache) { - // prepare computations - Scalar flux(0.0); - unsigned int i = 0; const auto& fluxVarsCache = elemFluxVarsCache[scvf]; const auto& tij = fluxVarsCache.heatConductionTij(); + const auto& Tj = fluxVarsCache.temperatures(); - // calculate Tij*tj - for (const auto volVarIdx : fluxVarsCache.heatConductionVolVarsStencil()) - flux += tij[i++]*elemVolVars[volVarIdx].temperature(); - - // add contributions from dirichlet BCs - for (const auto& d : fluxVarsCache.heatConductionDirichletData()) - flux += tij[i++]*elemVolVars[d.volVarIndex()].temperature(); + //! compute Tij*tj + Scalar flux = tij*Tj; + if (fluxVarsCache.heatConductionSwitchFluxSign()) + flux *= -1.0; - // return overall resulting flux - return fluxVarsCache.heatConductionSwitchFluxSign() ? -flux : flux; + return flux; } }; diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh index bb707613a506bbed39d7788007ca04cf8fe69fa0..fc4cdb07a16fd8f6b155a6daf26ac4b82ab8ff7c 100644 --- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh +++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh @@ -18,13 +18,15 @@ *****************************************************************************/ /*! * \file - * \brief Base class for a local finite volume geometry for cell-centered MPFA models + * \ingroup CCMpfaDiscretization + * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered mpfa models * This builds up the sub control volumes and sub control volume faces * for each element in the local scope we are restricting to, e.g. stencil or element. */ #ifndef DUMUX_DISCRETIZATION_CCMPFA_FV_ELEMENT_GEOMETRY_HH #define DUMUX_DISCRETIZATION_CCMPFA_FV_ELEMENT_GEOMETRY_HH +#include <dune/common/exceptions.hh> #include <dune/common/iteratorrange.hh> #include <dune/geometry/referenceelements.hh> @@ -35,30 +37,35 @@ namespace Dumux { /*! - * \ingroup ImplicitModel - * \brief Base class for the finite volume geometry vector for cell-centered MPFA models + * \ingroup CCMpfaDiscretization + * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered mpfa models * This builds up the sub control volumes and sub control volume faces - * for each element. + * for each element in the local scope we are restricting to, e.g. stencil or element. + * \note This class is specialized for versions with and without caching the fv geometries on the grid view */ template<class TypeTag, bool EnableFVGridGeometryCache> -class CCMpfaFVElementGeometry -{}; +class CCMpfaFVElementGeometry; -//! specialization in case the FVElementGeometries are stored globally -//! In this case we just forward internally to the global object +/*! + * \ingroup CCMpfaDiscretization + * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered mpfa models + * Specialization for grid caching enabled + * \note The finite volume geometries are stored in the corresponding FVGridGeometry + */ template<class TypeTag> class CCMpfaFVElementGeometry<TypeTag, true> { using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; + using Element = typename GridView::template Codim<0>::Entity; + + using GridIndexType = typename GridView::IndexSet::IndexType; using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - using Element = typename GridView::template Codim<0>::Entity; using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>; - using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>; + using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<GridIndexType>, ThisType>; + using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>; public: //! Constructor @@ -66,20 +73,20 @@ public: : fvGridGeometryPtr_(&fvGridGeometry) {} //! Get an element sub control volume with a global scv index - const SubControlVolume& scv(IndexType scvIdx) const + const SubControlVolume& scv(GridIndexType scvIdx) const { return fvGridGeometry().scv(scvIdx); } //! Get an element sub control volume face with a global scvf index - const SubControlVolumeFace& scvf(IndexType scvfIdx) const + const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const { return fvGridGeometry().scvf(scvfIdx); } - //! Get an element sub control volume face with a global scvf index - //! We separate element and neighbor scvfs to speed up mapping - const SubControlVolumeFace& flipScvf(IndexType scvfIdx, unsigned int outsideScvIdx = 0) const + //! Get the scvf on the same face but from the other side + //! Note that e.g. the normals might be different in the case of surface grids + const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvIdx = 0) const { return fvGridGeometry().flipScvf(scvfIdx, outsideScvIdx); } @@ -131,7 +138,7 @@ public: //! Bind only element-local void bindElement(const Element& element) { - scvIndices_ = std::vector<IndexType>({fvGridGeometry().elementMapper().index(element)}); + scvIndices_ = std::vector<GridIndexType>({fvGridGeometry().elementMapper().index(element)}); } //! The global finite volume geometry we are a restriction of @@ -140,11 +147,15 @@ public: private: - std::vector<IndexType> scvIndices_; + std::vector<GridIndexType> scvIndices_; const FVGridGeometry* fvGridGeometryPtr_; }; -//! specialization in case the FVElementGeometries are not stored +/*! + * \ingroup CCMpfaDiscretization + * \brief Stencil-local finite volume geometry (scvs and scvfs) for cell-centered TPFA models + * Specialization for grid caching disabled + */ template<class TypeTag> class CCMpfaFVElementGeometry<TypeTag, false> { @@ -152,15 +163,16 @@ class CCMpfaFVElementGeometry<TypeTag, false> using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; + using Element = typename GridView::template Codim<0>::Entity; + using GridIndexType = typename GridView::IndexSet::IndexType; + using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - using Element = typename GridView::template Codim<0>::Entity; using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>; - using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>; + using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<GridIndexType>, ThisType>; + using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>; static const int dim = GridView::dimension; static const int dimWorld = GridView::dimensionworld; @@ -175,7 +187,7 @@ public: //! Get an elment sub control volume with a global scv index //! We separate element and neighbor scvs to speed up mapping - const SubControlVolume& scv(IndexType scvIdx) const + const SubControlVolume& scv(GridIndexType scvIdx) const { if (scvIdx == scvIndices_[0]) return scvs_[0]; @@ -185,7 +197,7 @@ public: //! Get an element sub control volume face with a global scvf index //! We separate element and neighbor scvfs to speed up mapping - const SubControlVolumeFace& scvf(IndexType scvfIdx) const + const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const { auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx); if (it != scvfIndices_.end()) @@ -194,9 +206,9 @@ public: return neighborScvfs_[findLocalIndex(scvfIdx, neighborScvfIndices_)]; } - //! Get an element sub control volume face with a global scvf index - //! We separate element and neighbor scvfs to speed up mapping - const SubControlVolumeFace& flipScvf(IndexType scvfIdx, unsigned int outsideScvIdx = 0) const + //! Get the scvf on the same face but from the other side + //! Note that e.g. the normals might be different in the case of surface grids + const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvIdx = 0) const { auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx); if (it != scvfIndices_.end()) @@ -207,7 +219,7 @@ public: else { const auto neighborScvfIdxLocal = findLocalIndex(scvfIdx, neighborScvfIndices_); - auto&& scvf = neighborScvfs_[neighborScvfIdxLocal]; + const auto& scvf = neighborScvfs_[neighborScvfIdxLocal]; if (scvf.outsideScvIdx(outsideScvIdx) == scvIndices_[0]) return scvfs_[neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx]]; @@ -252,6 +264,7 @@ public: //! called by the local assembler to prepare element assembly void bind(const Element& element) { + // make inside geometries bindElement(element); // get some references for convenience @@ -308,11 +321,11 @@ public: private: //! Computes the number of neighboring scvfs that have to be prepared - template<class DataJ> - unsigned int numNeighborScvfs_(const std::vector<DataJ>& dataJVector) + template<class DataJContainer> + std::size_t numNeighborScvfs_(const DataJContainer& dataJContainer) { - unsigned int numNeighborScvfs = 0; - for (const auto& dataJ : dataJVector) + std::size_t numNeighborScvfs = 0; + for (const auto& dataJ : dataJContainer) numNeighborScvfs += dataJ.scvfsJ.size() + dataJ.additionalScvfs.size(); return numNeighborScvfs; } @@ -329,7 +342,7 @@ private: const auto& scvFaceIndices = fvGridGeometry().scvfIndicesOfScv(eIdx); const auto& neighborVolVarIndices = fvGridGeometry().neighborVolVarIndices(eIdx); - // the quadrature point to be used on the scvf + // the quadrature point parameterizaion to be used on scvfs static const Scalar q = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Mpfa.Q"); // reserve memory for the scv faces @@ -402,7 +415,7 @@ private: //! create the scv and necessary scvfs of a neighboring element template<typename IndexVector> void makeNeighborGeometries(const Element& element, - IndexType eIdxGlobal, + GridIndexType eIdxGlobal, const IndexVector& scvfIndices, const IndexVector& additionalScvfs) { @@ -414,7 +427,7 @@ private: const auto& scvFaceIndices = fvGridGeometry().scvfIndicesOfScv(eIdxGlobal); const auto& neighborVolVarIndices = fvGridGeometry().neighborVolVarIndices(eIdxGlobal); - // the quadrature point to be used on the scvf + // the quadrature point parameterizaion to be used on scvfs static const Scalar q = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Mpfa.Q"); // for network grids we only want to do one scvf per half facet @@ -583,14 +596,16 @@ private: } } - const IndexType findLocalIndex(const IndexType idx, - const std::vector<IndexType>& indices) const + //! map a global index to the local storage index + const unsigned int findLocalIndex(const GridIndexType idx, + const std::vector<GridIndexType>& indices) const { auto it = std::find(indices.begin(), indices.end(), idx); assert(it != indices.end() && "Could not find the scv/scvf! Make sure to properly bind this class!"); return std::distance(indices.begin(), it); } + //! clear all containers void clear() { scvIndices_.clear(); @@ -602,24 +617,27 @@ private: neighborScvfIndices_.clear(); neighborScvs_.clear(); neighborScvfs_.clear(); + + flipScvfIndices_.clear(); + neighborFlipScvfIndices_.clear(); } const FVGridGeometry* fvGridGeometryPtr_; // local storage after binding an element - std::vector<IndexType> scvIndices_; - std::vector<IndexType> scvfIndices_; + std::vector<GridIndexType> scvIndices_; + std::vector<GridIndexType> scvfIndices_; std::vector<SubControlVolume> scvs_; std::vector<SubControlVolumeFace> scvfs_; - std::vector<IndexType> neighborScvIndices_; - std::vector<IndexType> neighborScvfIndices_; + std::vector<GridIndexType> neighborScvIndices_; + std::vector<GridIndexType> neighborScvfIndices_; std::vector<SubControlVolume> neighborScvs_; std::vector<SubControlVolumeFace> neighborScvfs_; - // flip index sets - std::vector< std::vector<IndexType> > flipScvfIndices_; - std::vector< std::vector<IndexType> > neighborFlipScvfIndices_; + // flip scvf index sets + std::vector< std::vector<GridIndexType> > flipScvfIndices_; + std::vector< std::vector<GridIndexType> > neighborFlipScvfIndices_; }; } // end namespace diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh index 858b4fec5b3484794142284dd80309d22b457220..dc34d0f5c0454dbb10e86123d576d7ef68522419 100644 --- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh +++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh @@ -18,9 +18,10 @@ *****************************************************************************/ /*! * \file - * \brief Base class for the finite volume geometry vector for mpfa models + * \ingroup CCMpfaDiscretization + * \brief The finite volume geometry (scvs and scvfs) for cell-centered mpfa models on a grid view * This builds up the sub control volumes and sub control volume faces - * for each element. + * for each element of the grid partition. */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_HH #define DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_HH @@ -37,129 +38,106 @@ namespace Dumux { /*! - * \ingroup Mpfa - * \brief Base class for the finite volume geometry vector for mpfa models + * \ingroup CCMpfaDiscretization + * \brief The finite volume geometry (scvs and scvfs) for cell-centered mpfa models on a grid view * This builds up the sub control volumes and sub control volume faces - * for each element. + * \note This class is specialized for versions with and without caching the fv geometries on the grid view */ template<class TypeTag, bool EnableFVElementGeometryCache> -class CCMpfaFVGridGeometry -{}; +class CCMpfaFVGridGeometry; -// specialization in case the finite volume grid geometries are stored +/*! + * \ingroup CCMpfaDiscretization + * \brief The finite volume geometry (scvs and scvfs) for cell-centered mpfa models on a grid view + * This builds up the sub control volumes and sub control volume faces + * \note For caching enabled we store the fv geometries for the whole grid view which is memory intensive but faster + */ template<class TypeTag> class CCMpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag> { using ParentType = BaseFVGridGeometry<TypeTag>; using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper); + using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); + using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); + using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using GridIVIndexSets = CCMpfaGridInteractionVolumeIndexSets<TypeTag>; - using ConnectivityMap = CCMpfaConnectivityMap<TypeTag>; using GridView = typename GET_PROP_TYPE(TypeTag, GridView); static constexpr int dim = GridView::dimension; static constexpr int dimWorld = GridView::dimensionworld; + using Element = typename GridView::template Codim<0>::Entity; using Vertex = typename GridView::template Codim<dim>::Entity; using Intersection = typename GridView::Intersection; using CoordScalar = typename GridView::ctype; - using IndexType = typename GridView::IndexSet::IndexType; + using GridIndexType = typename GridView::IndexSet::IndexType; + using LocalIndexType = typename PrimaryInteractionVolume::Traits::LocalIndexType; + + using GridIVIndexSets = CCMpfaGridInteractionVolumeIndexSets<TypeTag>; + using ConnectivityMap = CCMpfaConnectivityMap<TypeTag>; - using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>; using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>; public: - using SecondaryIvIndicator = std::function<bool(const Element&, const Intersection&, bool)>; + using SecondaryIvIndicatorType = std::function<bool(const Element&, const Intersection&, bool)>; //! Constructor without indicator function for secondary interaction volumes //! Per default, we use the secondary IVs at branching points & boundaries - explicit CCMpfaFVGridGeometry(const GridView& gridView) - : ParentType(gridView) - , secondaryIvIndicator_([] (const Element& e, const Intersection& is, bool isBranching) - { return is.boundary() || isBranching; } ) - {} - - //! Constructor with user-defined indicator function for secondary interaction volumes - explicit CCMpfaFVGridGeometry(const GridView& gridView, const SecondaryIvIndicator& indicator) - : ParentType(gridView) - , secondaryIvIndicator_(indicator) - {} + CCMpfaFVGridGeometry(const GridView& gridView) + : ParentType(gridView) + , secondaryIvIndicator_([] (const Element& e, const Intersection& is, bool isBranching) + { return is.boundary() || isBranching; } ) + {} + + //! Constructor with user-defined indicator function for secondary interaction volumes + CCMpfaFVGridGeometry(const GridView& gridView, const SecondaryIvIndicatorType& indicator) + : ParentType(gridView) + , secondaryIvIndicator_(indicator) + {} //! the element mapper is the dofMapper //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa... - const ElementMapper& dofMapper() const - { return this->elementMapper(); } + const ElementMapper& dofMapper() const { return this->elementMapper(); } + + //! The total number of sub control volumes + std::size_t numScv() const { return scvs_.size(); } + + //! The total number of sub control volume faces + std::size_t numScvf() const { return scvfs_.size(); } + + //! The total number of boundary sub control volume faces + std::size_t numBoundaryScvf() const { return numBoundaryScvf_; } + + //! The total number of degrees of freedom + std::size_t numDofs() const { return this->gridView().size(0); } - /*! - * \brief Returns the total number of sub control volumes. - */ - std::size_t numScv() const - { return scvs_.size(); } - - /*! - * \brief Returns the total number of sub control volume faces. - */ - std::size_t numScvf() const - { return scvfs_.size(); } - - /*! - * \brief Returns the number of scvfs on the domain boundary. - */ - std::size_t numBoundaryScvf() const - { return numBoundaryScvf_; } - - /*! - * \brief Returns the total number of degrees of freedom. - */ - std::size_t numDofs() const - { return this->gridView().size(0); } - - /*! - * \brief Gets an element from a global element index. - */ - Element element(IndexType eIdx) const - { return this->elementMap()[eIdx]; } - - /*! - * \brief Gets an element from a sub control volume contained in it. - */ - Element element(const SubControlVolume& scv) const - { return this->elementMap()[scv.elementIndex()]; } - - /*! - * \brief Returns true if primary interaction volumes are used around a given vertex, - * false otherwise. - */ + //! Get an element from a global element index + Element element(GridIndexType eIdx) const { return this->elementMap()[eIdx]; } + + //! Get an element from a sub control volume contained in it + Element element(const SubControlVolume& scv) const { return this->elementMap()[scv.elementIndex()]; } + + //! Returns true if primary interaction volumes are used around a given vertex. bool vertexUsesPrimaryInteractionVolume(const Vertex& v) const { return primaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; } - /*! - * \brief Returns true if primary interaction volumes are used around a given vertex index, - * false otherwise. - */ - bool vertexUsesPrimaryInteractionVolume(IndexType vIdxGlobal) const + //!Returns true if primary interaction volumes are used around a vertex (index). + bool vertexUsesPrimaryInteractionVolume(GridIndexType vIdxGlobal) const { return primaryInteractionVolumeVertices_[vIdxGlobal]; } - /*! - * \brief Returns if primary interaction volumes are used around a given vertex. - */ + //! Returns if primary interaction volumes are used around a given vertex. bool vertexUsesSecondaryInteractionVolume(const Vertex& v) const { return secondaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; } - /*! - * \brief Returns true if primary interaction volumes are used around a given vertex index, - * false otherwise. - */ - bool vertexUsesSecondaryInteractionVolume(IndexType vIdxGlobal) const + //! Returns true if primary interaction volumes are used around a given vertex index. + bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const { return secondaryInteractionVolumeVertices_[vIdxGlobal]; } - /*! - * \brief Updates all finite volume geometries of the grid. Hhas to be called again after grid adaptation. - */ + //! update all fvElementGeometries (do this again after grid adaption) void update() { ParentType::update(); @@ -189,10 +167,10 @@ public: const auto isGhostVertex = MpfaHelper::findGhostVertices(this->gridView(), this->vertexMapper()); // instantiate the dual grid index set (to be used for construction of interaction volumes) - CCMpfaDualGridIndexSet<TypeTag> dualIdSet(this->gridView()); + CCMpfaDualGridIndexSet<GridIndexType, LocalIndexType, dim> dualIdSet(this->gridView()); // Build the SCVs and SCV faces - IndexType scvfIdx = 0; + GridIndexType scvfIdx = 0; numBoundaryScvf_ = 0; for (const auto& element : elements(this->gridView())) { @@ -202,12 +180,12 @@ public: auto elementGeometry = element.geometry(); // The local scvf index set - std::vector<IndexType> scvfIndexSet; + std::vector<GridIndexType> scvfIndexSet; scvfIndexSet.reserve(MpfaHelper::getNumLocalScvfs(elementGeometry.type())); // for network grids there might be multiple intersections with the same geometryInInside // we indentify those by the indexInInside for now (assumes conforming grids at branching facets) - std::vector<std::vector<IndexType>> outsideIndices; + std::vector<std::vector<GridIndexType>> outsideIndices; if (dim < dimWorld) { outsideIndices.resize(element.subEntities(1)); @@ -269,7 +247,7 @@ public: secondaryInteractionVolumeVertices_[vIdxGlobal] = true; } - // the quadrature point to be used on the scvfs + // the quadrature point parameterizarion to be used on scvfs static const Scalar q = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Mpfa.Q"); // make the scv face (for non-boundary scvfs on network grids, use precalculated outside indices) @@ -278,15 +256,15 @@ public: if (!boundary) { return dim == dimWorld ? - std::vector<IndexType>({this->elementMapper().index(is.outside())}) : + std::vector<GridIndexType>({this->elementMapper().index(is.outside())}) : outsideIndices[indexInInside]; } else { // compute boundary scv idx and increment counter - const IndexType bIdx = numScvs + numBoundaryScvf_; + const GridIndexType bIdx = numScvs + numBoundaryScvf_; numBoundaryScvf_++; - return std::vector<IndexType>(1, bIdx); + return std::vector<GridIndexType>(1, bIdx); } } (); @@ -368,43 +346,27 @@ public: std::cout << "Initializing of the connectivity map took " << timer.elapsed() << " seconds." << std::endl; } - /*! - * \brief Returns a sub control volume for a given global scv index. - */ - const SubControlVolume& scv(IndexType scvIdx) const - { return scvs_[scvIdx]; } - - /*! - * \brief Returns a sub control volume face for a given global scvf index. - */ - const SubControlVolumeFace& scvf(IndexType scvfIdx) const - { return scvfs_[scvfIdx]; } - - /*! - * \brief Returns the "flipped" scvf, i.e. the correspongin scvf in an outside element. - * The second argument is optional and only comes into play on network/surface grids. - */ - const SubControlVolumeFace& flipScvf(IndexType scvfIdx, unsigned int outsideScvfIdx = 0) const - { return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; } + //! Get a sub control volume with a global scv index + const SubControlVolume& scv(GridIndexType scvIdx) const { return scvs_[scvIdx]; } - /*! - * \brief Returns the sub control volume face indices of an scv by global index. - */ - const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const - { return scvfIndicesOfScv_[scvIdx]; } + //! Get a sub control volume face with a global scvf index + const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const { return scvfs_[scvfIdx]; } + + //! Returns the connectivity map of which dofs + //! have derivatives with respect to a given dof. + const ConnectivityMap& connectivityMap() const { return connectivityMap_; } + + //! Returns the grid interaction volume seeds class. + const GridIVIndexSets& gridInteractionVolumeIndexSets() const { return ivIndexSets_; } - /*! - * \brief Returns the connectivity map of which dofs have derivatives with respect - * to a given dof. - */ - const ConnectivityMap& connectivityMap() const - { return connectivityMap_; } + //! Get the scvf on the same face but from the other side + //! Note that e.g. the normals might be different in the case of surface grids + const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const + { return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; } - /*! - * \brief Returns the grit interaction volume seeds class. - */ - const GridIVIndexSets& gridInteractionVolumeIndexSets() const - { return ivIndexSets_; } + //! Get the sub control volume face indices of an scv by global index + const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const + { return scvfIndicesOfScv_[scvIdx]; } private: // connectivity map for efficient assembly @@ -415,45 +377,56 @@ private: std::vector<SubControlVolumeFace> scvfs_; // containers storing the global data - std::vector<std::vector<IndexType>> scvfIndicesOfScv_; + std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_; std::vector<bool> primaryInteractionVolumeVertices_; std::vector<bool> secondaryInteractionVolumeVertices_; std::size_t numBoundaryScvf_; // needed for embedded surface and network grids (dim < dimWorld) - std::vector<std::vector<IndexType>> flipScvfIndices_; + std::vector<std::vector<GridIndexType>> flipScvfIndices_; // The grid interaction volume index set GridIVIndexSets ivIndexSets_; // Indicator function on where to use the secondary IVs - SecondaryIvIndicator secondaryIvIndicator_; + SecondaryIvIndicatorType secondaryIvIndicator_; }; -// specialization in case the FVElementGeometries are not stored +/*! + * \ingroup CCMpfaDiscretization + * \brief The finite volume geometry (scvs and scvfs) for cell-centered mpfa models on a grid view + * This builds up the sub control volumes and sub control volume faces + * \note For caching disabled we store only some essential index maps to build up local systems on-demand in + * the corresponding FVElementGeometry + */ template<class TypeTag> class CCMpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag> { using ParentType = BaseFVGridGeometry<TypeTag>; using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper); + using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); + using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); + using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using GridIVIndexSets = CCMpfaGridInteractionVolumeIndexSets<TypeTag>; - using ConnectivityMap = CCMpfaConnectivityMap<TypeTag>; using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper); static constexpr int dim = GridView::dimension; static constexpr int dimWorld = GridView::dimensionworld; + using Element = typename GridView::template Codim<0>::Entity; using Vertex = typename GridView::template Codim<dim>::Entity; using Intersection = typename GridView::Intersection; using CoordScalar = typename GridView::ctype; - using IndexType = typename GridView::IndexSet::IndexType; + using GridIndexType = typename GridView::IndexSet::IndexType; + using LocalIndexType = typename PrimaryInteractionVolume::Traits::LocalIndexType; + + using GridIVIndexSets = CCMpfaGridInteractionVolumeIndexSets<TypeTag>; + using ConnectivityMap = CCMpfaConnectivityMap<TypeTag>; - using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>; using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>; public: @@ -461,102 +434,64 @@ public: //! Constructor without indicator function for secondary interaction volumes //! Per default, we use the secondary IVs at branching points & boundaries - explicit CCMpfaFVGridGeometry(const GridView& gridView) - : ParentType(gridView) - , secondaryIvIndicator_([] (const Element& e, const Intersection& is, bool isBranching) - { return is.boundary() || isBranching; } ) - {} - - //! Constructor with user-defined indicator function for secondary interaction volumes - explicit CCMpfaFVGridGeometry(const GridView& gridView, const SecondaryIvIndicator& indicator) - : ParentType(gridView) - , secondaryIvIndicator_(indicator) - {} + CCMpfaFVGridGeometry(const GridView& gridView) + : ParentType(gridView) + , secondaryIvIndicator_([] (const Element& e, const Intersection& is, bool isBranching) + { return is.boundary() || isBranching; } ) + {} + + //! Constructor with user-defined indicator function for secondary interaction volumes + CCMpfaFVGridGeometry(const GridView& gridView, const SecondaryIvIndicator& indicator) + : ParentType(gridView) + , secondaryIvIndicator_(indicator) + {} //! the element mapper is the dofMapper //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa... const ElementMapper& dofMapper() const { return this->elementMapper(); } - /*! - * \brief Returns the total number of sub control volumes. - */ - std::size_t numScv() const - { return numScvs_; } - - /*! - * \brief Returns the total number of sub control volume faces. - */ - std::size_t numScvf() const - { return numScvf_; } - - /*! - * \brief Returns the number of scvfs on the domain boundary. - */ - std::size_t numBoundaryScvf() const - { return numBoundaryScvf_; } - - /*! - * \brief Returns the total number of degrees of freedom. - */ - std::size_t numDofs() const - { return this->gridView().size(0); } - - /*! - * \brief Gets an element from a global element index. - */ - Element element(IndexType eIdx) const - { return this->elementMap()[eIdx]; } - - /*! - * \brief Gets an element from a sub control volume contained in it. - */ - Element element(const SubControlVolume& scv) const - { return this->elementMap()[scv.elementIndex()]; } - - /*! - * \brief Returns true if primary interaction volumes are used around a given vertex, - * false otherwise. - */ + //! Returns the total number of sub control volumes. + std::size_t numScv() const { return numScvs_; } + + //! Returns the total number of sub control volume faces. + std::size_t numScvf() const { return numScvf_; } + + //! Returns the number of scvfs on the domain boundary. + std::size_t numBoundaryScvf() const { return numBoundaryScvf_; } + + //! Returns the total number of degrees of freedom. + std::size_t numDofs() const { return this->gridView().size(0); } + + //! Gets an element from a global element index. + Element element(GridIndexType eIdx) const { return this->elementMap()[eIdx]; } + + //! Gets an element from a sub control volume contained in it. + Element element(const SubControlVolume& scv) const { return this->elementMap()[scv.elementIndex()]; } + + //! Returns true if primary interaction volumes are used around a given vertex. bool vertexUsesPrimaryInteractionVolume(const Vertex& v) const { return primaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; } - /*! - * \brief Returns true if primary interaction volumes are used around a given vertex index, - * false otherwise. - */ - bool vertexUsesPrimaryInteractionVolume(IndexType vIdxGlobal) const + //! Returns true if primary interaction volumes are used around a given vertex (index). + bool vertexUsesPrimaryInteractionVolume(GridIndexType vIdxGlobal) const { return primaryInteractionVolumeVertices_[vIdxGlobal]; } - /*! - * \brief Returns if primary interaction volumes are used around a given vertex. - */ + //! Returns if primary interaction volumes are used around a given vertex. bool vertexUsesSecondaryInteractionVolume(const Vertex& v) const { return secondaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; } - /*! - * \brief Returns true if primary interaction volumes are used around a given vertex index, - * false otherwise. - */ - bool vertexUsesSecondaryInteractionVolume(IndexType vIdxGlobal) const + //! Returns true if primary interaction volumes are used around a given vertex (index). + bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const { return secondaryInteractionVolumeVertices_[vIdxGlobal]; } - /*! - * \brief Returns true if a given vertex lies on a processor boundary inside a ghost element. - */ - bool isGhostVertex(const Vertex& v) const - { return isGhostVertex_[this->vertexMapper().index(v)]; } - - /*! - * \brief Returns true if the vertex corresponding to a given vertex index lies on a - * processor boundary inside a ghost element. - */ - bool isGhostVertex(IndexType vIdxGlobal) const - { return isGhostVertex_[vIdxGlobal]; } - - /*! - * \brief Updates all finite volume geometries of the grid. Has to be called again after grid adaption. - */ + //! Returns true if a given vertex lies on a processor boundary inside a ghost element. + bool isGhostVertex(const Vertex& v) const { return isGhostVertex_[this->vertexMapper().index(v)]; } + + //! Returns true if the vertex (index) lies on a processor boundary inside a ghost element. + bool isGhostVertex(GridIndexType vIdxGlobal) const { return isGhostVertex_[vIdxGlobal]; } + + //! Updates all finite volume geometries of the grid. Has to be called again after grid adaption. void update() { ParentType::update(); @@ -579,7 +514,7 @@ public: isGhostVertex_ = MpfaHelper::findGhostVertices(this->gridView(), this->vertexMapper()); // instantiate the dual grid index set (to be used for construction of interaction volumes) - CCMpfaDualGridIndexSet<TypeTag> dualIdSet(this->gridView()); + CCMpfaDualGridIndexSet<GridIndexType, LocalIndexType, dim> dualIdSet(this->gridView()); // Build the SCVs and SCV faces numScvf_ = 0; @@ -593,14 +528,14 @@ public: // the element-wise index sets for finite volume geometry const auto numLocalFaces = MpfaHelper::getNumLocalScvfs(elementGeometry.type()); - std::vector<IndexType> scvfsIndexSet; - std::vector< std::vector<IndexType> > neighborVolVarIndexSet; + std::vector<GridIndexType> scvfsIndexSet; + std::vector< std::vector<GridIndexType> > neighborVolVarIndexSet; scvfsIndexSet.reserve(numLocalFaces); neighborVolVarIndexSet.reserve(numLocalFaces); // for network grids there might be multiple intersections with the same geometryInInside // we indentify those by the indexInInside for now (assumes conforming grids at branching facets) - std::vector<std::vector<IndexType>> outsideIndices; + std::vector<std::vector<GridIndexType>> outsideIndices; if (dim < dimWorld) { outsideIndices.resize(element.subEntities(1)); @@ -661,15 +596,15 @@ public: if (!boundary) { return dim == dimWorld ? - std::vector<IndexType>({this->elementMapper().index(is.outside())}) : + std::vector<GridIndexType>({this->elementMapper().index(is.outside())}) : outsideIndices[indexInInside]; } else { // compute boundary scv idx and increment counter - const IndexType bIdx = numScvs_ + numBoundaryScvf_; + const GridIndexType bIdx = numScvs_ + numBoundaryScvf_; numBoundaryScvf_++; - return std::vector<IndexType>(1, bIdx); + return std::vector<GridIndexType>(1, bIdx); } } (); @@ -705,38 +640,28 @@ public: std::cout << "Initializing of the connectivity map took " << timer.elapsed() << " seconds." << std::endl; } - /*! - * \brief Returns the sub control volume face indices of an scv by global index. - */ - const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const + //! Returns the sub control volume face indices of an scv by global index. + const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const { return scvfIndicesOfScv_[scvIdx]; } - /*! - * \brief Returns the neighboring vol var indices for each scvf contained in an scv. - */ - const std::vector< std::vector<IndexType> >& neighborVolVarIndices(IndexType scvIdx) const + //! Returns the neighboring vol var indices for each scvf contained in an scv. + const std::vector< std::vector<GridIndexType> >& neighborVolVarIndices(GridIndexType scvIdx) const { return neighborVolVarIndices_[scvIdx]; } - /*! - * \brief Returns the connectivity map of which dofs have derivatives with respect - * to a given dof. - */ - const ConnectivityMap& connectivityMap() const - { return connectivityMap_; } + //! Returns the connectivity map of which dofs + //! have derivatives with respect to a given dof. + const ConnectivityMap& connectivityMap() const { return connectivityMap_; } - /*! - * \brief Returns the grit interaction volume seeds class. - */ - const GridIVIndexSets& gridInteractionVolumeIndexSets() const - { return ivIndexSets_; } + //! Returns the grid interaction volume seeds class. + const GridIVIndexSets& gridInteractionVolumeIndexSets() const { return ivIndexSets_; } private: // connectivity map for efficient assembly ConnectivityMap connectivityMap_; // containers storing the global data - std::vector<std::vector<IndexType>> scvfIndicesOfScv_; - std::vector< std::vector< std::vector<IndexType> > > neighborVolVarIndices_; + std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_; + std::vector< std::vector< std::vector<GridIndexType> > > neighborVolVarIndices_; std::vector<bool> primaryInteractionVolumeVertices_; std::vector<bool> secondaryInteractionVolumeVertices_; std::vector<bool> isGhostVertex_; @@ -751,6 +676,6 @@ private: SecondaryIvIndicator secondaryIvIndicator_; }; -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/generalconnectivitymap.hh b/dumux/discretization/cellcentered/mpfa/generalconnectivitymap.hh index 26beb8a932812cf6d21b8d26458627330fc212e7..24166a0539f63d0e3d2b472858c29df8a5cd50f2 100644 --- a/dumux/discretization/cellcentered/mpfa/generalconnectivitymap.hh +++ b/dumux/discretization/cellcentered/mpfa/generalconnectivitymap.hh @@ -18,6 +18,7 @@ *****************************************************************************/ /*! * \file + * \ingroup CCMpfaDiscretization * \brief Stores the face indices corresponding to the neighbors of an element * that contribute to the derivative calculation */ @@ -34,7 +35,7 @@ namespace Dumux { /*! - * \ingroup CellCentered + * \ingroup CCMpfaDiscretization * \brief General version of the assembly map for cellcentered schemes. To each * cell I we store a list of cells J that are needed to compute the fluxes * in these cells J that depend on cell I. Furthermore, we store for each cell J @@ -47,7 +48,7 @@ class CCMpfaGeneralConnectivityMap using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; + using GridIndexType = typename GridView::IndexSet::IndexType; using FluxStencil = Dumux::FluxStencil<TypeTag>; // To each cell "globalI" there will be a list of "globalJ", in which globalI is part @@ -55,9 +56,9 @@ class CCMpfaGeneralConnectivityMap // additional scvfs which are needed temporarily to set up the transmissibilities of the scvfsJ struct DataJ { - IndexType globalJ; - std::vector<IndexType> scvfsJ; - std::vector<IndexType> additionalScvfs; + GridIndexType globalJ; + std::vector<GridIndexType> scvfsJ; + std::vector<GridIndexType> additionalScvfs; }; using Map = std::vector<std::vector<DataJ>>; @@ -81,7 +82,7 @@ public: fvGeometry.bindElement(element); // obtain the data of J in elements I - std::vector<std::pair<IndexType, std::vector<DataJ>>> dataJForI; + std::vector<std::pair<GridIndexType, std::vector<DataJ>>> dataJForI; // loop over sub control faces for (auto&& scvf : scvfs(fvGeometry)) @@ -111,7 +112,7 @@ public: // land in the list of additional scvfs. Of that list we will delete those // that are already in the list of scvfsJ later... const auto scvfVectorAtVertex = MpfaHelper::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry); - std::vector<IndexType> scvfIndicesAtVertex(scvfVectorAtVertex.size()); + std::vector<GridIndexType> scvfIndicesAtVertex(scvfVectorAtVertex.size()); for (std::size_t i = 0; i < scvfVectorAtVertex.size(); ++i) scvfIndicesAtVertex[i] = scvfVectorAtVertex[i]->index(); globalJDataJ.additionalScvfs.insert(globalJDataJ.additionalScvfs.end(), @@ -131,7 +132,7 @@ public: // if entry for globalJ2 does not exist yet, add globalJ2 to the J-data of globalI // with an empty set of scvfs over which I and J are coupled (i.e. they aren't coupled) if (it2 == it->second.end()) - it->second.push_back(DataJ({globalJ2, std::vector<IndexType>()})); + it->second.push_back(DataJ({globalJ2, std::vector<GridIndexType>()})); } } else @@ -139,13 +140,13 @@ public: // No DataJ for globalI exists yet. Make it and insert data on the actual // global J as first entry in the vector of DataJs belonging to globalI dataJForI.emplace_back(std::make_pair(globalI, - std::vector<DataJ>({DataJ({globalJ, std::vector<IndexType>({scvf.index()})})}))); + std::vector<DataJ>({DataJ({globalJ, std::vector<GridIndexType>({scvf.index()})})}))); // Also, all scvfs connected to a vertex together with the actual scvf // land in the list of additional scvfs. Of that list we will delete those // that are already in the list of scvfsJ later... const auto scvfVectorAtVertex = MpfaHelper::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry); - std::vector<IndexType> scvfIndicesAtVertex(scvfVectorAtVertex.size()); + std::vector<GridIndexType> scvfIndicesAtVertex(scvfVectorAtVertex.size()); for (unsigned int i = 0; i < scvfVectorAtVertex.size(); ++i) scvfIndicesAtVertex[i] = scvfVectorAtVertex[i]->index(); dataJForI.back().second[0].additionalScvfs.insert(dataJForI.back().second[0].additionalScvfs.end(), @@ -155,7 +156,7 @@ public: // all the other dofs in the stencil will be "globalJ" to globalI as well for (auto globalJ2 : stencil) if (globalJ2 != globalJ && globalJ2 != globalI) - dataJForI.back().second.push_back(DataJ({globalJ2, std::vector<IndexType>()})); + dataJForI.back().second.push_back(DataJ({globalJ2, std::vector<GridIndexType>()})); } } } @@ -191,7 +192,8 @@ public: } } - const std::vector<DataJ>& operator[] (const IndexType globalI) const + //! Returns the assembly map of the element with given grid index + const std::vector<DataJ>& operator[] (const GridIndexType globalI) const { return map_[globalI]; } private: diff --git a/dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh index e56c7288accdeca7c5427270de2d2bcbe01c1f4d..1ac49768a65ae2b545caac995bf8006ce52451a6 100644 --- a/dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh +++ b/dumux/discretization/cellcentered/mpfa/gridfluxvariablescache.hh @@ -18,7 +18,8 @@ *****************************************************************************/ /*! * \file - * \brief The global object of flux var caches + * \ingroup CCMpfaDiscretization + * \brief Flux variable caches on a gridview */ #ifndef DUMUX_DISCRETIZATION_CCMPFA_GRID_FLUXVARSCACHE_HH #define DUMUX_DISCRETIZATION_CCMPFA_GRID_FLUXVARSCACHE_HH @@ -30,54 +31,56 @@ namespace Dumux { /*! - * \ingroup Mpfa - * \brief Base class for the flux variables cache vector, we store one cache per face + * \ingroup CCMpfaDiscretization + * \brief Flux variable caches on a gridview + * \note The class is specialized for a version with and without grid caching */ template<class TypeTag, bool EnableGridFluxVariablesCache> class CCMpfaGridFluxVariablesCache; - /*! - * \ingroup Mpfa - * \brief Spezialization when caching globally + * \ingroup CCMpfaDiscretization + * \brief Flux variable caches on a gridview with grid caching enabled + * \note The flux caches of the gridview are stored which is memory intensive but faster */ template<class TypeTag> class CCMpfaGridFluxVariablesCache<TypeTag, true> { - // the flux variables cache filler needs to be friend to fill - // the interaction volumes and data handles - friend CCMpfaFluxVariablesCacheFiller<TypeTag>; - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; + using IndexType = typename GridView::IndexSet::IndexType; + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); - using IndexType = typename GridView::IndexSet::IndexType; using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); + using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume); - using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle; + using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle; using FluxVariablesCacheFiller = CCMpfaFluxVariablesCacheFiller<TypeTag>; public: + //! The constructor CCMpfaGridFluxVariablesCache(const Problem& problem) : problemPtr_(&problem) {} - // When global caching is enabled, precompute transmissibilities for all scv faces + //! When global caching is enabled, precompute transmissibilities for all scv faces void update(const FVGridGeometry& fvGridGeometry, const GridVolumeVariables& gridVolVars, const SolutionVector& sol, bool forceUpdate = false) { - // only do the update if fluxes are solution dependent or if update is forced + // Update only if the filler puts solution-dependent + // stuff into the caches or if update is enforced if (FluxVariablesCacheFiller::isSolDependent || forceUpdate) { - // clear data if forced update is desired + // clear previous data if forced update is desired if (forceUpdate) { clear_(); @@ -89,17 +92,20 @@ public: secondaryInteractionVolumes_.reserve(numSecondaryIVs); primaryIvDataHandles_.reserve(numPrimaryIvs); secondaryIvDataHandles_.reserve(numSecondaryIVs); - } - // reserve memory estimate for caches, interaction volumes and corresponding data - fluxVarsCache_.resize(fvGridGeometry.numScvf()); + // reserve memory estimate for caches, interaction volumes and corresponding data + fluxVarsCache_.resize(fvGridGeometry.numScvf()); + } // instantiate helper class to fill the caches FluxVariablesCacheFiller filler(problem()); + // set all the caches to "outdated" + for (auto& cache : fluxVarsCache_) + cache.setUpdateStatus(false); + for (const auto& element : elements(fvGridGeometry.gridView())) { - // Prepare the geometries within the elements of the stencil auto fvGeometry = localView(fvGridGeometry); fvGeometry.bind(element); @@ -118,7 +124,8 @@ public: const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) { - // update only if transmissibilities are solution-dependent + // Update only if the filler puts + // solution-dependent stuff into the caches if (FluxVariablesCacheFiller::isSolDependent) { const auto& fvGridGeometry = fvGeometry.fvGridGeometry(); @@ -150,7 +157,6 @@ public: auto& scvfCache = fluxVarsCache_[scvfIdx]; if (!scvfCache.isUpdated()) { - // update cache const auto& scvf = fvGeometry.scvf(scvfIdx); filler.fill(*this, scvfCache, elementJ, fvGeometry, elemVolVars, scvf); } @@ -167,18 +173,52 @@ public: friend inline ElementFluxVariablesCache localView(const CCMpfaGridFluxVariablesCache& global) { return ElementFluxVariablesCache(global); } - // access operators in the case of caching + //! access operators in the case of caching const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const { return fluxVarsCache_[scvf.index()]; } + //! access operators in the case of caching FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) { return fluxVarsCache_[scvf.index()]; } + //! access to the stored interaction volumes + const std::vector<PrimaryInteractionVolume>& primaryInteractionVolumes() const + { return primaryInteractionVolumes_; } + + //! access to the stored interaction volumes + std::vector<PrimaryInteractionVolume>& primaryInteractionVolumes() + { return primaryInteractionVolumes_; } + + //! access to the stored data handles + const std::vector<PrimaryIvDataHandle>& primaryDataHandles() const + { return primaryIvDataHandles_; } + + //! access to the stored data handles + std::vector<PrimaryIvDataHandle>& primaryDataHandles() + { return primaryIvDataHandles_; } + + //! access to the stored interaction volumes + const std::vector<SecondaryInteractionVolume>& secondaryInteractionVolumes() const + { return secondaryInteractionVolumes_; } + + //! access to the stored interaction volumes + std::vector<SecondaryInteractionVolume>& secondaryInteractionVolumes() + { return secondaryInteractionVolumes_; } + + //! access to the stored data handles + const std::vector<SecondaryIvDataHandle>& secondaryDataHandles() const + { return secondaryIvDataHandles_; } + + //! access to the stored data handles + std::vector<SecondaryIvDataHandle>& secondaryDataHandles() + { return secondaryIvDataHandles_; } + const Problem& problem() const { return *problemPtr_; } private: + //! clear all containers void clear_() { fluxVarsCache_.clear(); @@ -194,13 +234,13 @@ private: // store the interaction volumes and handles std::vector<PrimaryInteractionVolume> primaryInteractionVolumes_; std::vector<SecondaryInteractionVolume> secondaryInteractionVolumes_; - std::vector<DataHandle> primaryIvDataHandles_; - std::vector<DataHandle> secondaryIvDataHandles_; + std::vector<PrimaryIvDataHandle> primaryIvDataHandles_; + std::vector<SecondaryIvDataHandle> secondaryIvDataHandles_; }; /*! - * \ingroup ImplicitModel - * \brief Spezialization when not using global caching + * \ingroup CCMpfaDiscretization + * \brief Flux variable caches on a gridview with grid caching disabled */ template<class TypeTag> class CCMpfaGridFluxVariablesCache<TypeTag, false> @@ -216,14 +256,16 @@ class CCMpfaGridFluxVariablesCache<TypeTag, false> using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); public: + //! The constructor CCMpfaGridFluxVariablesCache(const Problem& problem) : problemPtr_(&problem) {} - // When global flux variables caching is disabled, we don't need to update the cache + //! When global flux variables caching is disabled, we don't need to update the cache void update(const FVGridGeometry& fvGridGeometry, const GridVolumeVariables& gridVolVars, const SolutionVector& sol, bool forceUpdate = false) {} + //! When global flux variables caching is disabled, we don't need to update the cache void updateElement(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars) {} @@ -243,6 +285,6 @@ private: const Problem* problemPtr_; }; -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh index 8c34070b8d44502372abe23afbafbe5ed749fd97..0cc347bf1b8dbaddd29f667868b7dcdaf6918b4a 100644 --- a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh +++ b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh @@ -18,6 +18,7 @@ *****************************************************************************/ /*! * \file + * \ingroup CCMpfaDiscretization * \brief Class for the grid interaction volume index sets of mpfa schemes. */ #ifndef DUMUX_DISCRETIZATION_MPFA_O_GRIDINTERACTIONVOLUME_INDEXSETS_HH @@ -26,34 +27,39 @@ #include <memory> #include <dumux/common/properties.hh> +#include "dualgridindexset.hh" + namespace Dumux { -// forward declaration -template<class TypeTag> -class CCMpfaDualGridIndexSet; /*! - * \ingroup Mpfa - * \brief The grid interaction volume index sets class for the mpfa-o scheme. + * \ingroup CCMpfaDiscretization + * \brief Class that holds all interaction volume index sets on a grid view. */ template<class TypeTag> class CCMpfaGridInteractionVolumeIndexSets { using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using IndexType = typename GridView::IndexSet::IndexType; + using GridIndexType = typename GridView::IndexSet::IndexType; using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using PrimaryIV = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); using PrimaryIVIndexSet = typename PrimaryIV::Traits::IndexSet; using SecondaryIV = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume); using SecondaryIVIndexSet = typename SecondaryIV::Traits::IndexSet; - using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); - using DualGridIndexSet = CCMpfaDualGridIndexSet<TypeTag>; - static const int dim = GridView::dimension; + static constexpr int dim = GridView::dimension; + using LocalIndexType = typename PrimaryIV::Traits::LocalIndexType; + using DualGridIndexSet = CCMpfaDualGridIndexSet< GridIndexType, LocalIndexType, dim>; public: - + /*! + * \brief Construct all interaction volume index sets on the grid view + * + * \param fvGridGeometry The finite volume geometry on the grid view + * \param dualGridIdSet The index sets of the dual grid on the grid view + */ void update(FVGridGeometry& fvGridGeometry, DualGridIndexSet&& dualGridIdSet) { dualGridIndexSet_ = std::make_unique<DualGridIndexSet>(std::move(dualGridIdSet)); @@ -91,32 +97,37 @@ public: } } + //! Return the iv index set in which a given scvf is embedded in const PrimaryIVIndexSet& primaryIndexSet(const SubControlVolumeFace& scvf) const { return primaryIndexSet(scvf.index()); } - const PrimaryIVIndexSet& primaryIndexSet(const IndexType scvfIdx) const + //! Return the iv index set in which a given scvf (index) is embedded in + const PrimaryIVIndexSet& primaryIndexSet(const GridIndexType scvfIdx) const { return primaryIVIndexSets_[scvfIndexMap_[scvfIdx]]; } + //! Return the iv index set in which a given scvf is embedded in const SecondaryIVIndexSet& secondaryIndexSet(const SubControlVolumeFace& scvf) const { return secondaryIndexSet(scvf.index()); } - const SecondaryIVIndexSet& secondaryIndexSet(const IndexType scvfIdx) const + //! Return the iv index set in which a given scvf (index) is embedded in + const SecondaryIVIndexSet& secondaryIndexSet(const GridIndexType scvfIdx) const { return secondaryIVIndexSets_[scvfIndexMap_[scvfIdx]]; } + //! Returns number of primary/secondary interaction volumes on the grid view std::size_t numPrimaryInteractionVolumes() const { return numPrimaryIV_; } std::size_t numSecondaryInteractionVolumes() const { return numSecondaryIV_; } private: std::vector<PrimaryIVIndexSet> primaryIVIndexSets_; std::vector<SecondaryIVIndexSet> secondaryIVIndexSets_; - std::vector<IndexType> scvfIndexMap_; + std::vector<GridIndexType> scvfIndexMap_; std::size_t numPrimaryIV_; std::size_t numSecondaryIV_; std::unique_ptr<DualGridIndexSet> dualGridIndexSet_; }; -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh index 988abe331ed52f3e74aae3116f0ea7f37e388837..93d05eb5bce6846391d0414180ec7379db78e29b 100644 --- a/dumux/discretization/cellcentered/mpfa/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/helper.hh @@ -18,7 +18,8 @@ *****************************************************************************/ /*! * \file - * \brief Helper class to get information required for mpfa scheme. + * \ingroup CCMpfaDiscretization + * \brief Helper class to get data required for mpfa scheme. */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_HELPER_HH #define DUMUX_DISCRETIZATION_CC_MPFA_HELPER_HH @@ -37,19 +38,16 @@ namespace Dumux { /*! - * \brief Mpfa method-specific implementation of the helper class (dimension-dependent) - */ -template<class TypeTag, MpfaMethods Method, int dim, int dimWorld> -class MpfaMethodHelper; - -/*! - * \brief Dimension-specific implementation of the helper class (common for all methods) + * \ingroup CCMpfaDiscretization + * \brief Dimension-specific helper class to get data required for mpfa scheme. */ template<class TypeTag, int dim, int dimWorld> class MpfaDimensionHelper; /*! - * \brief Specialization for dim == 2 & dimWorld == 2 + * \ingroup CCMpfaDiscretization + * \brief Dimension-specific helper class to get data required for mpfa scheme. + * Specialization for dim == 2 & dimWorld == 2 */ template<class TypeTag> class MpfaDimensionHelper<TypeTag, /*dim*/2, /*dimWorld*/2> @@ -60,19 +58,20 @@ class MpfaDimensionHelper<TypeTag, /*dim*/2, /*dimWorld*/2> using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage; using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using ScvBasis = typename InteractionVolume::Traits::ScvBasis; + using LocalScvType = typename InteractionVolume::Traits::LocalScvType; + using ScvBasis = typename LocalScvType::LocalBasis; - //! We know that dim = 2 & dimworld = 2. However, the specialization for - //! dim = 2 & dimWorld = 3 reuses some methods of this class, thus, we - //! obtain the world dimension from the grid view here to get the - //! GlobalPosition vector right. Be picky about the dimension though. + // We know that dim = 2 & dimworld = 2. However, the specialization for + // dim = 2 & dimWorld = 3 reuses some methods of this class, thus, we + // obtain the world dimension from the grid view here to get the + // GlobalPosition vector right. Be picky about the dimension though. static_assert(GridView::dimension == 2, "The chosen mpfa helper expects a grid view with dim = 2"); static const int dim = 2; static const int dimWorld = GridView::dimensionworld; using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; - //! Container to store the positions of intersections required for scvf - //! corner computation. In 2d, these are the center plus the two corners + // Container to store the positions of intersections required for scvf + // corner computation. In 2d, these are the center plus the two corners using ScvfPositionsOnIntersection = std::array<GlobalPosition, 3>; public: @@ -104,7 +103,7 @@ public: * * \param scvBasis The basis of an scv */ - static Scalar calculateDetX(const ScvBasis& scvBasis) + static typename LocalScvType::ctype calculateDetX(const ScvBasis& scvBasis) { using std::abs; return abs(crossProduct<Scalar>(scvBasis[0], scvBasis[1])); @@ -227,15 +226,19 @@ public: }; /*! - * \brief Specialization for dim == 2 & dimWorld == 2. Reuses some - * functionality of the specialization for dim = dimWorld = 2 + * \ingroup CCMpfaDiscretization + * \brief Dimension-specific helper class to get data required for mpfa scheme. + * Specialization for dim == 2 & dimWorld == 2. Reuses some functionality + * of the specialization for dim = dimWorld = 2 */ template<class TypeTag> -class MpfaDimensionHelper<TypeTag, /*dim*/2, /*dimWorld*/3> : public MpfaDimensionHelper<TypeTag, 2, 2> +class MpfaDimensionHelper<TypeTag, /*dim*/2, /*dimWorld*/3> + : public MpfaDimensionHelper<TypeTag, 2, 2> { using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using ScvBasis = typename InteractionVolume::Traits::ScvBasis; + using LocalScvType = typename InteractionVolume::Traits::LocalScvType; + using ScvBasis = typename LocalScvType::LocalBasis; public: @@ -269,7 +272,7 @@ public: * * \param scvBasis The basis of an scv */ - static Scalar calculateDetX(const ScvBasis& scvBasis) + static typename LocalScvType::ctype calculateDetX(const ScvBasis& scvBasis) { using std::abs; return abs(crossProduct<Scalar>(scvBasis[0], scvBasis[1]).two_norm()); @@ -285,9 +288,10 @@ public: static bool isRightHandSystem(const ScvBasis& scvBasis) { return true; } }; - /*! - * \brief Specialization for dim == 3 & dimWorld == 3. + * \ingroup CCMpfaDiscretization + * \brief Dimension-specific helper class to get data required for mpfa scheme. + * Specialization for dim == 3 & dimWorld == 3. */ template<class TypeTag> class MpfaDimensionHelper<TypeTag, /*dim*/3, /*dimWorld*/3> @@ -298,7 +302,8 @@ class MpfaDimensionHelper<TypeTag, /*dim*/3, /*dimWorld*/3> using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage; using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using ScvBasis = typename InteractionVolume::Traits::ScvBasis; + using LocalScvType = typename InteractionVolume::Traits::LocalScvType; + using ScvBasis = typename LocalScvType::LocalBasis; // Be picky about the dimensions static_assert(GridView::dimension == 3 && GridView::dimensionworld == 3, @@ -343,7 +348,7 @@ public: * * \param scvBasis The basis of an scv */ - static Scalar calculateDetX(const ScvBasis& scvBasis) + static typename LocalScvType::ctype calculateDetX(const ScvBasis& scvBasis) { using std::abs; return abs(tripleProduct<Scalar>(scvBasis[0], scvBasis[1], scvBasis[2])); @@ -542,13 +547,12 @@ public: }; /*! - * \ingroup Mpfa + * \ingroup CCMpfaDiscretization * \brief Helper class to get the required information on an interaction volume. * Specializations depending on the method and dimension are provided. */ template<class TypeTag, MpfaMethods Method, int dim, int dimWorld> -class CCMpfaHelperImplementation : public MpfaDimensionHelper<TypeTag, dim, dimWorld>, - public MpfaMethodHelper<TypeTag, Method, dim, dimWorld> +class CCMpfaHelperImplementation : public MpfaDimensionHelper<TypeTag, dim, dimWorld> { using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); @@ -560,12 +564,11 @@ class CCMpfaHelperImplementation : public MpfaDimensionHelper<TypeTag, dim, dimW using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage; + using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType; - using ScvBasis = typename InteractionVolume::Traits::ScvBasis; + using GlobalPosition = typename InteractionVolume::Traits::GlobalPosition; - using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; - using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>; using CoordScalar = typename GridView::ctype; using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>; @@ -628,25 +631,11 @@ public: template<typename V1, typename V2> static bool vectorContainsValue(const std::vector<V1>& vector, const V2 value) { return std::find(vector.begin(), vector.end(), value) != vector.end(); } - - // calculates the product of a transposed vector n, a Matrix M and another vector v (n^T M v) - static Scalar nT_M_v(const GlobalPosition& n, const DimWorldMatrix& M, const GlobalPosition& v) - { - GlobalPosition tmp; - M.mv(v, tmp); - return n*tmp; - } - - // calculates the product of a transposed vector n, a Scalar M and another vector v (n^T M v) - static Scalar nT_M_v(const GlobalPosition& n, const Scalar M, const GlobalPosition& v) - { - return M*(n*v); - } }; /*! - * \ingroup Mpfa - * \brief Helper class for the mpfa methods for the construction of the interaction regions etc. + * \ingroup CCMpfaDiscretization + * \brief Helper class to obtain data required for mpfa methods. * It inherits from dimension-, dimensionworld- and method-specific implementations. */ template<class TypeTag> @@ -657,7 +646,4 @@ using CCMpfaHelper = CCMpfaHelperImplementation<TypeTag, } // end namespace Dumux -// The implemented helper classes need to be included here -#include <dumux/discretization/cellcentered/mpfa/omethod/helper.hh> - #endif diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/interactionvolume.hh index b3bd11aab4a55b42ce566d504672bafab1970dde..0923b59f2d0f9b28b92be31683fff12f81d854b4 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolume.hh @@ -18,7 +18,8 @@ *****************************************************************************/ /*! * \file - * \brief Base classes for interaction volume of mpfa methods. + * \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 @@ -28,18 +29,20 @@ namespace Dumux { -// forward declaration of the base class + +// forward declaration of the method-specific implementation template<class TypeTag, MpfaMethods MpfaMethod> class CCMpfaInteractionVolumeImplementation; /*! - * \ingroup Mpfa - * \brief Base class for the interaction volumes of the mpfa method + * \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 +} // 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> diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh index acc01de45ce209fa815826e6b1efdbdd13f7f6d8..b84a708764d11d695c0e0a7f6fc4aa0252f29e37 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh @@ -18,200 +18,126 @@ *****************************************************************************/ /*! * \file + * \ingroup CCMpfaDiscretization * \brief Base class for interaction volumes of mpfa methods. */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEBASE_HH #define DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEBASE_HH -#include <dune/common/dynmatrix.hh> - -#include <dumux/discretization/cellcentered/mpfa/methods.hh> -#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh> -#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh> +#include <dumux/common/properties.hh> namespace Dumux { -/*! - * \ingroup Mpfa - * \brief Base class for the interaction volume traits. The types stated here - * have to be defined in interaction volume traits. It is recommended - * that different implementations inherit from this class and overwrite the - * desired types or publicly state the types of this base class. +/* + * \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. + * + * \code + * //! export the interaction volume type to be used on boundaries etc. + * using SecondaryInteractionVolume = ...; + * //! 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) + * using LocalFaceDataContainer = ...; + * //! export the type used for iv-local matrices + * using Matrix = ...; + * //! export the type used for iv-local vectors + * using Vector = ...; + * //! export the data handle type for this iv + * using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >; + * \endcode */ -template<class TypeTag> -class CCMpfaInteractionVolumeTraitsBase -{ - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - - static const int dim = GridView::dimension; - static const int dimWorld = GridView::dimensionworld; - using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; - -public: - using LocalIndexType = std::uint8_t; - - //! The dynamic types are used e.g. by the mpfa-o method. - //! To be compatible with schemes using both dynamic and static - //! array types (e.g. L-method using mpfa-o interaction volumes - //! on the boudaries), other classes interacting with the interaction - //! volumes (e.g. flux vars cache) export the dynamic types. If your - //! scheme is fully static on the entire grid, overwrite these traits. - using DynamicLocalIndexContainer = std::vector<LocalIndexType>; - using DynamicGlobalIndexContainer = std::vector<typename GridView::IndexSet::IndexType>; - using DynamicMatrix = Dune::DynamicMatrix<Scalar>; - using DynamicVector = typename DynamicMatrix::row_type; - - //! The data handle type. Uses the dynamic types as well per default... - using DataHandle = InteractionVolumeDataHandle<TypeTag>; - - using ScvfVector = std::array<const SubControlVolumeFace*, dim>; - using ScvBasis = std::array<GlobalPosition, dim>; - - //! for network grids this means that we assume the tensors - //! to be given in world coordinates! If a transformation of - //! given data has to be performed, it has to be done in the - //! spatial parameters method where the tensor is returned or - //! in the volume variables where it is stored - using Tensor = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>; -}; +template< class InteractionVolume > +struct CCMpfaInteractionVolumeTraits {}; /*! - * \ingroup Mpfa + * \ingroup CCMpfaDiscretization * \brief Base class for the interaction volumes of mpfa methods. It defines * the interface and actual implementations should derive from this class. */ -template<class TypeTag, typename T> +template<class TypeTag, class Implementation> class CCMpfaInteractionVolumeBase { - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + // Curiously recurring template pattern + Implementation & asImp() { return static_cast<Implementation&>(*this); } + const Implementation & asImp() const { return static_cast<const Implementation&>(*this); } + using GridView = typename GET_PROP_TYPE(TypeTag, 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 SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); - using DualGridNodalIndexSet = typename CCMpfaDualGridIndexSet<TypeTag>::NodalIndexSet; - using DataHandle = typename T::DataHandle; - using IndexSet = typename T::IndexSet; - using LocalIndexContainer = typename T::DynamicLocalIndexContainer; - using LocalIndexType = typename LocalIndexContainer::value_type; - using GlobalIndexContainer = typename T::DynamicGlobalIndexContainer; - using GlobalIndexType = typename GlobalIndexContainer::value_type; + //! state the traits type publicly + using Traits = CCMpfaInteractionVolumeTraits< Implementation >; - static const int dim = GridView::dimension; - static const int dimWorld = GridView::dimensionworld; - using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; - -public: - // state the traits type publicly - using Traits = T; - - class LocalFaceData - { - LocalIndexType ivLocalScvfIndex_; //!< the iv-local scvf index this scvf maps to - LocalIndexType ivLocalInsideScvIndex_; //!< the iv-local index of the scvfs' inside scv - LocalIndexType ivLocalOutsideScvfIndex_; //!< the index of this scvf in the iv-local outside faces - LocalIndexType scvfLocalOutsideScvfIndex_; //!< the index of this scvf in the scvf-local outside faces - GlobalIndexType globalScvfIndex_; //!< the index of the corresponding global scvf - bool isOutside_; //!< indicates if this face maps to the iv-local index from "outside" - - public: - //! Constructor for "inside" faces - explicit LocalFaceData(LocalIndexType faceIndex, - LocalIndexType scvIndex, - GlobalIndexType globalScvfIndex) - : ivLocalScvfIndex_(faceIndex), - ivLocalInsideScvIndex_(scvIndex), - globalScvfIndex_(globalScvfIndex), - isOutside_(false) {} - - //! Constructor for "outside" faces - explicit LocalFaceData(LocalIndexType faceIndex, - LocalIndexType scvIndex, - LocalIndexType indexInIvOutsideFaces, - LocalIndexType indexInScvfOutsideFaces, - GlobalIndexType globalScvfIndex) - : ivLocalScvfIndex_(faceIndex), - ivLocalInsideScvIndex_(scvIndex), - ivLocalOutsideScvfIndex_(indexInIvOutsideFaces), - scvfLocalOutsideScvfIndex_(indexInScvfOutsideFaces), - globalScvfIndex_(globalScvfIndex), - isOutside_(true) {} - - //! The index of the scvf within the inside faces - LocalIndexType ivLocalScvfIndex() const { return ivLocalScvfIndex_; } - LocalIndexType ivLocalInsideScvIndex() const { return ivLocalInsideScvIndex_; } - LocalIndexType ivLocalOutsideScvfIndex() const { assert(isOutside_); return ivLocalOutsideScvfIndex_; } - LocalIndexType scvfLocalOutsideScvfIndex() const { assert(isOutside_); return scvfLocalOutsideScvfIndex_; } - GlobalIndexType globalScvfIndex() const { return globalScvfIndex_; } - bool isOutside() const { return isOutside_; } - }; - - struct DirichletData - { - GlobalIndexType volVarIndex_; - GlobalPosition ipGlobal_; - - public: - explicit DirichletData(const GlobalIndexType index, const GlobalPosition& ip) - : volVarIndex_(index) - , ipGlobal_(ip) - {} - - const GlobalPosition& ipGlobal() const { return ipGlobal_; } - GlobalIndexType volVarIndex() const { return volVarIndex_; } - }; - - using DirichletDataContainer = std::vector<DirichletData>; - using LocalFaceDataContainer = std::vector<LocalFaceData>; - - //! Sets up the local scope (geometries etc) for a given iv index set! - void setUpLocalScope(const IndexSet& indexSet, + //! Prepares everything for the assembly + void setUpLocalScope(const typename Traits::IndexSet& indexSet, const Problem& problem, const FVElementGeometry& fvGeometry) - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a setUpLocalScope() method."); } + { asImp().setUpLocalScope(); } + + //! returns the number of "primary" scvfs of this interaction volume + std::size_t numFaces() const { return asImp().numFaces(); } + + //! returns the number of intermediate unknowns within this interaction volume + std::size_t numUnknowns() const { return asImp().numUnknowns(); } + + //! returns the number of (in this context) known solution values within this interaction volume + std::size_t numKnowns() const { return asImp().numKnowns(); } + + //! returns the number of scvs embedded in this interaction volume + std::size_t numScvs() const { return asImp().numScvs(); } - //! sets the sizes of the corresponding matrices in the data handle - void prepareDataHandle(DataHandle& dataHandle) - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a prepareDataHandle() method."); } + //! returns the number of scvfs embedded in this interaction volume + std::size_t numScvfs() const { return asImp().numScvfs(); } - //! solves for the transmissibilities subject to a given tensor - template<typename GetTensorFunction> - void solveLocalSystem(const GetTensorFunction& getTensor, - const Problem& problem, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - DataHandle& dataHandle) - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a solveLocalSystem() method."); } + //! returns a reference to the container with the local face data + const typename Traits::LocalFaceDataContainer& localFaceData() const { asImp().localFaceData(); } - //! obtain the local data object for a given global scvf - const LocalFaceData& getLocalFaceData(const SubControlVolumeFace& scvf) const - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a getLocalFaceData() method."); } + //! returns the local scvf entity corresponding to a given iv-local scvf idx + const typename Traits::LocalScvfType& localScvf(typename Traits::LocalIndexType ivLocalScvfIdx) const + { return asImp().localScvf(ivLocalScvfIdx); } - //!returns a reference to the container with the local face data - const LocalFaceDataContainer& localFaceData() const - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a localFaceData() method."); } + //! returns the local scv entity corresponding to a given iv-local scv idx + const typename Traits::LocalScvType& localScv(typename Traits::LocalIndexType ivLocalScvIdx) const + { return asImp().localScv(ivLocalScvIdx); } - //! returns the indices of the volvars in the stencil of the interaction volume - const GlobalIndexContainer& volVarsStencil() const - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a volVarsStencil() method."); } + //! returns the element in which the scv with the given local idx is embedded in + const Element& element(typename Traits::LocalIndexType ivLocalScvIdx) const + { return asImp().element(); } //! returns the number of interaction volumes living around a vertex - static std::size_t numInteractionVolumesAtVertex(const DualGridNodalIndexSet& nodalIndexSet) - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a numInteractionVolumesAtVertex() method."); } + template< class NodalIndexSet > + static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) + { return Implementation::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 - template<class IvIndexSetContainer, class ScvfIndexMap> + template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet> static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer, ScvfIndexMap& scvfIndexMap, - const DualGridNodalIndexSet& nodalIndexSet) - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a numInteractionVolumesAtVertex() method."); } + const NodalIndexSet& nodalIndexSet) + { Implementation::addInteractionVolumeIndexSets(ivIndexSetContainer, scvfIndexMap, nodalIndexSet); } }; -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh index 520372856f35d807635adbe647eaf9db8588db1b..8bad1b63e8eb70938b9075ad5417a267d0a31972 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh @@ -24,383 +24,301 @@ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEDATAHANDLE_HH #define DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEDATAHANDLE_HH -#include <dune/common/exceptions.hh> +#include <vector> + #include <dumux/common/properties.hh> namespace Dumux { - //! Empty data handle class - template<class TypeTag> - class EmptyDataHandle - { - //! we use the dynamic types here to be compatible on the boundary - using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using DirichletDataContainer = typename InteractionVolume::DirichletDataContainer; - using GlobalIndexContainer = typename InteractionVolume::Traits::DynamicGlobalIndexContainer; - using Matrix = typename InteractionVolume::Traits::DynamicMatrix; - using Vector = typename InteractionVolume::Traits::DynamicVector; - - public: - //! diffusion caches need to set phase and component index - void setDiffusionContext(unsigned int phaseIdx, unsigned int compIdx) {} - - //! functions to set the size of the matrices - void resizeT(unsigned int n, unsigned int m) {} - void resizeAB(unsigned int n, unsigned int m) {} - void resizeOutsideTij(unsigned int n, unsigned int m) {} - - //! functions to set the pointers to the stencil - void setVolVarsStencilPointer(const GlobalIndexContainer& stencil) {} - - //! return functions for the stored data - const GlobalIndexContainer& volVarsStencil() const { return throw_<const GlobalIndexContainer&>(); } - const DirichletDataContainer& dirichletData() const { return throw_<DirichletDataContainer&>(); } - - const Matrix& T() const { return throw_<const Matrix&>(); } - Matrix& T() { return throw_<Matrix&>(); } - - const Matrix& AB() const { return throw_<const Matrix&>(); } - Matrix& AB() { return throw_<Matrix&>(); } - - const Matrix& outsideTij() const { return throw_<const Matrix&>(); } - Matrix& outsideTij() { return throw_<Matrix&>(); } - - private: - template<class ReturnType> - ReturnType throw_() const { DUNE_THROW(Dune::InvalidStateException, "Trying to access data for a deactivated physical process"); } - }; - - //! Data handle for quantities related to advection - template<class TypeTag, bool EnableAdvection> - class AdvectionDataHandle - { - //! we use the dynamic types here to be compatible on the boundary - using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using DirichletDataContainer = typename InteractionVolume::DirichletDataContainer; - using GlobalIndexContainer = typename InteractionVolume::Traits::DynamicGlobalIndexContainer; - using Matrix = typename InteractionVolume::Traits::DynamicMatrix; - using Vector = typename InteractionVolume::Traits::DynamicVector; - - public: - //! set the sizes of the matrices - void resizeT(unsigned int n, unsigned int m) { advectionT_.resize(n, m); } - void resizeAB(unsigned int n, unsigned int m) { advectionAB_.resize(n, m); } - void resizeOutsideTij(unsigned int n, unsigned int m) { advectionTout_.resize(n, m); } - - //! sets the pointer to the stencil - void setVolVarsStencilPointer(const GlobalIndexContainer& stencil) { advectionVolVarsStencil_ = &stencil; } - - //! return functions for the stored data - const GlobalIndexContainer& volVarsStencil() const { return *advectionVolVarsStencil_; } - - const Matrix& T() const { return advectionT_; } - Matrix& T() { return advectionT_; } - - const Matrix& AB() const { return advectionAB_; } - Matrix& AB() { return advectionAB_; } - - const Matrix& outsideTij() const { return advectionTout_; } - Matrix& outsideTij() { return advectionTout_; } - - private: - // advection-related variables - const GlobalIndexContainer* advectionVolVarsStencil_; //!< Pointer to the global volvar indices (stored in the interaction volume) - Matrix advectionT_; //!< The transmissibilities - Matrix advectionAB_; //!< Coefficients for gradient reconstruction - Matrix advectionTout_; //!< The transmissibilities associated with "outside" faces (only necessary on surface grids) - }; - - //! Data handle for quantities related to diffusion - template<class TypeTag, bool EnableDiffusion> - class DiffusionDataHandle - { - //! we use the dynamic types here to be compatible on the boundary - using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using DirichletDataContainer = typename InteractionVolume::DirichletDataContainer; - using GlobalIndexContainer = typename InteractionVolume::Traits::DynamicGlobalIndexContainer; - using Matrix = typename InteractionVolume::Traits::DynamicMatrix; - using Vector = typename InteractionVolume::Traits::DynamicVector; - - static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); - static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents); - - public: - //! diffusion caches need to set phase and component index - void setDiffusionContext(unsigned int phaseIdx, unsigned int compIdx) - { - contextPhaseIdx_ = phaseIdx; - contextCompIdx_ = compIdx; - } - - //! set the sizes of the matrices - void resizeT(unsigned int n, unsigned int m) - { - for (auto& array : diffusionT_) - for (auto& matrix : array) - matrix.resize(n, m); - } - - void resizeAB(unsigned int n, unsigned int m) - { - for (auto& array : diffusionAB_) - for (auto& matrix : array) - matrix.resize(n, m); - } - void resizeOutsideTij(unsigned int n, unsigned int m) - { - for (auto& array : diffusionTout_) - for (auto& matrix : array) - matrix.resize(n, m); - } - - //! sets the pointer to stencil - void setVolVarsStencilPointer(const GlobalIndexContainer& stencil) - { - diffusionVolVarsStencil_[contextPhaseIdx_][contextCompIdx_] = &stencil; - } - - //! return functions for the stored data - const GlobalIndexContainer& volVarsStencil() const - { return *diffusionVolVarsStencil_[contextPhaseIdx_][contextCompIdx_]; } +//! Empty data handle class +template<class InteractionVolume> +class EmptyDataHandle +{ +public: + void resize(const InteractionVolume& iv) {} +}; - const Matrix& T() const { return diffusionT_[contextPhaseIdx_][contextCompIdx_]; } - Matrix& T() { return diffusionT_[contextPhaseIdx_][contextCompIdx_]; } +//! Data handle for quantities related to advection +template<class TypeTag, class InteractionVolume, bool EnableAdvection> +class AdvectionDataHandle +{ + //! export matrix & vector types from interaction volume + using Matrix = typename InteractionVolume::Traits::Matrix; + using Vector = typename InteractionVolume::Traits::Vector; + using Scalar = typename Vector::value_type; + using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType; - const Matrix& AB() const { return diffusionAB_[contextPhaseIdx_][contextCompIdx_]; } - Matrix& AB() { return diffusionAB_[contextPhaseIdx_][contextCompIdx_]; } + using OutsideDataContainer = std::vector< std::vector<Vector> >; - const Matrix& outsideTij() const { return diffusionTout_[contextPhaseIdx_][contextCompIdx_]; } - Matrix& outsideTij() { return diffusionTout_[contextPhaseIdx_][contextCompIdx_]; } + static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension; + static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld; + static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); - private: - // diffusion-related variables (see comments in AdvectionDataHandle) - unsigned int contextPhaseIdx_; //!< The phase index set for the context - unsigned int contextCompIdx_; //!< The component index set for the context - std::array<std::array<const GlobalIndexContainer*, numComponents>, numPhases> diffusionVolVarsStencil_; - std::array<std::array<Matrix, numComponents>, numPhases> diffusionT_; - std::array<std::array<Matrix, numComponents>, numPhases> diffusionAB_; - std::array<std::array<Matrix, numComponents>, numPhases> diffusionTout_; - }; +public: - //! Data handle for quantities related to advection - template<class TypeTag, bool EnableHeatConduction> - class HeatConductionDataHandle - { - //! we use the dynamic types here to be compatible on the boundary - using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using DirichletDataContainer = typename InteractionVolume::DirichletDataContainer; - using GlobalIndexContainer = typename InteractionVolume::Traits::DynamicGlobalIndexContainer; - using Matrix = typename InteractionVolume::Traits::DynamicMatrix; - using Vector = typename InteractionVolume::Traits::DynamicVector; - - public: - //! set the sizes of the matrices - void resizeT(unsigned int n, unsigned int m) { heatConductionT_.resize(n, m); } - void resizeAB(unsigned int n, unsigned int m) { heatConductionAB_.resize(n, m); } - void resizeOutsideTij(unsigned int n, unsigned int m) { heatConductionTout_.resize(n, m); } - - //! sets the pointer to the stencil - void setVolVarsStencilPointer(const GlobalIndexContainer& stencil) { heatConductionVolVarsStencil_ = &stencil; } - - //! return functions for the stored data - const GlobalIndexContainer& volVarsStencil() const { return *heatConductionVolVarsStencil_; } - - const Matrix& T() const { return heatConductionT_; } - Matrix& T() { return heatConductionT_; } - - const Matrix& AB() const { return heatConductionAB_; } - Matrix& AB() { return heatConductionAB_; } - - const Matrix& outsideTij() const { return heatConductionTout_; } - Matrix& outsideTij() { return heatConductionTout_; } - - private: - // heat conduction-related variables - const GlobalIndexContainer* heatConductionVolVarsStencil_; //!< Pointer to the global volvar indices (stored in the interaction volume) - Matrix heatConductionT_; //!< The transmissibilities - Matrix heatConductionAB_; //!< Coefficients for gradient reconstruction - Matrix heatConductionTout_; //!< The transmissibilities associated with "outside" faces (only necessary on surface grids) - }; - - //! Process-dependet data handle when related process is disabled - template<class TypeTag> class AdvectionDataHandle<TypeTag, false> : public EmptyDataHandle<TypeTag> {}; - template<class TypeTag> class DiffusionDataHandle<TypeTag, false> : public EmptyDataHandle<TypeTag> {}; - template<class TypeTag> class HeatConductionDataHandle<TypeTag, false> : public EmptyDataHandle<TypeTag> {}; - - //! Interaction volume data handle class - template<class TypeTag> - class InteractionVolumeDataHandle : public AdvectionDataHandle<TypeTag, GET_PROP_VALUE(TypeTag, EnableAdvection)>, - public DiffusionDataHandle<TypeTag, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>, - public HeatConductionDataHandle<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)> + //! prepare data handle for subsequent fill (normal grids) + template< int d = dim, int dw = dimWorld, std::enable_if_t<(d==dw), int> = 0> + void resize(const InteractionVolume& iv) { - using AdvectionHandle = AdvectionDataHandle<TypeTag, GET_PROP_VALUE(TypeTag, EnableAdvection)>; - using DiffusionHandle = DiffusionDataHandle<TypeTag, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>; - using HeatConductionHandle = HeatConductionDataHandle<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>; - - //! we use the dynamic types here to be compatible on the boundary - using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using DirichletDataContainer = typename InteractionVolume::DirichletDataContainer; - using GlobalIndexContainer = typename InteractionVolume::Traits::DynamicGlobalIndexContainer; - using Matrix = typename InteractionVolume::Traits::DynamicMatrix; - using Vector = typename InteractionVolume::Traits::DynamicVector; - - public: - enum class Contexts : unsigned int + // resize transmissibility matrix & pressure vectors + T_.resize(iv.numFaces(), iv.numKnowns()); + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + p_[pIdx].resize(iv.numKnowns()); + + // maybe resize gravity container + static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"); + if (enableGravity) { - undefined, - advection, - diffusion, - heatConduction - }; - - //! The constructor - InteractionVolumeDataHandle() : context_(Contexts::undefined) {} - - //! set the context of the cache - void setAdvectionContext() { context_ = Contexts::advection; } - void setHeatConductionContext() { context_ = Contexts::heatConduction; } - void setDiffusionContext(unsigned int phaseIdx, unsigned int compIdx) - { - context_ = Contexts::diffusion; - DiffusionHandle::setDiffusionContext(phaseIdx, compIdx); + CA_.resize(iv.numFaces(), iv.numUnknowns()); + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + g_[pIdx].resize(iv.numFaces()); } + } - //! returns the current context - Contexts getContext() const { return context_; } - - //! set the sizes of the matrices - void resizeT(unsigned int n, unsigned int m) - { - AdvectionHandle::resizeT(n, m); - DiffusionHandle::resizeT(n, m); - HeatConductionHandle::resizeT(n, m); - } - void resizeAB(unsigned int n, unsigned int m) - { - AdvectionHandle::resizeAB(n, m); - DiffusionHandle::resizeAB(n, m); - HeatConductionHandle::resizeAB(n, m); - } + //! prepare data handle for subsequent fill (surface grids) + template< 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"); - void resizeOutsideTij(unsigned int n, unsigned int m) + if (!enableGravity) { - AdvectionHandle::resizeOutsideTij(n, m); - DiffusionHandle::resizeOutsideTij(n, m); - HeatConductionHandle::resizeOutsideTij(n, m); + T_.resize(iv.numFaces(), iv.numKnowns()); + outsideT_.resize(iv.numFaces()); + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + p_[pIdx].resize(iv.numKnowns()); + + for (LocalIndexType i = 0; i < iv.numFaces(); ++i) + { + const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1; + outsideT_[i].resize(numNeighbors); + for (LocalIndexType j = 0; j < numNeighbors; ++j) + outsideT_[i][j].resize(iv.numKnowns()); + } } - //! sets the pointer to the stencil - void setVolVarsStencilPointer(const GlobalIndexContainer& stencil) + else { - if (context_ == Contexts::advection) - AdvectionHandle::setVolVarsStencilPointer(stencil); - else if (context_ == Contexts::diffusion) - DiffusionHandle::setVolVarsStencilPointer(stencil); - else if (context_ == Contexts::heatConduction) - HeatConductionHandle::setVolVarsStencilPointer(stencil); - else - DUNE_THROW(Dune::InvalidStateException, "No valid context set!"); - } + T_.resize(iv.numFaces(), iv.numKnowns()); + CA_.resize(iv.numFaces(), iv.numUnknowns()); + A_.resize(iv.numUnknowns(), iv.numUnknowns()); + outsideT_.resize(iv.numFaces()); + + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + { + p_[pIdx].resize(iv.numKnowns()); + g_[pIdx].resize(iv.numFaces()); + outsideG_[pIdx].resize(iv.numFaces()); + } + + for (LocalIndexType i = 0; i < iv.numFaces(); ++i) + { + const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1; + outsideT_[i].resize(numNeighbors); + + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + outsideG_[pIdx][i].resize(numNeighbors); + + for (LocalIndexType j = 0; j < numNeighbors; ++j) + outsideT_[i][j].resize(iv.numKnowns()); + } - //! sets the dirichlet data container - void setDirichletData(DirichletDataContainer&& data) - { - dirichletData_ = std::move(data); } + } + + //! Access to the iv-wide pressure of one phase + const Vector& pressures(unsigned int pIdx) const { return p_[pIdx]; } + Vector& pressures(unsigned int pIdx) { return p_[pIdx]; } + + //! Access to the gravitational flux contributions for one phase + const Vector& gravity(unsigned int pIdx) const { return g_[pIdx]; } + Vector& gravity(unsigned int pIdx) { return g_[pIdx]; } + + //! Access to the gravitational flux contributions for all phases + const std::array< Vector, numPhases >& gravity() const { return g_; } + std::array< Vector, numPhases >& gravity() { return g_; } + + //! Projection matrix for gravitational acceleration + const Matrix& advectionCA() const { return CA_; } + Matrix& advectionCA() { return CA_; } + + //! Additional projection matrix needed on surface grids + const Matrix& advectionA() const { return A_; } + Matrix& advectionA() { return A_; } + + //! The transmissibilities associated with advective fluxes + const Matrix& advectionT() const { return T_; } + Matrix& advectionT() { return T_; } + + //! The transmissibilities for "outside" faces (used on surface grids) + const std::vector< std::vector<Vector> >& advectionTout() const { return outsideT_; } + std::vector< std::vector<Vector> >& advectionTout() { return outsideT_; } + + //! The gravitational acceleration for "outside" faces (used on surface grids) + const std::array< std::vector<Vector>, numPhases >& gravityOutside() const { return outsideG_; } + std::array< std::vector<Vector>, numPhases >& gravityOutside() { return outsideG_; } + + //! The gravitational acceleration for one phase on "outside" faces (used on surface grids) + const std::vector<Vector>& gravityOutside(unsigned int pIdx) const { return outsideG_[pIdx]; } + std::vector<Vector>& gravityOutside(unsigned int pIdx) { return outsideG_[pIdx]; } + +private: + //! advection-related variables + Matrix T_; //!< The transmissibilities such that f_i = T_ij*p_j + Matrix CA_; //!< Matrix to project gravitational acceleration to all scvfs + Matrix A_; //!< Matrix additionally needed for the projection on surface grids + std::array< Vector, numPhases > p_; //!< The interaction volume-wide phase pressures + std::array< Vector, numPhases > g_; //!< The gravitational acceleration at each scvf (only for enabled gravity) + std::vector< std::vector<Vector> > outsideT_; //!< The transmissibilities for "outside" faces (only on surface grids) + std::array< std::vector<Vector>, numPhases > outsideG_; //!< The gravitational acceleration on "outside" faces (only on surface grids) +}; + +//! Data handle for quantities related to diffusion +template<class TypeTag, class InteractionVolume, bool EnableDiffusion> +class DiffusionDataHandle +{ + //! export matrix & vector types from interaction volume + using Matrix = typename InteractionVolume::Traits::Matrix; + using Vector = typename InteractionVolume::Traits::Vector; + using OutsideTContainer = std::vector< std::vector<Vector> >; - //! return functions for the stored data - const DirichletDataContainer& dirichletData() const { return dirichletData_; } + static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension; + static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld; - const GlobalIndexContainer& volVarsStencil() const - { - if (context_ == Contexts::advection) - return AdvectionHandle::volVarsStencil(); - else if (context_ == Contexts::diffusion) - return DiffusionHandle::volVarsStencil(); - else if (context_ == Contexts::heatConduction) - return HeatConductionHandle::volVarsStencil(); - else - DUNE_THROW(Dune::InvalidStateException, "No valid context set!"); - } + static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); + static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents); - const Matrix& T() const - { - if (context_ == Contexts::advection) - return AdvectionHandle::T(); - else if (context_ == Contexts::diffusion) - return DiffusionHandle::T(); - else if (context_ == Contexts::heatConduction) - return HeatConductionHandle::T(); - else - DUNE_THROW(Dune::InvalidStateException, "No valid context set!"); - } +public: + //! diffusion caches need to set phase and component index + void setPhaseIndex(unsigned int phaseIdx) { phaseIdx_ = phaseIdx; } + void setComponentIndex(unsigned int compIdx) { compIdx_ = compIdx; } - Matrix& T() + //! prepare data handle for subsequent fill + void resize(const InteractionVolume& iv) + { + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) { - if (context_ == Contexts::advection) - return AdvectionHandle::T(); - else if (context_ == Contexts::diffusion) - return DiffusionHandle::T(); - else if (context_ == Contexts::heatConduction) - return HeatConductionHandle::T(); - else - DUNE_THROW(Dune::InvalidStateException, "No valid context set!"); + for (unsigned int cIdx = 0; cIdx < numComponents; ++cIdx) + { + if (pIdx == cIdx) + continue; + + //! resize transmissibility matrix & mole fraction vector + T_[pIdx][cIdx].resize(iv.numFaces(), iv.numKnowns()); + xj_[pIdx][cIdx].resize(iv.numKnowns()); + + //! resize outsideTij on surface grids + if (dim < dimWorld) + { + outsideT_[pIdx][cIdx].resize(iv.numFaces()); + + using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType; + for (LocalIndexType i = 0; i < iv.numFaces(); ++i) + outsideT_[pIdx][cIdx][i].resize(iv.localScvf(i).neighboringLocalScvIndices().size()-1); + } + } } + } + + //! Access to the iv-wide mole fractions of a component in one phase + const Vector& moleFractions(unsigned int pIdx, unsigned int compIdx) const { return xj_[pIdx][compIdx]; } + Vector& moleFractions(unsigned int pIdx, unsigned int compIdx) { return xj_[pIdx][compIdx]; } + + //! The transmissibilities associated with diffusive fluxes + const Matrix& diffusionT() const { return T_[phaseIdx_][compIdx_]; } + Matrix& diffusionT() { return T_[phaseIdx_][compIdx_]; } + + //! The transmissibilities for "outside" faces (used on surface grids) + const OutsideTContainer& diffusionTout() const { return outsideT_[phaseIdx_][compIdx_]; } + OutsideTContainer& diffusionTout() { return outsideT_[phaseIdx_][compIdx_]; } + +private: + //! diffusion-related variables + unsigned int phaseIdx_; //!< The phase index set for the context + unsigned int compIdx_; //!< The component index set for the context + std::array< std::array<Matrix, numComponents>, numPhases > T_; //!< The transmissibilities such that f_i = T_ij*x_j + std::array< std::array<Vector, numComponents>, numPhases > xj_; //!< The interaction volume-wide mole fractions + std::array< std::array<OutsideTContainer, numComponents>, numPhases> outsideT_; +}; + +//! Data handle for quantities related to heat conduction +template<class TypeTag, class InteractionVolume, bool EnableHeatConduction> +class HeatConductionDataHandle +{ + //! export matrix & vector types from interaction volume + using Matrix = typename InteractionVolume::Traits::Matrix; + using Vector = typename InteractionVolume::Traits::Vector; - const Matrix& AB() const - { - if (context_ == Contexts::advection) - return AdvectionHandle::AB(); - else if (context_ == Contexts::diffusion) - return DiffusionHandle::AB(); - else if (context_ == Contexts::heatConduction) - return HeatConductionHandle::AB(); - else - DUNE_THROW(Dune::InvalidStateException, "No valid context set!"); - } + static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension; + static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld; - Matrix& AB() - { - if (context_ == Contexts::advection) - return AdvectionHandle::AB(); - else if (context_ == Contexts::diffusion) - return DiffusionHandle::AB(); - else if (context_ == Contexts::heatConduction) - return HeatConductionHandle::AB(); - else - DUNE_THROW(Dune::InvalidStateException, "No valid context set!"); - } +public: + //! prepare data handle for subsequent fill + void resize(const InteractionVolume& iv) + { + //! resize transmissibility matrix & temperature vector + T_.resize(iv.numFaces(), iv.numKnowns()); + Tj_.resize(iv.numKnowns()); - const Matrix& outsideTij() const + //! resize outsideTij on surface grids + if (dim < dimWorld) { - if (context_ == Contexts::advection) - return AdvectionHandle::outsideTij(); - else if (context_ == Contexts::diffusion) - return DiffusionHandle::outsideTij(); - else if (context_ == Contexts::heatConduction) - return HeatConductionHandle::outsideTij(); - else - DUNE_THROW(Dune::InvalidStateException, "No valid context set!"); - } + outsideT_.resize(iv.numFaces()); - Matrix& outsideTij() - { - if (context_ == Contexts::advection) - return AdvectionHandle::outsideTij(); - else if (context_ == Contexts::diffusion) - return DiffusionHandle::outsideTij(); - else if (context_ == Contexts::heatConduction) - return HeatConductionHandle::outsideTij(); - else - DUNE_THROW(Dune::InvalidStateException, "No valid context set!"); + using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType; + for (LocalIndexType i = 0; i < iv.numFaces(); ++i) + outsideT_[i].resize(iv.localScvf(i).neighboringLocalScvIndices().size()-1); } + } + + //! Access to the iv-wide temperatures + const Vector& temperatures() const { return Tj_; } + Vector& temperatures() { return Tj_; } + + //! The transmissibilities associated with conductive fluxes + const Matrix& heatConductionT() const { return T_; } + Matrix& heatConductionT() { return T_; } + + //! The transmissibilities for "outside" faces (used on surface grids) + const std::vector< std::vector<Vector> >& heatConductionTout() const { return outsideT_; } + std::vector< std::vector<Vector> >& heatConductionTout() { return outsideT_; } + +private: + // heat conduction-related variables + Matrix T_; //!< The transmissibilities such that f_i = T_ij*T_j + Vector Tj_; //!< The interaction volume-wide temperatures + std::vector< std::vector<Vector> > outsideT_; //!< The transmissibilities for "outside" faces (only necessary on surface grids) +}; + +//! 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)> +{ + 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)>; - private: - Contexts context_; //!< The context variable - DirichletDataContainer dirichletData_; //!< The dirichlet data container of this iv - }; +public: + //! prepare data handles for subsequent fills + void resize(const InteractionVolume& iv) + { + AdvectionHandle::resize(iv); + DiffusionHandle::resize(iv); + HeatConductionHandle::resize(iv); + } +}; } // end namespace Dumux diff --git a/dumux/discretization/cellcentered/mpfa/localassembler.hh b/dumux/discretization/cellcentered/mpfa/localassembler.hh new file mode 100644 index 0000000000000000000000000000000000000000..ca270a907bf3aa24254362d843ae800aba135001 --- /dev/null +++ b/dumux/discretization/cellcentered/mpfa/localassembler.hh @@ -0,0 +1,261 @@ +// -*- 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 Defines the general interface of classes used for the assembly + * of the local systems of equations involved in the transmissibility + * computaion in mpfa schemes. + */ +#ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_HH +#define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_HH + +#include <dune/common/exceptions.hh> + +namespace Dumux +{ +//! Forward declaration of the implementation +template< class InteractionVolume > class InteractionVolumeAssemblerImpl; + +//! Alias to select the right implementation. +template< class InteractionVolume > +using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< InteractionVolume >; + +/*! + * \ingroup CCMpfaDiscretization + * \brief Defines the general interface of the local assembler + * classes for the assembly of the interaction volume-local + * transmissibility matrix. Specializations have to be provided + * for the available interaction volume implementations. these + * should derive from this base clases. + * + * \tparam Interaction Volume The interaction volume implementation + * used by the mpfa scheme + */ +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; + + public: + /*! + * \brief The constructor. + * Sets pointers to the objects required for a subsequent call to assemble(). + * + * \param problem The problem to be solved (boundary/initial conditions etc.) + * \param fvGeometry The local view on the finite volume grid geometry + * \param elemVolVars The local view on the primary/secondary variables + */ + InteractionVolumeAssemblerBase(const Problem& problem, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars) + { + problemPtr_ = &problem; + fvGeometryPtr_ = &fvGeometry; + elemVolVarsPtr_ = &elemVolVars; + } + + //! return functions to the local views & problem + const Problem& problem() const { return *problemPtr_; } + const FVElementGeometry& fvGeometry() const { return *fvGeometryPtr_; } + const ElementVolumeVariables& elemVolVars() const { return *elemVolVarsPtr_; } + + /*! + * \brief General interface of a function assembling the + * interaction volume-local transmissibility matrix. + * + * \tparam GetTensorFunction Lambda to obtain the tensor w.r.t. + * which the local system is to be solved + * + * \param T The transmissibility matrix to be assembled + * \param iv The interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class GetTensorFunction > + void assemble(Matrix& T, InteractionVolume& iv, const GetTensorFunction& getTensor) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function"); + } + + /*! + * \brief General interface of a function assembling the interaction + * volume-local transmissibilities matrix for surface grids. The + * transmissibility associated with "outside" faces are stored + * in a separate container. + * + * \param outsideTij tij on "outside" faces to be assembled + * \param T The transmissibility matrix tij to be assembled + * \param iv The mpfa-o interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class OutsideTijContainer, class GetTensorFunction > + void assemble(OutsideTijContainer& outsideTij, Matrix& T, InteractionVolume& iv, const GetTensorFunction& getTensor) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function to be used on surface grids"); + } + + /*! + * \brief General interface of a function assembling the interaction + * volume-local transmissibility matrix in the case that gravity + * is to be considered in the local system of equations. + * + * \tparam GravityContainer The type of container used to store the + * gravitational acceleration per scvf & phase + * \tparam GetTensorFunction Lambda to obtain the tensor w.r.t. + * which the local system is to be solved + * + * \param T The transmissibility matrix to be assembled + * \param g Container to assemble gravity per scvf & phase + * \param CA Matrix to store matrix product C*A^-1 + * \param iv The interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class GravityContainer, class GetTensorFunction > + void assembleWithGravity(Matrix& T, + GravityContainer& g, + Matrix& CA, + InteractionVolume& iv, + const GetTensorFunction& getTensor) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleWithGravity() function"); + } + + /*! + * \brief General interface of a function assembling the interaction + * volume-local transmissibility matrix in the case that gravity + * is to be considered in the local system of equations. This + * specialization is to be used on surface grids, where the + * gravitational flux contributions on "outside" faces are stored + * in a separate container. + * + * \param outsideTij tij on "outside" faces to be assembled + * \param T The transmissibility matrix to be assembled + * \param outsideG Container to assemble gravity on "outside" faces + * \param g Container to assemble gravity per scvf & phase + * \param CA Matrix to store matrix product C*A^-1 + * \param A Matrix to store the inverse A^-1 + * \param iv The interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class GravityContainer, + class OutsideGravityContainer, + class OutsideTijContainer, + class GetTensorFunction > + void assembleWithGravity(OutsideTijContainer& outsideTij, + Matrix& T, + OutsideGravityContainer& outsideG, + GravityContainer& g, + Matrix& CA, + Matrix& A, + InteractionVolume& iv, + const GetTensorFunction& getTensor) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleWithGravity() function to be used on surface grids"); + } + + /*! + * \brief General interface for the assembly of the vector of + * primary (cell) unknowns and (maybe) Dirichlet boundary + * conditions within an interaction volume. + * + * \tparam GetU Lambda to obtain the cell unknowns from grid indices + * + * \param u The vector to be filled with the cell unknowns + * \param iv The interaction volume + * \param getU Lambda to obtain the desired cell value from grid indices + */ + template< class GetU > + void assemble(Vector& u, const InteractionVolume& iv, const GetU& getU) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function for the cell unknowns"); + } + + /*! + * \brief General interface for the assembly of the gravitational + * flux contributions on the scvfs within an interaction volume. + * + * \note For each face, the gravity term in the form of \f$\rho \mathbf{n K g}\f$ is + * evaluated. Thus, make sure to only call this with a lambda that returns the + * hydraulic conductivity. + * + * \param g Container to assemble gravity per scvf & phase + * \param iv The interaction volume + * \param CA Projection matrix transforming the gravity terms in the local system of + * equations to the entire set of faces within the interaction volume + * \param getTensor Lambda to evaluate scv-wise hydraulic conductivities + */ + template< class GravityContainer, class GetTensorFunction > + void assembleGravity(GravityContainer& g, + const InteractionVolume& iv, + const Matrix& CA, + const GetTensorFunction& getTensor) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleGravity() function"); + } + + /*! + * \brief General interface for the assembly of the gravitational + * flux contributions on the scvfs within an interaction volume. + * This specialization is to be used on surface grids, where the gravitational + * flux contributions on "outside" faces are stored in a separate container. + * + * \note For each face, the gravity term in the form of \f$\rho \mathbf{n K g}\f$ is + * evaluated. Thus, make sure to only call this with a lambda that returns the + * hydraulic conductivity. + * + * \param g Container to store gravity per scvf & phase + * \param outsideG Container to store gravity per "outside" scvf & phase + * \param iv The mpfa-o interaction volume + * \param CA Projection matrix transforming the gravity terms in the local system of + * equations to the entire set of faces within the interaction volume + * \param A Matrix needed for the "reconstruction" of face unknowns as a function of gravity + * \param getTensor Lambda to evaluate scv-wise hydraulic conductivities + */ + template< class GravityContainer, + class OutsideGravityContainer, + class GetTensorFunction > + void assembleGravity(GravityContainer& g, + OutsideGravityContainer& outsideG, + const InteractionVolume& iv, + const Matrix& CA, + const Matrix& A, + const GetTensorFunction& getTensor) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleGravity() function to be used on surface grids"); + } + + private: + // pointers to the data required for assembly + const Problem* problemPtr_; + const FVElementGeometry* fvGeometryPtr_; + const ElementVolumeVariables* elemVolVarsPtr_; +}; + +} // end namespace Dumux + +//! include all specializations for different mpfa schemes +#include <dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh> + +#endif diff --git a/dumux/discretization/cellcentered/mpfa/localfacedata.hh b/dumux/discretization/cellcentered/mpfa/localfacedata.hh new file mode 100644 index 0000000000000000000000000000000000000000..1693f214c998d2e0441eb835c187b3b80164952b --- /dev/null +++ b/dumux/discretization/cellcentered/mpfa/localfacedata.hh @@ -0,0 +1,80 @@ +// -*- 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 Data structure holding interaction volume-local information for + * a grid subb-control volume face embedded in it. + */ +#ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_FACE_DATA_HH +#define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_FACE_DATA_HH + +namespace Dumux +{ +/*! + * \ingroup CCMpfaDiscretization + * \brief General implementation of a data structure holding interaction + * volume-local information for a grid subb-control volume face embedded in it. + * + * \tparam GridIndexType The type used for indices on the grid + * \tparam LocalIndexType The type used for indices inside interaction volumes + */ +template< class GridIndexType, class LocalIndexType > +class InteractionVolumeLocalFaceData +{ + LocalIndexType ivLocalScvfIndex_; //!< the iv-local scvf index this scvf maps to + LocalIndexType ivLocalInsideScvIndex_; //!< the iv-local index of the scvfs' inside scv + LocalIndexType scvfLocalOutsideScvfIndex_; //!< the index of this scvf in the scvf-local outside faces + GridIndexType globalScvfIndex_; //!< the index of the corresponding global scvf + bool isOutside_; //!< indicates if this face maps to the iv-local index from "outside" + +public: + //! Constructor + InteractionVolumeLocalFaceData(LocalIndexType faceIndex, + LocalIndexType scvIndex, + GridIndexType globalScvfIndex) + : ivLocalScvfIndex_(faceIndex) + , ivLocalInsideScvIndex_(scvIndex) + , globalScvfIndex_(globalScvfIndex) + , isOutside_(false) + {} + + //! Constructor for "outside" faces + InteractionVolumeLocalFaceData(LocalIndexType faceIndex, + LocalIndexType scvIndex, + LocalIndexType indexInScvfOutsideFaces, + GridIndexType globalScvfIndex) + : ivLocalScvfIndex_(faceIndex) + , ivLocalInsideScvIndex_(scvIndex) + , scvfLocalOutsideScvfIndex_(indexInScvfOutsideFaces) + , globalScvfIndex_(globalScvfIndex) + , isOutside_(true) + {} + + //! Functions to return stored data + LocalIndexType ivLocalScvfIndex() const { return ivLocalScvfIndex_; } + LocalIndexType ivLocalInsideScvIndex() const { return ivLocalInsideScvIndex_; } + LocalIndexType scvfLocalOutsideScvfIndex() const { assert(isOutside_); return scvfLocalOutsideScvfIndex_; } + GridIndexType globalScvfIndex() const { return globalScvfIndex_; } + bool isOutside() const { return isOutside_; } +}; + +} // end namespace + +#endif diff --git a/dumux/discretization/cellcentered/mpfa/methods.hh b/dumux/discretization/cellcentered/mpfa/methods.hh index 0c9c57d24a744c632476fd4b04967b7cdc563484..1f4e3244d0cd384f12badb2fa099f21c03f8a377 100644 --- a/dumux/discretization/cellcentered/mpfa/methods.hh +++ b/dumux/discretization/cellcentered/mpfa/methods.hh @@ -18,17 +18,23 @@ *****************************************************************************/ /*! * \file - * \brief Available implemented mpfa schemes + * \ingroup CCMpfaDiscretization + * \brief The available mpfa schemes in Dumux */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_METHODS_HH #define DUMUX_DISCRETIZATION_CC_MPFA_METHODS_HH namespace Dumux { + /*! + * \brief The available mpfa schemes in Dumux + * \ingroup CCMpfaDiscretization + */ enum class MpfaMethods : unsigned int { oMethod }; -} // end namespace + +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/omethod/CMakeLists.txt b/dumux/discretization/cellcentered/mpfa/omethod/CMakeLists.txt index 8a45cecd5c6bb656f2980953de3803e5caebbe9c..d2c61c2a7b45d6ba942dfb77a7933b45875f6e52 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/CMakeLists.txt +++ b/dumux/discretization/cellcentered/mpfa/omethod/CMakeLists.txt @@ -1,6 +1,6 @@ install(FILES -helper.hh interactionvolume.hh interactionvolumeindexset.hh +localassembler.hh localsubcontrolentities.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/cellcentered/mpfa/omethod) diff --git a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh deleted file mode 100644 index 8c9cb79fe4615088ff75a32e47defce0be008184..0000000000000000000000000000000000000000 --- a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh +++ /dev/null @@ -1,66 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- -// vi: set et ts=4 sw=4 sts=4: -/***************************************************************************** - * See the file COPYING for full copying permissions. * - * * - * This program is free software: you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation, either version 2 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program. If not, see <http://www.gnu.org/licenses/>. * - *****************************************************************************/ - -/*! - * \file - * \brief Helper class for interaction volumes using the mpfa-o scheme. - */ -#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_HELPER_HH -#define DUMUX_DISCRETIZATION_CC_MPFA_O_HELPER_HH - -#include <dumux/common/math.hh> -#include <dumux/discretization/cellcentered/mpfa/methods.hh> - -namespace Dumux -{ - -// forward declaration -template<class TypeTag, MpfaMethods m, int dim, int dimWorld> -class MpfaMethodHelper; - -/*! - * \ingroup Mpfa - * \brief Helper class to get the required information on an interaction volume. - * Specialization for the Mpfa-O method in two dimensions. This scheme does - * not require any additional functionality to be implemented. - */ -template<class TypeTag> -class MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, /*dim*/2, /*dimWorld*/2> {}; - -/*! - * \ingroup Mpfa - * \brief Helper class to get the required information on an interaction volume. - * Specialization for the Mpfa-O method in two dimensions embedded in a 3d world. - * This scheme does not require any additional functionality to be implemented. - */ -template<class TypeTag> -class MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, /*dim*/2, /*dimWorld*/3> {}; - -/*! - * \ingroup Mpfa - * \brief Helper class to get the required information on an interaction volume. - * Specialization for the Mpfa-O method in three dimensions. - * This scheme does not require any additional functionality to be implemented. - */ -template<class TypeTag> -class MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, /*dim*/3, /*dimWorld*/3> {}; - -} // end namespace - -#endif diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh index 02cf2ec7ac09177641e3f43043e2782926727562..7bec27dadd139dbe82c6ce2e7607198be4e944ec 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh @@ -18,128 +18,144 @@ *****************************************************************************/ /*! * \file - * \brief Base classes for interaction volume of mpfa methods. + * \ingroup CCMpfaDiscretization + * \brief Classes for the interaction volume of the mpfa-o scheme. */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH #define DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH +#include <dune/common/dynmatrix.hh> +#include <dune/common/dynvector.hh> +#include <dune/common/fvector.hh> + #include <dumux/common/math.hh> +#include <dumux/common/properties.hh> +#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh> #include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh> -#include <dumux/discretization/cellcentered/mpfa/methods.hh> #include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh> +#include <dumux/discretization/cellcentered/mpfa/localfacedata.hh> +#include <dumux/discretization/cellcentered/mpfa/methods.hh> #include "localsubcontrolentities.hh" #include "interactionvolumeindexset.hh" namespace Dumux { +//! Forward declaration of the interaction volume specialization for the mpfa-o scheme +template< class TypeTag > +class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>; + //! Specialization of the interaction volume traits class for the mpfa-o method -template<class TypeTag> -class CCMpfaOInteractionVolumeTraits : public CCMpfaInteractionVolumeTraitsBase<TypeTag> +template< class TypeTag > +struct CCMpfaInteractionVolumeTraits< CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> > { - using BaseTraits = CCMpfaInteractionVolumeTraitsBase<TypeTag>; - using NodalIndexSet = typename CCMpfaDualGridIndexSet<TypeTag>::NodalIndexSet; - -public: - using SecondaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>; +private: + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using typename BaseTraits::DynamicLocalIndexContainer; - using typename BaseTraits::DynamicGlobalIndexContainer; - using IndexSet = CCMpfaOInteractionVolumeIndexSet<NodalIndexSet, DynamicGlobalIndexContainer, DynamicLocalIndexContainer>; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; - // In the o-scheme, matrix & vector types are dynamic types - using typename BaseTraits::DynamicVector; - using typename BaseTraits::DynamicMatrix; - using Vector = DynamicVector; - using Matrix = DynamicMatrix; + using InteractionVolumeType = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>; - using LocalScvType = CCMpfaOInteractionVolumeLocalScv<TypeTag, IndexSet>; - using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf<TypeTag>; +public: + //! Per default, we use the same interaction volume everywhere (also on boundaries etc...) + using SecondaryInteractionVolume = InteractionVolumeType; + + //! 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 >; + //! export the type of interaction-volume local scvfs + using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf< IndexSet >; + //! export the type of used for the iv-local face data + using LocalFaceData = InteractionVolumeLocalFaceData<GridIndexType, LocalIndexType>; + //! 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 >; + //! export the type used for iv-local vectors + using Vector = Dune::DynamicVector< Scalar >; + //! export the data handle type for this iv + using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >; }; -//! Forward declaration of the mpfa-o interaction volume -template<class TypeTag, class Traits> -class CCMpfaOInteractionVolume; - /*! - * \ingroup Mpfa - * \brief Base class for the interaction volumes of the mpfa-o method. - * We introduce one more level of inheritance here because the o-method with - * full pressure support uses the mpfa-o interaction volume but uses a different - * traits class. + * \ingroup CCMpfaDiscretization + * \brief Class for the interaction volume of the mpfa-o method. */ template<class TypeTag> -class CCMpfaInteractionVolumeImplementation< TypeTag, MpfaMethods::oMethod> - : public CCMpfaOInteractionVolume< TypeTag, CCMpfaOInteractionVolumeTraits<TypeTag> > +class CCMpfaInteractionVolumeImplementation< TypeTag, MpfaMethods::oMethod > + : public CCMpfaInteractionVolumeBase< TypeTag, + CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> > { - using TraitsType = CCMpfaOInteractionVolumeTraits<TypeTag>; -public: - // state the traits class type - using Traits = TraitsType; -}; + using ThisType = CCMpfaInteractionVolumeImplementation< TypeTag, MpfaMethods::oMethod >; + using ParentType = CCMpfaInteractionVolumeBase< TypeTag, ThisType >; -template<class TypeTag, class Traits> -class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Traits> -{ - using ParentType = CCMpfaInteractionVolumeBase<TypeTag, Traits>; - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); 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 FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); - using DualGridNodalIndexSet = typename CCMpfaDualGridIndexSet<TypeTag>::NodalIndexSet; - - static const int dim = GridView::dimension; - static const int dimWorld = GridView::dimensionworld; using Element = typename GridView::template Codim<0>::Entity; + + static constexpr int dim = GridView::dimension; using DimVector = Dune::FieldVector<Scalar, dim>; - using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; - using Vector = typename Traits::DynamicVector; - using Matrix = typename Traits::DynamicMatrix; - using Tensor = typename Traits::Tensor; + 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 LocalScvType = typename Traits::LocalScvType; - using LocalScvfType = typename Traits::LocalScvfType; + using IndexSet = typename TraitsType::IndexSet; + using GridIndexType = typename TraitsType::GridIndexType; + using LocalIndexType = typename TraitsType::LocalIndexType; - using IndexSet = typename Traits::IndexSet; - using LocalIndexContainer = typename Traits::DynamicLocalIndexContainer; - using LocalIndexType = typename LocalIndexContainer::value_type; - using GlobalIndexContainer = typename Traits::DynamicGlobalIndexContainer; - using DataHandle = typename Traits::DataHandle; -public: + //! Data attached to scvf touching Dirichlet boundaries. + //! For the default o-scheme, we only store the corresponding vol vars index. + class DirichletData + { + GridIndexType volVarIndex_; + + public: + //! Constructor + DirichletData(const GridIndexType index) : volVarIndex_(index) {} + + //! Return corresponding vol var index + GridIndexType volVarIndex() const { return volVarIndex_; } + }; +public: //! publicly state the mpfa-scheme this interaction volume is associated with static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod; - using typename ParentType::LocalFaceData; - using typename ParentType::DirichletDataContainer; - using typename ParentType::LocalFaceDataContainer; + //! 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); - //! Sets up the local scope for a given iv index set! + //! Sets up the local scope for a given iv index set void setUpLocalScope(const IndexSet& indexSet, const Problem& problem, const FVElementGeometry& fvGeometry) { - //! store a pointer to the index set - indexSetPtr_ = &indexSet; - - //! clear previous data - clear_(); - - //! number of interaction-volume-local faces + // number of interaction-volume-local (= node-local for o-scheme) scvs/scvf numFaces_ = indexSet.numFaces(); - - //! number of interaction-volume-local (and node-local) scvs - const auto& scvIndices = indexSet.globalScvIndices(); const auto numLocalScvs = indexSet.numScvs(); - - //! number of global scvfs appearing in this interaction volume const auto numGlobalScvfs = indexSet.nodalIndexSet().numScvfs(); - //! reserve memory for local entities + // reserve memory for local entities + elements_.clear(); + scvs_.clear(); + scvfs_.clear(); + localFaceData_.clear(); + dirichletData_.clear(); elements_.reserve(numLocalScvs); scvs_.reserve(numLocalScvs); scvfs_.reserve(numFaces_); @@ -147,142 +163,101 @@ public: localFaceData_.reserve(numGlobalScvfs); // set up quantities related to sub-control volumes + const auto& scvIndices = indexSet.globalScvIndices(); for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numLocalScvs; scvIdxLocal++) { - const auto scvIdxGlobal = scvIndices[scvIdxLocal]; - scvs_.emplace_back(fvGeometry, fvGeometry.scv(scvIdxGlobal), scvIdxLocal, indexSet); - elements_.emplace_back(fvGeometry.fvGridGeometry().element(scvIdxGlobal)); + scvs_.emplace_back(Helper(), + fvGeometry, + fvGeometry.scv( scvIndices[scvIdxLocal] ), + scvIdxLocal, + indexSet); + elements_.emplace_back(fvGeometry.fvGridGeometry().element( scvIndices[scvIdxLocal] )); } // keep track of the number of unknowns etc numUnknowns_ = 0; - numOutsideFaces_ = 0; - numPotentials_ = numLocalScvs; + numKnowns_ = numLocalScvs; + + // resize omega storage container + wijk_.resize(numFaces_); // set up quantitites related to sub-control volume faces for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numFaces_; ++faceIdxLocal) { - const auto scvfIdxGlobal = indexSet.scvfIdxGlobal(faceIdxLocal); - const auto& neighborScvIndicesLocal = indexSet.neighboringLocalScvIndices(faceIdxLocal); - const auto insideLocalScvIdx = neighborScvIndicesLocal[0]; + // get corresponding grid scvf + const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(faceIdxLocal)); - // we have to use the "inside" scv face here - const auto& scvf = fvGeometry.scvf(scvfIdxGlobal); + // the neighboring scvs in local indices (order: 0 - inside scv, 1..n - outside scvs) + const auto& neighborScvIndicesLocal = indexSet.neighboringLocalScvIndices(faceIdxLocal); // create local face data object for this face - localFaceData_.emplace_back(faceIdxLocal, insideLocalScvIdx, scvf.index()); + localFaceData_.emplace_back(faceIdxLocal, neighborScvIndicesLocal[0], scvf.index()); // create iv-local scvf object if (scvf.boundary()) { - const auto insideElement = elements_[insideLocalScvIdx]; - const auto bcTypes = problem.boundaryTypes(insideElement, scvf); + const auto bcTypes = problem.boundaryTypes(elements_[neighborScvIndicesLocal[0]], scvf); if (bcTypes.hasOnlyDirichlet()) { - scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/true, numPotentials_++); - dirichletData_.emplace_back(scvf.outsideScvIdx(), scvf.ipGlobal()); + scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numKnowns_++, /*isDirichlet*/true); + dirichletData_.emplace_back(scvf.outsideScvIdx()); } else - scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/false, numUnknowns_++); + scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numUnknowns_++, /*isDirichlet*/false); + + // on boundary faces we will only need one inside omega + wijk_[faceIdxLocal].resize(1); } else { - scvfs_.emplace_back(scvf, neighborScvIndicesLocal, /*isDirichlet*/false, numUnknowns_++); + scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numUnknowns_++, /*isDirichlet*/false); - // add local face data object for the outside faces - for (LocalIndexType i = 1; i < neighborScvIndicesLocal.size(); ++i) - { - const auto outsideLocalScvIdx = neighborScvIndicesLocal[i]; + // we will need as many omegas as scvs around the face + const auto numNeighborScvs = neighborScvIndicesLocal.size(); + wijk_[faceIdxLocal].resize(numNeighborScvs); + // add local face data objects for the outside faces + for (LocalIndexType i = 1; i < numNeighborScvs; ++i) + { // loop over scvfs in outside scv until we find the one coinciding with current scvf + const auto outsideLocalScvIdx = neighborScvIndicesLocal[i]; for (int coord = 0; coord < dim; ++coord) { if (indexSet.scvfIdxLocal(outsideLocalScvIdx, coord) == faceIdxLocal) { const auto globalScvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord); const auto& flipScvf = fvGeometry.scvf(globalScvfIdx); - localFaceData_.emplace_back(faceIdxLocal, //! iv-local scvf idx - outsideLocalScvIdx, //! iv-local scv index - numOutsideFaces_++, //! iv-local index in outside faces - i-1, //! scvf-local index in outside faces - flipScvf.index()); //!< global scvf index + localFaceData_.emplace_back(faceIdxLocal, // iv-local scvf idx + outsideLocalScvIdx, // iv-local scv index + i-1, // scvf-local index in outside faces + flipScvf.index()); // global scvf index } } } } } - // resize the local matrices + // Resize local matrices A_.resize(numUnknowns_, numUnknowns_); - B_.resize(numUnknowns_, numPotentials_); + B_.resize(numUnknowns_, numKnowns_); C_.resize(numFaces_, numUnknowns_); - D_.resize(numFaces_, numPotentials_); } - //! sets the sizes of the corresponding matrices in the data handle - void prepareDataHandle(DataHandle& dataHandle) - { - // resize the transmissibility matrix in the data handle - dataHandle.resizeT(numFaces_, numPotentials_); + //! returns the number of primary scvfs of this interaction volume + std::size_t numFaces() const { return numFaces_; } - // move the dirichlet data into the handle - dataHandle.setDirichletData(std::move(dirichletData_)); - - // resize possible additional containers in the data handle - if (requireABMatrix_()) - dataHandle.resizeAB(numUnknowns_, numPotentials_); - if (dim < dimWorld) - dataHandle.resizeOutsideTij(numOutsideFaces_, numPotentials_); - } - - //! solves for the transmissibilities subject to a given tensor - template<typename GetTensorFunction> - void solveLocalSystem(const GetTensorFunction& getTensor, - const Problem& problem, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - DataHandle& dataHandle) - { - // if only dirichlet faces are present, assemble T_ directly - if (numUnknowns_ == 0) - assemblePureDirichletT_(getTensor, problem, fvGeometry, elemVolVars, dataHandle.T()); - else - { - // assemble - assembleLocalMatrices_(getTensor, problem, fvGeometry, elemVolVars); - - // solve - A_.invert(); - - // T = C*A^-1*B + D - dataHandle.T() = multiplyMatrices(C_.rightmultiply(A_), B_); - dataHandle.T() += D_; - - // store A-1B only when gradient reconstruction is necessary - if (requireABMatrix_()) - dataHandle.AB() = B_.leftmultiply(A_); - } + //! returns the number of intermediate unknowns within this interaction volume + std::size_t numUnknowns() const { return numUnknowns_; } - // set vol vars stencil pointer in handle - dataHandle.setVolVarsStencilPointer(indexSet().globalScvIndices()); + //! returns the number of (in this context) known solution values within this interaction volume + std::size_t numKnowns() const { return numKnowns_; } - // on surface grids, additionally prepare the outside transmissibilities - if (dim < dimWorld) - computeOutsideTransmissibilities_(dataHandle); - } + //! returns the number of scvs embedded in this interaction volume + std::size_t numScvs() const { return scvs_.size(); } - //! obtain the local data object for a given global scvf - const LocalFaceData& getLocalFaceData(const SubControlVolumeFace& scvf) const - { - //! find corresponding entry in the local face data container - const auto scvfIdxGlobal = scvf.index(); - auto it = std::find_if(localFaceData_.begin(), - localFaceData_.end(), - [scvfIdxGlobal] (const LocalFaceData& d) { return d.globalScvfIndex() == scvfIdxGlobal; }); - assert(it != localFaceData_.end() && "Could not find the local face data corresponding to the given scvf"); - return localFaceData_[std::distance(localFaceData_.begin(), it)]; - } + //! returns the number of scvfs embedded in this interaction volume + std::size_t numScvfs() const { return scvfs_.size(); } //! returns the grid element corresponding to a given iv-local scv idx const Element& element(const LocalIndexType ivLocalScvIdx) const { return elements_[ivLocalScvIdx]; } @@ -296,19 +271,33 @@ public: //! returns a reference to the container with the local face data const std::vector<LocalFaceData>& localFaceData() const { return localFaceData_; } - //! returns a reference to the index set of this iv - const IndexSet& indexSet() const { return *indexSetPtr_; } + //! returns a reference to the information container on Dirichlet BCs within this iv + const std::vector<DirichletData>& dirichletData() const { return dirichletData_; } + + //! return functions for matrices involved in local system of equations + const Matrix& A() const { return A_; } + Matrix& A() { return A_; } + + const Matrix& B() const { return B_; } + Matrix& B() { return B_; } + + const Matrix& C() const { return C_; } + Matrix& C() { return C_; } + + const std::vector< std::vector< DimVector > >& omegas() const { return wijk_; } + std::vector< std::vector< DimVector > >& omegas() { return wijk_; } //! returns the number of interaction volumes living around a vertex //! the mpfa-o scheme always constructs one iv per vertex - static std::size_t numInteractionVolumesAtVertex(const DualGridNodalIndexSet& nodalIndexSet) { return 1; } + template< class NodalIndexSet > + static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) { return 1; } //! 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 - template<class IvIndexSetContainer, class ScvfIndexMap> + template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet> static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer, ScvfIndexMap& scvfIndexMap, - const DualGridNodalIndexSet& nodalIndexSet) + const NodalIndexSet& nodalIndexSet) { // the global index of the iv index set that is about to be created const auto curGlobalIndex = ivIndexSetContainer.size(); @@ -322,301 +311,25 @@ public: } private: - //! returns a boolean whether or not the AB matrix has to be passed to the handles - static bool requireABMatrix_() - { - static const bool requireAB = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Vtk.AddVelocity") || dim < dimWorld; - return requireAB; - } - - //! clears all the containers - void clear_() - { - elements_.clear(); - scvs_.clear(); - scvfs_.clear(); - localFaceData_.clear(); - dirichletData_.clear(); - } - - //! Assembles the local matrices that define the local system of equations and flux expressions - template<typename GetTensorFunction> - void assembleLocalMatrices_(const GetTensorFunction& getTensor, - const Problem& problem, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars) - { - // reset matrices - A_ = 0.0; - B_ = 0.0; - C_ = 0.0; - D_ = 0.0; - - // reserve space for the omegas - wijk_.resize(scvfs_.size()); - - // loop over the local faces - for (LocalIndexType faceIdx = 0; faceIdx < numFaces_; ++faceIdx) - { - const auto& curLocalScvf = localScvf(faceIdx); - const auto& curGlobalScvf = fvGeometry.scvf(curLocalScvf.globalScvfIndex()); - const auto curIsDirichlet = curLocalScvf.isDirichlet(); - const auto curLocalDofIdx = curLocalScvf.localDofIndex(); - - // get diffusion tensor in "positive" sub volume - const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); - const auto posLocalScvIdx = neighborScvIndices[0]; - const auto& posLocalScv = localScv(posLocalScvIdx); - const auto& posGlobalScv = fvGeometry.scv(posLocalScv.globalScvIndex()); - const auto& posVolVars = elemVolVars[posGlobalScv]; - const auto& posElement = element(posLocalScvIdx); - const auto tensor = getTensor(problem, posElement, posVolVars, fvGeometry, posGlobalScv); - - // the omega factors of the "positive" sub volume - auto posWijk = calculateOmegas_(posLocalScv, curGlobalScvf.unitOuterNormal(), curGlobalScvf.area(), tensor); - posWijk *= posVolVars.extrusionFactor(); - - // go over the coordinate directions in the positive sub volume - for (unsigned int localDir = 0; localDir < dim; localDir++) - { - const auto otherLocalScvfIdx = posLocalScv.scvfIdxLocal(localDir); - const auto& otherLocalScvf = localScvf(otherLocalScvfIdx); - const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); - - // if we are not on a Dirichlet face, add entries associated with unknown face pressures - // i.e. in matrix C and maybe A (if current face is not a Dirichlet face) - if (!otherLocalScvf.isDirichlet()) - { - C_[faceIdx][otherLocalDofIdx] -= posWijk[localDir]; - if (!curIsDirichlet) - A_[curLocalDofIdx][otherLocalDofIdx] -= posWijk[localDir]; - } - // the current face is a Dirichlet face and creates entries in D & maybe B - else - { - D_[faceIdx][otherLocalDofIdx] -= posWijk[localDir]; - if (!curIsDirichlet) - B_[curLocalDofIdx][otherLocalDofIdx] += posWijk[localDir]; - } - - // add entries related to pressures at the scv centers (dofs) - const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); - D_[faceIdx][posScvLocalDofIdx] += posWijk[localDir]; - - if (!curIsDirichlet) - B_[curLocalDofIdx][posScvLocalDofIdx] -= posWijk[localDir]; - } - - // store the omegas - wijk_[faceIdx].emplace_back(std::move(posWijk)); - - // If we are on an interior face, add values from negative sub volume - if (!curGlobalScvf.boundary()) - { - // loop over all the outside neighbors of this face and add entries - for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) - { - const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1]; - const auto& negLocalScv = localScv(negLocalScvIdx); - const auto& negGlobalScv = fvGeometry.scv(negLocalScv.globalScvIndex()); - const auto& negVolVars = elemVolVars[negGlobalScv]; - const auto& negElement = element(negLocalScvIdx); - const auto negTensor = getTensor(problem, negElement, negVolVars, fvGeometry, negGlobalScv); - - // the omega factors of the "negative" sub volume - DimVector negWijk; - - // if dim < dimWorld, use outside normal vector - if (dim < dimWorld) - { - const auto& flipScvf = fvGeometry.flipScvf(curGlobalScvf.index(), idxInOutside); - auto negNormal = flipScvf.unitOuterNormal(); - negNormal *= -1.0; - negWijk = calculateOmegas_(negLocalScv, negNormal, curGlobalScvf.area(), negTensor); - } - else - negWijk = calculateOmegas_(negLocalScv, curGlobalScvf.unitOuterNormal(), curGlobalScvf.area(), negTensor); - - // scale by extrusion factpr - negWijk *= negVolVars.extrusionFactor(); - - // go over the coordinate directions in the positive sub volume - for (int localDir = 0; localDir < dim; localDir++) - { - const auto otherLocalScvfIdx = negLocalScv.scvfIdxLocal(localDir); - const auto& otherLocalScvf = localScvf(otherLocalScvfIdx); - const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); - - if (!otherLocalScvf.isDirichlet()) - A_[curLocalDofIdx][otherLocalDofIdx] += negWijk[localDir]; - else - B_[curLocalDofIdx][otherLocalDofIdx] -= negWijk[localDir]; - - // add entries to matrix B - B_[curLocalDofIdx][negLocalScv.localDofIndex()] += negWijk[localDir]; - } - - // store the omegas - wijk_[faceIdx].emplace_back(std::move(negWijk)); - } - } - } - } - - //! for interaction volumes that have only dirichlet scvfs, - //! the transmissibility matrix can be assembled directly - template<typename GetTensorFunction> - void assemblePureDirichletT_(const GetTensorFunction& getTensor, - const Problem& problem, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - Matrix& T) - { - // reset the transmissibility matrix beforehand - T = 0.0; - - // Loop over all the faces, in this case these are all dirichlet boundaries - for (unsigned int faceIdx = 0; faceIdx < numFaces_; ++faceIdx) - { - const auto& curLocalScvf = localScvf(faceIdx); - const auto& curGlobalScvf = fvGeometry.scvf(curLocalScvf.globalScvfIndex()); - - // get diffusion tensor in "positive" sub volume - const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); - const auto posLocalScvIdx = neighborScvIndices[0]; - const auto& posLocalScv = localScv(posLocalScvIdx); - const auto& posGlobalScv = fvGeometry.scv(posLocalScv.globalScvIndex()); - const auto& posVolVars = elemVolVars[posGlobalScv]; - const auto& posElement = element(posLocalScvIdx); - const auto tensor = getTensor(problem, posElement, posVolVars, fvGeometry, posGlobalScv); - - // the omega factors of the "positive" sub volume - auto posWijk = calculateOmegas_(posLocalScv, curGlobalScvf.unitOuterNormal(), curGlobalScvf.area(), tensor); - posWijk *= posVolVars.extrusionFactor(); - - const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); - for (LocalIndexType localDir = 0; localDir < dim; localDir++) - { - const auto otherLocalScvfIdx = posLocalScv.scvfIdxLocal(localDir); - const auto& otherLocalScvf = localScvf(otherLocalScvfIdx); - const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); - T[faceIdx][otherLocalDofIdx] -= posWijk[localDir]; - T[faceIdx][posScvLocalDofIdx] += posWijk[localDir]; - } - } - } - - //! computes the transmissibilities associated with "outside" faces on surface grids - void computeOutsideTransmissibilities_(DataHandle& dataHandle) const - { - assert(dim < dimWorld && "only for dim < dimWorld the outside transmissiblity container has the right size"); - - for (const auto& localFaceData : localFaceData_) - { - //! continue only for "outside" faces - if (!localFaceData.isOutside()) continue; - - const auto localScvIdx = localFaceData.ivLocalInsideScvIndex(); - const auto localScvfIdx = localFaceData.ivLocalScvfIndex(); - const auto& posLocalScv = localScv(localScvIdx); - const auto& wijk = wijk_[localScvfIdx][localFaceData.scvfLocalOutsideScvfIndex() + 1]; - - //! store the calculated transmissibilities in the data handle - auto& tij = dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()]; - - //! reset transmissibility vector - tij = 0.0; - - //! add contributions from all local directions - for (LocalIndexType localDir = 0; localDir < dim; localDir++) - { - //! the scvf corresponding to this local direction in the scv - const auto& curLocalScvf = localScvf(posLocalScv.scvfIdxLocal(localDir)); - - //! on interior faces the coefficients of the AB matrix come into play - if (!curLocalScvf.isDirichlet()) - { - auto tmp = dataHandle.AB()[curLocalScvf.localDofIndex()]; - tmp *= wijk[localDir]; - tij -= tmp; - } - else - tij[curLocalScvf.localDofIndex()] -= wijk[localDir]; - - // add entry from the scv unknown - tij[localScvIdx] += wijk[localDir]; - } - } - } - - // calculates n_i^T*K_j*nu_k - DimVector calculateOmegas_(const LocalScvType& localScv, - const GlobalPosition& normal, - const Scalar area, - const Tensor& T) const - { - // make sure we have positive definite diffsion tensors - assert(this->tensorIsPositiveDefinite(T) && "only positive definite tensors can be handled by mpfa methods"); - - DimVector wijk; - GlobalPosition tmp; - for (LocalIndexType dir = 0; dir < dim; ++dir) - { - T.mv(localScv.innerNormal(dir), tmp); - wijk[dir] = tmp*normal; - } - wijk *= area; - wijk /= localScv.detX(); - - return wijk; - } - - // calculates n_i^T*K_j*nu_k - DimVector calculateOmegas_(const LocalScvType& localScv, - const GlobalPosition& normal, - const Scalar area, - const Scalar t) const - { - // make sure we have positive diffusion coefficients - assert(t > 0.0 && "non-positive diffusion coefficients cannot be handled by mpfa methods"); - - DimVector wijk; - GlobalPosition tmp(normal); - tmp *= t; - - for (LocalIndexType dir = 0; dir < dim; ++dir) - wijk[dir] = tmp*localScv.innerNormal(dir); - wijk *= area; - wijk /= localScv.detX(); - - return wijk; - } - - const IndexSet* indexSetPtr_; - // Variables defining the local scope std::vector<Element> elements_; std::vector<LocalScvType> scvs_; std::vector<LocalScvfType> scvfs_; std::vector<LocalFaceData> localFaceData_; - DirichletDataContainer dirichletData_; + std::vector<DirichletData> dirichletData_; - // sizes involved in the local matrices - unsigned int numFaces_; - unsigned int numUnknowns_; - unsigned int numPotentials_; - unsigned int numOutsideFaces_; - - // The omega factors and the matrix Aâ»1*B are stored - // in order to recover the transmissibilities of outside faces on network grids - std::vector< std::vector< DimVector > > wijk_; - Matrix CAinv_; - - // Matrices involved in transmissibility calculations + // Matrices needed for computation of transmissibilities Matrix A_; Matrix B_; Matrix C_; - Matrix D_; + + // The omega factors are stored during assemble of local system + std::vector< std::vector< DimVector > > wijk_; + + // sizes involved in the local system equations + std::size_t numFaces_; + std::size_t numUnknowns_; + std::size_t numKnowns_; }; } // end namespace diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh index 2645b214ad35139af3c0ec4cbf17a164edfb311a..84a26f7a4e4a253f905b1ca6405f6b1e2a2bdd65 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh @@ -18,6 +18,7 @@ *****************************************************************************/ /*! * \file + * \ingroup CCMpfaDiscretization * \brief Class for the index set within an interaction volume of the mpfa-o scheme. */ #ifndef DUMUX_DISCRETIZATION_MPFA_O_INTERACTIONVOLUME_INDEXSET_HH @@ -28,20 +29,42 @@ namespace Dumux { /*! - * \ingroup Mpfa + * \ingroup CCMpfaDiscretization * \brief The interaction volume index set class for the mpfa-o scheme. + * + * \tparam DualGridNodalIndexSet The type used for the nodal index set in the dual grid. */ -template<class DualGridNodalIndexSet, class GlobalIndexContainer, class LocalIndexContainer> +template< class DualGridNodalIndexSet > class CCMpfaOInteractionVolumeIndexSet { - using LocalIndexType = typename LocalIndexContainer::value_type; - using GlobalIndexType = typename GlobalIndexContainer::value_type; - public: - explicit CCMpfaOInteractionVolumeIndexSet(const DualGridNodalIndexSet& nodalIndexSet) : nodalIndexSet_(nodalIndexSet) + using LocalIndexType = typename DualGridNodalIndexSet::LocalIndexType; + using GridIndexType = typename DualGridNodalIndexSet::GridIndexType; + + using LocalIndexContainer = typename DualGridNodalIndexSet::LocalIndexContainer; + using GridIndexContainer = typename DualGridNodalIndexSet::GridIndexContainer; + + /*! + * \brief The constructor + * \note The actual type used for the nodal index sets might be different, as maybe + * a different type for the local indexes is used. We therefore template this + * constructor. However, a static assertion enforces you to use the same LocalIndexType + * in the traits for both the secondary and the primary interaction volume traits. + * + * \tparam NodalIndexSet Possibly differing type for the DualGridNodalIndexSet + */ + template< class NodalIndexSet > + CCMpfaOInteractionVolumeIndexSet(const NodalIndexSet& nodalIndexSet) + : nodalIndexSet_( static_cast<const DualGridNodalIndexSet&>(nodalIndexSet) ) { - //! determine the number of iv-local faces for memory reservation - //! note that this might be a vast overestimation on surface grids! + // make sure the index types used are the same in order to avoid any losses due to type conversion + static_assert(std::is_same<GridIndexType, typename NodalIndexSet::GridIndexType>::value, + "Provided nodal index set does not use the same type for grid indices as the given template argument"); + static_assert(std::is_same<LocalIndexType, typename NodalIndexSet::LocalIndexType>::value, + "Provided nodal index set does not use the same type for local indices as the given template argument"); + + // determine the number of iv-local faces for memory reservation + // note that this might be a vast overestimation on surface grids! const auto numNodalScvfs = nodalIndexSet.numScvfs(); const auto numBoundaryScvfs = nodalIndexSet.numBoundaryScvfs(); const std::size_t numFaceEstimate = numBoundaryScvfs + (numNodalScvfs-numBoundaryScvfs)/2; @@ -49,26 +72,26 @@ public: // make sure we found a reasonable number of faces assert((numNodalScvfs-numBoundaryScvfs)%2 == 0); - //! index transformation from interaction-volume-local to node-local + // index transformation from interaction-volume-local to node-local ivToNodeScvf_.reserve(numFaceEstimate); nodeToIvScvf_.resize(numNodalScvfs); - //! the local neighboring scv indices of the faces + // the local neighboring scv indices of the faces scvfNeighborScvLocalIndices_.reserve(numFaceEstimate); - //! keeps track of which nodal scvfs have been handled already + // keeps track of which nodal scvfs have been handled already std::vector<bool> isHandled(numNodalScvfs, false); - //! go over faces in nodal index set, check if iv-local face has been - //! inserted already for this scvf and if not, insert index mapping + // go over faces in nodal index set, check if iv-local face has been + // inserted already for this scvf and if not, insert index mapping numFaces_ = 0; for (LocalIndexType i = 0; i < numNodalScvfs; ++i) { - //! check if the nodal scvf still has to be handled + // check if the nodal scvf still has to be handled if (isHandled[i]) continue; - //! for scvfs touching the boundary there are no "outside" scvfs + // for scvfs touching the boundary there are no "outside" scvfs if (nodalIndexSet.scvfIsOnBoundary(i)) { nodeToIvScvf_[i] = ivToNodeScvf_.size(); @@ -78,8 +101,8 @@ public: continue; } - //! We insert a new iv-local face and find all "outside" scvfs that map - //! to this face as well by comparing the set of neighboring scv indices. + // We insert a new iv-local face and find all "outside" scvfs that map + // to this face as well by comparing the set of neighboring scv indices. const auto scvIndices = [&nodalIndexSet, i] () { auto tmp = nodalIndexSet.neighboringScvIndices(i); @@ -91,7 +114,7 @@ public: std::vector<LocalIndexType> outsideScvfs; for (LocalIndexType j = i+1; j < numNodalScvfs; ++j) { - //! a face that has been handled already cannot be an "outside" face here + // a face that has been handled already cannot be an "outside" face here if (!isHandled[j]) { const auto scvIndices2 = [&nodalIndexSet, j] () @@ -125,8 +148,8 @@ public: } // compute local neighboring scv indices for the iv-local scvfs - scvfNeighborScvLocalIndices_.resize(numFaces()); - for (unsigned int i = 0; i < numFaces(); ++i) + scvfNeighborScvLocalIndices_.resize(numFaces_); + for (unsigned int i = 0; i < numFaces_; ++i) { const auto& neighborsGlobal = nodalIndexSet_.neighboringScvIndices(ivToNodeScvf_[i]); const auto numNeighbors = nodalIndexSet_.scvfIsOnBoundary(ivToNodeScvf_[i]) ? 1 : neighborsGlobal.size(); @@ -141,16 +164,12 @@ public: const DualGridNodalIndexSet& nodalIndexSet() const { return nodalIndexSet_; } //! returns a global scvf idx for a given iv_local scvf index - GlobalIndexType scvfIdxGlobal(LocalIndexType ivLocalScvfIdx) const - { return nodalIndexSet_.scvfIdxGlobal(ivToNodeScvf_[ivLocalScvfIdx]); } + GridIndexType scvfIdxGlobal(LocalIndexType ivLocalScvfIdx) const + { return nodalIndexSet_.scvfIdxGlobal( ivToNodeScvf_[ivLocalScvfIdx] ); } //! returns the iv-local scvf idx of the i-th scvfs embedded in a local scv LocalIndexType scvfIdxLocal(LocalIndexType scvIdxLocal, unsigned int i) const - { return nodeToIvScvf_[nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i)]; } - - //! returns whether or not an scvf touches the boundary - bool scvfTouchesBoundary(LocalIndexType ivLocalScvfIdx) const - { return nodalIndexSet_.scvfTouchesBoundary(ivToNodeScvf_[ivLocalScvfIdx]); } + { return nodeToIvScvf_[ nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) ]; } //! returns the local indices of the neighboring scvs of an scvf const LocalIndexContainer& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const @@ -163,14 +182,14 @@ public: std::size_t numScvs() const { return nodalIndexSet_.numScvs(); } //! returns the global scv indices connected to this dual grid node - const GlobalIndexContainer& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); } + const GridIndexContainer& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); } //! returns the global scvf indices connected to this dual grid node - const GlobalIndexContainer& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); } + const GridIndexContainer& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); } private: //! returns the local scv index to a given global scv index - unsigned int findLocalScvIdx_(GlobalIndexType globalScvIdx) const + unsigned int findLocalScvIdx_(GridIndexType globalScvIdx) const { auto it = std::find( nodalIndexSet_.globalScvIndices().begin(), nodalIndexSet_.globalScvIndices().end(), globalScvIdx ); assert(it != nodalIndexSet_.globalScvIndices().end() && "Global scv index not found in local container!"); @@ -180,13 +199,13 @@ private: const DualGridNodalIndexSet& nodalIndexSet_; std::size_t numFaces_; - LocalIndexContainer ivToNodeScvf_; - LocalIndexContainer nodeToIvScvf_; - //! maps to each scvf a list of neighbouring scv indices - //! ordering: 0 - inside scv idx; 1..n - outside scv indices - std::vector<LocalIndexContainer> scvfNeighborScvLocalIndices_; + std::vector<LocalIndexType> ivToNodeScvf_; + std::vector<LocalIndexType> nodeToIvScvf_; + // maps to each scvf a list of neighbouring scv indices + // ordering: 0 - inside scv idx; 1..n - outside scv indices + std::vector< LocalIndexContainer > scvfNeighborScvLocalIndices_; }; -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh new file mode 100644 index 0000000000000000000000000000000000000000..3b245375437f940c52fd20a053aea9cdeeb2c534 --- /dev/null +++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh @@ -0,0 +1,740 @@ +// -*- 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 for the assembly of the local systems of equations + * involved in the transmissibility computaion in the mpfa-o scheme. + */ +#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH +#define DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH + +#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> + +namespace Dumux +{ + +/*! + * \ingroup CCMpfaDiscretization + * \brief Specialization of the interaction volume-local + * assembler class for the mpfa-o scheme. + */ +template< class TypeTag > +class InteractionVolumeAssemblerImpl< CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> > + : public InteractionVolumeAssemblerBase< CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod> > +{ + using InteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethod>; + using ParentType = InteractionVolumeAssemblerBase< InteractionVolume >; + + using Traits = typename InteractionVolume::Traits; + using LocalIndexType = typename Traits::LocalIndexType; + using Matrix = typename Traits::Matrix; + using Vector = typename Traits::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); + +public: + //! Use the constructor of the base class + using ParentType::ParentType; + + /*! + * \brief Assembles the transmissibility matrix + * within an interaction volume for the mpfa-o scheme. + * + * \param T The transmissibility matrix to be assembled + * \param iv The mpfa-o interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class GetTensorFunction > + void assemble(Matrix& T, InteractionVolume& iv, const GetTensorFunction& getTensor) + { + //! assemble D into T directly + assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor); + + //! maybe solve the local system + if (iv.numUnknowns() > 0) + { + //! T = C*A^-1*B + D + iv.A().invert(); + iv.C().rightmultiply(iv.A()); + T += multiplyMatrices(iv.C(), iv.B()); + } + } + + /*! + * \brief Assembles the interaction volume-local transmissibility + * matrix for surface grids. The transmissibilities associated + * with "outside" faces are stored in a separate container. + * + * \param outsideTij tij on "outside" faces to be assembled + * \param T The transmissibility matrix tij to be assembled + * \param iv The mpfa-o interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class OutsideTijContainer, class GetTensorFunction > + void assemble(OutsideTijContainer& outsideTij, Matrix& T, InteractionVolume& iv, const GetTensorFunction& getTensor) + { + //! assemble D into T directly + assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor); + + //! maybe solve the local system + if (iv.numUnknowns() > 0) + { + // T = C*A^-1*B + D + iv.A().invert(); + iv.B().leftmultiply(iv.A()); + T += multiplyMatrices(iv.C(), iv.B()); + + //! compute outside transmissibilities + for (const auto& localFaceData : iv.localFaceData()) + { + // continue only for "outside" faces + if (!localFaceData.isOutside()) + continue; + + const auto localScvIdx = localFaceData.ivLocalInsideScvIndex(); + const auto localScvfIdx = localFaceData.ivLocalScvfIndex(); + const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex(); + const auto& posLocalScv = iv.localScv(localScvIdx); + const auto& wijk = iv.omegas()[localScvfIdx][idxInOutside+1]; + + // store the calculated transmissibilities in the data handle + auto& tij = outsideTij[localScvfIdx][idxInOutside]; + tij = 0.0; + + // add contributions from all local directions + for (LocalIndexType localDir = 0; localDir < dim; localDir++) + { + // the scvf corresponding to this local direction in the scv + const auto& curLocalScvf = iv.localScvf(posLocalScv.scvfIdxLocal(localDir)); + + // on interior faces the coefficients of the AB matrix come into play + if (!curLocalScvf.isDirichlet()) + { + auto tmp = iv.B()[curLocalScvf.localDofIndex()]; + tmp *= wijk[localDir]; + tij -= tmp; + } + else + tij[curLocalScvf.localDofIndex()] -= wijk[localDir]; + + // add entry from the scv unknown + tij[localScvIdx] += wijk[localDir]; + } + } + } + } + + /*! + * \brief Assemble the transmissibility matrix within an interaction + * volume for the mpfa-o scheme to be used for advective flux + * computation in the case that gravity is to be considered in + * the local system of equations. + * + * \param T The transmissibility matrix to be assembled + * \param g Container to assemble gravity per scvf & phase + * \param CA Matrix to store matrix product C*A^-1 + * \param iv The mpfa-o interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class GravityContainer, class GetTensorFunction > + void assembleWithGravity(Matrix& T, + GravityContainer& g, + Matrix& CA, + InteractionVolume& iv, + const GetTensorFunction& getTensor) + { + //! assemble D into T & C into CA directly + assembleLocalMatrices_(iv.A(), iv.B(), CA, T, iv, getTensor); + + //! maybe solve the local system + if (iv.numUnknowns() > 0) + { + //! T = C*A^-1*B + D + iv.A().invert(); + CA.rightmultiply(iv.A()); + T += multiplyMatrices(CA, iv.B()); + } + + //! assemble gravitational acceleration container (enforce usage of mpfa-o type version) + assembleGravity(g, iv, CA, getTensor); + } + + /*! + * \brief Assembles the interaction volume-local transmissibility + * matrix in the case that gravity is to be considered in the + * local system of equations. This specialization is to be used + * on surface grids, where the gravitational flux contributions + * on "outside" faces are stored in a separate container. + * + * \param outsideTij tij on "outside" faces to be assembled + * \param T The transmissibility matrix to be assembled + * \param outsideG Container to assemble gravity on "outside" faces + * \param g Container to assemble gravity per scvf & phase + * \param CA Matrix to store matrix product C*A^-1 + * \param A Matrix to store the inverse A^-1 + * \param iv The mpfa-o interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class GravityContainer, + class OutsideGravityContainer, + class OutsideTijContainer, + class GetTensorFunction > + void assembleWithGravity(OutsideTijContainer& outsideTij, + Matrix& T, + OutsideGravityContainer& outsideG, + GravityContainer& g, + Matrix& CA, + Matrix& A, + InteractionVolume& iv, + const GetTensorFunction& getTensor) + { + //! assemble D into T directly + assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor); + + //! maybe solve the local system + if (iv.numUnknowns() > 0) + { + // T = C*A^-1*B + D + iv.A().invert(); + iv.B().leftmultiply(iv.A()); + T += multiplyMatrices(iv.C(), iv.B()); + A = iv.A(); + CA = iv.C().rightmultiply(A); + + //! compute outside transmissibilities + for (const auto& localFaceData : iv.localFaceData()) + { + // continue only for "outside" faces + if (!localFaceData.isOutside()) + continue; + + const auto localScvIdx = localFaceData.ivLocalInsideScvIndex(); + const auto localScvfIdx = localFaceData.ivLocalScvfIndex(); + const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex(); + const auto& posLocalScv = iv.localScv(localScvIdx); + const auto& wijk = iv.omegas()[localScvfIdx][idxInOutside+1]; + + // store the calculated transmissibilities in the data handle + auto& tij = outsideTij[localScvfIdx][idxInOutside]; + tij = 0.0; + + // add contributions from all local directions + for (LocalIndexType localDir = 0; localDir < dim; localDir++) + { + // the scvf corresponding to this local direction in the scv + const auto& curLocalScvf = iv.localScvf(posLocalScv.scvfIdxLocal(localDir)); + + // on interior faces the coefficients of the AB matrix come into play + if (!curLocalScvf.isDirichlet()) + { + auto tmp = iv.B()[curLocalScvf.localDofIndex()]; + tmp *= wijk[localDir]; + tij -= tmp; + } + else + tij[curLocalScvf.localDofIndex()] -= wijk[localDir]; + + // add entry from the scv unknown + tij[localScvIdx] += wijk[localDir]; + } + } + } + + assembleGravity(g, outsideG, iv, CA, A, getTensor); + } + + /*! + * \brief Assembles the vector of primary (cell) unknowns and (maybe) + * Dirichlet boundary conditions within an interaction volume. + * + * \param u The vector to be filled with the cell unknowns + * \param iv The mpfa-o interaction volume + * \param getU Lambda to obtain the desired cell/Dirichlet value from grid index + */ + template< class GetU > + void assemble(Vector& u, const InteractionVolume& iv, const GetU& getU) + { + //! resize given container + u.resize(iv.numKnowns()); + + //! put the cell pressures first + + for (LocalIndexType i = 0; i < iv.numScvs(); ++i) + u[i] = getU( iv.localScv(i).globalScvIndex() ); + + //! Dirichlet BCs come afterwards + unsigned int i = iv.numScvs(); + for (const auto& data : iv.dirichletData()) + u[i++] = getU( data.volVarIndex() ); + } + + /*! + * \brief Assemble the gravitational flux contributions on the scvfs within + * an mpfa-o interaction volume. + * + * \note For each face, the gravity term in the form of \f$\rho \mathbf{n K g}\f$ is + * evaluated. Thus, make sure to only call this with a lambda that returns the + * hydraulic conductivity. + * + * \param g Container to assemble gravity per scvf & phase + * \param iv The mpfa-o interaction volume + * \param CA Projection matrix transforming the gravity terms in the local system of + * equations to the entire set of faces within the interaction volume + * \param getTensor Lambda to evaluate scv-wise hydraulic conductivities + */ + template< class GravityContainer, class GetTensorFunction > + void assembleGravity(GravityContainer& g, + const InteractionVolume& iv, + 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"); + assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns() && "Matrix CA does not have the correct size"); + + //! For each face, we... + //! - arithmetically average the phase densities + //! - compute the term alpha:= A*rho*n^T*K*g in each neighboring cell + //! - compute sum_alphas = alpha_outside - alpha_inside + 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); + } + + for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) + { + //! gravitational acceleration on this face + const auto& curLocalScvf = iv.localScvf(faceIdx); + const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.globalScvfIndex()); + const auto gravity = this->problem().gravityAtPos(curGlobalScvf.ipGlobal()); + + //! get permeability tensor in "positive" sub volume + const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); + const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); + const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.globalScvIndex()); + const auto& posVolVars = this->elemVolVars()[posGlobalScv]; + const auto& posElement = iv.element(neighborScvIndices[0]); + const auto tensor = getTensor(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv); + + //! This function should never be called for surface grids, + //! 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; + const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), tensor, gravity); + if (!curLocalScvf.isDirichlet()) + { + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + rho[pIdx] = posVolVars.density(pIdx); + + if (!curGlobalScvf.boundary()) + { + //! obtain outside tensor + const auto& negLocalScv = iv.localScv( neighborScvIndices[1] ); + const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.globalScvIndex()); + const auto& negVolVars = this->elemVolVars()[negGlobalScv]; + const auto& negElement = iv.element( neighborScvIndices[1] ); + const auto negTensor = getTensor(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv); + + const auto sum_alpha = negVolVars.extrusionFactor() + * vtmv(curGlobalScvf.unitOuterNormal(), negTensor, gravity) + - alpha_inside; + + const auto localDofIdx = curLocalScvf.localDofIndex(); + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + { + rho[pIdx] = 0.5*( rho[pIdx] + negVolVars.density(pIdx) ); + sum_alphas[pIdx][localDofIdx] = sum_alpha*rho[pIdx]*curGlobalScvf.area(); + } + } + else + { + const auto localDofIdx = curLocalScvf.localDofIndex(); + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + sum_alphas[pIdx][localDofIdx] -= alpha_inside*rho[pIdx]*curGlobalScvf.area(); + } + } + //! use Dirichlet BC densities + else + { + const auto& dirichletVolVars = this->elemVolVars()[curGlobalScvf.outsideScvIdx()]; + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + rho[pIdx] = dirichletVolVars.density(pIdx); + } + + //! add "inside" alpha to gravity container + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + g[pIdx][faceIdx] += alpha_inside*rho[pIdx]*curGlobalScvf.area(); + } + + //! g += CA*sum_alphas + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + CA.umv(sum_alphas[pIdx], g[pIdx]); + } + + /*! + * \brief Assembles the gravitational flux contributions on the scvfs within an mpfa-o + * interaction volume. This specialization is to be used on surface grids, where the + * gravitational flux contributions on "outside" faces are stored in a separate container. + * + * \note For each face, the gravity term in the form of \f$\rho \mathbf{n K g}\f$ is + * evaluated. Thus, make sure to only call this with a lambda that returns the + * hydraulic conductivity. + * + * \param g Container to store gravity per scvf & phase + * \param outsideG Container to store gravity per "outside" scvf & phase + * \param iv The mpfa-o interaction volume + * \param CA Projection matrix transforming the gravity terms in the local system of + * equations to the entire set of faces within the interaction volume + * \param A Matrix needed for the "reconstruction" of face unknowns as a function of gravity + * \param getTensor Lambda to evaluate scv-wise hydraulic conductivities + */ + template< class GravityContainer, + class OutsideGravityContainer, + class GetTensorFunction > + void assembleGravity(GravityContainer& g, + OutsideGravityContainer& outsideG, + const InteractionVolume& iv, + const Matrix& CA, + const Matrix& A, + const GetTensorFunction& getTensor) + { + //! 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"); + + //! For each face, we... + //! - arithmetically average the phase densities + //! - compute the term alpha:= A*rho*n^T*K*g in each neighboring cell + //! - compute sum_alphas = alpha_outside - alpha_inside + using Scalar = typename Vector::value_type; + + // reset everything to zero + std::array< Vector, numPhases > sum_alphas; + 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) + { + //! gravitational acceleration on this face + const auto& curLocalScvf = iv.localScvf(faceIdx); + const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.globalScvfIndex()); + const auto gravity = this->problem().gravityAtPos(curGlobalScvf.ipGlobal()); + + //! get permeability tensor in "positive" sub volume + const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); + const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); + const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.globalScvIndex()); + const auto& posVolVars = this->elemVolVars()[posGlobalScv]; + const auto& posElement = iv.element(neighborScvIndices[0]); + const auto tensor = getTensor(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv); + + 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; + + if (!curLocalScvf.isDirichlet()) + { + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + rho[pIdx] = posVolVars.density(pIdx); + + // arithmetically average density on inside faces + const auto localDofIdx = curLocalScvf.localDofIndex(); + if (!curGlobalScvf.boundary()) + { + for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) + { + //! obtain outside tensor + const auto& negLocalScv = iv.localScv( neighborScvIndices[idxInOutside] ); + const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.globalScvIndex()); + const auto& negVolVars = this->elemVolVars()[negGlobalScv]; + const auto& negElement = iv.element( neighborScvIndices[idxInOutside] ); + const auto negTensor = getTensor(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv); + + const auto& flipScvf = this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside); + alpha_outside[idxInOutside] = negVolVars.extrusionFactor() + * vtmv(flipScvf.unitOuterNormal(), negTensor, gravity); + + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + { + rho[pIdx] += negVolVars.density(pIdx); + sum_alphas[pIdx][localDofIdx] -= alpha_outside[idxInOutside]; + } + } + } + + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + { + rho[pIdx] /= numOutsideFaces + 1; + sum_alphas[pIdx][localDofIdx] -= alpha_inside; + sum_alphas[pIdx][localDofIdx] *= rho[pIdx]*curGlobalScvf.area(); + } + } + //! use Dirichlet BC densities + else + { + const auto& dirichletVolVars = this->elemVolVars()[curGlobalScvf.outsideScvIdx()]; + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + rho[pIdx] = dirichletVolVars.density(pIdx); + } + + //! add "inside" & "outside" alphas to gravity containers + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + { + g[pIdx][faceIdx] += alpha_inside*rho[pIdx]*curGlobalScvf.area(); + unsigned int i = 0; + for (const auto& alpha : alpha_outside) + outsideG[pIdx][faceIdx][i++] -= alpha*rho[pIdx]*curGlobalScvf.area(); + } + } + + //! g += CA*sum_alphas + //! outsideG = wikj*A^-1*sum_alphas + outsideG + for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) + { + CA.umv(sum_alphas[pIdx], g[pIdx]); + + Vector AG(iv.numUnknowns()); + A.mv(sum_alphas[pIdx], AG); + + //! compute gravitational accelerations + for (const auto& localFaceData : iv.localFaceData()) + { + // continue only for "outside" faces + if (!localFaceData.isOutside()) + continue; + + const auto localScvIdx = localFaceData.ivLocalInsideScvIndex(); + const auto localScvfIdx = localFaceData.ivLocalScvfIndex(); + const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex(); + const auto& posLocalScv = iv.localScv(localScvIdx); + const auto& wijk = iv.omegas()[localScvfIdx][idxInOutside+1]; + + // add contributions from all local directions + for (LocalIndexType localDir = 0; localDir < dim; localDir++) + { + // the scvf corresponding to this local direction in the scv + const auto& curLocalScvf = iv.localScvf(posLocalScv.scvfIdxLocal(localDir)); + + // on interior faces the coefficients of the AB matrix come into play + if (!curLocalScvf.isDirichlet()) + outsideG[pIdx][localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()]; + } + } + } + } + +private: + /*! + * \brief Assemble the matrices involved in the flux expressions + * across the scvfs inside an interaction volume as well as those involved + * in the interaction volume-local system of equations resulting from flux + * and solution continuity across the scvfs. + * + * Flux expressions: + * \f$\mathbf{f} = \mathbf{C} \bar{matbf{u}} + \mathbf{D} matbf{u}\f$. + * Continuity equations + * \f$\mathbf{A} \, \bar{matbf{u}} = \mathbf{B} \, matbf{u}\f$, + * + * \note The matrices are expected to have been resized beforehand. + * + * \param A The A matrix of the iv-local equation system + * \param B The B matrix of the iv-local equation system + * \param C The C matrix of the iv-local flux expressions + * \param D The D matrix of the iv-local flux expressions + * \param iv The mpfa-o interaction volume + * \param getTensor Lambda to evaluate the scv-wise tensors + */ + template< class GetTensorFunction > + void assembleLocalMatrices_(Matrix& A, Matrix& B, Matrix& C, Matrix& D, + InteractionVolume& iv, + const GetTensorFunction& getTensor) + { + //! Matrix D is assumed to have the right size already + assert(D.rows() == iv.numFaces() && D.cols() == iv.numKnowns() && "Matrix D does not have the correct size"); + + //! if only Dirichlet faces are present in the iv, + //! the matrices A, B & C are undefined and D = T + if (iv.numUnknowns() == 0) + { + //! reset matrix beforehand + D = 0.0; + + //! Loop over all the faces, in this case these are all dirichlet boundaries + for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) + { + const auto& curLocalScvf = iv.localScvf(faceIdx); + const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.globalScvfIndex()); + const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); + + //! get tensor in "positive" sub volume + const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); + const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.globalScvIndex()); + const auto& posVolVars = this->elemVolVars()[posGlobalScv]; + const auto& posElement = iv.element(neighborScvIndices[0]); + const auto tensor = getTensor(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv); + + //! the omega factors of the "positive" sub volume + const auto wijk = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); + + const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); + for (LocalIndexType localDir = 0; localDir < dim; localDir++) + { + const auto& otherLocalScvf = iv.localScvf( posLocalScv.scvfIdxLocal(localDir) ); + const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); + D[faceIdx][otherLocalDofIdx] -= wijk[localDir]; + D[faceIdx][posScvLocalDofIdx] += wijk[localDir]; + } + } + } + else + { + //! we require the matrices A,B,C to have the correct size already + assert(A.rows() == iv.numUnknowns() && A.cols() == iv.numUnknowns() && "Matrix A does not have the correct size"); + assert(B.rows() == iv.numUnknowns() && B.cols() == iv.numKnowns() && "Matrix B does not have the correct size"); + assert(C.rows() == iv.numFaces() && C.cols() == iv.numKnowns() && "Matrix C does not have the correct size"); + + //! reset matrices + A = 0.0; + B = 0.0; + C = 0.0; + D = 0.0; + + auto& wijk = iv.omegas(); + for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) + { + const auto& curLocalScvf = iv.localScvf(faceIdx); + const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.globalScvfIndex()); + const auto curIsDirichlet = curLocalScvf.isDirichlet(); + const auto curLocalDofIdx = curLocalScvf.localDofIndex(); + + //! get tensor in "positive" sub volume + const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); + const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); + const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.globalScvIndex()); + const auto& posVolVars = this->elemVolVars()[posGlobalScv]; + const auto& posElement = iv.element(neighborScvIndices[0]); + const auto tensor = getTensor(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv); + + //! the omega factors of the "positive" sub volume + wijk[faceIdx][0] = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); + + //! go over the coordinate directions in the positive sub volume + for (unsigned int localDir = 0; localDir < dim; localDir++) + { + const auto& otherLocalScvf = iv.localScvf( posLocalScv.scvfIdxLocal(localDir) ); + const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); + + //! if we are not on a Dirichlet face, add entries associated with unknown face pressures + //! i.e. in matrix C and maybe A (if current face is not a Dirichlet face) + if (!otherLocalScvf.isDirichlet()) + { + C[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; + if (!curIsDirichlet) + A[curLocalDofIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; + } + //! the current face is a Dirichlet face and creates entries in D & maybe B + else + { + D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; + if (!curIsDirichlet) + B[curLocalDofIdx][otherLocalDofIdx] += wijk[faceIdx][0][localDir]; + } + + //! add entries related to pressures at the scv centers (dofs) + const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); + D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; + + if (!curIsDirichlet) + B[curLocalDofIdx][posScvLocalDofIdx] -= wijk[faceIdx][0][localDir]; + } + + // If we are on an interior face, add values from negative sub volume + if (!curGlobalScvf.boundary()) + { + // loop over all the outside neighbors of this face and add entries + for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) + { + const auto idxOnScvf = idxInOutside+1; + const auto& negLocalScv = iv.localScv( neighborScvIndices[idxOnScvf] ); + const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.globalScvIndex()); + const auto& negVolVars = this->elemVolVars()[negGlobalScv]; + const auto& negElement = iv.element( neighborScvIndices[idxOnScvf] ); + const auto negTensor = getTensor(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv); + + // On surface grids, use outside face for "negative" transmissibility calculation + const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside) + : curGlobalScvf; + wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility(negLocalScv, scvf, negTensor, negVolVars.extrusionFactor()); + + // flip sign on surface grids (since we used the "outside" normal) + if (dim < dimWorld) + wijk[faceIdx][idxOnScvf] *= -1.0; + + // go over the coordinate directions in the positive sub volume + for (int localDir = 0; localDir < dim; localDir++) + { + const auto otherLocalScvfIdx = negLocalScv.scvfIdxLocal(localDir); + const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx); + const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); + + if (!otherLocalScvf.isDirichlet()) + A[curLocalDofIdx][otherLocalDofIdx] += wijk[faceIdx][idxOnScvf][localDir]; + else + B[curLocalDofIdx][otherLocalDofIdx] -= wijk[faceIdx][idxOnScvf][localDir]; + + // add entries to matrix B + B[curLocalDofIdx][negLocalScv.localDofIndex()] += wijk[faceIdx][idxOnScvf][localDir]; + } + } + } + } + } + } +}; + +} // end namespace + +#endif diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh index 4336c1017b798cf9da41f26f7a378a841305ca88..68432b87fc5856213e3d75f3e4241d0b38a1989e 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh @@ -18,50 +18,67 @@ *****************************************************************************/ /*! * \file - * \brief Base class for sub control entities of the mpfa-o method. + * \ingroup CCMpfaDiscretization + * \brief Classes for sub control entities of the mpfa-o method. */ -#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_LOCALSUBCONTROLENTITIES_HH -#define DUMUX_DISCRETIZATION_CC_MPFA_O_LOCALSUBCONTROLENTITIES_HH +#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_SUBCONTROLENTITIES_HH +#define DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_SUBCONTROLENTITIES_HH #include <dune/common/fvector.hh> #include <dumux/common/properties.hh> namespace Dumux { -template<class TypeTag, class IvIndexSet> + +/*! + * \ingroup CCMpfaDiscretization + * \brief Class for the interaction volume-local sub-control volume used + * 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 + */ +template< class IvIndexSet, class GC, int dim > class CCMpfaOInteractionVolumeLocalScv { - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using GlobalIndexType = typename GridView::IndexSet::IndexType; - using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); - - //! The mpfa-o method always uses the dynamic types - using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType; - using LocalBasis = typename InteractionVolume::Traits::ScvBasis; - - static const int dim = GridView::dimension; - static const int dimWorld = GridView::dimensionworld; - using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; public: - explicit CCMpfaOInteractionVolumeLocalScv(const FVElementGeometry& fvGeometry, - const SubControlVolume& scv, - const LocalIndexType localIndex, - const IvIndexSet& indexSet) - : indexSet_(indexSet) - , globalScvIndex_(scv.dofIndex()) - , localDofIndex_(localIndex) + //! publicly state some types + using GridIndexType = typename IvIndexSet::GridIndexType; + using LocalIndexType = typename IvIndexSet::LocalIndexType; + using GlobalCoordinate = GC; + using ctype = typename GlobalCoordinate::value_type; + using LocalBasis = std::array< GlobalCoordinate, dim >; + + static constexpr int myDimension = dim; + static constexpr int worldDimension = GlobalCoordinate::dimension; + + /*! + * \brief The constructor + * + * \param helper Helper class for mpfa schemes + * \param fvGeometry The element finite volume geometry + * \param scv The grid sub-control volume + * \param localIndex The iv-local index of this scvIdx + * \param indexSet The interaction volume index set + */ + template<class MpfaHelper, class FVElementGeometry, class SubControlVolume> + CCMpfaOInteractionVolumeLocalScv(const MpfaHelper& helper, + const FVElementGeometry& fvGeometry, + const SubControlVolume& scv, + const LocalIndexType localIndex, + const IvIndexSet& indexSet) + : indexSet_(indexSet) + , globalScvIndex_(scv.dofIndex()) + , localDofIndex_(localIndex) { // center of the global scv - const auto center = scv.center(); + const auto& center = scv.center(); // set up local basis LocalBasis localBasis; - for (unsigned int coordIdx = 0; coordIdx < dim; ++coordIdx) + for (unsigned int coordIdx = 0; coordIdx < myDimension; ++coordIdx) { const auto scvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(localDofIndex_, coordIdx); const auto& scvf = fvGeometry.scvf(scvfIdx); @@ -69,67 +86,83 @@ public: localBasis[coordIdx] -= center; } - innerNormals_ = Helper::calculateInnerNormals(localBasis); - detX_ = Helper::calculateDetX(localBasis); + nus_ = helper.calculateInnerNormals(localBasis); + detX_ = helper.calculateDetX(localBasis); } //! detX is needed for setting up the omegas in the interaction volumes - Scalar detX() const { return detX_; } - - //! the inner normals are needed for setting up the omegas in the interaction volumes - const GlobalPosition& innerNormal(const LocalIndexType coordDir) const - { - assert(coordDir < dim); - return innerNormals_[coordDir]; - } + ctype detX() const { return detX_; } //! grid view-global index related to this scv - GlobalIndexType globalScvIndex() const { return globalScvIndex_; } + GridIndexType globalScvIndex() const { return globalScvIndex_; } //! returns the index in the set of cell unknowns of the iv LocalIndexType localDofIndex() const { return localDofIndex_; } //! iv-local index of the coordir's scvf in this scv - LocalIndexType scvfIdxLocal(const unsigned int coordDir) const + LocalIndexType scvfIdxLocal(unsigned int coordDir) const { - assert(coordDir < dim); + assert(coordDir < myDimension); return indexSet_.scvfIdxLocal(localDofIndex_, coordDir); } + //! the nu vectors are needed for setting up the omegas of the iv + const GlobalCoordinate& nu(unsigned int coordDir) const + { + assert(coordDir < myDimension); + return nus_[coordDir]; + } + private: const IvIndexSet& indexSet_; - GlobalIndexType globalScvIndex_; + GridIndexType globalScvIndex_; LocalIndexType localDofIndex_; - LocalBasis innerNormals_; - Scalar detX_; + LocalBasis nus_; + ctype detX_; }; - -template<class TypeTag> +/*! + * \ingroup CCMpfaDiscretization + * \brief Class for the interaction volume-local sub-control volume face + * used in the mpfa-o scheme. + * + * \tparam IvIndexSet The type used for index sets within interaction volumes + */ +template< class IvIndexSet > struct CCMpfaOInteractionVolumeLocalScvf { - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - - //! The mpfa-o method always uses the dynamic types - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using GlobalIndexType = typename GridView::IndexSet::IndexType; - using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume); - using LocalIndexContainer = typename InteractionVolume::Traits::DynamicLocalIndexContainer; - using LocalIndexType = typename LocalIndexContainer::value_type; + using LocalIndexContainer = typename IvIndexSet::LocalIndexContainer; public: + //! export index types + using GridIndexType = typename IvIndexSet::GridIndexType; + using LocalIndexType = typename IvIndexSet::LocalIndexType; + + /*! + * \brief The constructor + * + * \param scvf The grid sub-control volume face + * \param localScvIndices The iv-local neighboring scv indices + * \param localDofIdx This scvf's interaction volume-local dof index + * \param isDirichlet Specifies if this scv is on a Dirichlet boundary + */ + template< class SubControlVolumeFace > CCMpfaOInteractionVolumeLocalScvf(const SubControlVolumeFace& scvf, const LocalIndexContainer& localScvIndices, - const bool isDirichlet, - const LocalIndexType localDofIdx) - : scvfIdxGlobal_(scvf.index()) - , neighborScvIndicesLocal_(localScvIndices) - , isDirichlet_(isDirichlet) + const LocalIndexType localDofIdx, + const bool isDirichlet) + : isDirichlet_(isDirichlet) + , scvfIdxGlobal_(scvf.index()) , localDofIndex_(localDofIdx) + , neighborScvIndicesLocal_(localScvIndices) {} + //! This is either the iv-local index of the intermediate unknown (interior/Neumann face) + //! or the index of the Dirichlet boundary within the vol vars (Dirichlet faces) + LocalIndexType localDofIndex() const { return localDofIndex_; } + //! returns the grid view-global index of this scvf - GlobalIndexType globalScvfIndex() const { return scvfIdxGlobal_; } + GridIndexType globalScvfIndex() const { return scvfIdxGlobal_; } //! Returns the local indices of the scvs neighboring this scvf const LocalIndexContainer& neighboringLocalScvIndices() const { return neighborScvIndicesLocal_; } @@ -137,17 +170,13 @@ public: //! states if this is scvf is on a Dirichlet boundary bool isDirichlet() const { return isDirichlet_; } - //! This is either the iv-local index of the intermediate unknown (interior/Neumann face) - //! or the index of the Dirichlet boundary within the vol vars (Dirichlet faces) - LocalIndexType localDofIndex() const { return localDofIndex_; } - private: - GlobalIndexType scvfIdxGlobal_; - const LocalIndexContainer& neighborScvIndicesLocal_; - bool isDirichlet_; + GridIndexType scvfIdxGlobal_; LocalIndexType localDofIndex_; + const LocalIndexContainer& neighborScvIndicesLocal_; }; -} // end namespace + +} // end namespace Dumux #endif diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh index 78f21d0b7a378efa509191c803e67b2715649801..f306b74653c3fce598c07726c7bf08cd074f55e6 100644 --- a/dumux/discretization/cellcentered/mpfa/properties.hh +++ b/dumux/discretization/cellcentered/mpfa/properties.hh @@ -17,13 +17,11 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. * *****************************************************************************/ /*! - * \ingroup Properties * \file - * - * \brief Defines a type tag and some properties for models using - * a cell-centered scheme with multi-point flux approximation. + * \ingroup CCMpfaDiscretization + * \brief Properties for all models using cell-centered finite volume scheme with mpfa + * \note Inherit from these properties to use a cell-centered finite volume scheme with mpfa */ - #ifndef DUMUX_CC_MPFA_PROPERTIES_HH #define DUMUX_CC_MPFA_PROPERTIES_HH diff --git a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh index 9830e8991de30d28f037fd15c10658f834d15cf2..7d9d112b3226999c848f8b7c73e9907f62b6db74 100644 --- a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh +++ b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh @@ -18,7 +18,8 @@ *****************************************************************************/ /*! * \file - * \brief Class for the sub-control volume face in mpfa schemes + * \ingroup CCMpfaDiscretization + * \brief The sub control volume face */ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_SUBCONTROLVOLUMEFACE_HH #define DUMUX_DISCRETIZATION_CC_MPFA_SUBCONTROLVOLUMEFACE_HH @@ -32,7 +33,7 @@ namespace Dumux { /*! - * \ingroup Mpfa + * \ingroup CCMpfaDiscretization * \brief Default implementation of the class for a sub-control volume face in mpfa methods. */ template<class ScvfGeometryTraits> @@ -51,7 +52,7 @@ public: /*! * \brief Constructor * - * \param geomHelper The mpfa geometry helper + * \param helper The helper class for mpfa schemes * \param corners The corners of the scv face * \param unitOuterNormal The unit outer normal vector of the scvf * \param vIdxGlobal The global vertex index the scvf is connected to @@ -63,35 +64,35 @@ public: * \param boundary Boolean to specify whether or not the scvf is on a boundary */ template<class MpfaHelper> - explicit 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) - : boundary_(boundary), - vertexIndex_(vIdxGlobal), - scvfIndex_(scvfIndex), - insideScvIdx_(insideScvIdx), - outsideScvIndices_(outsideScvIndices), - vIdxInElement_(vIdxLocal), - corners_(std::move(corners)), - center_(0.0), - unitOuterNormal_(std::move(unitOuterNormal)) - { - // compute the center of the scvf - for (const auto& corner : corners_) - center_ += corner; - center_ /= corners_.size(); - - // use helper class to obtain area & integration point - ipGlobal_ = helper.getScvfIntegrationPoint(corners_, q); - area_ = helper.getScvfArea(corners_); - } + 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) + : boundary_(boundary) + , vertexIndex_(vIdxGlobal) + , scvfIndex_(scvfIndex) + , insideScvIdx_(insideScvIdx) + , outsideScvIndices_(outsideScvIndices) + , vIdxInElement_(vIdxLocal) + , corners_(std::move(corners)) + , center_(0.0) + , unitOuterNormal_(std::move(unitOuterNormal)) + { + // compute the center of the scvf + for (const auto& corner : corners_) + center_ += corner; + center_ /= corners_.size(); + + // use helper class to obtain area & integration point + ipGlobal_ = helper.getScvfIntegrationPoint(corners_, q); + area_ = helper.getScvfArea(corners_); + } //! The area of the sub control volume face Scalar area() const { return area_; } @@ -169,7 +170,7 @@ private: }; /*! - * \ingroup Mpfa + * \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, diff --git a/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh b/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh index 30a39cbf6d1881cdf3e32652e6b76b55c020a278..ea0682ee25ee90c9eece9363985ed4f26d7e060b 100644 --- a/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh +++ b/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh @@ -18,6 +18,7 @@ *****************************************************************************/ /*! * \file + * \ingroup CCMpfaDiscretization * \brief Helper class to be used to obtain lambda functions for the tensors * involved in the laws that describe the different kind of fluxes that * occur in DuMuX models (i.e. advective, diffusive and heat conduction fluxes). @@ -37,7 +38,7 @@ namespace Dumux { /*! - * \ingroup MpfaModel + * \ingroup CCMpfaDiscretization * \brief Helper class to be used to obtain lambda functions for the tensors * involved in the laws that describe the different kind of fluxes that * occur in DuMuX models (i.e. advective, diffusive and heat conduction fluxes). @@ -55,8 +56,9 @@ public: //! We return zero scalars here in the functions below. //! We have to return something as the local systems expect a type - //! to perform actions on. We return 0.0 as this call should never happen - //! for a tensor which is not treated by an mpfa method anyway. + //! to perform actions on, thus a compiler error will occur. + //! We return 0.0 as this call should never happen for a tensor + //! which is not treated by an mpfa method anyway. //! lambda for the law describing the advective term static auto getAdvectionLambda()