diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt index e7d5753277724e689b5578ffdf0f26dd363267b6..abda2600c744aa263385ddf6150a83b753082577 100644 --- a/dumux/common/CMakeLists.txt +++ b/dumux/common/CMakeLists.txt @@ -20,7 +20,6 @@ intersectionmapper.hh intrange.hh loggingparametertree.hh math.hh -matrixvectorhelper.hh numericdifferentiation.hh optional.hh parameters.hh diff --git a/dumux/common/matrixvectorhelper.hh b/dumux/common/matrixvectorhelper.hh deleted file mode 100644 index 44769ee23f09a9b64508bec4eb4ea8d4b14c874b..0000000000000000000000000000000000000000 --- a/dumux/common/matrixvectorhelper.hh +++ /dev/null @@ -1,97 +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 . * - *****************************************************************************/ -/*! - * \file - * \ingroup Common - * \brief Helper functions to be used for memory allocation operations on matrices & vectors - * that can be called for both static and dynamic types. This is useful wherever it - * is not always known if the matrix & vector types at hand are static or dynamic. - */ -#ifndef DUMUX_COMMON_MATRIX_VECTOR_HELPER_HH -#define DUMUX_COMMON_MATRIX_VECTOR_HELPER_HH - -#include - -namespace Dumux { - -//! Determines whether or not a matrix has a resize() function -template -struct matrix_has_resize_method -{ -private: - typedef std::true_type yes; - typedef std::false_type no; - - // resize function is called with two indices for matrices - template static auto test(int) -> decltype(std::declval().resize(0, 0), yes()); - template static no test(...); - -public: - static constexpr bool value = std::is_same(0)), yes>::value; -}; - -//! determines whether or not a vector has a resize() function -template -struct vector_has_resize_method -{ -private: - typedef std::true_type yes; - typedef std::false_type no; - - // resize function is called with one index for vectors - template static auto test(int) -> decltype(std::declval().resize(0), yes()); - template static no test(...); - -public: - static constexpr bool value = std::is_same(0)), yes>::value; -}; - -//! resizes a matrix to the given sizes (specialization for dynamic matrix type) -template< class Matrix, - class size_type, - std::enable_if_t::value, int> = 0 > -void resizeMatrix(Matrix& M, size_type rows, size_type cols) -{ - M.resize(rows, cols); -} - -//! resizes a matrix to the given sizes (specialization for static matrix type - do nothing) -template< class Matrix, - class size_type, - std::enable_if_t::value, int> = 0 > -void resizeMatrix(Matrix& M, size_type rows, size_type cols) {} - -//! resizes a vector to the given size (specialization for dynamic matrix type) -template< class Vector, - class size_type, - std::enable_if_t::value, int> = 0 > -void resizeVector(Vector& v, size_type size) -{ - v.resize(size); -} - -//! resizes a vector to the given size (specialization for static vector type - do nothing) -template< class Vector, - class size_type, - std::enable_if_t::value, int> = 0 > -void resizeVector(Vector& v, size_type rows) {} - -} // end namespace Dumux - -#endif diff --git a/dumux/discretization/cellcentered/mpfa/CMakeLists.txt b/dumux/discretization/cellcentered/mpfa/CMakeLists.txt index d1940659e201b9a8ba3bc6d6796fb655800bc45e..51527f79f90d7a6cb5c54a78d5cc0ceb6fd90bd7 100644 --- a/dumux/discretization/cellcentered/mpfa/CMakeLists.txt +++ b/dumux/discretization/cellcentered/mpfa/CMakeLists.txt @@ -19,7 +19,7 @@ gridvolumevariables.hh helper.hh interactionvolumebase.hh interactionvolumedatahandle.hh -localassembler.hh +localassemblerbase.hh localfacedata.hh methods.hh properties.hh diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh index fd18b34195c4f965928a8051d2bb57d69178917d..10d186ec216e63484e71341380260ed3c1eb9839 100644 --- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh @@ -50,6 +50,9 @@ class DarcysLawImplementation using GridView = GetPropType; using Element = typename GridView::template Codim<0>::Entity; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + using FVGridGeometry = GetPropType; using FVElementGeometry = typename FVGridGeometry::LocalView; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; @@ -74,44 +77,34 @@ class DarcysLawImplementation { // get interaction volume related data from the filler class & upate the cache if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) - scvfFluxVarsCache.updateAdvection(problem, - fluxVarsCacheFiller.secondaryInteractionVolume(), + scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.secondaryInteractionVolume(), fluxVarsCacheFiller.secondaryIvLocalFaceData(), - fluxVarsCacheFiller.secondaryIvDataHandle(), - scvf); + fluxVarsCacheFiller.secondaryIvDataHandle()); else - scvfFluxVarsCache.updateAdvection(problem, - fluxVarsCacheFiller.primaryInteractionVolume(), + scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.primaryInteractionVolume(), fluxVarsCacheFiller.primaryIvLocalFaceData(), - fluxVarsCacheFiller.primaryIvDataHandle(), - scvf); + fluxVarsCacheFiller.primaryIvDataHandle()); } }; //! The cache used in conjunction with the mpfa Darcy's Law class MpfaDarcysLawCache { - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - static constexpr int numPhases = GetPropType::numPhases(); - using DualGridNodalIndexSet = GetPropType; using Stencil = typename DualGridNodalIndexSet::NodalGridStencilType; - using MpfaHelper = typename FVGridGeometry::MpfaHelper; - static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs(); + static constexpr bool considerSecondaryIVs = FVGridGeometry::MpfaHelper::considerSecondaryIVs(); + using PrimaryDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle::AdvectionHandle; + using SecondaryDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle::AdvectionHandle; - using PrimaryInteractionVolume = GetPropType; - using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData; - using PrimaryIvDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle; - using PrimaryIvCellVector = typename PrimaryInteractionVolume::Traits::MatVecTraits::CellVector; - using PrimaryIvTij = typename PrimaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type; + //! sets the pointer to the data handle (overload for secondary data handles) + template< bool doSecondary = considerSecondaryIVs, std::enable_if_t = 0 > + void setHandlePointer_(const SecondaryDataHandle& dataHandle) + { secondaryHandlePtr_ = &dataHandle; } - using SecondaryInteractionVolume = GetPropType; - using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData; - using SecondaryIvDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle; - using SecondaryIvCellVector = typename SecondaryInteractionVolume::Traits::MatVecTraits::CellVector; - using SecondaryIvTij = typename SecondaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type; + //! sets the pointer to the data handle (overload for primary data handles) + void setHandlePointer_(const PrimaryDataHandle& dataHandle) + { primaryHandlePtr_ = &dataHandle; } public: // export the filler type @@ -121,152 +114,39 @@ class DarcysLawImplementation * \brief Update cached objects (transmissibilities and gravity). * This is used for updates with primary interaction volumes. * - * \param problem The problem * \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 - void updateAdvection(const Problem& problem, - const PrimaryInteractionVolume& iv, - const PrimaryIvLocalFaceData& localFaceData, - const PrimaryIvDataHandle& dataHandle, - const SubControlVolumeFace &scvf) + template + void updateAdvection(const IV& iv, + const LocalFaceData& localFaceData, + const DataHandle& dataHandle) { - switchFluxSign_ = localFaceData.isOutside(); + switchFluxSign_ = localFaceData.isOutsideFace(); stencil_ = &iv.stencil(); - - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - primaryPj_[pIdx] = &dataHandle.pressures(pIdx); - - static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity"); - const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); - - // standard grids - if (dim == dimWorld) - { - primaryTij_ = &dataHandle.advectionT()[ivLocalIdx]; - if (enableGravity) - for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx]; - } - // surface grids - else - { - if (!localFaceData.isOutside()) - { - primaryTij_ = &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(); - primaryTij_ = &dataHandle.advectionTout()[ivLocalIdx][idxInOutsideFaces]; - if (enableGravity) - for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - g_[phaseIdx] = dataHandle.gravityOutside(phaseIdx)[ivLocalIdx][idxInOutsideFaces]; - } - } - } - - /*! - * \brief Update cached objects (transmissibilities and gravity). - * This is used for updates with secondary interaction volumes. - * - * \param problem The problem - * \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 = 0 > - void updateAdvection(const Problem& problem, - const SecondaryInteractionVolume& iv, - const SecondaryIvLocalFaceData& localFaceData, - const SecondaryIvDataHandle& dataHandle, - const SubControlVolumeFace &scvf) - { - switchFluxSign_ = localFaceData.isOutside(); - stencil_ = &iv.stencil(); - - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - secondaryPj_[pIdx] = &dataHandle.pressures(pIdx); - - static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity"); - const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); - - // standard grids - if (dim == dimWorld) - { - secondaryTij_ = &dataHandle.advectionT()[ivLocalIdx]; - if (enableGravity) - for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx]; - } - // surface grids - else - { - if (!localFaceData.isOutside()) - { - secondaryTij_ = &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(); - secondaryTij_ = &dataHandle.advectionTout()[ivLocalIdx][idxInOutsideFaces]; - if (enableGravity) - for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - g_[phaseIdx] = dataHandle.gravityOutside(phaseIdx)[ivLocalIdx][idxInOutsideFaces]; - } - } + setHandlePointer_(dataHandle.advectionHandle()); } //! The stencil corresponding to the transmissibilities (primary type) const Stencil& advectionStencil() const { return *stencil_; } - //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (primary type) - const PrimaryIvTij& advectionTijPrimaryIv() const { return *primaryTij_; } - - //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (secondary type) - const SecondaryIvTij& advectionTijSecondaryIv() const { return *secondaryTij_; } - - //! The cell (& maybe Dirichlet) pressures within this interaction volume (primary type) - const PrimaryIvCellVector& pressuresPrimaryIv(unsigned int phaseIdx) const { return *primaryPj_[phaseIdx]; } - - //! The cell (& maybe Dirichlet) pressures within this interaction volume (secondary type) - const SecondaryIvCellVector& pressuresSecondaryIv(unsigned int phaseIdx) const { return *secondaryPj_[phaseIdx]; } + //! The corresponding data handles + const PrimaryDataHandle& advectionPrimaryDataHandle() const { return *primaryHandlePtr_; } + const SecondaryDataHandle& advectionSecondaryDataHandle() const { return *secondaryHandlePtr_; } - //! The gravitational acceleration for a phase on this scvf - Scalar gravity(unsigned int phaseIdx) const { return g_[phaseIdx]; } - - //! 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. + //! Returns whether or not this scvf is an "outside" face in the scope of the iv. bool advectionSwitchFluxSign() const { return switchFluxSign_; } private: bool switchFluxSign_; + //! pointers to the corresponding iv-data handles + const PrimaryDataHandle* primaryHandlePtr_; + const SecondaryDataHandle* secondaryHandlePtr_; + //! The stencil, i.e. the grid indices j const Stencil* stencil_; - - //! The transmissibilities such that f = Tij*pj - const PrimaryIvTij* primaryTij_; - const SecondaryIvTij* secondaryTij_; - - //! The interaction-volume wide phase pressures pj - std::array primaryPj_; - std::array secondaryPj_; - - //! Gravitational flux contribution on this face - std::array< Scalar, numPhases > g_; }; public: @@ -287,28 +167,41 @@ public: { const auto& fluxVarsCache = elemFluxVarsCache[scvf]; - // compute t_ij*p_j - Scalar scvfFlux; + // forward to the private function taking the iv data handle if (fluxVarsCache.usesSecondaryIv()) - { - const auto& tij = fluxVarsCache.advectionTijSecondaryIv(); - const auto& pj = fluxVarsCache.pressuresSecondaryIv(phaseIdx); - scvfFlux = tij*pj; - } + return flux_(problem, fluxVarsCache, fluxVarsCache.advectionSecondaryDataHandle(), phaseIdx); else - { - const auto& tij = fluxVarsCache.advectionTijPrimaryIv(); - const auto& pj = fluxVarsCache.pressuresPrimaryIv(phaseIdx); - scvfFlux = tij*pj; - } + return flux_(problem, fluxVarsCache, fluxVarsCache.advectionPrimaryDataHandle(), phaseIdx); + } + +private: + template< class Problem, class FluxVarsCache, class DataHandle > + static Scalar flux_(const Problem& problem, + const FluxVarsCache& cache, + const DataHandle& dataHandle, + int phaseIdx) + { + dataHandle.setPhaseIndex(phaseIdx); + + const bool switchSign = cache.advectionSwitchFluxSign(); + + const auto localFaceIdx = cache.ivLocalFaceIndex(); + const auto idxInOutside = cache.indexInOutsideFaces(); + const auto& pj = dataHandle.uj(); + const auto& tij = dim == dimWorld ? dataHandle.T()[localFaceIdx] + : (!switchSign ? dataHandle.T()[localFaceIdx] + : dataHandle.tijOutside()[localFaceIdx][idxInOutside]); + Scalar scvfFlux = tij*pj; // maybe add gravitational acceleration static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity"); if (enableGravity) - scvfFlux += fluxVarsCache.gravity(phaseIdx); + scvfFlux += dim == dimWorld ? dataHandle.g()[localFaceIdx] + : (!switchSign ? dataHandle.g()[localFaceIdx] + : dataHandle.gOutside()[localFaceIdx][idxInOutside]); // switch the sign if necessary - if (fluxVarsCache.advectionSwitchFluxSign()) + if (cache.advectionSwitchFluxSign()) scvfFlux *= -1.0; return scvfFlux; diff --git a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh index 64802261dc9a22e9bdfe47a4d434c71a52b47120..fcd012caafc51d75c8085b00967ddbe15a960d75 100644 --- a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh +++ b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh @@ -70,8 +70,8 @@ class CCMpfaDualGridNodalIndexSet using LI = typename T::LocalIndexType; using GI = typename T::GridIndexType; - using DimLocalIndexVector = Dune::ReservedVector; - using ScvfIndicesInScvStorage = typename T::template NodalScvDataStorage< DimLocalIndexVector >; + using DimIndexVector = Dune::ReservedVector; + using ScvfIndicesInScvStorage = typename T::template NodalScvDataStorage< DimIndexVector >; public: //! Export the traits type @@ -96,9 +96,7 @@ public: template void insert(const SubControlVolumeFace& scvf) { - insert(scvf.index(), - scvf.insideScvIdx(), - scvf.boundary()); + insert(scvf.index(), scvf.insideScvIdx(), scvf.boundary()); } //! Inserts scvf data @@ -148,11 +146,11 @@ public: { return numBoundaryScvfs_; } //! returns the grid scv indices connected to this dual grid node - const NodalGridStencilType& globalScvIndices() const + const NodalGridStencilType& gridScvIndices() const { return scvIndices_; } //! returns the grid scvf indices connected to this dual grid node - const NodalGridScvfStencilType& globalScvfIndices() const + const NodalGridScvfStencilType& gridScvfIndices() const { return scvfIndices_; } //! returns whether or not the i-th scvf is on a domain boundary @@ -163,21 +161,21 @@ public: } //! returns the grid scv idx of the i-th scv - GridIndexType scvIdxGlobal(unsigned int i) const + GridIndexType gridScvIndex(unsigned int i) const { assert(i < numScvs()); return scvIndices_[i]; } //! returns the index of the i-th scvf - GridIndexType scvfIdxGlobal(unsigned int i) const + GridIndexType gridScvfIndex(unsigned int i) const { assert(i < numScvfs()); return scvfIndices_[i]; } //! returns the grid index of the j-th scvf embedded in the i-th scv - GridIndexType scvfIdxGlobal(unsigned int i, unsigned int j) const + GridIndexType gridScvfIndex(unsigned int i, unsigned int j) const { assert(i < numScvs()); assert(j < localScvfIndicesInScv_[i].size()); @@ -185,7 +183,7 @@ public: } //! 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 + LocalIndexType localScvfIndex(unsigned int i, unsigned int j) const { assert(i < numScvs()); assert(j < localScvfIndicesInScv_[i].size()); @@ -193,7 +191,7 @@ public: } //! returns the node-local index of the inside scv of the i-th scvf - LocalIndexType insideScvIdxLocal(unsigned int i) const + LocalIndexType insideScvLocalIndex(unsigned int i) const { assert(i < numScvfs()); return scvfInsideScvIndices_[i]; diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index 675980bf32850b70718c5f1f70b8803fb99f7284..7f8635d08aa24ea8cfe711daff5db1ed7adfbc17 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -317,7 +317,7 @@ private: void addBoundaryVolVars_(const Problem& problem, const FVElementGeometry& fvGeometry, const NodalIndexSet& nodalIndexSet) { // check each scvf in the index set for boundary presence - for (auto scvfIdx : nodalIndexSet.globalScvfIndices()) + for (auto scvfIdx : nodalIndexSet.gridScvfIndices()) { const auto& ivScvf = fvGeometry.scvf(scvfIdx); diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh index d9afded35c0da64b34acb164e3cf66c40058a9fb..3d5e759d4d95babc40e4bd7a5de687110319cde0 100644 --- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh @@ -46,6 +46,9 @@ class FicksLawImplementation using GridView = GetPropType; using Element = typename GridView::template Codim<0>::Entity; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + using FVGridGeometry = GetPropType; using FVElementGeometry = typename FVGridGeometry::LocalView; using FluidSystem = GetPropType; @@ -79,12 +82,12 @@ class FicksLawImplementation scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.secondaryInteractionVolume(), fluxVarsCacheFiller.secondaryIvLocalFaceData(), fluxVarsCacheFiller.secondaryIvDataHandle(), - scvf, phaseIdx, compIdx); + phaseIdx, compIdx); else scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.primaryInteractionVolume(), fluxVarsCacheFiller.primaryIvLocalFaceData(), fluxVarsCacheFiller.primaryIvDataHandle(), - scvf, phaseIdx, compIdx); + phaseIdx, compIdx); } }; @@ -94,24 +97,19 @@ class FicksLawImplementation using DualGridNodalIndexSet = GetPropType; using Stencil = typename DualGridNodalIndexSet::NodalGridStencilType; - using MpfaHelper = typename FVGridGeometry::MpfaHelper; - static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs(); - - using PrimaryInteractionVolume = GetPropType; - using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData; - using PrimaryIvDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle; - using PrimaryIvCellVector = typename PrimaryInteractionVolume::Traits::MatVecTraits::CellVector; - using PrimaryIvTij = typename PrimaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type; + static constexpr int numPhases = GetPropType::numPhases(); + static constexpr bool considerSecondaryIVs = FVGridGeometry::MpfaHelper::considerSecondaryIVs(); + using PrimaryDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle::DiffusionHandle; + using SecondaryDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle::DiffusionHandle; - using SecondaryInteractionVolume = GetPropType; - using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData; - using SecondaryIvDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle; - using SecondaryIvCellVector = typename SecondaryInteractionVolume::Traits::MatVecTraits::CellVector; - using SecondaryIvTij = typename SecondaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type; + //! sets the pointer to the data handle (overload for secondary data handles) + template< bool doSecondary = considerSecondaryIVs, std::enable_if_t = 0 > + void setHandlePointer_(const SecondaryDataHandle& dataHandle) + { secondaryHandlePtr_ = &dataHandle; } - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - static constexpr int numPhases = GetPropType::numPhases(); + //! sets the pointer to the data handle (overload for primary data handles) + void setHandlePointer_(const PrimaryDataHandle& dataHandle) + { primaryHandlePtr_ = &dataHandle; } public: // export filler type @@ -124,98 +122,39 @@ class FicksLawImplementation * \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 */ - void updateDiffusion(const PrimaryInteractionVolume& iv, - const PrimaryIvLocalFaceData& localFaceData, - const PrimaryIvDataHandle& dataHandle, - const SubControlVolumeFace &scvf, + template + void updateDiffusion(const IV& iv, + const LocalFaceData& localFaceData, + const DataHandle& dataHandle, unsigned int phaseIdx, unsigned int compIdx) { stencil_[phaseIdx][compIdx] = &iv.stencil(); - switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside(); - - // store pointer to the mole fraction vector of this iv - primaryXj_[phaseIdx][compIdx] = &dataHandle.moleFractions(phaseIdx, compIdx); - - const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); - if (dim == dimWorld) - primaryTij_[phaseIdx][compIdx] = &dataHandle.diffusionT()[ivLocalIdx]; - else - primaryTij_[phaseIdx][compIdx] = localFaceData.isOutside() ? &dataHandle.diffusionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()] - : &dataHandle.diffusionT()[ivLocalIdx]; + switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutsideFace(); + setHandlePointer_(dataHandle.diffusionHandle()); } - /*! - * \brief Update cached objects (transmissibilities). - * This is used for updates with secondary interaction volumes. - * - * \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< bool doSecondary = considerSecondaryIVs, std::enable_if_t = 0 > - void updateDiffusion(const SecondaryInteractionVolume& iv, - const SecondaryIvLocalFaceData& localFaceData, - const SecondaryIvDataHandle& dataHandle, - const SubControlVolumeFace &scvf, - unsigned int phaseIdx, unsigned int compIdx) - { - stencil_[phaseIdx][compIdx] = &iv.stencil(); - switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside(); - - // store pointer to the mole fraction vector of this iv - secondaryXj_[phaseIdx][compIdx] = &dataHandle.moleFractions(phaseIdx, compIdx); + //! The stencils corresponding to the transmissibilities + const Stencil& diffusionStencil(unsigned int phaseIdx, unsigned int compIdx) const + { return *stencil_[phaseIdx][compIdx]; } - const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); - if (dim == dimWorld) - secondaryTij_[phaseIdx][compIdx] = &dataHandle.diffusionT()[ivLocalIdx]; - else - secondaryTij_[phaseIdx][compIdx] = localFaceData.isOutside() ? &dataHandle.diffusionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()] - : &dataHandle.diffusionT()[ivLocalIdx]; - } + //! The corresponding data handles + const PrimaryDataHandle& diffusionPrimaryDataHandle() const { return *primaryHandlePtr_; } + const SecondaryDataHandle& diffusionSecondaryDataHandle() const { return *secondaryHandlePtr_; } - //! 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. + //! Returns whether or not this scvf is an "outside" face in the scope of the iv. bool diffusionSwitchFluxSign(unsigned int phaseIdx, unsigned int compIdx) const { return switchFluxSign_[phaseIdx][compIdx]; } - //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (primary type) - const PrimaryIvTij& diffusionTijPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const - { return *primaryTij_[phaseIdx][compIdx]; } - - //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (secondary type) - const SecondaryIvTij& diffusionTijSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const - { return *secondaryTij_[phaseIdx][compIdx]; } - - //! The cell (& maybe Dirichlet) mole fractions within this interaction volume (primary type) - const PrimaryIvCellVector& moleFractionsPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const - { return *primaryXj_[phaseIdx][compIdx]; } - - //! The cell (& maybe Dirichlet) mole fractions within this interaction volume (secondary type) - const SecondaryIvCellVector& moleFractionsSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const - { return *secondaryXj_[phaseIdx][compIdx]; } - - //! The stencils corresponding to the transmissibilities - const Stencil& diffusionStencil(unsigned int phaseIdx, unsigned int compIdx) const - { return *stencil_[phaseIdx][compIdx]; } private: + //! phase-/component- specific data std::array< std::array, numPhases > switchFluxSign_; - - //! The stencils, i.e. the grid indices j std::array< std::array, numPhases > stencil_; - //! The transmissibilities such that f = Tij*xj - std::array< std::array, numPhases > primaryTij_; - std::array< std::array, numPhases > secondaryTij_; - - //! The interaction-volume wide mole fractions xj - std::array< std::array, numPhases > primaryXj_; - std::array< std::array, numPhases > secondaryXj_; + //! pointers to the corresponding iv-data handles + const PrimaryDataHandle* primaryHandlePtr_; + const SecondaryDataHandle* secondaryHandlePtr_; }; public: @@ -253,25 +192,17 @@ public: // calculate the density at the interface const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx); - // calculate Tij*xj - Scalar flux; + // compute the flux if (fluxVarsCache.usesSecondaryIv()) - { - const auto& tij = fluxVarsCache.diffusionTijSecondaryIv(phaseIdx, compIdx); - const auto& xj = fluxVarsCache.moleFractionsSecondaryIv(phaseIdx, compIdx); - flux = tij*xj; - } + componentFlux[compIdx] = rho*effFactor*computeVolumeFlux(problem, + fluxVarsCache, + fluxVarsCache.diffusionSecondaryDataHandle(), + phaseIdx, compIdx); else - { - const auto& tij = fluxVarsCache.diffusionTijPrimaryIv(phaseIdx, compIdx); - const auto& xj = fluxVarsCache.moleFractionsPrimaryIv(phaseIdx, compIdx); - flux = tij*xj; - } - - if (fluxVarsCache.diffusionSwitchFluxSign(phaseIdx, compIdx)) - flux *= -1.0; - - componentFlux[compIdx] = flux*rho*effFactor; + componentFlux[compIdx] = rho*effFactor*computeVolumeFlux(problem, + fluxVarsCache, + fluxVarsCache.diffusionPrimaryDataHandle(), + phaseIdx, compIdx); } // accumulate the phase component flux @@ -283,6 +214,32 @@ public: } private: + template< class Problem, class FluxVarsCache, class DataHandle > + static Scalar computeVolumeFlux(const Problem& problem, + const FluxVarsCache& cache, + const DataHandle& dataHandle, + int phaseIdx, int compIdx) + { + dataHandle.setPhaseIndex(phaseIdx); + dataHandle.setComponentIndex(compIdx); + + const bool switchSign = cache.diffusionSwitchFluxSign(phaseIdx, compIdx); + + const auto localFaceIdx = cache.ivLocalFaceIndex(); + const auto idxInOutside = cache.indexInOutsideFaces(); + const auto& xj = dataHandle.uj(); + const auto& tij = dim == dimWorld ? dataHandle.T()[localFaceIdx] + : (!switchSign ? dataHandle.T()[localFaceIdx] + : dataHandle.tijOutside()[localFaceIdx][idxInOutside]); + Scalar scvfFlux = tij*xj; + + // switch the sign if necessary + if (cache.advectionSwitchFluxSign()) + scvfFlux *= -1.0; + + return scvfFlux; + } + //! compute the density at branching facets for network grids as arithmetic mean static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf, diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh index 8d8b4968a906578f3dbfee41c7ee425f738c96c8..30ccbd1282166edaca392bf3ca1e400c9bb0b284 100644 --- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh +++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh @@ -29,7 +29,6 @@ #include #include -#include namespace Dumux { @@ -120,12 +119,11 @@ public: const auto& indexSet = fvGridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf); fluxVarsCacheContainer.secondaryInteractionVolumes().emplace_back(); secondaryIv_ = &fluxVarsCacheContainer.secondaryInteractionVolumes().back(); - secondaryIv_->setUpLocalScope(indexSet, problem(), fvGeometry); + secondaryIv_->bind(indexSet, problem(), fvGeometry); // prepare the corresponding data handle 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_, *secondaryIvDataHandle_, ivIndexInContainer, true); @@ -151,12 +149,11 @@ public: const auto& indexSet = fvGridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf); fluxVarsCacheContainer.primaryInteractionVolumes().emplace_back(); primaryIv_ = &fluxVarsCacheContainer.primaryInteractionVolumes().back(); - primaryIv_->setUpLocalScope(indexSet, problem(), fvGeometry); + primaryIv_->bind(indexSet, problem(), fvGeometry); // prepare the corresponding data handle 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_, *primaryIvDataHandle_, ivIndexInContainer, true); @@ -227,12 +224,15 @@ private: for (const auto& d : iv.localFaceData()) { // obtain the scvf - const auto& scvfJ = fvGeometry().scvf(d.globalScvfIndex()); + const auto& scvfJ = fvGeometry().scvf(d.gridScvfIndex()); ivScvfs[i] = &scvfJ; ivFluxVarCaches[i] = &fluxVarsCacheContainer[scvfJ]; ivFluxVarCaches[i]->setIvIndexInContainer(ivIndexInContainer); ivFluxVarCaches[i]->setUpdateStatus(true); ivFluxVarCaches[i]->setSecondaryIvUsage(isSecondary); + ivFluxVarCaches[i]->setIvLocalFaceIndex(d.ivLocalScvfIndex()); + if (dim < dimWorld) + ivFluxVarCaches[i]->setIndexInOutsideFaces(d.scvfLocalOutsideScvfIndex()); i++; } @@ -316,6 +316,7 @@ private: for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + handle.diffusionHandle().setPhaseIndex(phaseIdx); for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx) { using FluidSystem = GetPropType; @@ -323,6 +324,7 @@ private: continue; // fill data in the handle + handle.diffusionHandle().setComponentIndex(compIdx); fillDiffusionHandle(iv, handle, forceUpdateAll, phaseIdx, compIdx); // fill diffusion caches @@ -434,75 +436,29 @@ private: using LambdaFactory = TensorLambdaFactory; // get instance of the interaction volume-local assembler - static constexpr MpfaMethods M = InteractionVolume::MpfaMethod; - using IvLocalAssembler = InteractionVolumeAssembler< Problem, FVElementGeometry, ElementVolumeVariables, M >; + using Traits = typename InteractionVolume::Traits; + using IvLocalAssembler = typename Traits::template LocalAssembler; IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars()); - // Use different assembly if gravity is enabled - static const bool enableGravity = getParamFromGroup(problem().paramGroup(), "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() ); - } + localAssembler.assembleMatrices(handle.advectionHandle(), iv, LambdaFactory::getAdvectionLambda()); // assemble pressure vectors for (unsigned int pIdx = 0; pIdx < ModelTraits::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); + // set context in handle + handle.advectionHandle().setPhaseIndex(pIdx); + + // maybe (re-)assemble gravity contribution vector + auto getRho = [pIdx] (const auto& volVars) { return volVars.density(pIdx); }; + static const bool enableGravity = getParamFromGroup(problem().paramGroup(), "Problem.EnableGravity"); + if (enableGravity) + localAssembler.assembleGravity(handle.advectionHandle(), iv, getRho); + + // reassemble pressure vector + auto getPressure = [pIdx] (const auto& volVars) { return volVars.pressure(pIdx); }; + localAssembler.assembleU(handle.advectionHandle(), iv, getPressure); } } @@ -519,33 +475,19 @@ private: using LambdaFactory = TensorLambdaFactory; // get instance of the interaction volume-local assembler - static constexpr MpfaMethods M = InteractionVolume::MpfaMethod; - using IvLocalAssembler = InteractionVolumeAssembler< Problem, FVElementGeometry, ElementVolumeVariables, M >; + using Traits = typename InteractionVolume::Traits; + using IvLocalAssembler = typename Traits::template LocalAssembler; IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars()); - // solve the local system subject to the tensor and update the handle - handle.setPhaseIndex(phaseIdx); - handle.setComponentIndex(compIdx); - - // assemble T + // maybe (re-)assemble matrices 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) ); - } + localAssembler.assembleMatrices(handle.diffusionHandle(), + 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); + auto getMoleFraction = [phaseIdx, compIdx] (const auto& volVars) { return volVars.moleFraction(phaseIdx, compIdx); }; + localAssembler.assembleU(handle.diffusionHandle(), iv, getMoleFraction); } //! prepares the quantities necessary for conductive fluxes in the handle @@ -556,29 +498,22 @@ private: void fillHeatConductionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll) { using LambdaFactory = TensorLambdaFactory; + using ThermCondModel = GetPropType; // get instance of the interaction volume-local assembler - static constexpr MpfaMethods M = InteractionVolume::MpfaMethod; - using IvLocalAssembler = InteractionVolumeAssembler< Problem, FVElementGeometry, ElementVolumeVariables, M >; + using Traits = typename InteractionVolume::Traits; + using IvLocalAssembler = typename Traits::template LocalAssembler; IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars()); - if (forceUpdateAll || soldependentAdvection) - { - if (dim < dimWorld) - localAssembler.assemble( handle.heatConductionTout(), - handle.heatConductionT(), - iv, - LambdaFactory::template getHeatConductionLambda>() ); - else - localAssembler.assemble( handle.heatConductionT(), - iv, - LambdaFactory::template getHeatConductionLambda>() ); - } + // maybe (re-)assemble matrices + if (forceUpdateAll || soldependentHeatConduction) + localAssembler.assembleMatrices(handle.heatConductionHandle(), + iv, + LambdaFactory::template 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); + auto getMoleFraction = [] (const auto& volVars) { return volVars.temperature(); }; + localAssembler.assembleU(handle.heatConductionHandle(), iv, getMoleFraction); } //! fill handle only when advection uses mpfa diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh index db27fa6202ed1c534e10289f50a3f6a700c5f2c3..7b56dfb9e543021a68d34a1a465c57c0e92224fb 100644 --- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh @@ -48,6 +48,9 @@ class FouriersLawImplementation using GridView = GetPropType; using Element = typename GridView::template Codim<0>::Entity; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + using FVGridGeometry = GetPropType; using FVElementGeometry = typename FVGridGeometry::LocalView; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; @@ -75,13 +78,11 @@ class FouriersLawImplementation if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.secondaryInteractionVolume(), fluxVarsCacheFiller.secondaryIvLocalFaceData(), - fluxVarsCacheFiller.secondaryIvDataHandle(), - scvf); + fluxVarsCacheFiller.secondaryIvDataHandle()); else scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.primaryInteractionVolume(), fluxVarsCacheFiller.primaryIvLocalFaceData(), - fluxVarsCacheFiller.primaryIvDataHandle(), - scvf); + fluxVarsCacheFiller.primaryIvDataHandle()); } }; @@ -91,23 +92,18 @@ class FouriersLawImplementation using DualGridNodalIndexSet = GetPropType; using Stencil = typename DualGridNodalIndexSet::NodalGridStencilType; - using MpfaHelper = typename FVGridGeometry::MpfaHelper; - static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs(); - - using PrimaryInteractionVolume = GetPropType; - using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData; - using PrimaryIvDataHandle = typename ElementFluxVarsCache::PrimaryIvDataHandle; - using PrimaryIvCellVector = typename PrimaryInteractionVolume::Traits::MatVecTraits::CellVector; - using PrimaryIvTij = typename PrimaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type; + static constexpr bool considerSecondaryIVs = FVGridGeometry::MpfaHelper::considerSecondaryIVs(); + using PrimaryDataHandle = typename ElementFluxVarsCache::PrimaryIvDataHandle::HeatConductionHandle; + using SecondaryDataHandle = typename ElementFluxVarsCache::SecondaryIvDataHandle::HeatConductionHandle; - using SecondaryInteractionVolume = GetPropType; - using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData; - using SecondaryIvDataHandle = typename ElementFluxVarsCache::SecondaryIvDataHandle; - using SecondaryIvCellVector = typename SecondaryInteractionVolume::Traits::MatVecTraits::CellVector; - using SecondaryIvTij = typename SecondaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type; + //! sets the pointer to the data handle (overload for secondary data handles) + template< bool doSecondary = considerSecondaryIVs, std::enable_if_t = 0 > + void setHandlePointer_(const SecondaryDataHandle& dataHandle) + { secondaryHandlePtr_ = &dataHandle; } - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; + //! sets the pointer to the data handle (overload for primary data handles) + void setHandlePointer_(const PrimaryDataHandle& dataHandle) + { primaryHandlePtr_ = &dataHandle; } public: // export filler type @@ -120,90 +116,36 @@ class FouriersLawImplementation * \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 */ - void updateHeatConduction(const PrimaryInteractionVolume& iv, - const PrimaryIvLocalFaceData& localFaceData, - const PrimaryIvDataHandle& dataHandle, - const SubControlVolumeFace &scvf) + template + void updateHeatConduction(const IV& iv, + const LocalFaceData& localFaceData, + const DataHandle& dataHandle) { + switchFluxSign_ = localFaceData.isOutsideFace(); stencil_ = &iv.stencil(); - switchFluxSign_ = localFaceData.isOutside(); - - // store pointer to the temperature vector of this iv - primaryTj_ = &dataHandle.temperatures(); - - const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); - if (dim == dimWorld) - primaryTij_ = &dataHandle.heatConductionT()[ivLocalIdx]; - else - primaryTij_ = localFaceData.isOutside() ? &dataHandle.heatConductionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()] - : &dataHandle.heatConductionT()[ivLocalIdx]; + setHandlePointer_(dataHandle.heatConductionHandle()); } - /*! - * \brief Update cached objects (transmissibilities). - * This is used for updates with secondary interaction volumes. - * - * \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< bool doSecondary = considerSecondaryIVs, std::enable_if_t = 0 > - void updateHeatConduction(const SecondaryInteractionVolume& iv, - const SecondaryIvLocalFaceData& localFaceData, - const SecondaryIvDataHandle& dataHandle, - const SubControlVolumeFace &scvf) - { - stencil_ = &iv.stencil(); - switchFluxSign_ = localFaceData.isOutside(); - - // store pointer to the temperature vector of this iv - secondaryTj_ = &dataHandle.temperatures(); - - const auto ivLocalIdx = localFaceData.ivLocalScvfIndex(); - if (dim == dimWorld) - secondaryTij_ = &dataHandle.heatConductionT()[ivLocalIdx]; - else - secondaryTij_ = localFaceData.isOutside() ? &dataHandle.heatConductionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()] - : &dataHandle.heatConductionT()[ivLocalIdx]; - } - - //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (primary type) - const PrimaryIvTij& heatConductionTijPrimaryIv() const { return *primaryTij_; } - - //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (secondary type) - const SecondaryIvTij& heatConductionTijSecondaryIv() const { return *secondaryTij_; } - //! The stencil corresponding to the transmissibilities (primary type) const Stencil& heatConductionStencil() const { return *stencil_; } - //! The cell (& Dirichlet) temperatures within this interaction volume (primary type) - const PrimaryIvCellVector& temperaturesPrimaryIv() const { return *primaryTj_; } - - //! The cell (& Dirichlet) temperatures within this interaction volume (secondary type) - const SecondaryIvCellVector& temperaturesSecondaryIv() const { return *secondaryTj_; } + //! The corresponding data handles + const PrimaryDataHandle& heatConductionPrimaryDataHandle() const { return *primaryHandlePtr_; } + const SecondaryDataHandle& heatConductionSecondaryDataHandle() const { return *secondaryHandlePtr_; } - //! 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. + //! Returns whether or not this scvf is an "outside" face in the scope of the iv. bool heatConductionSwitchFluxSign() const { return switchFluxSign_; } private: bool switchFluxSign_; + //! pointers to the corresponding iv-data handles + const PrimaryDataHandle* primaryHandlePtr_; + const SecondaryDataHandle* secondaryHandlePtr_; + //! The stencil, i.e. the grid indices j const Stencil* stencil_; - - //! The transmissibilities such that f = Tij*Tj - const PrimaryIvTij* primaryTij_; - const SecondaryIvTij* secondaryTij_; - - //! The interaction-volume wide temperature Tj - const PrimaryIvCellVector* primaryTj_; - const SecondaryIvCellVector* secondaryTj_; }; public: @@ -223,25 +165,34 @@ public: { const auto& fluxVarsCache = elemFluxVarsCache[scvf]; - // compute Tij*tj - Scalar flux; + // forward to the private function taking the iv data handle if (fluxVarsCache.usesSecondaryIv()) - { - const auto& tij = fluxVarsCache.heatConductionTijSecondaryIv(); - const auto& Tj = fluxVarsCache.temperaturesSecondaryIv(); - flux = tij*Tj; - } + return flux_(problem, fluxVarsCache, fluxVarsCache.heatConductionSecondaryDataHandle()); else - { - const auto& tij = fluxVarsCache.heatConductionTijPrimaryIv(); - const auto& Tj = fluxVarsCache.temperaturesPrimaryIv(); - flux = tij*Tj; - } + return flux_(problem, fluxVarsCache, fluxVarsCache.heatConductionPrimaryDataHandle()); + } + +private: + template< class Problem, class FluxVarsCache, class DataHandle > + static Scalar flux_(const Problem& problem, + const FluxVarsCache& cache, + const DataHandle& dataHandle) + { + const bool switchSign = cache.advectionSwitchFluxSign(); + + const auto localFaceIdx = cache.ivLocalFaceIndex(); + const auto idxInOutside = cache.indexInOutsideFaces(); + const auto& Tj = dataHandle.uj(); + const auto& tij = dim == dimWorld ? dataHandle.T()[localFaceIdx] + : (!switchSign ? dataHandle.T()[localFaceIdx] + : dataHandle.tijOutside()[localFaceIdx][idxInOutside]); + Scalar scvfFlux = tij*Tj; - if (fluxVarsCache.heatConductionSwitchFluxSign()) - flux *= -1.0; + // switch the sign if necessary + if (cache.advectionSwitchFluxSign()) + scvfFlux *= -1.0; - return flux; + return scvfFlux; } }; diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh index a1f315aec636f5cf901abd20f2f102139e5bec44..8c73955c03d2f4c2c0baccdad89980e1b1cf3880 100644 --- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh +++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh @@ -166,7 +166,8 @@ public: //! Returns true if secondary interaction volumes are used around a given vertex (index). //! If the use of secondary interaction volumes is disabled, this can be evaluated at compile time. template = 0> - constexpr bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const { return false; } + constexpr bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const + { return false; } //! update all fvElementGeometries (do this again after grid adaption) void update() @@ -354,7 +355,7 @@ public: // building the geometries has finished std::cout << "Initializing of the grid finite volume geometry took " << timer.elapsed() << " seconds." << std::endl; - // Initialize the grid interaction volume seeds + // Initialize the grid interaction volume index sets timer.reset(); ivIndexSets_.update(*this, std::move(dualIdSet)); std::cout << "Initializing of the grid interaction volume index sets took " << timer.elapsed() << " seconds." << std::endl; @@ -533,7 +534,8 @@ public: //! Returns true if secondary interaction volumes are used around a given vertex (index). //! If the use of secondary interaction volumes is disabled, this can be evaluated at compile time. template = 0> - constexpr bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const { return false; } + constexpr bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const + { return false; } //! Returns true if a given vertex lies on a processor boundary inside a ghost element. bool isGhostVertex(const Vertex& v) const @@ -715,12 +717,12 @@ public: // building the geometries has finished std::cout << "Initializing of the grid finite volume geometry took " << timer.elapsed() << " seconds." << std::endl; - // Initialize the grid interaction volume seeds + // Initialize the grid interaction volume index sets timer.reset(); ivIndexSets_.update(*this, std::move(dualIdSet)); std::cout << "Initializing of the grid interaction volume index sets took " << timer.elapsed() << " seconds." << std::endl; - // build the connectivity map for an effecient assembly + // build the connectivity map for an efficient assembly timer.reset(); connectivityMap_.update(*this); std::cout << "Initializing of the connectivity map took " << timer.elapsed() << " seconds." << std::endl; diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh index 318df5512a1ead029d2c06b953f33b37f69396b2..63debf94c6b5c0181818d5a774a73f33e87a1d5f 100644 --- a/dumux/discretization/cellcentered/mpfa/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/helper.hh @@ -65,7 +65,6 @@ class MpfaDimensionHelper public: /*! * \brief Calculates the inner normal vectors to a given scv basis. - * * \param scvBasis The basis of an scv */ template< class ScvBasis > @@ -89,7 +88,6 @@ public: /*! * \brief Calculates the determinant of an scv basis. * This is equal to the cross product for dim = dimWorld = 2 - * * \param scvBasis The basis of an scv */ template< class ScvBasis > @@ -102,7 +100,6 @@ public: /*! * \brief Returns the global number of scvfs in the grid. Assumes the grid to be made up of only * basic geometry types. Overlad this function if you want to use different geometry types. - * * \param gridView The grid view to be checked */ static std::size_t getGlobalNumScvf(const GridView& gridView) @@ -118,7 +115,6 @@ public: /*! * \brief Checks whether or not a given scv basis forms a right hand system. - * * \param scvBasis The basis of an scv */ template< class ScvBasis > @@ -172,15 +168,12 @@ public: assert(cornerIdx < 2 && "provided index exceeds the number of corners of facets in 2d"); // create & return the scvf corner vector - if (cornerIdx == 0) - return ScvfCornerVector({p[0], p[1]}); - else - return ScvfCornerVector({p[0], p[2]}); + return cornerIdx == 0 ? ScvfCornerVector({p[0], p[1]}) + : ScvfCornerVector({p[0], p[2]}); } /*! * \brief Calculates the area of an scvf. - * * \param scvfCorners Container with the corners of the scvf */ static CoordScalar getScvfArea(const ScvfCornerVector& scvfCorners) @@ -188,7 +181,6 @@ public: /*! * \brief Calculates the number of scvfs in a given element geometry type. - * * \param gt The element geometry type */ static std::size_t getNumLocalScvfs(const Dune::GeometryType gt) @@ -209,7 +201,7 @@ public: */ template class MpfaDimensionHelper - : public MpfaDimensionHelper +: public MpfaDimensionHelper { using GridView = typename FVGridGeometry::GridView; using CoordScalar = typename GridView::ctype; @@ -217,18 +209,18 @@ public: /*! * \brief Calculates the inner normal vectors to a given scv basis. - * * \param scvBasis The basis of an scv */ template< class ScvBasis > static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis) { // compute vector normal to the basis plane - const auto normal = [&] () { - auto n = crossProduct(scvBasis[0], scvBasis[1]); - n /= n.two_norm(); - return n; - } (); + const auto normal = [&scvBasis] () + { + auto n = crossProduct(scvBasis[0], scvBasis[1]); + n /= n.two_norm(); + return n; + } (); // compute inner normals using the normal vector ScvBasis innerNormals; @@ -302,11 +294,7 @@ public: innerNormals[2] = crossProduct(scvBasis[0], scvBasis[1]); if (!isRightHandSystem(scvBasis)) - { - innerNormals[0] *= -1.0; - innerNormals[1] *= -1.0; - innerNormals[2] *= -1.0; - } + std::for_each(innerNormals.begin(), innerNormals.end(), [] (auto& n) { n *= -1.0; }); return innerNormals; } @@ -349,7 +337,6 @@ public: /*! * \brief Checks whether or not a given scv basis forms a right hand system. - * * \param scvBasis The basis of an scv */ template< class ScvBasis > @@ -479,7 +466,6 @@ public: /*! * \brief Calculates the area of an scvf. - * * \param scvfCorners Container with the corners of the scvf */ static CoordScalar getScvfArea(const ScvfCornerVector& scvfCorners) diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh index fc14cadc8a0d00dabbdfb24a103aeae7fc61bfb8..f399ac51dff9d05395ff2433b06839a572c37b41 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh @@ -51,6 +51,8 @@ namespace Dumux * using LocalFaceData = ...; * //! export the matrix/vector type traits to be used by the iv * using MatVecTraits = ...; + * //! export the type used for the assembly of the iv's local eq system + * using LocalAssembler = ...; * \endcode */ @@ -59,16 +61,11 @@ namespace Dumux * \brief Base class for the interaction volumes of mpfa methods. It defines * the interface and actual implementations should derive from this class. * - * \tparam Impl The actual implementation of the interaction volume * \tparam T The traits class to be used */ -template< class Impl, class T > +template< class T > class CCMpfaInteractionVolumeBase { - // Curiously recurring template pattern - Impl& asImp() { return static_cast(*this); } - const Impl& asImp() const { return static_cast(*this); } - using GridView = typename T::GridView; using Element = typename GridView::template Codim<0>::Entity; @@ -83,53 +80,53 @@ public: //! Prepares everything for the assembly template< class Problem, class FVElementGeometry > - void setUpLocalScope(const typename Traits::IndexSet& indexSet, - const Problem& problem, - const FVElementGeometry& fvGeometry) - { asImp().setUpLocalScope(); } + void bind(const typename Traits::IndexSet& indexSet, + const Problem& problem, + const FVElementGeometry& fvGeometry) + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a bind() function"); } //! returns the number of "primary" scvfs of this interaction volume std::size_t numFaces() const - { return asImp().numFaces(); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a numFaces() function"); } //! returns the number of intermediate unknowns within this interaction volume std::size_t numUnknowns() const - { return asImp().numUnknowns(); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a numUnknowns() function"); } //! returns the number of (in this context) known solution values within this interaction volume std::size_t numKnowns() const - { return asImp().numKnowns(); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a numKnowns() function"); } //! returns the number of scvs embedded in this interaction volume std::size_t numScvs() const - { return asImp().numScvs(); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a numScvs() function"); } //! Returns a reference to the container with the local face data. The actual type of //! the container depends on the interaction volume implementation. At this point we throw //! an exception and force the implementation to overload this function. const std::vector& localFaceData() const - { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a localFaceData() funtion"); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a localFaceData() function"); } //! returns the cell-stencil of this interaction volume const NodalStencilType& stencil() const - { return asImp().stencil(); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a stencil() function"); } //! returns the local scvf entity corresponding to a given iv-local scvf idx const LocalScvfType& localScvf(LocalIndexType ivLocalScvfIdx) const - { return asImp().localScvf(ivLocalScvfIdx); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a localScvf() function"); } //! returns the local scv entity corresponding to a given iv-local scv idx const LocalScvType& localScv(LocalIndexType ivLocalScvIdx) const - { return asImp().localScv(ivLocalScvIdx); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a localScv() function"); } //! returns the element in which the scv with the given local idx is embedded in const Element& element(LocalIndexType ivLocalScvIdx) const - { return asImp().element(); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide an element() function"); } //! returns the number of interaction volumes living around a vertex template< class NodalIndexSet > static std::size_t numIVAtVertex(const NodalIndexSet& nodalIndexSet) - { return Impl::numIVAtVertex(nodalIndexSet); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide a numIVAtVertex() function"); } //! 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 @@ -141,7 +138,7 @@ public: ScvfIndexMap& scvfIndexMap, const NodalIndexSet& nodalIndexSet, const FlipScvfIndexSet& flipScvfIndexSet) - { Impl::addIVIndexSets(ivIndexSetContainer, scvfIndexMap, nodalIndexSet, flipScvfIndexSet); } + { DUNE_THROW(Dune::NotImplemented, "Interaction volume does not provide an addIVIndexSets() function"); } }; } // end namespace Dumux diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh index 72e2ddc632184a61f1083d35e6d09501218890b4..ebcebfd5a5f364a21a7bd2a1e1dbac0f261bc381 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh @@ -24,278 +24,189 @@ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEDATAHANDLE_HH #define DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEDATAHANDLE_HH +#include #include #include #include -#include -namespace Dumux +namespace Dumux { +namespace Detail { + +/*! + * \ingroup CCMpfaDiscretization + * \brief Common base class to all handles. Stores arrays of the + * matrices involved in the interaction volume-local systems + * of equations. Apart from the transmissibility matrix we + * store those matrices that are needed e.g. for later face + * pressure reconstruction. + * The fluxes as well as the local systems of equations can + * be expressed as functions of the intermediate unknown face + * face values \f$\bar{\mathbf{u}}\f$ and the known cell/Dirichlet + * values \f$\mathbf{u}\f$ using the matrices \f$\mathbf{A}\f$, + * \f$\mathbf{B}\f$, \f$\mathbf{C}\f$, \f$\mathbf{D}\f$ and the + * vector of Neumann fluxes \f$\mathbf{N}\f$ as follows: + * + * Fluxes: \f$\mathbf{f} = \mathbf{C}\bar{\mathbf{u}} + \mathbf{D}\mathbf{u}\f$ + * Local eq system: \f$\mathbf{A}\bar{\mathbf{u}} = \mathbf{B}\mathbf{u} + \mathbf{N}\f$ + * + * \tparam MVT The matrix/vector traits collecting type information used by the iv + * \tparam size1 first size specifier for the arrays + * \tparam size2 second size specifier for the arrays + */ +template +class MatrixDataHandleBase { + using AMatrix = typename MVT::AMatrix; + using BMatrix = typename MVT::BMatrix; + using CMatrix = typename MVT::CMatrix; + using TMatrix = typename MVT::TMatrix; + using CellVector = typename MVT::CellVector; + using OutsideTij = std::vector< std::vector >; + using OmegaStorage = typename MVT::OmegaStorage; -//! Empty data handle class -class EmptyDataHandle +public: + //! Access functions to context-dependent data + const CMatrix& CA() const { return CA_[contextIdx1_][contextIdx2_]; } + CMatrix& CA() { return CA_[contextIdx1_][contextIdx2_]; } + + const AMatrix& A() const { return A_[contextIdx1_][contextIdx2_]; } + AMatrix& A() { return A_[contextIdx1_][contextIdx2_]; } + + const BMatrix& AB() const { return AB_[contextIdx1_][contextIdx2_]; } + BMatrix& AB() { return AB_[contextIdx1_][contextIdx2_]; } + + const TMatrix& T() const { return T_[contextIdx1_][contextIdx2_]; } + TMatrix& T() { return T_[contextIdx1_][contextIdx2_]; } + + const OutsideTij& tijOutside() const { return tijOutside_[contextIdx1_][contextIdx2_]; } + OutsideTij& tijOutside() { return tijOutside_[contextIdx1_][contextIdx2_]; } + + const OmegaStorage& omegas() const { return wijk_[contextIdx1_][contextIdx2_]; } + OmegaStorage& omegas() { return wijk_[contextIdx1_][contextIdx2_]; } + + //! functionality to set the context indices + void setContextIndex1(unsigned int idx) const { assert(idx < size1); contextIdx1_ = idx; } + void setContextIndex2(unsigned int idx) const { assert(idx < size2); contextIdx2_ = idx; } + +protected: + //! indices to be set before accessing data + mutable unsigned int contextIdx1_{0}; + mutable unsigned int contextIdx2_{0}; + + std::array< std::array, size1 > wijk_; //!< The omega factors that form the entries of the matrices + + std::array< std::array, size1 > T_; //!< The transmissibility matrix + std::array< std::array, size1 > A_; //!< Inverse of the iv-local system matrix + std::array< std::array, size1 > AB_; //!< A_ left multiplied to B + std::array< std::array, size1 > CA_; //!< A_ right multiplied to C + std::array< std::array, size1 > tijOutside_; //!< The transmissibilities for "outside" faces (on surface grids) +}; + +/*! + * \ingroup CCMpfaDiscretization + * \brief Common base class to all handles. Stores arrays of the + * vectors of known cell/Dirichlet values within the interaction volume. + * + * \tparam MVT The matrix/vector traits collecting type information used by the iv + * \tparam size1 first size specifier for the arrays + * \tparam size2 second size specifier for the arrays + */ +template +class VectorDataHandleBase { + using CellVector = typename MVT::CellVector; + public: - template< class InteractionVolume > - void resize(const InteractionVolume& iv) {} + //! Access to the iv-wide known cell/Dirichlet values + const CellVector& uj() const { return u_[contextIdx1_][contextIdx2_]; } + CellVector& uj() { return u_[contextIdx1_][contextIdx2_]; } + +protected: + //! functionality to set the context indices + void setContextIndex1(unsigned int idx) const { assert(idx < size1); contextIdx1_ = idx; } + void setContextIndex2(unsigned int idx) const { assert(idx < size2); contextIdx2_ = idx; } + + //! indices to be set before accessing data + mutable unsigned int contextIdx1_{0}; + mutable unsigned int contextIdx2_{0}; + + //! The interaction volume-local known values + std::array< std::array, size1 > u_; }; +} // end namespace Detail + +//! Empty data handle class +class EmptyDataHandle {}; + //! Data handle for quantities related to advection template class AdvectionDataHandle +: public Detail::MatrixDataHandleBase +, public Detail::VectorDataHandleBase { - // export matrix & vector types from interaction volume - using AMatrix = typename MatVecTraits::AMatrix; - using CMatrix = typename MatVecTraits::CMatrix; - using TMatrix = typename MatVecTraits::TMatrix; - using CellVector = typename MatVecTraits::CellVector; - using FaceVector = typename MatVecTraits::FaceVector; - using FaceScalar = typename FaceVector::value_type; - using OutsideGravityStorage = std::vector< Dune::DynamicVector >; + // we only have one local system for all phases since we + // solve them w.r.t. permeability tensor (unique for all phases) + using Base1 = Detail::MatrixDataHandleBase; + // we do have cell/Dirichlet values for all phases though! static constexpr int numPhases = PhysicsTraits::numPhases; + using Base2 = Detail::VectorDataHandleBase; + + using UnknownVector = typename MatVecTraits::AMatrix::row_type; + using FaceVector = typename MatVecTraits::FaceVector; + using FaceScalar = typename FaceVector::value_type; + using OutsideGravityStorage = std::vector< std::vector >; public: + //! Set the phase index of the context + void setPhaseIndex(unsigned int phaseIdx) const { Base2::setContextIndex1(phaseIdx); } + + //! The gravitational flux contributions for a phase on all faces + const FaceVector& g() const { return g_[Base2::contextIdx1_]; } + FaceVector& g() { return g_[Base2::contextIdx1_]; } - //! prepare data handle for subsequent fill (normal grids) - template< class IV, - std::enable_if_t<(IV::Traits::GridView::dimension==IV::Traits::GridView::dimensionworld), int> = 0> - void resize(const IV& iv) - { - // resize transmissibility matrix & pressure vectors - resizeMatrix(T_, iv.numFaces(), iv.numKnowns()); - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - resizeVector(p_[pIdx], iv.numKnowns()); - - // maybe resize gravity container - static const bool enableGravity = getParam("Problem.EnableGravity"); - if (enableGravity) - { - resizeMatrix(CA_, iv.numFaces(), iv.numUnknowns()); - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - resizeVector(g_[pIdx], iv.numFaces()); - } - } - - - //! prepare data handle for subsequent fill (surface grids) - template< class IV, - std::enable_if_t<(IV::Traits::GridView::dimension = 0> - void resize(const IV& iv) - { - using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; - - static const bool enableGravity = getParam("Problem.EnableGravity"); - if (!enableGravity) - { - resizeMatrix(T_, iv.numFaces(), iv.numKnowns()); - outsideT_.resize(iv.numFaces()); - - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - resizeVector(p_[pIdx], 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) - resizeVector(outsideT_[i][j], iv.numKnowns()); - } - } - - else - { - resizeMatrix(T_, iv.numFaces(), iv.numKnowns()); - resizeMatrix(CA_, iv.numFaces(), iv.numUnknowns()); - resizeMatrix(A_, iv.numUnknowns(), iv.numUnknowns()); - outsideT_.resize(iv.numFaces()); - - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - { - resizeVector(p_[pIdx], iv.numKnowns()); - resizeVector(g_[pIdx], 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) - resizeVector(outsideG_[pIdx][i], numNeighbors); - for (LocalIndexType j = 0; j < numNeighbors; ++j) - resizeVector(outsideT_[i][j], iv.numKnowns()); - } - } - } - - //! Access to the iv-wide pressure of one phase - const CellVector& pressures(unsigned int pIdx) const { return p_[pIdx]; } - CellVector& pressures(unsigned int pIdx) { return p_[pIdx]; } - - //! Access to the gravitational flux contributions for one phase - const FaceVector& gravity(unsigned int pIdx) const { return g_[pIdx]; } - FaceVector& gravity(unsigned int pIdx) { return g_[pIdx]; } - - //! Access to the gravitational flux contributions for all phases - const std::array< FaceVector, numPhases >& gravity() const { return g_; } - std::array< FaceVector, numPhases >& gravity() { return g_; } - - //! Projection matrix for gravitational acceleration - const CMatrix& advectionCA() const { return CA_; } - CMatrix& advectionCA() { return CA_; } - - //! Additional projection matrix needed on surface grids - const AMatrix& advectionA() const { return A_; } - AMatrix& advectionA() { return A_; } - - //! The transmissibilities associated with advective fluxes - const TMatrix& advectionT() const { return T_; } - TMatrix& advectionT() { return T_; } - - //! The transmissibilities for "outside" faces (used on surface grids) - const std::vector< std::vector >& advectionTout() const { return outsideT_; } - std::vector< std::vector >& advectionTout() { return outsideT_; } - - //! The gravitational acceleration for "outside" faces (used on surface grids) - const std::array< OutsideGravityStorage, numPhases >& gravityOutside() const { return outsideG_; } - std::array< OutsideGravityStorage, numPhases >& gravityOutside() { return outsideG_; } + //! The deltaG vector for gravity within the iv-local eq-system + const UnknownVector& deltaG() const { return deltaG_[Base2::contextIdx1_]; } + UnknownVector& deltaG() { return deltaG_[Base2::contextIdx1_]; } //! The gravitational acceleration for one phase on "outside" faces (used on surface grids) - const OutsideGravityStorage& gravityOutside(unsigned int pIdx) const { return outsideG_[pIdx]; } - OutsideGravityStorage& gravityOutside(unsigned int pIdx) { return outsideG_[pIdx]; } + const OutsideGravityStorage& gOutside() const { return outsideG_[Base2::contextIdx1_]; } + OutsideGravityStorage& gOutside() { return outsideG_[Base2::contextIdx1_]; } private: - TMatrix T_; //!< The transmissibilities such that f_i = T_ij*p_j - CMatrix CA_; //!< Matrix to project gravitational acceleration to all scvfs - AMatrix A_; //!< Matrix additionally needed for the projection on surface grids - std::array< CellVector, numPhases > p_; //!< The interaction volume-wide phase pressures - std::array< FaceVector, numPhases > g_; //!< The gravitational acceleration at each scvf (only for enabled gravity) - std::vector< std::vector > outsideT_; //!< The transmissibilities for "outside" faces (only on surface grids) + std::array< FaceVector, numPhases > g_; //!< The gravitational acceleration at each scvf (only for enabled gravity) + std::array< UnknownVector, numPhases > deltaG_; //!< The gravity coefficients forming part of iv-local eq-system std::array< OutsideGravityStorage, numPhases > outsideG_; //!< The gravitational acceleration on "outside" faces (only on surface grids) }; //! Data handle for quantities related to diffusion template class DiffusionDataHandle +: public Detail::MatrixDataHandleBase +, public Detail::VectorDataHandleBase { - using TMatrix = typename MatVecTraits::TMatrix; - using CellVector = typename MatVecTraits::CellVector; - using OutsideTContainer = std::vector< std::vector >; - static constexpr int numPhases = PhysicsTraits::numPhases; static constexpr int numComponents = PhysicsTraits::numComponents; + using Base1 = Detail::MatrixDataHandleBase; + using Base2 = Detail::VectorDataHandleBase; public: //! diffusion caches need to set phase and component index - void setPhaseIndex(unsigned int phaseIdx) { phaseIdx_ = phaseIdx; } - void setComponentIndex(unsigned int compIdx) { compIdx_ = compIdx; } - - //! prepare data handle for subsequent fill - template< class InteractionVolume > - void resize(const InteractionVolume& iv) - { - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - { - for (unsigned int cIdx = 0; cIdx < numComponents; ++cIdx) - { - // resize transmissibility matrix & mole fraction vector - resizeMatrix(T_[pIdx][cIdx], iv.numFaces(), iv.numKnowns()); - resizeVector(xj_[pIdx][cIdx], iv.numKnowns()); - - // resize outsideTij on surface grids - using GridView = typename InteractionVolume::Traits::GridView; - if (int(GridView::dimension) < int(GridView::dimensionworld)) - { - outsideT_[pIdx][cIdx].resize(iv.numFaces()); - - using LocalIndexType = typename InteractionVolume::Traits::IndexSet::LocalIndexType; - for (LocalIndexType i = 0; i < iv.numFaces(); ++i) - { - const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1; - outsideT_[pIdx][cIdx][i].resize(numNeighbors); - for (LocalIndexType j = 0; j < numNeighbors; ++j) - resizeVector(outsideT_[pIdx][cIdx][i][j], iv.numKnowns()); - } - } - } - } - } - - //! Access to the iv-wide mole fractions of a component in one phase - const CellVector& moleFractions(unsigned int pIdx, unsigned int compIdx) const { return xj_[pIdx][compIdx]; } - CellVector& moleFractions(unsigned int pIdx, unsigned int compIdx) { return xj_[pIdx][compIdx]; } - - //! The transmissibilities associated with diffusive fluxes - const TMatrix& diffusionT() const { return T_[phaseIdx_][compIdx_]; } - TMatrix& 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, numPhases > T_; //!< The transmissibilities such that f_i = T_ij*x_j - std::array< std::array, numPhases > xj_; //!< The interaction volume-wide mole fractions - std::array< std::array, numPhases> outsideT_; + void setPhaseIndex(unsigned int phaseIdx) const + { Base1::setContextIndex1(phaseIdx); Base2::setContextIndex1(phaseIdx); } + void setComponentIndex(unsigned int compIdx) const + { Base1::setContextIndex2(compIdx); Base2::setContextIndex2(compIdx); } }; //! Data handle for quantities related to heat conduction template class HeatConductionDataHandle -{ - using TMatrix = typename MatVecTraits::TMatrix; - using CellVector = typename MatVecTraits::CellVector; - -public: - //! prepare data handle for subsequent fill - template< class InteractionVolume > - void resize(const InteractionVolume& iv) - { - //! resize transmissibility matrix & temperature vector - resizeMatrix(T_, iv.numFaces(), iv.numKnowns()); - resizeVector(Tj_, iv.numKnowns()); - - //! resize outsideTij on surface grids - using GridView = typename InteractionVolume::Traits::GridView; - if (int(GridView::dimension) < int(GridView::dimensionworld)) - { - outsideT_.resize(iv.numFaces()); - - using LocalIndexType = typename InteractionVolume::Traits::IndexSet::LocalIndexType; - 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) - resizeVector(outsideT_[i][j], iv.numKnowns()); - } - } - } - - //! Access to the iv-wide temperatures - const CellVector& temperatures() const { return Tj_; } - CellVector& temperatures() { return Tj_; } - - //! The transmissibilities associated with conductive fluxes - const TMatrix& heatConductionT() const { return T_; } - TMatrix& heatConductionT() { return T_; } - - //! The transmissibilities for "outside" faces (used on surface grids) - const std::vector< std::vector >& heatConductionTout() const { return outsideT_; } - std::vector< std::vector >& heatConductionTout() { return outsideT_; } - -private: - // heat conduction-related variables - TMatrix T_; //!< The transmissibilities such that f_i = T_ij*T_j - CellVector Tj_; //!< The interaction volume-wide temperatures - std::vector< std::vector > outsideT_; //!< The transmissibilities for "outside" faces (only necessary on surface grids) -}; +: public Detail::MatrixDataHandleBase +, public Detail::VectorDataHandleBase +{}; //! Process-dependent data handles when related process is disabled template @@ -313,23 +224,31 @@ class HeatConductionDataHandle : public Empt * \tparam PT The physics traits collecting data on the physical processes to be considered */ template -class InteractionVolumeDataHandle : public AdvectionDataHandle, - public DiffusionDataHandle, - public HeatConductionDataHandle +class InteractionVolumeDataHandle { + +public: + //! export the underlying process-specific handle types using AdvectionHandle = AdvectionDataHandle; using DiffusionHandle = DiffusionDataHandle; using HeatConductionHandle = HeatConductionDataHandle; -public: - //! prepare data handles for subsequent fills - template< class InteractionVolume > - void resize(const InteractionVolume& iv) - { - AdvectionHandle::resize(iv); - DiffusionHandle::resize(iv); - HeatConductionHandle::resize(iv); - } + //! return references to the handle containing data related to advection + const AdvectionHandle& advectionHandle() const { return advectionHandle_; } + AdvectionHandle& advectionHandle() { return advectionHandle_; } + + //! return references to the handle containing data related to diffusion + const DiffusionHandle& diffusionHandle() const { return diffusionHandle_; } + DiffusionHandle& diffusionHandle() { return diffusionHandle_; } + + //! return references to the handle containing data related to heat conduction + const HeatConductionHandle& heatConductionHandle() const { return heatConductionHandle_; } + HeatConductionHandle& heatConductionHandle() { return heatConductionHandle_; } + +private: + AdvectionHandle advectionHandle_; + DiffusionHandle diffusionHandle_; + HeatConductionHandle heatConductionHandle_; }; } // end namespace Dumux diff --git a/dumux/discretization/cellcentered/mpfa/localassembler.hh b/dumux/discretization/cellcentered/mpfa/localassembler.hh deleted file mode 100644 index ad06867aee3fb0c3b0c2ca9b89eb32e9ad6e3e91..0000000000000000000000000000000000000000 --- a/dumux/discretization/cellcentered/mpfa/localassembler.hh +++ /dev/null @@ -1,280 +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 . * - *****************************************************************************/ -/*! - * \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 - -#include - -namespace Dumux -{ -//! Forward declaration of the implementation -template< class P, class EG, class EV, MpfaMethods M > class InteractionVolumeAssemblerImpl; - -//! Alias to select the right implementation. -template< class P, class EG, class EV, MpfaMethods M > -using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< P, EG, EV, M >; - -/*! - * \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 P The problem type - * \tparam EG The element finite volume geometry - * \tparam EV The element volume variables type - */ -template< class P, class EG, class EV > -class InteractionVolumeAssemblerBase -{ - using Problem = P; - using FVElementGeometry = EG; - using ElementVolumeVariables = EV; - - 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 IV Interaction volume type implementation - * \tparam TensorFunc 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 getT Lambda to evaluate the scv-wise tensors - */ - template< class IV, class TensorFunc > - void assemble(typename IV::Traits::MatVecTraits::MatVecTraits::TMatrix& T, IV& iv, const TensorFunc& getT) - { - 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 - * transmissibilities associated with "outside" faces are stored - * in a separate container. - * - * \tparam TOutside Container to store the "outside" transmissibilities - * \tparam IV Interaction volume type implementation - * \tparam TensorFunc Lambda to obtain the tensor w.r.t. - * which the local system is to be solved - * - * \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 getT Lambda to evaluate the scv-wise tensors - */ - template< class TOutside, class IV, class TensorFunc > - void assemble(TOutside& outsideTij, typename IV::Traits::MatVecTraits::TMatrix& T, IV& iv, const TensorFunc& getT) - { - 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 GC The type of container used to store the - * gravitational acceleration per scvf & phase - * \tparam IV The interaction volume type implementation - * \tparam TensorFunc 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 getT Lambda to evaluate the scv-wise tensors - */ - template< class GC, class IV, class TensorFunc > - void assembleWithGravity(typename IV::Traits::MatVecTraits::TMatrix& T, - GC& g, - typename IV::Traits::MatVecTraits::CMatrix& CA, - IV& iv, - const TensorFunc& getT) - { - 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. - * - * \tparam GC The type of container used to store the - * gravitational acceleration per scvf & phase - * \tparam GOut Type of container used to store gravity on "outside" faces - * \tparam TOutside Container to store the "outside" transmissibilities - * \tparam IV The interaction volume type implementation - * \tparam TensorFunc Lambda to obtain the tensor w.r.t. - * which the local system is to be solved - * - * \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 getT Lambda to evaluate the scv-wise tensors - */ - template< class GC, class GOut, class TOutside, class IV, class TensorFunc > - void assembleWithGravity(TOutside& outsideTij, - typename IV::Traits::MatVecTraits::TMatrix& T, - GOut& outsideG, - GC& g, - typename IV::Traits::MatVecTraits::CMatrix& CA, - typename IV::Traits::MatVecTraits::AMatrix& A, - IV& iv, - const TensorFunc& getT) - { - 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 IV The interaction volume type implementation - * \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 IV, class GetU > - void assemble(typename IV::Traits::MatVecTraits::CellVector& u, const IV& 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. - * - * \tparam GC The type of container used to store the - * gravitational acceleration per scvf & phase - * \tparam TensorFunc Lambda to obtain the tensor w.r.t. - * which the local system is to be solved - * - * \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 getT Lambda to evaluate scv-wise hydraulic conductivities - */ - template< class GC, class IV, class TensorFunc > - void assembleGravity(GC& g, const IV& iv, const typename IV::Traits::MatVecTraits::CMatrix& CA, const TensorFunc& getT) - { - 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. - * - * \tparam GC The type of container used to store the - * gravitational acceleration per scvf & phase - * \tparam GOut Type of container used to store gravity on "outside" faces - * \tparam IV The interaction volume type implementation - * \tparam TensorFunc Lambda to obtain the tensor w.r.t. - * which the local system is to be solved - * - * \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 getT Lambda to evaluate scv-wise hydraulic conductivities - */ - template< class GC, class GOut, class IV, class TensorFunc > - void assembleGravity(GC& g, - GOut& outsideG, - const IV& iv, - const typename IV::Traits::MatVecTraits::CMatrix& CA, - const typename IV::Traits::MatVecTraits::AMatrix& A, - const TensorFunc& getT) - { - 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 - -#endif diff --git a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh new file mode 100644 index 0000000000000000000000000000000000000000..56c350025deb6a37dc9c334205dd14db955ac4e1 --- /dev/null +++ b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh @@ -0,0 +1,325 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup CCMpfaDiscretization + * \brief Defines the general interface of classes used for the assembly + * of the local systems of equations involved in the transmissibility + * computation in mpfa schemes. + */ +#ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_BASE_HH +#define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_BASE_HH + +#include + +#include + +namespace Dumux +{ +/*! + * \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 P The problem type + * \tparam EG The element finite volume geometry + * \tparam EV The element volume variables type + */ +template< class P, class EG, class EV > +class InteractionVolumeAssemblerBase +{ + using Problem = P; + using FVElementGeometry = EG; + using ElementVolumeVariables = EV; + + typedef std::true_type yes; + typedef std::false_type no; + + //! Determines whether or not a matrix has a resize() function + template + struct matrix_has_resize_method + { + private: + // resize function is called with two indices for matrices + template static auto test(int) -> decltype(std::declval().resize(0, 0), yes()); + template static no test(...); + public: + static constexpr bool value = std::is_same(0)), yes>::value; + }; + + //! determines whether or not a vector has a resize() function + template + struct vector_has_resize_method + { + private: + // resize function is called with one index for vectors + template static auto test(int) -> decltype(std::declval().resize(0), yes()); + template static no test(...); + public: + static constexpr bool value = std::is_same(0)), yes>::value; + }; + + 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 Assembles the matrices involved in the flux + * expressions and the local system of equations + * within an mpfa interaction volume. + * + * \tparam IV The interaction volume type implementation + * \tparam TensorFunc Lambda to obtain the tensor w.r.t. + * which the local system is to be solved + * + * \param handle The data handle in which the matrices are stored + * \param iv The interaction volume + * \param getT Lambda to evaluate the scv-wise tensors + */ + template< class DataHandle, class IV, class TensorFunc > + void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleMatrices() function"); + } + + /*! + * \brief Assembles the vector of primary (cell) unknowns and (maybe) + * Dirichlet boundary conditions within an interaction volume. + * + * \tparam IV The interaction volume type implementation + * \tparam GetU Lambda to obtain the cell unknowns from grid indices + * + * \param handle The data handle in which the vector is stored + * \param iv The mpfa-o interaction volume + * \param getU Lambda to obtain the desired cell/Dirichlet value from vol vars + */ + template< class DataHandle, class IV, class GetU > + void assembleU(DataHandle& handle, const IV& iv, const GetU& getU) + { + DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function for the cell/Dirichlet unknowns"); + } + + /*! + * \brief Assembles the gravitational flux contributions on the scvfs within an + * interaction volume. + * + * \param handle The data handle in which the vector is stored + * \param iv The mpfa-o interaction volume + * \param getRho Lambda to obtain the density from volume variables + */ + template< class DataHandle, class IV, class GetRho > + void assembleGravity(DataHandle& handle, const IV& iv, const GetRho& getRho) + { + using GridView = typename IV::Traits::GridView; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + static constexpr bool isSurfaceGrid = dim < dimWorld; + + // resize the gravity vectors + auto& g = handle.g(); + auto& deltaG = handle.deltaG(); + auto& outsideG = handle.gOutside(); + resizeVector_(g, iv.numFaces()); + resizeVector_(deltaG, iv.numUnknowns()); + if (isSurfaceGrid) + resizeVector_(outsideG, iv.numFaces()); + + // we require the CA matrix to have the correct size already + assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns()); + + //! For each face, we... + //! - arithmetically average the phase densities + //! - compute the term \f$ \alpha := \mathbf{A} \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell + //! - compute \f$ \alpha^* = \sum{\alpha_{outside, i}} - \alpha_{inside} \f$ + using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type; + using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; + + for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) + { + // gravitational acceleration on this face + const auto& curLocalScvf = iv.localScvf(faceIdx); + const auto& curGlobalScvf = fvGeometry().scvf(curLocalScvf.gridScvfIndex()); + const auto& gravity = problem().gravityAtPos(curGlobalScvf.ipGlobal()); + + // get permeability tensor in "positive" sub volume + const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); + const auto& posGlobalScv = fvGeometry().scv(iv.localScv(neighborScvIndices[0]).gridScvIndex()); + const auto& posVolVars = elemVolVars()[posGlobalScv]; + const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), + posVolVars.permeability(), + gravity); + + const auto numOutsideFaces = !curGlobalScvf.boundary() ? curGlobalScvf.numOutsideScvs() : 0; + using OutsideAlphaStorage = std::conditional_t< isSurfaceGrid, + std::vector, + Dune::ReservedVector >; + OutsideAlphaStorage alpha_outside; alpha_outside.resize(numOutsideFaces); + std::fill(alpha_outside.begin(), alpha_outside.end(), 0.0); + Scalar rho; + + if (isSurfaceGrid) + resizeVector_(outsideG[faceIdx], numOutsideFaces); + + if (!curLocalScvf.isDirichlet()) + { + const auto localDofIdx = curLocalScvf.localDofIndex(); + + rho = getRho(posVolVars); + deltaG[localDofIdx] = 0.0; + + if (!curGlobalScvf.boundary()) + { + for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) + { + // obtain outside tensor + const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1]; + const auto& negGlobalScv = fvGeometry().scv(iv.localScv(negLocalScvIdx).gridScvIndex()); + const auto& negVolVars = elemVolVars()[negGlobalScv]; + const auto& flipScvf = !isSurfaceGrid ? curGlobalScvf + : fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside); + + alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*vtmv(flipScvf.unitOuterNormal(), + negVolVars.permeability(), + gravity); + if (isSurfaceGrid) + alpha_outside[idxInOutside] *= -1.0; + + rho += getRho(negVolVars); + deltaG[localDofIdx] += alpha_outside[idxInOutside]; + } + } + + rho /= numOutsideFaces + 1; + deltaG[localDofIdx] -= alpha_inside; + deltaG[localDofIdx] *= rho*curGlobalScvf.area(); + } + // use density resulting from Dirichlet BCs + else + rho = getRho(elemVolVars()[curGlobalScvf.outsideScvIdx()]); + + // add "inside" & "outside" alphas to gravity containers + g[faceIdx] = alpha_inside*rho*curGlobalScvf.area(); + + if (isSurfaceGrid) + { + unsigned int i = 0; + for (const auto& alpha : alpha_outside) + outsideG[faceIdx][i++] = alpha*rho*curGlobalScvf.area(); + } + } + + // add iv-wide contributions to gravity vectors + handle.CA().umv(deltaG, g); + if (isSurfaceGrid) + { + using FaceVector = typename IV::Traits::MatVecTraits::FaceVector; + FaceVector AG; + resizeVector_(AG, iv.numUnknowns()); + handle.A().mv(deltaG, AG); + + // compute gravitational accelerations + for (const auto& localFaceData : iv.localFaceData()) + { + // continue only for "outside" faces + if (!localFaceData.isOutsideFace()) + 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 = handle.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.localScvfIndex(localDir)); + if (!curLocalScvf.isDirichlet()) + outsideG[localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()]; + } + } + } + } + +protected: + //! resizes a matrix to the given sizes (specialization for dynamic matrix type) + template< class Matrix, + class size_type, + std::enable_if_t::value, int> = 0 > + void resizeMatrix_(Matrix& M, size_type rows, size_type cols) + { + M.resize(rows, cols); + } + + //! resizes a matrix to the given sizes (specialization for static matrix type - do nothing) + template< class Matrix, + class size_type, + std::enable_if_t::value, int> = 0 > + void resizeMatrix_(Matrix& M, size_type rows, size_type cols) + {} + + //! resizes a vector to the given size (specialization for dynamic matrix type) + template< class Vector, + class size_type, + std::enable_if_t::value, int> = 0 > + void resizeVector_(Vector& v, size_type size) + { + v.resize(size); + } + + //! resizes a vector to the given size (specialization for static vector type - do nothing) + template< class Vector, + class size_type, + std::enable_if_t::value, int> = 0 > + void resizeVector_(Vector& v, size_type rows) + {} + +private: + // pointers to the data required for assembly + const Problem* problemPtr_; + const FVElementGeometry* fvGeometryPtr_; + const ElementVolumeVariables* elemVolVarsPtr_; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/discretization/cellcentered/mpfa/localfacedata.hh b/dumux/discretization/cellcentered/mpfa/localfacedata.hh index 1fe6489b2deb3729fc1da0f0f70f76f8f8811a9c..8cc387d95b9afabbfd8b258e09fe138125a069df 100644 --- a/dumux/discretization/cellcentered/mpfa/localfacedata.hh +++ b/dumux/discretization/cellcentered/mpfa/localfacedata.hh @@ -32,7 +32,7 @@ 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. + * volume-local information for a grid sub-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 @@ -43,8 +43,8 @@ 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" + GridIndexType gridScvfIndex_; //!< the index of the corresponding global scvf + bool isOutside_; //!< indicates if this face is an "outside" face in the iv-local system public: //! Default constructor @@ -53,10 +53,10 @@ public: //! Constructor InteractionVolumeLocalFaceData(LocalIndexType faceIndex, LocalIndexType scvIndex, - GridIndexType globalScvfIndex) + GridIndexType gridScvfIndex) : ivLocalScvfIndex_(faceIndex) , ivLocalInsideScvIndex_(scvIndex) - , globalScvfIndex_(globalScvfIndex) + , gridScvfIndex_(gridScvfIndex) , isOutside_(false) {} @@ -64,11 +64,11 @@ public: InteractionVolumeLocalFaceData(LocalIndexType faceIndex, LocalIndexType scvIndex, LocalIndexType indexInScvfOutsideFaces, - GridIndexType globalScvfIndex) + GridIndexType gridScvfIndex) : ivLocalScvfIndex_(faceIndex) , ivLocalInsideScvIndex_(scvIndex) , scvfLocalOutsideScvfIndex_(indexInScvfOutsideFaces) - , globalScvfIndex_(globalScvfIndex) + , gridScvfIndex_(gridScvfIndex) , isOutside_(true) {} @@ -76,8 +76,8 @@ public: 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_; } + GridIndexType gridScvfIndex() const { return gridScvfIndex_; } + bool isOutsideFace() const { return isOutside_; } }; } // end namespace diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh index d5fa3020cf860ae94e957638c064228520873a11..a81919eaa3809ef8e71cb335bac1bc6d51e0d2cc 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh @@ -32,13 +32,13 @@ #include #include -#include #include #include #include #include +#include "localassembler.hh" #include "localsubcontrolentities.hh" #include "interactionvolumeindexset.hh" @@ -67,9 +67,16 @@ private: static constexpr int dim = NodalIndexSet::Traits::GridView::dimension; static constexpr int dimWorld = NodalIndexSet::Traits::GridView::dimensionworld; + using DimVector = Dune::FieldVector; + using FaceOmegas = typename std::conditional< (dim, + Dune::ReservedVector >::type; + //! Matrix/Vector traits to be used by the data handle struct MVTraits { + using OmegaStorage = std::vector< FaceOmegas >; + using AMatrix = Dune::DynamicMatrix< Scalar >; using BMatrix = Dune::DynamicMatrix< Scalar >; using CMatrix = Dune::DynamicMatrix< Scalar >; @@ -92,6 +99,10 @@ public: using LocalFaceData = InteractionVolumeLocalFaceData; //! export the matrix/vector traits to be used by the iv using MatVecTraits = MVTraits; + + //! the type of assembler used for the o-method's iv-local eq systems + template + using LocalAssembler = MpfaOInteractionVolumeAssembler; }; /*! @@ -102,7 +113,7 @@ public: */ template< class Traits > class CCMpfaOInteractionVolume - : public CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume, Traits > +: public CCMpfaInteractionVolumeBase< Traits > { using GridView = typename Traits::GridView; using Element = typename GridView::template Codim<0>::Entity; @@ -112,20 +123,6 @@ class CCMpfaOInteractionVolume using LocalIndexType = typename IndexSet::LocalIndexType; using Stencil = typename IndexSet::NodalGridStencilType; - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - - //! export scalar type from T matrix and define omegas - using Scalar = typename Traits::MatVecTraits::TMatrix::value_type; - using DimVector = Dune::FieldVector; - using FaceOmegas = typename std::conditional< (dim, - Dune::ReservedVector >::type; - - using AMatrix = typename Traits::MatVecTraits::AMatrix; - using BMatrix = typename Traits::MatVecTraits::BMatrix; - using CMatrix = typename Traits::MatVecTraits::CMatrix; - using LocalScvType = typename Traits::LocalScvType; using LocalScvfType = typename Traits::LocalScvfType; using LocalFaceData = typename Traits::LocalFaceData; @@ -144,21 +141,18 @@ public: GridIndexType volVarIndex() const { return volVarIndex_; } }; - //! export the type used for transmissibility storage - using TransmissibilityStorage = std::vector< FaceOmegas >; - //! publicly state the mpfa-scheme this interaction volume is associated with static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod; //! Sets up the local scope for a given iv index set template< class Problem, class FVElementGeometry > - void setUpLocalScope(const IndexSet& indexSet, - const Problem& problem, - const FVElementGeometry& fvGeometry) + void bind(const IndexSet& indexSet, + const Problem& problem, + const FVElementGeometry& fvGeometry) { // for the o-scheme, the stencil is equal to the scv // index set of the dual grid's nodal index set - stencil_ = &indexSet.nodalIndexSet().globalScvIndices(); + stencil_ = &indexSet.nodalIndexSet().gridScvIndices(); // number of interaction-volume-local scvs(=node-local for o-scheme) and scvfs numFaces_ = indexSet.numFaces(); @@ -166,25 +160,19 @@ public: const auto numGlobalScvfs = indexSet.nodalIndexSet().numScvfs(); // reserve memory for local entities - elements_.clear(); - scvs_.clear(); - scvfs_.clear(); - localFaceData_.clear(); - dirichletData_.clear(); - elements_.reserve(numLocalScvs); - scvs_.reserve(numLocalScvs); - scvfs_.reserve(numFaces_); - dirichletData_.reserve(numFaces_); - localFaceData_.reserve(numGlobalScvfs); + elements_.clear(); elements_.reserve(numLocalScvs); + scvs_.clear(); scvs_.reserve(numLocalScvs); + scvfs_.clear(); scvfs_.reserve(numFaces_); + localFaceData_.clear(); localFaceData_.reserve(numGlobalScvfs); + dirichletData_.clear(); dirichletData_.reserve(numFaces_); // set up stuff related to sub-control volumes - const auto& scvIndices = indexSet.globalScvIndices(); for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numLocalScvs; scvIdxLocal++) { - elements_.emplace_back(fvGeometry.fvGridGeometry().element( scvIndices[scvIdxLocal] )); + elements_.emplace_back(fvGeometry.fvGridGeometry().element( stencil()[scvIdxLocal] )); scvs_.emplace_back(fvGeometry.fvGridGeometry().mpfaHelper(), fvGeometry, - fvGeometry.scv( scvIndices[scvIdxLocal] ), + fvGeometry.scv( stencil()[scvIdxLocal] ), scvIdxLocal, indexSet); } @@ -194,27 +182,19 @@ public: numKnowns_ = numLocalScvs; // set up quantitites related to sub-control volume faces - wijk_.resize(numFaces_); for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numFaces_; ++faceIdxLocal) { - const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(faceIdxLocal)); + const auto& scvf = fvGeometry.scvf(indexSet.gridScvfIndex(faceIdxLocal)); // 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, neighborScvIndicesLocal[0], scvf.index()); - - // we will need as many omegas as scvs around the face const auto numNeighborScvs = neighborScvIndicesLocal.size(); - wijk_[faceIdxLocal].resize(numNeighborScvs); + localFaceData_.emplace_back(faceIdxLocal, neighborScvIndicesLocal[0], scvf.index()); // create iv-local scvf object if (scvf.boundary()) { - const auto bcTypes = problem.boundaryTypes(elements_[neighborScvIndicesLocal[0]], scvf); - - if (bcTypes.hasOnlyDirichlet()) + if (problem.boundaryTypes(elements_[neighborScvIndicesLocal[0]], scvf).hasOnlyDirichlet()) { scvfs_.emplace_back(scvf, neighborScvIndicesLocal, numKnowns_++, /*isDirichlet*/true); dirichletData_.emplace_back(scvf.outsideScvIdx()); @@ -231,11 +211,11 @@ public: { // 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) + for (int coord = 0; coord < GridView::dimension; ++coord) { - if (indexSet.scvfIdxLocal(outsideLocalScvIdx, coord) == faceIdxLocal) + if (indexSet.localScvfIndex(outsideLocalScvIdx, coord) == faceIdxLocal) { - const auto globalScvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord); + const auto globalScvfIdx = indexSet.nodalIndexSet().gridScvfIndex(outsideLocalScvIdx, coord); const auto& flipScvf = fvGeometry.scvf(globalScvfIdx); localFaceData_.emplace_back(faceIdxLocal, // iv-local scvf idx outsideLocalScvIdx, // iv-local scv index @@ -250,11 +230,6 @@ public: } } } - - // Maybe resize local matrices if dynamic types are used - resizeMatrix(A_, numUnknowns_, numUnknowns_); - resizeMatrix(B_, numUnknowns_, numKnowns_); - resizeMatrix(C_, numFaces_, numUnknowns_); } //! returns the number of primary scvfs of this interaction volume @@ -297,25 +272,10 @@ public: const std::vector& dirichletData() const { return dirichletData_; } - //! returns the matrix associated with face unknowns in local equation system - const AMatrix& A() const { return A_; } - AMatrix& A() { return A_; } - - //! returns the matrix associated with cell unknowns in local equation system - const BMatrix& B() const { return B_; } - BMatrix& B() { return B_; } - - //! returns the matrix associated with face unknowns in flux expressions - const CMatrix& C() const { return C_; } - CMatrix& C() { return C_; } - - //! returns container storing the transmissibilities for each face & coordinate - const TransmissibilityStorage& omegas() const { return wijk_; } - TransmissibilityStorage& omegas() { return wijk_; } - //! returns the number of interaction volumes living around a vertex template< class NI > - static constexpr std::size_t numIVAtVertex(const NI& nodalIndexSet) { return 1; } + static constexpr std::size_t numIVAtVertex(const NI& 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 @@ -335,7 +295,7 @@ public: ivIndexSetContainer.emplace_back(nodalIndexSet, flipScvfIndexSet); // store the index mapping - for (const auto scvfIdx : nodalIndexSet.globalScvfIndices()) + for (const auto scvfIdx : nodalIndexSet.gridScvfIndices()) scvfIndexMap[scvfIdx] = curGlobalIndex; } @@ -350,14 +310,6 @@ private: std::vector localFaceData_; std::vector dirichletData_; - // Matrices needed for computation of transmissibilities - AMatrix A_; - BMatrix B_; - CMatrix C_; - - // The omega factors are stored during assembly of local system - TransmissibilityStorage wijk_; - // sizes involved in the local system equations std::size_t numFaces_; std::size_t numUnknowns_; diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh index 2ab9082e1da9bf531957191ab3eab76267a16830..ca930c1651df6e9cddf154d34eeb3e48f7fc0bcf 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh @@ -78,7 +78,7 @@ public: // for scvfs touching the boundary there are no "outside" scvfs if (nodalIndexSet.scvfIsOnBoundary(i)) { - scvfNeighborScvLocalIndices_.push_back({nodalIndexSet.insideScvIdxLocal(i)}); + scvfNeighborScvLocalIndices_.push_back({nodalIndexSet.insideScvLocalIndex(i)}); nodeToIvScvf_[i] = ivToNodeScvf_.size(); isHandled[i] = true; ivToNodeScvf_.push_back(i); @@ -92,12 +92,12 @@ public: isHandled[i] = true; // construct local index sets - const auto& flipScvfIndices = flipScvfIndexSet[nodalIndexSet.scvfIdxGlobal(i)]; + const auto& flipScvfIndices = flipScvfIndexSet[nodalIndexSet.gridScvfIndex(i)]; const auto numFlipIndices = flipScvfIndices.size(); ScvfNeighborLocalIndexSet neighborsLocal; neighborsLocal.resize(numFlipIndices + 1); - neighborsLocal[0] = nodalIndexSet.insideScvIdxLocal(i); + neighborsLocal[0] = nodalIndexSet.insideScvLocalIndex(i); // mappings for all flip scvf for (unsigned int j = 0; j < numFlipIndices; ++j) @@ -105,9 +105,9 @@ public: const auto outsideScvfIdx = flipScvfIndices[j]; for (unsigned int nodeLocalIdx = 0; nodeLocalIdx < nodalIndexSet.numScvfs(); ++nodeLocalIdx) { - if (nodalIndexSet.scvfIdxGlobal(nodeLocalIdx) == outsideScvfIdx) + if (nodalIndexSet.gridScvfIndex(nodeLocalIdx) == outsideScvfIdx) { - neighborsLocal[j+1] = nodalIndexSet.insideScvIdxLocal(nodeLocalIdx); + neighborsLocal[j+1] = nodalIndexSet.insideScvLocalIndex(nodeLocalIdx); nodeToIvScvf_[nodeLocalIdx] = curIvLocalScvfIdx; isHandled[nodeLocalIdx] = true; break; // go to next outside scvf @@ -126,8 +126,8 @@ public: { return nodalIndexSet_; } //! returns the global scv indices connected to this dual grid node - const NodalGridStencilType& globalScvIndices() const - { return nodalIndexSet_.globalScvIndices(); } + const NodalGridStencilType& gridScvIndices() const + { return nodalIndexSet_.gridScvIndices(); } //! returns the number of faces in the interaction volume std::size_t numFaces() const @@ -137,18 +137,18 @@ public: std::size_t numScvs() const { return nodalIndexSet_.numScvs(); } - //! returns a global scvf idx for a given iv-local scvf index - GridIndexType scvfIdxGlobal(LocalIndexType ivLocalScvfIdx) const + //! returns a grid scvf idx for a given iv-local scvf index + GridIndexType gridScvfIndex(LocalIndexType ivLocalScvfIdx) const { assert(ivLocalScvfIdx < numFaces()); - return nodalIndexSet_.scvfIdxGlobal( ivToNodeScvf_[ivLocalScvfIdx] ); + return nodalIndexSet_.gridScvfIndex( ivToNodeScvf_[ivLocalScvfIdx] ); } //! returns the iv-local scvf idx of the i-th scvf embedded in a local scv - LocalIndexType scvfIdxLocal(LocalIndexType scvIdxLocal, unsigned int i) const + LocalIndexType localScvfIndex(LocalIndexType scvIdxLocal, unsigned int i) const { - assert(nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) < nodeToIvScvf_.size()); - return nodeToIvScvf_[ nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) ]; + assert(nodalIndexSet_.localScvfIndex(scvIdxLocal, i) < nodeToIvScvf_.size()); + return nodeToIvScvf_[ nodalIndexSet_.localScvfIndex(scvIdxLocal, i) ]; } //! returns the local indices of the neighboring scvs of an scvf diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh index c0c35e8d7dd2a026870e9aff6b50d870d54f8422..a4039645784954e00004fc1caf71c3dc756fedc6 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh @@ -25,11 +25,11 @@ #ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH #define DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH -#include -#include +#include +#include #include -#include +#include #include namespace Dumux @@ -45,242 +45,99 @@ namespace Dumux * \tparam EV The element volume variables type */ template< class P, class EG, class EV > -class InteractionVolumeAssemblerImpl< P, EG, EV, MpfaMethods::oMethod > - : public InteractionVolumeAssemblerBase< P, EG, EV > +class MpfaOInteractionVolumeAssembler +: public InteractionVolumeAssemblerBase< P, EG, EV > { using ParentType = InteractionVolumeAssemblerBase< P, EG, EV >; public: - //! Use the constructor of the base class + //! Pull up constructor of the base class using ParentType::ParentType; /*! - * \brief Assembles the transmissibility matrix within an - * interaction volume in an mpfa-o type way. + * \brief Assembles the matrices involved in the flux + * expressions and the local system of equations + * within an interaction volume in an mpfa-o type way. * * \tparam IV The interaction volume type implementation * \tparam TensorFunc 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 getT Lambda to evaluate the scv-wise tensors - */ - template< class IV, class TensorFunc > - void assemble(typename IV::Traits::MatVecTraits::TMatrix& T, IV& iv, const TensorFunc& getT) - { - // assemble D into T directly - assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getT); - - // 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. - * - * \tparam TOutside Container to store the "outside" transmissibilities - * \tparam IV The interaction volume type implementation - * \tparam TensorFunc Lambda to obtain the tensor w.r.t. * which the local system is to be solved * - * \param outsideTij tij on "outside" faces to be assembled - * \param T The transmissibility matrix tij to be assembled + * \param handle The data handle in which the matrices are stored * \param iv The interaction volume * \param getT Lambda to evaluate the scv-wise tensors */ - template< class TOutside, class IV, class TensorFunc > - void assemble(TOutside& outsideTij, typename IV::Traits::MatVecTraits::TMatrix& T, IV& iv, const TensorFunc& getT) + template< class DataHandle, class IV, class TensorFunc > + void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT) { - // assemble D into T directly - assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getT); + assembleLocalMatrices_(handle.A(), handle.AB(), handle.CA(), handle.T(), handle.omegas(), iv, getT); // 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()) + // T = C*(A^-1)*B + D + handle.A().invert(); + handle.CA().rightmultiply(handle.A()); + handle.T() += multiplyMatrices(handle.CA(), handle.AB()); + handle.AB().leftmultiply(handle.A()); + + // On surface grids, compute the "outside" transmissibilities + using GridView = typename IV::Traits::GridView; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + if (dim < dimWorld) { - // 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 + // bring outside tij container to the right size + auto& tijOut = handle.tijOutside(); + tijOut.resize(iv.numFaces()); using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; - for (LocalIndexType localDir = 0; localDir < IV::Traits::GridView::dimension; localDir++) + for (LocalIndexType fIdx = 0; fIdx < iv.numFaces(); ++fIdx) { - // 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]; + const auto& curGlobalScvf = this->fvGeometry().scvf(iv.localScvf(fIdx).gridScvfIndex()); + const auto numOutsideFaces = curGlobalScvf.boundary() ? 0 : curGlobalScvf.numOutsideScvs(); + // resize each face entry to the right number of outside faces + tijOut[fIdx].resize(numOutsideFaces); + std::for_each(tijOut[fIdx].begin(), + tijOut[fIdx].end(), + [&](auto& v) { this->resizeVector_(v, iv.numKnowns()); }); } - } - } - } - - /*! - * \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. - * - * \tparam GC The type of container used to store the - * gravitational acceleration per scvf & phase - * \tparam IV The interaction volume type implementation - * \tparam TensorFunc 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 mpfa-o interaction volume - * \param getT Lambda to evaluate the scv-wise tensors - */ - template< class GC, class IV, class TensorFunc > - void assembleWithGravity(typename IV::Traits::MatVecTraits::TMatrix& T, - GC& g, - typename IV::Traits::MatVecTraits::CMatrix& CA, - IV& iv, - const TensorFunc& getT) - { - // assemble D into T & C into CA directly - assembleLocalMatrices_(iv.A(), iv.B(), CA, T, iv, getT); - // 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()); - } + // compute outside transmissibilities + for (const auto& localFaceData : iv.localFaceData()) + { + // continue only for "outside" faces + if (!localFaceData.isOutsideFace()) + continue; - // assemble gravitational acceleration container (enforce usage of mpfa-o type version) - assembleGravity(g, iv, CA, getT); - } + const auto scvIdx = localFaceData.ivLocalInsideScvIndex(); + const auto scvfIdx = localFaceData.ivLocalScvfIndex(); + const auto idxInOut = localFaceData.scvfLocalOutsideScvfIndex(); - /*! - * \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. - * - * \tparam GC The type of container used to store the - * gravitational acceleration per scvf & phase - * \tparam GOut Type of container used to store gravity on "outside" faces - * \tparam TOutside Container to store the "outside" transmissibilities - * \tparam IV The interaction volume type implementation - * \tparam TensorFunc Lambda to obtain the tensor w.r.t. - * which the local system is to be solved - * - * \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 getT Lambda to evaluate the scv-wise tensors - */ - template< class GC, class GOut, class TOutside, class IV, class TensorFunc > - void assembleWithGravity(TOutside& outsideTij, - typename IV::Traits::MatVecTraits::TMatrix& T, - GOut& outsideG, - GC& g, - typename IV::Traits::MatVecTraits::CMatrix& CA, - typename IV::Traits::MatVecTraits::AMatrix& A, - IV& iv, - const TensorFunc& getT) - { - // assemble D into T directly - assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getT); + const auto& wijk = handle.omegas()[scvfIdx][idxInOut+1]; + auto& tij = tijOut[scvfIdx][idxInOut]; - // 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; + tij = 0.0; + for (unsigned int dir = 0; dir < dim; dir++) + { + // the scvf corresponding to this local direction in the scv + const auto& scvf = iv.localScvf(iv.localScv(scvIdx).localScvfIndex(dir)); - // add contributions from all local directions - using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; - for (LocalIndexType localDir = 0; localDir < IV::Traits::GridView::dimension; 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 (!scvf.isDirichlet()) + { + auto tmp = handle.AB()[scvf.localDofIndex()]; + tmp *= wijk[dir]; + tij -= tmp; + } + else + tij[scvf.localDofIndex()] -= wijk[dir]; - // 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; + // add entry from the scv unknown + tij[scvIdx] += wijk[dir]; } - else - tij[curLocalScvf.localDofIndex()] -= wijk[localDir]; - - // add entry from the scv unknown - tij[localScvIdx] += wijk[localDir]; } } } - - assembleGravity(g, outsideG, iv, CA, A, getT); } /*! @@ -290,320 +147,22 @@ public: * \tparam IV The interaction volume type implementation * \tparam GetU Lambda to obtain the cell unknowns from grid indices * - * \param u The vector to be filled with the cell unknowns + * \param handle The data handle in which the vector is stored * \param iv The mpfa-o interaction volume - * \param getU Lambda to obtain the desired cell/Dirichlet value from grid index + * \param getU Lambda to obtain the desired cell/Dirichlet value from vol vars */ - template< class IV, class GetU > - void assemble(typename IV::Traits::MatVecTraits::CellVector& u, const IV& iv, const GetU& getU) + template< class DataHandle, class IV, class GetU > + void assembleU(DataHandle& handle, const IV& iv, const GetU& getU) { - // put the cell pressures first - using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; - for (LocalIndexType i = 0; i < iv.numScvs(); ++i) - u[i] = getU( iv.localScv(i).globalScvIndex() ); + auto& u = handle.uj(); + this->resizeVector_(u, iv.numKnowns()); - // Dirichlet BCs come afterwards - LocalIndexType i = iv.numScvs(); + // put the cell unknowns first, then Dirichlet values + typename IV::Traits::IndexSet::LocalIndexType i = 0; + for (; i < iv.numScvs(); i++) + u[i] = getU( this->elemVolVars()[iv.localScv(i).gridScvIndex()] ); 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. - * - * \tparam GC The type of container used to store the - * gravitational acceleration per scvf & phase - * \tparam TensorFunc Lambda to obtain the tensor w.r.t. - * which the local system is to be solved - * - * \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 getT Lambda to evaluate scv-wise hydraulic conductivities - */ - template< class GC, class IV, class TensorFunc > - void assembleGravity(GC& g, - const IV& iv, - const typename IV::Traits::MatVecTraits::CMatrix& CA, - const TensorFunc& getT) - { - //! We require the gravity container to be a two-dimensional vector/array type, structured as follows: - //! - first index adresses the respective phases - //! - second index adresses the face within the interaction volume - - // make sure g vector and CA matrix have the correct sizes already - assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); })); - assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns()); - - //! For each face, we... - //! - arithmetically average the phase densities - //! - compute the term \f$ \alpha := A \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell - //! - compute \f$ \alpha^* = \alpha_{outside} - \alpha_{inside} \f$ - using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type; - using FaceVector = typename IV::Traits::MatVecTraits::FaceVector; - using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; - - // reset gravity containers to zero - const auto numPhases = g.size(); - std::vector< FaceVector > sum_alphas(numPhases); - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - { - resizeVector(sum_alphas[pIdx], iv.numUnknowns()); - sum_alphas[pIdx] = 0.0; - g[pIdx] = 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 = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv); - - // On surface grids one should use the function specialization below - assert(neighborScvIndices.size() <= 2); - - std::vector< Scalar > rho(numPhases); - 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 = getT(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. - * - * \tparam GC The type of container used to store the - * gravitational acceleration per scvf & phase - * \tparam GOut Type of container used to store gravity on "outside" faces - * \tparam IV The interaction volume type implementation - * \tparam TensorFunc Lambda to obtain the tensor w.r.t. - * which the local system is to be solved - * - * \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 getT Lambda to evaluate scv-wise hydraulic conductivities - */ - template< class GC, class GOut, class IV, class TensorFunc > - void assembleGravity(GC& g, - GOut& outsideG, - const IV& iv, - const typename IV::Traits::MatVecTraits::CMatrix& CA, - const typename IV::Traits::MatVecTraits::AMatrix& A, - const TensorFunc& getT) - { - //! We require the gravity container to be a two-dimensional vector/array type, structured as follows: - //! - first index adresses the respective phases - //! - second index adresses the face within the interaction volume - //! We require the outside gravity container to be a three-dimensional vector/array type, structured as follows: - //! - first index adresses the respective phases - //! - second index adresses the face within the interaction volume - //! - third index adresses the i-th "outside" face of the current face - - // we require the CA matrix and the gravity containers to have the correct size already - assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns()); - assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); })); - assert(std::all_of(outsideG.begin(), outsideG.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); })); - - //! For each face, we... - //! - arithmetically average the phase densities - //! - compute the term \f$ \alpha := \mathbf{A} \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell - //! - compute \f$ \alpha^* = \sum{\alpha_{outside, i}} - \alpha_{inside} \f$ - using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type; - using FaceVector = typename IV::Traits::MatVecTraits::FaceVector; - using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; - - // reset everything to zero - const auto numPhases = g.size(); - std::vector< FaceVector > sum_alphas(numPhases); - for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx) - { - resizeVector(sum_alphas[pIdx], iv.numUnknowns()); - sum_alphas[pIdx] = 0.0; - g[pIdx] = 0.0; - std::for_each(outsideG[pIdx].begin(), outsideG[pIdx].end(), [] (auto& v) { v = 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 = getT(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::vector< Scalar > rho(numPhases); - - 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 = getT(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]); - - FaceVector 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]; - - // make sure the given outside gravity container has the right size - assert(outsideG[pIdx][localScvfIdx].size() == iv.localScvf(localScvfIdx).neighboringLocalScvIndices().size()-1); - - // add contributions from all local directions - for (LocalIndexType localDir = 0; localDir < IV::Traits::GridView::dimension; 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()]; - } - } - } + u[i++] = getU( this->elemVolVars()[data.volVarIndex()] ); } private: @@ -635,85 +194,82 @@ private: typename IV::Traits::MatVecTraits::BMatrix& B, typename IV::Traits::MatVecTraits::CMatrix& C, typename IV::Traits::MatVecTraits::DMatrix& D, + typename IV::Traits::MatVecTraits::OmegaStorage& wijk, IV& iv, const TensorFunc& getT) { using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; static constexpr int dim = IV::Traits::GridView::dimension; static constexpr int dimWorld = IV::Traits::GridView::dimensionworld; - // Matrix D is assumed to have the right size already - assert(D.rows() == iv.numFaces() && D.cols() == iv.numKnowns()); + // resize omegas + this->resizeVector_(wijk, iv.numFaces()); // 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; + // resize & reset D matrix + this->resizeMatrix_(D, iv.numFaces(), iv.numKnowns()); 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& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); 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& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex()); const auto& posVolVars = this->elemVolVars()[posGlobalScv]; const auto& posElement = iv.element(neighborScvIndices[0]); const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv); // the omega factors of the "positive" sub volume - const auto wijk = computeMpfaTransmissibility(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); + this->resizeVector_(wijk[faceIdx], /*no outside scvs present*/1); + wijk[faceIdx][0] = 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& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) ); const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); - D[faceIdx][otherLocalDofIdx] -= wijk[localDir]; - D[faceIdx][posScvLocalDofIdx] += wijk[localDir]; + D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; + D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; } } } else { - // we require the matrices A,B,C to have the correct size already - assert(A.rows() == iv.numUnknowns() && A.cols() == iv.numUnknowns()); - assert(B.rows() == iv.numUnknowns() && B.cols() == iv.numKnowns()); - assert(C.rows() == iv.numFaces() && C.cols() == iv.numUnknowns()); - - // reset matrices - A = 0.0; - B = 0.0; - C = 0.0; - D = 0.0; - - auto& wijk = iv.omegas(); + // resize & reset matrices + this->resizeMatrix_(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0; + this->resizeMatrix_(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0; + this->resizeMatrix_(C, iv.numFaces(), iv.numUnknowns()); C = 0.0; + this->resizeMatrix_(D, iv.numFaces(), iv.numKnowns()); D = 0.0; + for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) { const auto& curLocalScvf = iv.localScvf(faceIdx); - const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.globalScvfIndex()); + const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); 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& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex()); const auto& posVolVars = this->elemVolVars()[posGlobalScv]; const auto& posElement = iv.element(neighborScvIndices[0]); const auto tensor = getT(this->problem(), posElement, posVolVars, this->fvGeometry(), posGlobalScv); // the omega factors of the "positive" sub volume + this->resizeVector_(wijk[faceIdx], neighborScvIndices.size()); 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& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) ); const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); // if we are not on a Dirichlet face, add entries associated with unknown face pressures @@ -748,7 +304,7 @@ private: { const auto idxOnScvf = idxInOutside+1; const auto& negLocalScv = iv.localScv( neighborScvIndices[idxOnScvf] ); - const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.globalScvIndex()); + const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.gridScvIndex()); const auto& negVolVars = this->elemVolVars()[negGlobalScv]; const auto& negElement = iv.element( neighborScvIndices[idxOnScvf] ); const auto negTensor = getT(this->problem(), negElement, negVolVars, this->fvGeometry(), negGlobalScv); @@ -762,10 +318,10 @@ private: if (dim < dimWorld) wijk[faceIdx][idxOnScvf] *= -1.0; - // go over the coordinate directions in the positive sub volume + // go over the coordinate directions in the negative sub volume for (int localDir = 0; localDir < dim; localDir++) { - const auto otherLocalScvfIdx = negLocalScv.scvfIdxLocal(localDir); + const auto otherLocalScvfIdx = negLocalScv.localScvfIndex(localDir); const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx); const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh index 32f2775f43c9ee69e5e298ece62bcfdcb7870716..6c492a60a4cf4b965d692274f24bb90f135182e1 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh @@ -82,7 +82,7 @@ public: LocalBasis localBasis; for (unsigned int coordIdx = 0; coordIdx < myDimension; ++coordIdx) { - const auto scvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(localDofIndex_, coordIdx); + const auto scvfIdx = indexSet.nodalIndexSet().gridScvfIndex(localDofIndex_, coordIdx); const auto& scvf = fvGeometry.scvf(scvfIdx); localBasis[coordIdx] = scvf.ipGlobal(); localBasis[coordIdx] -= center; @@ -96,8 +96,8 @@ public: ctype detX() const { return detX_; } - //! grid view-global index related to this scv - GridIndexType globalScvIndex() const + //! grid index related to this scv + GridIndexType gridScvIndex() const { return globalScvIndex_; } //! returns the index in the set of cell unknowns of the iv @@ -105,10 +105,10 @@ public: { return localDofIndex_; } //! iv-local index of the coordir's scvf in this scv - LocalIndexType scvfIdxLocal(unsigned int coordDir) const + LocalIndexType localScvfIndex(unsigned int coordDir) const { assert(coordDir < myDimension); - return indexSet_->scvfIdxLocal(localDofIndex_, coordDir); + return indexSet_->localScvfIndex(localDofIndex_, coordDir); } //! the nu vectors are needed for setting up the omegas of the iv @@ -170,7 +170,7 @@ public: LocalIndexType localDofIndex() const { return localDofIndex_; } //! returns the grid view-global index of this scvf - GridIndexType globalScvfIndex() const { return scvfIdxGlobal_; } + GridIndexType gridScvfIndex() const { return scvfIdxGlobal_; } //! Returns the local indices of the scvs neighboring this scvf const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices() const { return *neighborScvIndicesLocal_; } diff --git a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh index dec48f6ba8af3e9cb05617c066071b5bf5cb5e16..14703fc1524f65268500bb5747544ce2879968c2 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh @@ -30,7 +30,6 @@ #include #include -#include #include #include @@ -38,6 +37,7 @@ #include #include +#include "localassembler.hh" #include "localsubcontrolentities.hh" #include "interactionvolumeindexset.hh" @@ -68,9 +68,14 @@ private: static constexpr int dim = NI::Traits::GridView::dimension; static constexpr int dimWorld = NI::Traits::GridView::dimensionworld; + using DimVector = Dune::FieldVector; + using FaceOmegas = Dune::ReservedVector; + //! Matrix/Vector traits to be used by the data handle struct MVTraits { + using OmegaStorage = std::array< FaceOmegas, F >; + using AMatrix = Dune::FieldMatrix< S, F, F >; using BMatrix = Dune::FieldMatrix< S, F, C >; using CMatrix = Dune::FieldMatrix< S, F, F >; @@ -97,6 +102,10 @@ public: static constexpr int numScvs = C; //! export the number of scvfs in the interaction volumes static constexpr int numScvfs = F; + + //! the type of assembler used for the o-method's iv-local eq systems + template + using LocalAssembler = MpfaOInteractionVolumeAssembler; }; /*! @@ -111,7 +120,7 @@ public: */ template< class Traits > class CCMpfaOStaticInteractionVolume - : public CCMpfaInteractionVolumeBase< CCMpfaOStaticInteractionVolume, Traits > +: public CCMpfaInteractionVolumeBase< Traits > { using GridView = typename Traits::GridView; using Element = typename GridView::template Codim<0>::Entity; @@ -121,18 +130,6 @@ class CCMpfaOStaticInteractionVolume using LocalIndexType = typename IndexSet::LocalIndexType; using Stencil = typename IndexSet::NodalGridStencilType; - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - - //! export scalar type from T matrix and define omegas - using Scalar = typename Traits::MatVecTraits::TMatrix::value_type; - using DimVector = Dune::FieldVector; - using FaceOmegas = typename Dune::ReservedVector; - - using AMatrix = typename Traits::MatVecTraits::AMatrix; - using BMatrix = typename Traits::MatVecTraits::BMatrix; - using CMatrix = typename Traits::MatVecTraits::CMatrix; - using LocalScvType = typename Traits::LocalScvType; using LocalScvfType = typename Traits::LocalScvfType; using LocalFaceData = typename Traits::LocalFaceData; @@ -141,34 +138,33 @@ class CCMpfaOStaticInteractionVolume static constexpr int numScv = Traits::numScvs; public: + //! This does not work on surface grids + static_assert(int(GridView::dimension)==int(GridView::dimensionworld), "static iv does not work on surface grids"); + //! export the standard o-methods dirichlet data using DirichletData = typename CCMpfaOInteractionVolume< Traits >::DirichletData; - //! export the type used for transmissibility storage - using TransmissibilityStorage = std::array< FaceOmegas, numScvf >; - //! publicly state the mpfa-scheme this interaction volume is associated with static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod; //! Sets up the local scope for a given iv index set template< class Problem, class FVElementGeometry > - void setUpLocalScope(const IndexSet& indexSet, - const Problem& problem, - const FVElementGeometry& fvGeometry) + void bind(const IndexSet& indexSet, + const Problem& problem, + const FVElementGeometry& fvGeometry) { // for the o-scheme, the stencil is equal to the scv // index set of the dual grid's nodal index set assert(indexSet.numScvs() == numScv); - stencil_ = &indexSet.nodalIndexSet().globalScvIndices(); + stencil_ = &indexSet.nodalIndexSet().gridScvIndices(); // set up stuff related to sub-control volumes - const auto& scvIndices = indexSet.globalScvIndices(); for (LocalIndexType scvIdxLocal = 0; scvIdxLocal < numScv; scvIdxLocal++) { - elements_[scvIdxLocal] = fvGeometry.fvGridGeometry().element( scvIndices[scvIdxLocal] ); + elements_[scvIdxLocal] = fvGeometry.fvGridGeometry().element( stencil()[scvIdxLocal] ); scvs_[scvIdxLocal] = LocalScvType(fvGeometry.fvGridGeometry().mpfaHelper(), fvGeometry, - fvGeometry.scv( scvIndices[scvIdxLocal] ), + fvGeometry.scv( stencil()[scvIdxLocal] ), scvIdxLocal, indexSet); } @@ -176,26 +172,23 @@ public: // set up quantitites related to sub-control volume faces for (LocalIndexType faceIdxLocal = 0; faceIdxLocal < numScvf; ++faceIdxLocal) { - const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(faceIdxLocal)); + const auto& scvf = fvGeometry.scvf(indexSet.gridScvfIndex(faceIdxLocal)); + assert(!scvf.boundary()); // the neighboring scvs in local indices (order: 0 - inside scv, 1..n - outside scvs) const auto& neighborScvIndicesLocal = indexSet.neighboringLocalScvIndices(faceIdxLocal); - // this does not work for network grids or on boundaries - assert(!scvf.boundary()); - assert(neighborScvIndicesLocal.size() == 2); - // create iv-local scvf objects scvfs_[faceIdxLocal] = LocalScvfType(scvf, neighborScvIndicesLocal, faceIdxLocal, /*isDirichlet*/false); localFaceData_[faceIdxLocal*2] = LocalFaceData(faceIdxLocal, neighborScvIndicesLocal[0], scvf.index()); // add local face data objects for the outside face const auto outsideLocalScvIdx = neighborScvIndicesLocal[1]; - for (int coord = 0; coord < dim; ++coord) + for (int coord = 0; coord < GridView::dimension; ++coord) { - if (indexSet.scvfIdxLocal(outsideLocalScvIdx, coord) == faceIdxLocal) + if (indexSet.localScvfIndex(outsideLocalScvIdx, coord) == faceIdxLocal) { - const auto globalScvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord); + const auto globalScvfIdx = indexSet.nodalIndexSet().gridScvfIndex(outsideLocalScvIdx, coord); const auto& flipScvf = fvGeometry.scvf(globalScvfIdx); localFaceData_[faceIdxLocal*2+1] = LocalFaceData(faceIdxLocal, // iv-local scvf idx outsideLocalScvIdx, // iv-local scv index @@ -208,11 +201,6 @@ public: // make sure we found it assert(localFaceData_[faceIdxLocal*2+1].ivLocalInsideScvIndex() == outsideLocalScvIdx); } - - // Maybe resize local matrices if dynamic types are used - resizeMatrix(A_, numScvf, numScvf); - resizeMatrix(B_, numScvf, numScv); - resizeMatrix(C_, numScvf, numScv); } //! returns the number of primary scvfs of this interaction volume @@ -256,22 +244,6 @@ public: const std::array& dirichletData() const { return dirichletData_; } - //! returns the matrix associated with face unknowns in local equation system - const AMatrix& A() const { return A_; } - AMatrix& A() { return A_; } - - //! returns the matrix associated with cell unknowns in local equation system - const BMatrix& B() const { return B_; } - BMatrix& B() { return B_; } - - //! returns the matrix associated with face unknowns in flux expressions - const CMatrix& C() const { return C_; } - CMatrix& C() { return C_; } - - //! returns container storing the transmissibilities for each face & coordinate - const TransmissibilityStorage& omegas() const { return wijk_; } - TransmissibilityStorage& omegas() { return wijk_; } - //! returns the number of interaction volumes living around a vertex template< class NI > static constexpr std::size_t numIVAtVertex(const NI& nodalIndexSet) @@ -288,15 +260,11 @@ public: const NodalIndexSet& nodalIndexSet, const FlipScvfIndexSet& flipScvfIndexSet) { - // the global index of the iv index set that is about to be created - const auto curGlobalIndex = ivIndexSetContainer.size(); - - // make the one index set for this node - ivIndexSetContainer.emplace_back(nodalIndexSet, flipScvfIndexSet); - - // store the index mapping - for (const auto scvfIdx : nodalIndexSet.globalScvfIndices()) - scvfIndexMap[scvfIdx] = curGlobalIndex; + // reuse the standard o-method's implementation of this + CCMpfaOInteractionVolume::addIVIndexSets(ivIndexSetContainer, + scvfIndexMap, + nodalIndexSet, + flipScvfIndexSet); } private: @@ -309,14 +277,6 @@ private: std::array scvfs_; std::array localFaceData_; - // Matrices needed for computation of transmissibilities - AMatrix A_; - BMatrix B_; - CMatrix C_; - - // The omega factors are stored during assembly of local system - TransmissibilityStorage wijk_; - // Dummy dirichlet data container (compatibility with dynamic o-iv) std::array dirichletData_; }; diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh index 9741f71c65fdcb923d4505d11d578be6adc1f46a..44180aba17a4cf5b986ddce8404b2b8dd7e53066 100644 --- a/dumux/discretization/cellcentered/mpfa/properties.hh +++ b/dumux/discretization/cellcentered/mpfa/properties.hh @@ -73,8 +73,14 @@ public: template struct PrimaryInteractionVolume { +private: + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet); + + // use the default traits + using Traits = CCMpfaODefaultInteractionVolumeTraits< NodalIndexSet, Scalar >; public: - using type = GetPropType; + using type = CCMpfaOInteractionVolume< Traits >; }; //! Per default, we use the dynamic mpfa-o interaction volume on boundaries diff --git a/dumux/discretization/fluxstencil.hh b/dumux/discretization/fluxstencil.hh index fb3e859234172184a6eccd4608e361c135ddb99c..c601e6f1cc761e8d28ca9a70e2192255d8decad3 100644 --- a/dumux/discretization/fluxstencil.hh +++ b/dumux/discretization/fluxstencil.hh @@ -117,9 +117,9 @@ public: // return the scv (element) indices in the interaction region if (fvGridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) - return fvGridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf).globalScvIndices(); + return fvGridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf).gridScvIndices(); else - return fvGridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf).globalScvIndices(); + return fvGridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf).gridScvIndices(); } }; diff --git a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh index 219394d476c1a4815dc8841278597dcb48f28e72..808936346297de3921b7e688303b9f45e902e380 100644 --- a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh +++ b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh @@ -126,41 +126,20 @@ public: / curElemVolVars[scvf.insideScvIdx()].viscosity(); const auto& fluxVarsCache = elemFluxVarsCache[scvf]; - const auto& stencil = fluxVarsCache.advectionStencil(); - if (fluxVarsCache.usesSecondaryIv()) - { - const auto& tij = fluxVarsCache.advectionTijSecondaryIv(); - - // We assume same the tij are order as the stencil up to stencil.size() - // any contribution of Dirichlet BCs is assumed to be placed afterwards - assert(stencil.size() <= tij.size()); + const auto localFaceIdx = fluxVarsCache.ivLocalFaceIndex(); + const auto usesSecondary = fluxVarsCache.usesSecondaryIv(); + const auto switchSign = fluxVarsCache.advectionSwitchFluxSign(); - // add partial derivatives to the respective given matrices - for (unsigned int i = 0; i < stencil.size();++i) - { - if (fluxVarsCache.advectionSwitchFluxSign()) - derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] -= tij[i]*up; - else - derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] += tij[i]*up; - } - } - else - { - const auto& tij = fluxVarsCache.advectionTijPrimaryIv(); - - // We assume same the tij are order as the stencil up to stencil.size() - // any contribution of Dirichlet BCs is assumed to be placed afterwards - assert(stencil.size() <= tij.size()); - - // add partial derivatives to the respective given matrices - for (unsigned int i = 0; i < stencil.size();++i) - { - if (fluxVarsCache.advectionSwitchFluxSign()) - derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] -= tij[i]*up; - else - derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] += tij[i]*up; - } - } + const auto& stencil = fluxVarsCache.advectionStencil(); + const auto& tij = usesSecondary ? fluxVarsCache.advectionSecondaryDataHandle().T()[localFaceIdx] + : fluxVarsCache.advectionPrimaryDataHandle().T()[localFaceIdx]; + + // We assume same the tij are order as the stencil up to stencil.size() + // any contribution of Dirichlet BCs is assumed to be placed afterwards + assert(stencil.size() <= tij.size()); + for (unsigned int i = 0; i < stencil.size();++i) + derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] += switchSign ? -tij[i]*up + : tij[i]*up; } //! flux derivatives for the box scheme diff --git a/dumux/porousmediumflow/fluxvariablescache.hh b/dumux/porousmediumflow/fluxvariablescache.hh index f35e0e691ca53c0d6b6375bec36d088d8c016450..82170bc1482833beaacd2e52605c333714dc676b 100644 --- a/dumux/porousmediumflow/fluxvariablescache.hh +++ b/dumux/porousmediumflow/fluxvariablescache.hh @@ -103,27 +103,32 @@ public: template< bool doSecondary = considerSecondary, std::enable_if_t = 0> bool usesSecondaryIv() const { return usesSecondaryIv_; } + //! Returns the index of the iv (this scvf is embedded in) in its container + GridIndexType ivIndexInContainer() const { return ivIndexInContainer_; } + //! Returns interaction volume-local face index + unsigned int ivLocalFaceIndex() const { return ivLocalFaceIdx_; } + //! Returns index of the face among "outside" faces of iv-local "positive" face + unsigned int indexInOutsideFaces() const { return idxInOutsideFaces_; } + //! Sets the update status. When set to true, consecutive updates will be skipped void setUpdateStatus(bool status) { isUpdated_ = status; } - - //! Sets if this cache is associated witha secondary iv + //! Sets if this cache is associated with a secondary iv void setSecondaryIvUsage(bool status) { usesSecondaryIv_ = status; } - //! Sets the index of the iv (this scvf is embedded in) in its container void setIvIndexInContainer(GridIndexType ivIndex) { ivIndexInContainer_ = ivIndex; } - - //! Returns the index of the iv (this scvf is embedded in) in its container - GridIndexType ivIndexInContainer() const { return ivIndexInContainer_; } + //! Sets the iv-local face index + void setIvLocalFaceIndex(unsigned int idx) { ivLocalFaceIdx_ = idx; } + //! Sets the index of the face among the "positive" face's outside scvfs + void setIndexInOutsideFaces(unsigned int idx) { idxInOutsideFaces_ = idx; } private: - //! indicates if cache has been fully updated - bool isUpdated_ = false; - //! indicates if cache is associated with secondary interaction volume - bool usesSecondaryIv_ = false; + bool isUpdated_ = false; //!< returns true if cache has been fully updated + bool usesSecondaryIv_ = false; //!< returns true if scvf is embedded in secondary interaction volume - //! the index of the iv (this scvf is embedded in) in its container - GridIndexType ivIndexInContainer_; + GridIndexType ivIndexInContainer_; //!< index of the iv (this scvf is embedded in) in its container + unsigned int ivLocalFaceIdx_; //!< the interaction volume-local face index of this scvf + unsigned int idxInOutsideFaces_; //!< index of scvf among outside scvfs of iv-local "positive" face (only surface grids) }; } // end namespace Dumux diff --git a/test/porousmediumflow/2p/implicit/fracture/CMakeLists.txt b/test/porousmediumflow/2p/implicit/fracture/CMakeLists.txt index 2ad42d586bc44dcb5d21dd5b784a196a26da7fd8..4724143f3eb1bd1c22d07ad880c0c3d813ad0b4b 100644 --- a/test/porousmediumflow/2p/implicit/fracture/CMakeLists.txt +++ b/test/porousmediumflow/2p/implicit/fracture/CMakeLists.txt @@ -27,7 +27,7 @@ dune_add_test(NAME test_2p_fracture_mpfa COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py CMD_ARGS --script fuzzy --files ${CMAKE_SOURCE_DIR}/test/references/test_2p_fracture_mpfa-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_mpfa-00031.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_mpfa-00030.vtu --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_fracture_mpfa params.input -Problem.Name test_2p_fracture_mpfa") # tests with gravity diff --git a/test/porousmediumflow/2p2c/implicit/injection/problem.hh b/test/porousmediumflow/2p2c/implicit/injection/problem.hh index cf5466effbd6a24a4140a4ed53ad34f53b278766..9f94940af5df5c54305ba98770a5d2d302a36277 100644 --- a/test/porousmediumflow/2p2c/implicit/injection/problem.hh +++ b/test/porousmediumflow/2p2c/implicit/injection/problem.hh @@ -26,6 +26,7 @@ #include +#include #include #include #include @@ -90,6 +91,24 @@ template struct EnableGridVolumeVariablesCache { static constexpr bool value = ENABLECACHING; }; template struct EnableGridFluxVariablesCache { static constexpr bool value = ENABLECACHING; }; + +// use the static interaction volume around interior vertices in the mpfa test +template +struct PrimaryInteractionVolume +{ +private: + using Scalar = GetPropType; + using NodalIndexSet = GetPropType; + + // structured two-d grid + static constexpr int numIvScvs = 4; + static constexpr int numIvScvfs = 4; + + // use the default traits + using Traits = CCMpfaODefaultStaticInteractionVolumeTraits< NodalIndexSet, Scalar, numIvScvs, numIvScvfs >; +public: + using type = CCMpfaOStaticInteractionVolume< Traits >; +}; } // end namespace Properties /*!