diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt
index 9743710a3e7e12c825c9032a76af3e3b3416acf6..445c4932df251ead8bd512527abc8b4a9c52093c 100644
--- a/dumux/common/CMakeLists.txt
+++ b/dumux/common/CMakeLists.txt
@@ -20,6 +20,7 @@ intersectionmapper.hh
intrange.hh
loggingparametertree.hh
math.hh
+matrixvectorhelper.hh
optional.hh
parameters.hh
pointsource.hh
diff --git a/dumux/common/matrixvectorhelper.hh b/dumux/common/matrixvectorhelper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..44769ee23f09a9b64508bec4eb4ea8d4b14c874b
--- /dev/null
+++ b/dumux/common/matrixvectorhelper.hh
@@ -0,0 +1,97 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup 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 <type_traits>
+
+namespace Dumux {
+
+//! Determines whether or not a matrix has a resize() function
+template<typename M>
+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<typename U> static auto test(int) -> decltype(std::declval<U>().resize(0, 0), yes());
+ template<typename> static no test(...);
+
+public:
+ static constexpr bool value = std::is_same<decltype(test<M>(0)), yes>::value;
+};
+
+//! determines whether or not a vector has a resize() function
+template<typename V>
+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<typename U> static auto test(int) -> decltype(std::declval<U>().resize(0), yes());
+ template<typename> static no test(...);
+
+public:
+ static constexpr bool value = std::is_same<decltype(test<V>(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<matrix_has_resize_method<Matrix>::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<!matrix_has_resize_method<Matrix>::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<vector_has_resize_method<Vector>::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<!vector_has_resize_method<Vector>::value, int> = 0 >
+void resizeVector(Vector& v, size_type rows) {}
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 9f4943bff4f4c77fda79aef01c08391addbbf696..0d7f2a7be9c275cea7faf65ef00b5d8562e03ea2 100644
--- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
@@ -44,7 +44,7 @@ class DarcysLawImplementation;
* \ingroup CCMpfaDiscretization
* \brief Darcy's law for cell-centered finite volume schemes with multi-point flux approximation.
*/
-template <class TypeTag>
+template<class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
@@ -93,73 +93,55 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- using GridIndexType = typename GridView::IndexSet::IndexType;
- // In the current implementation of the flux variables cache we cannot
- // make a disctinction between dynamic (mpfa-o) and static (mpfa-l)
- // matrix and vector types, as currently the cache class can only be templated
- // by a type tag (and there can only be one). We use a dynamic vector here to
- // make sure it works in case one of the two used interaction volume types uses
- // dynamic types performance is thus lowered for schemes using static types.
- // TODO: this has to be checked thoroughly as soon as a scheme using static types
- // is implemented. One idea to overcome the performance drop could be only
- // storing the iv-local index here and obtain tij always from the datahandle
- // of the fluxVarsCacheContainer
- using Vector = Dune::DynamicVector< Scalar >;
- using Matrix = Dune::DynamicMatrix< Scalar >;
- using Stencil = std::vector< GridIndexType >;
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
+ // In the current implementation of the flux variables cache we cannot make a
+ // disctinction between dynamic (e.g. mpfa-o unstructured) and static (e.g.mpfa-l)
+ // matrix and vector types, as currently the cache class can only be templated
+ // by a type tag (and there can only be one). Currently, pointers to both the
+ // primary and secondary iv data is stored. Before accessing it has to be checked
+ // whether or not the scvf is embedded in a secondary interaction volume.
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+ using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
+ using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
- using PrimaryStencil = typename PrimaryInteractionVolume::Traits::Stencil;
-
- static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
- "The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
- static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
- "The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );
- static_assert( std::is_convertible<PrimaryStencil*, Stencil*>::value,
- "The stencil type used in primary interaction volumes is not convertible to std::vector<GridIndexType>!" );
+ using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
+ using PrimaryIvStencil = typename PrimaryInteractionVolume::Traits::Stencil;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
+ using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
+ using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
- using SecondaryStencil = typename SecondaryInteractionVolume::Traits::Stencil;
-
- static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
- "The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
- static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
- "The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
- static_assert( std::is_convertible<SecondaryStencil*, Stencil*>::value,
- "The stencil type used in secondary interaction volumes is not convertible to std::vector<GridIndexType>!" );
+ using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
+ using SecondaryIvStencil = typename SecondaryInteractionVolume::Traits::Stencil;
public:
// export the filler type
using Filler = MpfaDarcysLawCacheFiller;
/*!
- * \brief Update cached objects (transmissibilities and gravity)
- *
- * \tparam InteractionVolume The (mpfa scheme-specific) interaction volume
- * \tparam LocalFaceData The object used to store iv-local info on an scvf
- * \tparam DataHandle The object used to store transmissibility matrices etc.
+ * \brief Update cached objects (transmissibilities and gravity).
+ * This is used for updates with primary 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< class InteractionVolume, class LocalFaceData, class DataHandle >
- void updateAdvection(const InteractionVolume& iv,
- const LocalFaceData& localFaceData,
- const DataHandle& dataHandle,
+ void updateAdvection(const PrimaryInteractionVolume& iv,
+ const PrimaryIvLocalFaceData& localFaceData,
+ const PrimaryIvDataHandle& dataHandle,
const SubControlVolumeFace &scvf)
{
- stencil_ = &iv.stencil();
switchFluxSign_ = localFaceData.isOutside();
+ primaryStencil_ = &iv.stencil();
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
- pj_[pIdx] = &dataHandle.pressures(pIdx);
+ primaryPj_[pIdx] = &dataHandle.pressures(pIdx);
static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
@@ -167,20 +149,17 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
// standard grids
if (dim == dimWorld)
{
- Tij_ = &dataHandle.advectionT()[ivLocalIdx];
-
+ 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())
{
- Tij_ = &dataHandle.advectionT()[ivLocalIdx];
-
+ primaryTij_ = &dataHandle.advectionT()[ivLocalIdx];
if (enableGravity)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx];
@@ -188,8 +167,60 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
else
{
const auto idxInOutsideFaces = localFaceData.scvfLocalOutsideScvfIndex();
- Tij_ = &dataHandle.advectionTout()[ivLocalIdx][idxInOutsideFaces];
+ 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 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<doSecondary, int > = 0 >
+ void updateAdvection(const SecondaryInteractionVolume& iv,
+ const SecondaryIvLocalFaceData& localFaceData,
+ const SecondaryIvDataHandle& dataHandle,
+ const SubControlVolumeFace &scvf)
+ {
+ switchFluxSign_ = localFaceData.isOutside();
+ secondaryStencil_ = &iv.stencil();
+
+ for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
+ secondaryPj_[pIdx] = &dataHandle.pressures(pIdx);
+
+ static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "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];
@@ -197,14 +228,23 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
}
}
- //! The stencil corresponding to the transmissibilities
- const Stencil& advectionStencil() const { return *stencil_; }
+ //! The stencil corresponding to the transmissibilities (primary type)
+ const PrimaryIvStencil& advectionStencilPrimaryIv() const { return *primaryStencil_; }
- //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
- const Vector& advectionTij() const { return *Tij_; }
+ //! The stencil corresponding to the transmissibilities (secondary type)
+ const SecondaryIvStencil& advectionStencilSecondaryIv() const { return *secondaryStencil_; }
- //! The cell (& Dirichlet) pressures within this interaction volume
- const Vector& pressures(unsigned int phaseIdx) const { return *pj_[phaseIdx]; }
+ //! 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 (& Dirichlet) pressures within this interaction volume (primary type)
+ const PrimaryIvVector& pressuresPrimaryIv(unsigned int phaseIdx) const { return *primaryPj_[phaseIdx]; }
+
+ //! The cell (& Dirichlet) pressures within this interaction volume (secondary type)
+ const SecondaryIvVector& pressuresSecondaryIv(unsigned int phaseIdx) const { return *secondaryPj_[phaseIdx]; }
//! The gravitational acceleration for a phase on this scvf
Scalar gravity(unsigned int phaseIdx) const { return g_[phaseIdx]; }
@@ -217,10 +257,21 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
private:
bool switchFluxSign_;
- const Stencil* stencil_; //!< The stencil, i.e. the grid indices j
- const Vector* Tij_; //!< The transmissibilities such that f = Tij*pj
- std::array< Scalar, numPhases > g_; //!< Gravitational flux contribution on this face
- std::array<const Vector*, numPhases> pj_; //!< The interaction-volume wide phase pressures pj
+
+ //! The stencil, i.e. the grid indices j
+ const PrimaryIvStencil* primaryStencil_;
+ const SecondaryIvStencil* secondaryStencil_;
+
+ //! The transmissibilities such that f = Tij*pj
+ const PrimaryIvTij* primaryTij_;
+ const SecondaryIvTij* secondaryTij_;
+
+ //! The interaction-volume wide phase pressures pj
+ std::array<const PrimaryIvVector*, numPhases> primaryPj_;
+ std::array<const SecondaryIvVector*, numPhases> secondaryPj_;
+
+ //! Gravitational flux contribution on this face
+ std::array< Scalar, numPhases > g_;
};
public:
@@ -240,11 +291,21 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache)
{
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
- const auto& tij = fluxVarsCache.advectionTij();
- const auto& pj = fluxVarsCache.pressures(phaseIdx);
// compute t_ij*p_j
- Scalar scvfFlux = tij*pj;
+ Scalar scvfFlux;
+ if (fluxVarsCache.usesSecondaryIv())
+ {
+ const auto& tij = fluxVarsCache.advectionTijSecondaryIv();
+ const auto& pj = fluxVarsCache.pressuresSecondaryIv(phaseIdx);
+ scvfFlux = tij*pj;
+ }
+ else
+ {
+ const auto& tij = fluxVarsCache.advectionTijPrimaryIv();
+ const auto& pj = fluxVarsCache.pressuresPrimaryIv(phaseIdx);
+ scvfFlux = tij*pj;
+ }
// maybe add gravitational acceleration
static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
index 313cf34659b7ecd614d60f39fd8e517a13672483..52d4a4568db07b73574042055b9865e479e404be 100644
--- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
@@ -123,7 +123,7 @@ class CCMpfaElementFluxVariablesCache<TypeTag, false>
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
- using SecondaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
+ using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
public:
CCMpfaElementFluxVariablesCache(const GridFluxVariablesCache& global)
diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
index 6269dd5ca43d74fe8bfef2bc84ffc33fd6722a98..cc4ae3707f2d2a6d314aeecaaf38a0a0b731eb1c 100644
--- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
@@ -194,7 +194,7 @@ public:
// Update boundary volume variables in the neighbors
for (const auto& scvf : scvfs(fvGeometry))
{
- if (fvGridGeometry.vertexUsesPrimaryInteractionVolume(scvf.vertexIndex()))
+ if (!fvGridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
{
const auto& nodalIndexSet = gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet();
// if present, insert boundary vol vars
@@ -298,7 +298,7 @@ private:
std::size_t numBoundaryVolVars = 0;
for (const auto& scvf : scvfs(fvGeometry))
{
- if (fvGridGeometry.vertexUsesPrimaryInteractionVolume(scvf.vertexIndex()))
+ if (!fvGridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
numBoundaryVolVars += gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet().numBoundaryScvfs();
else
numBoundaryVolVars += gridIvIndexSets.secondaryIndexSet(scvf).nodalIndexSet().numBoundaryScvfs();
diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
index d3b2d0cc857850c97ccaf277d159ecaf6ccce5d5..be10d54f81548e4e79bf883b6cffde8bf49ad9c2 100644
--- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
@@ -90,44 +90,30 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
//! The cache used in conjunction with the mpfa Fick's Law
class MpfaFicksLawCache
{
- // In the current implementation of the flux variables cache we cannot
- // make a disctinction between dynamic (mpfa-o) and static (mpfa-l)
- // matrix and vector types, as currently the cache class can only be templated
- // by a type tag (and there can only be one). We use a dynamic vector here to
- // make sure it works in case one of the two used interaction volume types uses
- // dynamic types performance is thus lowered for schemes using static types.
- // TODO: this has to be checked thoroughly as soon as a scheme using static types
- // is implemented. One idea to overcome the performance drop could be only
- // storing the iv-local index here and obtain tij always from the datahandle
- // of the fluxVarsCacheContainer
- using GridIndexType = typename GridView::IndexSet::IndexType;
- using Vector = Dune::DynamicVector< Scalar >;
- using Matrix = Dune::DynamicMatrix< Scalar >;
- using Stencil = std::vector< GridIndexType >;
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
+ // In the current implementation of the flux variables cache we cannot make a
+ // disctinction between dynamic (e.g. mpfa-o unstructured) and static (e.g.mpfa-l)
+ // matrix and vector types, as currently the cache class can only be templated
+ // by a type tag (and there can only be one). Currently, pointers to both the
+ // primary and secondary iv data is stored. Before accessing it has to be checked
+ // whether or not the scvf is embedded in a secondary interaction volume.
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+ using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
+ using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
- using PrimaryStencil = typename PrimaryInteractionVolume::Traits::Stencil;
-
- static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
- "The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
- static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
- "The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );
- static_assert( std::is_convertible<PrimaryStencil*, Stencil*>::value,
- "The stencil type used in primary interaction volumes is not convertible to std::vector<GridIndexType>!" );
+ using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
+ using PrimaryIvStencil = typename PrimaryInteractionVolume::Traits::Stencil;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
+ using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
+ using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
- using SecondaryStencil = typename SecondaryInteractionVolume::Traits::Stencil;
-
- static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
- "The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
- static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
- "The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
- static_assert( std::is_convertible<SecondaryStencil*, Stencil*>::value,
- "The stencil type used in secondary interaction volumes is not convertible to std::vector<GridIndexType>!" );
+ using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
+ using SecondaryIvStencil = typename SecondaryInteractionVolume::Traits::Stencil;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
@@ -138,36 +124,62 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using Filler = MpfaFicksLawCacheFiller;
/*!
- * \brief Update cached objects (transmissibilities)
+ * \brief Update cached objects (transmissibilities).
+ * This is used for updates with primary interaction volumes.
*
- * \tparam InteractionVolume The (mpfa scheme-specific) interaction volume
- * \tparam LocalFaceData The object used to store iv-local info on an scvf
- * \tparam DataHandle The object used to store transmissibility matrices etc.
+ * \param iv The interaction volume this scvf is embedded in
+ * \param localFaceData iv-local info on this scvf
+ * \param dataHandle Transmissibility matrix & gravity data of this iv
+ * \param scvf The sub-control volume face
+ */
+ void updateDiffusion(const PrimaryInteractionVolume& iv,
+ const PrimaryIvLocalFaceData& localFaceData,
+ const PrimaryIvDataHandle& dataHandle,
+ const SubControlVolumeFace &scvf,
+ unsigned int phaseIdx, unsigned int compIdx)
+ {
+ primaryStencil_[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];
+ }
+
+ /*!
+ * \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<class InteractionVolume, class LocalFaceData, class DataHandle>
- void updateDiffusion(const InteractionVolume& iv,
- const LocalFaceData& localFaceData,
- const DataHandle& dataHandle,
+ template< bool doSecondary = considerSecondaryIVs, std::enable_if_t<doSecondary, int > = 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();
+ secondaryStencil_[phaseIdx][compIdx] = &iv.stencil();
switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside();
// store pointer to the mole fraction vector of this iv
- xj_[phaseIdx][compIdx] = &dataHandle.moleFractions(phaseIdx, compIdx);
+ secondaryXj_[phaseIdx][compIdx] = &dataHandle.moleFractions(phaseIdx, compIdx);
const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
if (dim == dimWorld)
- Tij_[phaseIdx][compIdx] = &dataHandle.diffusionT()[ivLocalIdx];
+ secondaryTij_[phaseIdx][compIdx] = &dataHandle.diffusionT()[ivLocalIdx];
else
- Tij_[phaseIdx][compIdx] = localFaceData.isOutside() ? &dataHandle.diffusionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()]
- : &dataHandle.diffusionT()[ivLocalIdx];
+ secondaryTij_[phaseIdx][compIdx] = localFaceData.isOutside() ? &dataHandle.diffusionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()]
+ : &dataHandle.diffusionT()[ivLocalIdx];
}
//! In the interaction volume-local system of eq we have one unknown per face.
@@ -177,23 +189,44 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
bool diffusionSwitchFluxSign(unsigned int phaseIdx, unsigned int compIdx) const
{ return switchFluxSign_[phaseIdx][compIdx]; }
- //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
- const Vector& diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const
- { return *Tij_[phaseIdx][compIdx]; }
+ //! 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 (& Dirichlet) mole fractions within this interaction volume
- const Vector& moleFractions(unsigned int phaseIdx, unsigned int compIdx) const
- { return *xj_[phaseIdx][compIdx]; }
+ //! The cell (& Dirichlet) mole fractions within this interaction volume (primary type)
+ const PrimaryIvVector& moleFractionsPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const
+ { return *primaryXj_[phaseIdx][compIdx]; }
- //! The stencils corresponding to the transmissibilities
- const Stencil& diffusionStencil(unsigned int phaseIdx, unsigned int compIdx) const
- { return *stencil_[phaseIdx][compIdx]; }
+ //! The cell (& Dirichlet) mole fractions within this interaction volume (secondary type)
+ const SecondaryIvVector& moleFractionsSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
+ { return *secondaryXj_[phaseIdx][compIdx]; }
+
+ //! The stencils corresponding to the transmissibilities (primary type)
+ const PrimaryIvStencil& diffusionStencilPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const
+ { return *primaryStencil_[phaseIdx][compIdx]; }
+
+ //! The stencils corresponding to the transmissibilities (secondary type)
+ const SecondaryIvStencil& diffusionStencilSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
+ { return *secondaryStencil_[phaseIdx][compIdx]; }
private:
std::array< std::array<bool, numComponents>, numPhases > switchFluxSign_;
- std::array< std::array<const Stencil*, numComponents>, numPhases > stencil_; //!< The stencils, i.e. the grid indices j
- std::array< std::array<const Vector*, numComponents>, numPhases > Tij_; //!< The transmissibilities such that f = Tij*xj
- std::array< std::array<const Vector*, numComponents>, numPhases > xj_; //!< The interaction-volume wide mole fractions xj
+
+ //! The stencils, i.e. the grid indices j
+ std::array< std::array<const PrimaryIvStencil*, numComponents>, numPhases > primaryStencil_;
+ std::array< std::array<const SecondaryIvStencil*, numComponents>, numPhases > secondaryStencil_;
+
+ //! The transmissibilities such that f = Tij*xj
+ std::array< std::array<const PrimaryIvVector*, numComponents>, numPhases > primaryTij_;
+ std::array< std::array<const SecondaryIvVector*, numComponents>, numPhases > secondaryTij_;
+
+ //! The interaction-volume wide mole fractions xj
+ std::array< std::array<const PrimaryIvVector*, numComponents>, numPhases > primaryXj_;
+ std::array< std::array<const SecondaryIvVector*, numComponents>, numPhases > secondaryXj_;
};
public:
@@ -212,6 +245,9 @@ public:
const int phaseIdx,
const ElementFluxVariablesCache& elemFluxVarsCache)
{
+ // obtain this scvf's cache
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+
ComponentFluxVector componentFlux(0.0);
for (int compIdx = 0; compIdx < numComponents; compIdx++)
{
@@ -228,13 +264,21 @@ public:
// calculate the density at the interface
const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
- // prepare computations
- const auto& fluxVarsCache = elemFluxVarsCache[scvf];
- const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx);
- const auto& xj = fluxVarsCache.moleFractions(phaseIdx, compIdx);
-
// calculate Tij*xj
- Scalar flux = tij*xj;
+ Scalar flux;
+ if (fluxVarsCache.usesSecondaryIv())
+ {
+ const auto& tij = fluxVarsCache.diffusionTijSecondaryIv(phaseIdx, compIdx);
+ const auto& xj = fluxVarsCache.moleFractionsSecondaryIv(phaseIdx, compIdx);
+ flux = tij*xj;
+ }
+ 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;
diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index 452d0f6883b512869aea2cd05b1e4239bc03f622..f7172553cee3c69c95a8db8a3047dc3a4ff5de9a 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -43,6 +43,7 @@ template<class TypeTag>
class CCMpfaFluxVariablesCacheFiller
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
@@ -208,6 +209,10 @@ private:
unsigned int ivIndexInContainer,
bool forceUpdateAll = false)
{
+ // determine if secondary interaction volumes are used here
+ static constexpr bool isSecondary = MpfaHelper::considerSecondaryIVs()
+ && std::is_same<InteractionVolume, SecondaryInteractionVolume>::value;
+
// First we upate data which are not dependent on the physical processes.
// We store pointers to the other flux var caches, so that we have to obtain
// this data only once and can use it again in the sub-cache fillers.
@@ -224,6 +229,7 @@ private:
ivFluxVarCaches[i] = &fluxVarsCacheContainer[scvfJ];
ivFluxVarCaches[i]->setIvIndexInContainer(ivIndexInContainer);
ivFluxVarCaches[i]->setUpdateStatus(true);
+ ivFluxVarCaches[i]->setSecondaryIvUsage(isSecondary);
i++;
}
@@ -253,7 +259,8 @@ private:
if (AdvectionMethod == DiscretizationMethods::CCMpfa)
{
// get instance of the interaction volume-local assembler
- using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >;
+ using IVTraits = typename InteractionVolume::Traits;
+ using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
// Use different assembly if gravity is enabled
@@ -384,7 +391,8 @@ private:
static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
// get instance of the interaction volume-local assembler
- using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >;
+ using IVTraits = typename InteractionVolume::Traits;
+ using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
@@ -486,7 +494,8 @@ private:
if (HeatConductionMethod == DiscretizationMethods::CCMpfa)
{
// get instance of the interaction volume-local assembler
- using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >;
+ using IVTraits = typename InteractionVolume::Traits;
+ using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
if (forceUpdateAll || soldependentAdvection)
diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
index 1c769cf1461ee4dfa385e1d64a2471cd58b96ecf..11e2a4accd03b1fa82106cc6efae2618ee8650fa 100644
--- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
@@ -88,44 +88,30 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
//! The cache used in conjunction with the mpfa Fourier's Law
class MpfaFouriersLawCache
{
- // In the current implementation of the flux variables cache we cannot
- // make a disctinction between dynamic (mpfa-o) and static (mpfa-l)
- // matrix and vector types, as currently the cache class can only be templated
- // by a type tag (and there can only be one). We use a dynamic vector here to
- // make sure it works in case one of the two used interaction volume types uses
- // dynamic types performance is thus lowered for schemes using static types.
- // TODO: this has to be checked thoroughly as soon as a scheme using static types
- // is implemented. One idea to overcome the performance drop could be only
- // storing the iv-local index here and obtain tij always from the datahandle
- // of the fluxVarsCacheContainer
- using GridIndexType = typename GridView::IndexSet::IndexType;
- using Vector = Dune::DynamicVector< Scalar >;
- using Matrix = Dune::DynamicMatrix< Scalar >;
- using Stencil = std::vector< GridIndexType >;
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
+ // In the current implementation of the flux variables cache we cannot make a
+ // disctinction between dynamic (e.g. mpfa-o unstructured) and static (e.g.mpfa-l)
+ // matrix and vector types, as currently the cache class can only be templated
+ // by a type tag (and there can only be one). Currently, pointers to both the
+ // primary and secondary iv data is stored. Before accessing it has to be checked
+ // whether or not the scvf is embedded in a secondary interaction volume.
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+ using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
+ using PrimaryIvDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
- using PrimaryStencil = typename PrimaryInteractionVolume::Traits::Stencil;
-
- static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
- "The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
- static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
- "The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );
- static_assert( std::is_convertible<PrimaryStencil*, Stencil*>::value,
- "The stencil type used in primary interaction volumes is not convertible to std::vector<GridIndexType>!" );
+ using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
+ using PrimaryIvStencil = typename PrimaryInteractionVolume::Traits::Stencil;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
+ using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
+ using SecondaryIvDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
- using SecondaryStencil = typename SecondaryInteractionVolume::Traits::Stencil;
-
- static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
- "The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
- static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
- "The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
- static_assert( std::is_convertible<SecondaryStencil*, Stencil*>::value,
- "The stencil type used in secondary interaction volumes is not convertible to std::vector<GridIndexType>!" );
+ using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
+ using SecondaryIvStencil = typename SecondaryInteractionVolume::Traits::Stencil;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
@@ -135,42 +121,79 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using Filler = MpfaFouriersLawCacheFiller;
/*!
- * \brief Update cached objects (transmissibilities)
+ * \brief Update cached objects (transmissibilities).
+ * This is used for updates with primary interaction volumes.
*
- * \tparam InteractionVolume The (mpfa scheme-specific) interaction volume
- * \tparam LocalFaceData The object used to store iv-local info on an scvf
- * \tparam DataHandle The object used to store transmissibility matrices etc.
+ * \param iv The interaction volume this scvf is embedded in
+ * \param localFaceData iv-local info on this scvf
+ * \param dataHandle Transmissibility matrix & gravity data of this iv
+ * \param scvf The sub-control volume face
+ */
+ void updateHeatConduction(const PrimaryInteractionVolume& iv,
+ const PrimaryIvLocalFaceData& localFaceData,
+ const PrimaryIvDataHandle& dataHandle,
+ const SubControlVolumeFace &scvf)
+ {
+ primaryStencil_ = &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];
+ }
+
+ /*!
+ * \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<class InteractionVolume, class LocalFaceData, class DataHandle>
- void updateHeatConduction(const InteractionVolume& iv,
- const LocalFaceData& localFaceData,
- const DataHandle& dataHandle,
+ template< bool doSecondary = considerSecondaryIVs, std::enable_if_t<doSecondary, int > = 0 >
+ void updateHeatConduction(const SecondaryInteractionVolume& iv,
+ const SecondaryIvLocalFaceData& localFaceData,
+ const SecondaryIvDataHandle& dataHandle,
const SubControlVolumeFace &scvf)
{
- stencil_ = &iv.stencil();
+ secondaryStencil_ = &iv.stencil();
switchFluxSign_ = localFaceData.isOutside();
// store pointer to the temperature vector of this iv
- Tj_ = &dataHandle.temperatures();
+ secondaryTj_ = &dataHandle.temperatures();
const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
if (dim == dimWorld)
- Tij_ = &dataHandle.heatConductionT()[ivLocalIdx];
+ secondaryTij_ = &dataHandle.heatConductionT()[ivLocalIdx];
else
- Tij_ = localFaceData.isOutside() ? &dataHandle.heatConductionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()]
- : &dataHandle.heatConductionT()[ivLocalIdx];
+ secondaryTij_ = localFaceData.isOutside() ? &dataHandle.heatConductionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()]
+ : &dataHandle.heatConductionT()[ivLocalIdx];
}
- //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
- const Vector& heatConductionTij() const { return *Tij_; }
+ //! 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
- const Stencil& heatConductionStencil() const { return *stencil_; }
+ //! The stencil corresponding to the transmissibilities (primary type)
+ const PrimaryIvStencil& heatConductionStencilPrimaryIv() const { return *primaryStencil_; }
+
+ //! The stencil corresponding to the transmissibilities (secondary type)
+ const SecondaryIvStencil& heatConductionStencilSecondaryIv() const { return *secondaryStencil_; }
+
+ //! The cell (& Dirichlet) temperatures within this interaction volume (primary type)
+ const PrimaryIvVector& temperaturesPrimaryIv() const { return *primaryTj_; }
+
+ //! The cell (& Dirichlet) temperatures within this interaction volume (secondary type)
+ const SecondaryIvVector& temperaturesSecondaryIv() const { return *secondaryTj_; }
//! 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
@@ -178,14 +201,20 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
//! This function returns whether or not this scvf is an "outside" face in the iv.
bool heatConductionSwitchFluxSign() const { return switchFluxSign_; }
- //! The cell (& Dirichlet) temperatures within this interaction volume
- const Vector& temperatures() const { return *Tj_; }
-
private:
bool switchFluxSign_;
- const Stencil* stencil_; //!< The stencil, i.e. the grid indices j
- const Vector* Tij_; //!< The transmissibilities such that f = Tij*Tj
- const Vector* Tj_; //!< The interaction-volume wide temperature Tj
+
+ //! The stencil, i.e. the grid indices j
+ const PrimaryIvStencil* primaryStencil_;
+ const SecondaryIvStencil* secondaryStencil_;
+
+ //! The transmissibilities such that f = Tij*Tj
+ const PrimaryIvTij* primaryTij_;
+ const SecondaryIvTij* secondaryTij_;
+
+ //! The interaction-volume wide temperature Tj
+ const PrimaryIvVector* primaryTj_;
+ const SecondaryIvVector* secondaryTj_;
};
public:
@@ -204,11 +233,22 @@ public:
const ElementFluxVarsCache& elemFluxVarsCache)
{
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
- const auto& tij = fluxVarsCache.heatConductionTij();
- const auto& Tj = fluxVarsCache.temperatures();
// compute Tij*tj
- Scalar flux = tij*Tj;
+ Scalar flux;
+ if (fluxVarsCache.usesSecondaryIv())
+ {
+ const auto& tij = fluxVarsCache.heatConductionTijSecondaryIv();
+ const auto& Tj = fluxVarsCache.temperaturesSecondaryIv();
+ flux = tij*Tj;
+ }
+ else
+ {
+ const auto& tij = fluxVarsCache.heatConductionTijPrimaryIv();
+ const auto& Tj = fluxVarsCache.temperaturesPrimaryIv();
+ flux = tij*Tj;
+ }
+
if (fluxVarsCache.heatConductionSwitchFluxSign())
flux *= -1.0;
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index 6570faf63e0fefc9bdc36783619684bc6250d49a..1900565f2e1773126ff75095ce3117ee49bd9f11 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -124,22 +124,17 @@ public:
//! Get an element from a sub control volume contained in it
Element element(const SubControlVolume& scv) const { return this->elementMap()[scv.elementIndex()]; }
- //! Returns true if primary interaction volumes are used around a given vertex.
- bool vertexUsesPrimaryInteractionVolume(const Vertex& v) const
- { return primaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
-
- //!Returns true if primary interaction volumes are used around a vertex (index).
- bool vertexUsesPrimaryInteractionVolume(GridIndexType vIdxGlobal) const
- { return primaryInteractionVolumeVertices_[vIdxGlobal]; }
-
- //! Returns if primary interaction volumes are used around a given vertex.
- bool vertexUsesSecondaryInteractionVolume(const Vertex& v) const
- { return secondaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
-
- //! Returns true if primary interaction volumes are used around a given vertex index.
+ //! Returns true if secondary interaction volumes are used around a given vertex (index).
+ //! This specialization is enabled if the use of secondary interaction volumes is active.
+ template<bool useSecondary = MpfaHelper::considerSecondaryIVs(), std::enable_if_t<useSecondary, int> = 0>
bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const
{ return secondaryInteractionVolumeVertices_[vIdxGlobal]; }
+ //! 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<bool useSecondary = MpfaHelper::considerSecondaryIVs(), std::enable_if_t<!useSecondary, int> = 0>
+ constexpr bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const { return false; }
+
//! update all fvElementGeometries (do this again after grid adaption)
void update()
{
@@ -163,7 +158,6 @@ public:
// Some methods require to use a second type of interaction volume, e.g.
// around vertices on the boundary or branching points (surface grids)
- primaryInteractionVolumeVertices_.resize(numVert, true);
secondaryInteractionVolumeVertices_.resize(numVert, false);
// find vertices on processor boundaries
@@ -245,10 +239,7 @@ public:
// if this vertex is tagged to use the secondary IVs, store info
if (usesSecondaryIV)
- {
- primaryInteractionVolumeVertices_[vIdxGlobal] = false;
secondaryInteractionVolumeVertices_[vIdxGlobal] = true;
- }
// the quadrature point parameterizarion to be used on scvfs
static const Scalar q = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Mpfa.Q");
@@ -381,7 +372,6 @@ private:
// containers storing the global data
std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
- std::vector<bool> primaryInteractionVolumeVertices_;
std::vector<bool> secondaryInteractionVolumeVertices_;
std::size_t numBoundaryScvf_;
@@ -472,22 +462,17 @@ public:
//! Gets an element from a sub control volume contained in it.
Element element(const SubControlVolume& scv) const { return this->elementMap()[scv.elementIndex()]; }
- //! Returns true if primary interaction volumes are used around a given vertex.
- bool vertexUsesPrimaryInteractionVolume(const Vertex& v) const
- { return primaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
-
- //! Returns true if primary interaction volumes are used around a given vertex (index).
- bool vertexUsesPrimaryInteractionVolume(GridIndexType vIdxGlobal) const
- { return primaryInteractionVolumeVertices_[vIdxGlobal]; }
-
- //! Returns if primary interaction volumes are used around a given vertex.
- bool vertexUsesSecondaryInteractionVolume(const Vertex& v) const
- { return secondaryInteractionVolumeVertices_[this->vertexMapper().index(v)]; }
-
- //! Returns true if primary interaction volumes are used around a given vertex (index).
+ //! Returns true if secondary interaction volumes are used around a given vertex (index).
+ //! This specialization is enabled if the use of secondary interaction volumes is active.
+ template<bool useSecondary = MpfaHelper::considerSecondaryIVs(), std::enable_if_t<useSecondary, int> = 0>
bool vertexUsesSecondaryInteractionVolume(GridIndexType vIdxGlobal) const
{ return secondaryInteractionVolumeVertices_[vIdxGlobal]; }
+ //! 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<bool useSecondary = MpfaHelper::considerSecondaryIVs(), std::enable_if_t<!useSecondary, int> = 0>
+ 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 { return isGhostVertex_[this->vertexMapper().index(v)]; }
@@ -510,7 +495,6 @@ public:
// Some methods require to use a second type of interaction volume, e.g.
// around vertices on the boundary or branching points (surface grids)
const auto numVert = this->gridView().size(dim);
- primaryInteractionVolumeVertices_.resize(numVert, true);
secondaryInteractionVolumeVertices_.resize(numVert, false);
// find vertices on processor boundaries HERE!!
@@ -588,10 +572,7 @@ public:
// if this vertex is tagged to use the secondary IVs, store info
if (usesSecondaryIV)
- {
- primaryInteractionVolumeVertices_[vIdxGlobal] = false;
secondaryInteractionVolumeVertices_[vIdxGlobal] = true;
- }
// make the scv face (for non-boundary scvfs on network grids, use precalculated outside indices)
const auto& outsideScvIndices = [&] ()
@@ -665,7 +646,6 @@ private:
// containers storing the global data
std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
std::vector< std::vector< std::vector<GridIndexType> > > neighborVolVarIndices_;
- std::vector<bool> primaryInteractionVolumeVertices_;
std::vector<bool> secondaryInteractionVolumeVertices_;
std::vector<bool> isGhostVertex_;
std::size_t numScvs_;
diff --git a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
index 0cc347bf1b8dbaddd29f667868b7dcdaf6918b4a..5a7a7882f0ddfcfe07e9b4ed0005244d9e0ef336 100644
--- a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
+++ b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
@@ -75,7 +75,7 @@ public:
for (const auto& vertex : vertices(fvGridGeometry.gridView()))
{
const auto vIdxGlobal = fvGridGeometry.vertexMapper().index(vertex);
- if (fvGridGeometry.vertexUsesPrimaryInteractionVolume(vIdxGlobal))
+ if (!fvGridGeometry.vertexUsesSecondaryInteractionVolume(vIdxGlobal))
numPrimaryIV_ += PrimaryIV::numInteractionVolumesAtVertex((*dualGridIndexSet_)[vIdxGlobal]);
else
numSecondaryIV_ += SecondaryIV::numInteractionVolumesAtVertex((*dualGridIndexSet_)[vIdxGlobal]);
@@ -90,7 +90,7 @@ public:
for (const auto& vertex : vertices(fvGridGeometry.gridView()))
{
const auto vIdxGlobal = fvGridGeometry.vertexMapper().index(vertex);
- if (fvGridGeometry.vertexUsesPrimaryInteractionVolume(vIdxGlobal))
+ if (!fvGridGeometry.vertexUsesSecondaryInteractionVolume(vIdxGlobal))
PrimaryIV::addInteractionVolumeIndexSets(primaryIVIndexSets_, scvfIndexMap_, (*dualGridIndexSet_)[vIdxGlobal]);
else
SecondaryIV::addInteractionVolumeIndexSets(secondaryIVIndexSets_, scvfIndexMap_, (*dualGridIndexSet_)[vIdxGlobal]);
diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh
index e4ce6f1a1ca66cc0d0fa9cd95b2902638fc2c50c..c103d99aea74775dc5bdc26f2fdf5033a53f4fb8 100644
--- a/dumux/discretization/cellcentered/mpfa/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/helper.hh
@@ -560,16 +560,16 @@ class CCMpfaHelperImplementation : public MpfaDimensionHelper<TypeTag, dim, dimW
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
+ using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+ using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
+
using VertexMapper = typename GET_PROP_TYPE(TypeTag, VertexMapper);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage;
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
- using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType;
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
-
using CoordScalar = typename GridView::ctype;
+ using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
public:
@@ -627,6 +627,11 @@ public:
return ghostVertices;
}
+ //! Returns whether or not secondary interaction volumes have to be considered in the model.
+ //! This is always the case when the specified types for the interaction volumes differ.
+ static constexpr bool considerSecondaryIVs()
+ { return !std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value; }
+
//! returns whether or not a value exists in a vector
template<typename V1, typename V2>
static bool vectorContainsValue(const std::vector<V1>& vector, const V2 value)
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
index 381d501f304321864e68fbea9a29092590eee16d..51ac0be10eebe29b0c12ce5c8df5480725f1c059 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
@@ -24,6 +24,8 @@
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEBASE_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEBASE_HH
+#include <dune/common/exceptions.hh>
+
#include <dumux/common/properties.hh>
namespace Dumux
@@ -59,8 +61,6 @@ namespace Dumux
* using LocalScvfType = ...;
* //! export the type of used for the iv-local face data
* using LocalFaceData = ...;
- * //! export the type of face data container
- * using LocalFaceDataContainer = ...;
* //! export the type used for iv-local matrices
* using Matrix = ...;
* //! export the type used for iv-local vectors
@@ -84,8 +84,8 @@ template< class Impl, class T>
class CCMpfaInteractionVolumeBase
{
// Curiously recurring template pattern
- Impl & asImp() { return static_cast<Impl&>(*this); }
- const Impl & asImp() const { return static_cast<const Impl&>(*this); }
+ Impl& asImp() { return static_cast<Impl&>(*this); }
+ const Impl& asImp() const { return static_cast<const Impl&>(*this); }
using Problem = typename T::Problem;
using FVElementGeometry = typename T::FVElementGeometry;
@@ -116,11 +116,11 @@ public:
//! returns the number of scvs embedded in this interaction volume
std::size_t numScvs() const { return asImp().numScvs(); }
- //! returns the number of scvfs embedded in this interaction volume
- std::size_t numScvfs() const { return asImp().numScvfs(); }
-
- //! returns a reference to the container with the local face data
- const typename Traits::LocalFaceDataContainer& localFaceData() const { asImp().localFaceData(); }
+ //! 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<typename Traits::LocalFaceData>& localFaceData() const
+ { DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a localFaceData() funtion"); }
//! returns the cell-stencil of this interaction volume
const typename Traits::Stencil& stencil() const { return asImp().stencil(); }
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
index f91f5e9e97b9bff3a61488fed65abd01eeab9be3..7a41428485d51205c0db03158f21ed971346ac07 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
@@ -28,6 +28,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
+#include <dumux/common/matrixvectorhelper.hh>
namespace Dumux
{
@@ -63,17 +64,17 @@ public:
void resize(const InteractionVolume& iv)
{
// resize transmissibility matrix & pressure vectors
- T_.resize(iv.numFaces(), iv.numKnowns());
+ resizeMatrix(T_, iv.numFaces(), iv.numKnowns());
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
- p_[pIdx].resize(iv.numKnowns());
+ resizeVector(p_[pIdx], iv.numKnowns());
// maybe resize gravity container
static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
if (enableGravity)
{
- CA_.resize(iv.numFaces(), iv.numUnknowns());
+ resizeMatrix(CA_, iv.numFaces(), iv.numUnknowns());
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
- g_[pIdx].resize(iv.numFaces());
+ resizeVector(g_[pIdx], iv.numFaces());
}
}
@@ -86,31 +87,31 @@ public:
if (!enableGravity)
{
- T_.resize(iv.numFaces(), iv.numKnowns());
+ resizeMatrix(T_, iv.numFaces(), iv.numKnowns());
outsideT_.resize(iv.numFaces());
- for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
- p_[pIdx].resize(iv.numKnowns());
+ 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)
- outsideT_[i][j].resize(iv.numKnowns());
+ resizeVector(outsideT_[i][j], iv.numKnowns());
}
}
else
{
- T_.resize(iv.numFaces(), iv.numKnowns());
- CA_.resize(iv.numFaces(), iv.numUnknowns());
- A_.resize(iv.numUnknowns(), iv.numUnknowns());
+ 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)
{
- p_[pIdx].resize(iv.numKnowns());
- g_[pIdx].resize(iv.numFaces());
+ resizeVector(p_[pIdx], iv.numKnowns());
+ resizeVector(g_[pIdx], iv.numFaces());
outsideG_[pIdx].resize(iv.numFaces());
}
@@ -120,12 +121,10 @@ public:
outsideT_[i].resize(numNeighbors);
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
- outsideG_[pIdx][i].resize(numNeighbors);
-
+ resizeVector(outsideG_[pIdx][i], numNeighbors);
for (LocalIndexType j = 0; j < numNeighbors; ++j)
- outsideT_[i][j].resize(iv.numKnowns());
+ resizeVector(outsideT_[i][j], iv.numKnowns());
}
-
}
}
@@ -208,8 +207,8 @@ public:
continue;
// resize transmissibility matrix & mole fraction vector
- T_[pIdx][cIdx].resize(iv.numFaces(), iv.numKnowns());
- xj_[pIdx][cIdx].resize(iv.numKnowns());
+ resizeMatrix(T_[pIdx][cIdx], iv.numFaces(), iv.numKnowns());
+ resizeVector(xj_[pIdx][cIdx], iv.numKnowns());
// resize outsideTij on surface grids
if (dim < dimWorld)
@@ -218,7 +217,12 @@ public:
using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType;
for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
- outsideT_[pIdx][cIdx][i].resize(iv.localScvf(i).neighboringLocalScvIndices().size()-1);
+ {
+ 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());
+ }
}
}
}
@@ -261,8 +265,8 @@ public:
void resize(const InteractionVolume& iv)
{
//! resize transmissibility matrix & temperature vector
- T_.resize(iv.numFaces(), iv.numKnowns());
- Tj_.resize(iv.numKnowns());
+ resizeMatrix(T_, iv.numFaces(), iv.numKnowns());
+ resizeVector(Tj_, iv.numKnowns());
//! resize outsideTij on surface grids
if (dim < dimWorld)
@@ -271,7 +275,12 @@ public:
using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType;
for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
- outsideT_[i].resize(iv.localScvf(i).neighboringLocalScvIndices().size()-1);
+ {
+ 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());
+ }
}
}
diff --git a/dumux/discretization/cellcentered/mpfa/localassembler.hh b/dumux/discretization/cellcentered/mpfa/localassembler.hh
index 2c4aec915797e28ea9dc8e1381deebf22c33f506..d5eb6bf8289f43f07a475a4cf63aaf24546e5b77 100644
--- a/dumux/discretization/cellcentered/mpfa/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/localassembler.hh
@@ -28,14 +28,16 @@
#include <dune/common/exceptions.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
namespace Dumux
{
//! Forward declaration of the implementation
-template< class InteractionVolume > class InteractionVolumeAssemblerImpl;
+template< class IVTraits, MpfaMethods M > class InteractionVolumeAssemblerImpl;
//! Alias to select the right implementation.
-template< class InteractionVolume >
-using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< InteractionVolume >;
+template< class IVTraits, MpfaMethods M >
+using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< IVTraits, M >;
/*!
* \ingroup CCMpfaDiscretization
@@ -45,19 +47,17 @@ using InteractionVolumeAssembler = InteractionVolumeAssemblerImpl< InteractionVo
* for the available interaction volume implementations. these
* should derive from this base clases.
*
- * \tparam Interaction Volume The interaction volume implementation
- * used by the mpfa scheme
+ * \tparam IVTraits The Traits class used by the interaction volume
*/
-template< class InteractionVolume >
+template< class IVTraits >
class InteractionVolumeAssemblerBase
{
- using Traits = typename InteractionVolume::Traits;
- using Matrix = typename Traits::Matrix;
- using Vector = typename Traits::Vector;
+ using Matrix = typename IVTraits::Matrix;
+ using Vector = typename IVTraits::Vector;
- using Problem = typename Traits::Problem;
- using FVElementGeometry = typename Traits::FVElementGeometry;
- using ElementVolumeVariables = typename Traits::ElementVolumeVariables;
+ using Problem = typename IVTraits::Problem;
+ using FVElementGeometry = typename IVTraits::FVElementGeometry;
+ using ElementVolumeVariables = typename IVTraits::ElementVolumeVariables;
public:
/*!
@@ -86,15 +86,16 @@ class InteractionVolumeAssemblerBase
* \brief General interface of a function assembling the
* interaction volume-local transmissibility matrix.
*
- * \tparam GetTensorFunction Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam IV Interaction volume type implementation
+ * \tparam GetTensor Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
*
* \param T The transmissibility matrix to be assembled
* \param iv The interaction volume
* \param getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class GetTensorFunction >
- void assemble(Matrix& T, InteractionVolume& iv, const GetTensorFunction& getTensor)
+ template< class IV, class GetTensor >
+ void assemble(Matrix& T, IV& iv, const GetTensor& getTensor)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function");
}
@@ -102,16 +103,21 @@ class InteractionVolumeAssemblerBase
/*!
* \brief General interface of a function assembling the interaction
* volume-local transmissibilities matrix for surface grids. The
- * transmissibility associated with "outside" faces are stored
+ * 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 GetTensor 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 getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class OutsideTijContainer, class GetTensorFunction >
- void assemble(OutsideTijContainer& outsideTij, Matrix& T, InteractionVolume& iv, const GetTensorFunction& getTensor)
+ template< class TOutside, class IV, class GetTensor >
+ void assemble(TOutside& outsideTij, Matrix& T, IV& iv, const GetTensor& getTensor)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function to be used on surface grids");
}
@@ -121,10 +127,11 @@ class InteractionVolumeAssemblerBase
* volume-local transmissibility matrix in the case that gravity
* is to be considered in the local system of equations.
*
- * \tparam GravityContainer The type of container used to store the
- * gravitational acceleration per scvf & phase
- * \tparam GetTensorFunction Lambda to obtain the tensor w.r.t.
- * which the local system is to be solved
+ * \tparam GC The type of container used to store the
+ * gravitational acceleration per scvf & phase
+ * \tparam IV The interaction volume type implementation
+ * \tparam GetTensor 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
@@ -132,12 +139,8 @@ class InteractionVolumeAssemblerBase
* \param iv The interaction volume
* \param getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class GravityContainer, class GetTensorFunction >
- void assembleWithGravity(Matrix& T,
- GravityContainer& g,
- Matrix& CA,
- InteractionVolume& iv,
- const GetTensorFunction& getTensor)
+ template< class GC, class IV, class GetTensor >
+ void assembleWithGravity(Matrix& T, GC& g, Matrix& CA, IV& iv, const GetTensor& getTensor)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleWithGravity() function");
}
@@ -150,6 +153,14 @@ class InteractionVolumeAssemblerBase
* 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 GetTensor 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
@@ -159,18 +170,15 @@ class InteractionVolumeAssemblerBase
* \param iv The interaction volume
* \param getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class GravityContainer,
- class OutsideGravityContainer,
- class OutsideTijContainer,
- class GetTensorFunction >
- void assembleWithGravity(OutsideTijContainer& outsideTij,
+ template< class GC, class GOut, class TOutside, class IV, class GetTensor >
+ void assembleWithGravity(TOutside& outsideTij,
Matrix& T,
- OutsideGravityContainer& outsideG,
- GravityContainer& g,
+ GOut& outsideG,
+ GC& g,
Matrix& CA,
Matrix& A,
- InteractionVolume& iv,
- const GetTensorFunction& getTensor)
+ IV& iv,
+ const GetTensor& getTensor)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleWithGravity() function to be used on surface grids");
}
@@ -180,14 +188,15 @@ class InteractionVolumeAssemblerBase
* 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 GetU >
- void assemble(Vector& u, const InteractionVolume& iv, const GetU& getU)
+ template< class IV, class GetU >
+ void assemble(Vector& u, const IV& iv, const GetU& getU)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assemble() function for the cell unknowns");
}
@@ -200,17 +209,19 @@ class InteractionVolumeAssemblerBase
* 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 GetTensor 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 getTensor Lambda to evaluate scv-wise hydraulic conductivities
*/
- template< class GravityContainer, class GetTensorFunction >
- void assembleGravity(GravityContainer& g,
- const InteractionVolume& iv,
- const Matrix& CA,
- const GetTensorFunction& getTensor)
+ template< class GC, class IV, class GetTensor >
+ void assembleGravity(GC& g, const IV& iv, const Matrix& CA, const GetTensor& getTensor)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleGravity() function");
}
@@ -225,6 +236,13 @@ class InteractionVolumeAssemblerBase
* 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 GetTensor 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
@@ -233,15 +251,13 @@ class InteractionVolumeAssemblerBase
* \param A Matrix needed for the "reconstruction" of face unknowns as a function of gravity
* \param getTensor Lambda to evaluate scv-wise hydraulic conductivities
*/
- template< class GravityContainer,
- class OutsideGravityContainer,
- class GetTensorFunction >
- void assembleGravity(GravityContainer& g,
- OutsideGravityContainer& outsideG,
- const InteractionVolume& iv,
+ template< class GC, class GOut, class IV, class GetTensor >
+ void assembleGravity(GC& g,
+ GOut& outsideG,
+ const IV& iv,
const Matrix& CA,
const Matrix& A,
- const GetTensorFunction& getTensor)
+ const GetTensor& getTensor)
{
DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleGravity() function to be used on surface grids");
}
diff --git a/dumux/discretization/cellcentered/mpfa/localfacedata.hh b/dumux/discretization/cellcentered/mpfa/localfacedata.hh
index 886677085a732e72ee3531954facd04e3ecd5eac..1fe6489b2deb3729fc1da0f0f70f76f8f8811a9c 100644
--- a/dumux/discretization/cellcentered/mpfa/localfacedata.hh
+++ b/dumux/discretization/cellcentered/mpfa/localfacedata.hh
@@ -47,6 +47,9 @@ class InteractionVolumeLocalFaceData
bool isOutside_; //!< indicates if this face maps to the iv-local index from "outside"
public:
+ //! Default constructor
+ InteractionVolumeLocalFaceData() = default;
+
//! Constructor
InteractionVolumeLocalFaceData(LocalIndexType faceIndex,
LocalIndexType scvIndex,
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/CMakeLists.txt b/dumux/discretization/cellcentered/mpfa/omethod/CMakeLists.txt
index d2c61c2a7b45d6ba942dfb77a7933b45875f6e52..9b738e874ef3af320d2723311d20062f1b75d976 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/CMakeLists.txt
+++ b/dumux/discretization/cellcentered/mpfa/omethod/CMakeLists.txt
@@ -3,4 +3,5 @@ interactionvolume.hh
interactionvolumeindexset.hh
localassembler.hh
localsubcontrolentities.hh
+staticinteractionvolume.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/cellcentered/mpfa/omethod)
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 0e8cd035849e40762e10acbbd009997d8317a4f5..298ec506bd6461955c8b68fe6ee4c25d3417d8ab 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -30,6 +30,7 @@
#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
+#include <dumux/common/matrixvectorhelper.hh>
#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh>
#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
@@ -77,8 +78,6 @@ public:
using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf< IndexSet >;
//! export the type of used for the iv-local face data
using LocalFaceData = InteractionVolumeLocalFaceData<GridIndexType, LocalIndexType>;
- //! export the type of face data container (use dynamic container here)
- using LocalFaceDataContainer = std::vector< LocalFaceData >;
//! export the type used for iv-local matrices
using Matrix = Dune::DynamicMatrix< ScalarType >;
//! export the type used for iv-local vectors
@@ -92,6 +91,8 @@ public:
/*!
* \ingroup CCMpfaDiscretization
* \brief Class for the interaction volume of the mpfa-o method.
+ * This implementation creates dynamic objects of the local geometries
+ * and can be used at boundaries and on unstructured grids.
*/
template< class Traits >
class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >
@@ -120,6 +121,14 @@ class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInte
using LocalIndexType = typename Traits::LocalIndexType;
using Stencil = typename Traits::Stencil;
+ //! For the o method, the interaction volume stencil can be taken directly
+ //! from the nodal index set, which always uses dynamic types to be compatible
+ //! on the boundaries and unstructured grids. Thus, we have to make sure that
+ //! the type set for the stencils in the traits is castable.
+ static_assert( std::is_convertible<Stencil*, typename IndexSet::GridIndexContainer*>::value,
+ "The o-method uses the (dynamic) nodal index set's stencil as the interaction volume stencil. "
+ "Using a different type is not permissive here." );
+
//! Data attached to scvf touching Dirichlet boundaries.
//! For the default o-scheme, we only store the corresponding vol vars index.
class DirichletData
@@ -240,10 +249,10 @@ public:
}
}
- // Resize local matrices
- A_.resize(numUnknowns_, numUnknowns_);
- B_.resize(numUnknowns_, numKnowns_);
- C_.resize(numFaces_, numUnknowns_);
+ // 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
@@ -258,9 +267,6 @@ public:
//! returns the number of scvs embedded in this interaction volume
std::size_t numScvs() const { return scvs_.size(); }
- //! returns the number of scvfs embedded in this interaction volume
- std::size_t numScvfs() const { return scvfs_.size(); }
-
//! returns the cell-stencil of this interaction volume
const Stencil& stencil() const { return *stencil_; }
@@ -271,7 +277,7 @@ public:
const LocalScvfType& localScvf(const LocalIndexType ivLocalScvfIdx) const { return scvfs_[ivLocalScvfIdx]; }
//! returns the local scv entity corresponding to a given iv-local scv idx
- const LocalScvType& localScv(const LocalIndexType ivLocalScvfIdx) const { return scvs_[ivLocalScvfIdx]; }
+ const LocalScvType& localScv(const LocalIndexType ivLocalScvIdx) const { return scvs_[ivLocalScvIdx]; }
//! returns a reference to the container with the local face data
const std::vector<LocalFaceData>& localFaceData() const { return localFaceData_; }
@@ -334,7 +340,7 @@ private:
Matrix B_;
Matrix C_;
- // The omega factors are stored during assemble of local system
+ // The omega factors are stored during assembly of local system
std::vector< std::vector< DimVector > > wijk_;
// sizes involved in the local system equations
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
index aabba86066a1ab1dd05ddb45f1a62f3006210c49..4b5f11aa27cf41a64c3a2cb8b0124f6db21bfb4d 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
@@ -20,17 +20,18 @@
* \file
* \ingroup CCMpfaDiscretization
* \brief Class for the assembly of the local systems of equations
- * involved in the transmissibility computaion in the mpfa-o scheme.
+ * involved in the transmissibility computaion in an mpfa-o type manner.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH
#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
+#include <dumux/common/matrixvectorhelper.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
#include <dumux/discretization/cellcentered/mpfa/localassembler.hh>
#include <dumux/discretization/cellcentered/mpfa/computetransmissibility.hh>
-#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
namespace Dumux
{
@@ -38,16 +39,15 @@ namespace Dumux
/*!
* \ingroup CCMpfaDiscretization
* \brief Specialization of the interaction volume-local
- * assembler class for the mpfa-o scheme.
+ * assembler class for the schemes using an mpfa-o type assembly.
*
- * \tparam IVTraits The traits class of the interaction volume
+ * \tparam IVTraits The traits class used by the interaction volume implemetation
*/
template< class IVTraits >
-class InteractionVolumeAssemblerImpl< CCMpfaOInteractionVolume<IVTraits> >
- : public InteractionVolumeAssemblerBase< CCMpfaOInteractionVolume<IVTraits> >
+class InteractionVolumeAssemblerImpl< IVTraits, MpfaMethods::oMethod >
+ : public InteractionVolumeAssemblerBase< IVTraits >
{
- using InteractionVolume = CCMpfaOInteractionVolume< IVTraits >;
- using ParentType = InteractionVolumeAssemblerBase< InteractionVolume >;
+ using ParentType = InteractionVolumeAssemblerBase< IVTraits >;
using LocalIndexType = typename IVTraits::LocalIndexType;
using Matrix = typename IVTraits::Matrix;
@@ -61,15 +61,19 @@ public:
using ParentType::ParentType;
/*!
- * \brief Assembles the transmissibility matrix
- * within an interaction volume for the mpfa-o scheme.
+ * \brief Assembles the transmissibility matrix within an
+ * interaction volume in an mpfa-o type way.
+ *
+ * \tparam IV The interaction volume type implementation
+ * \tparam GetTensor 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 mpfa-o interaction volume
+ * \param iv The interaction volume
* \param getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class GetTensorFunction >
- void assemble(Matrix& T, InteractionVolume& iv, const GetTensorFunction& getTensor)
+ template< class IV, class GetTensor >
+ void assemble(Matrix& T, IV& iv, const GetTensor& getTensor)
{
// assemble D into T directly
assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor);
@@ -89,13 +93,18 @@ public:
* 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 GetTensor 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 iv The interaction volume
* \param getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class OutsideTijContainer, class GetTensorFunction >
- void assemble(OutsideTijContainer& outsideTij, Matrix& T, InteractionVolume& iv, const GetTensorFunction& getTensor)
+ template< class TOutside, class IV, class GetTensor >
+ void assemble(TOutside& outsideTij, Matrix& T, IV& iv, const GetTensor& getTensor)
{
// assemble D into T directly
assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor);
@@ -154,18 +163,20 @@ public:
* 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 GetTensor 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 getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class GravityContainer, class GetTensorFunction >
- void assembleWithGravity(Matrix& T,
- GravityContainer& g,
- Matrix& CA,
- InteractionVolume& iv,
- const GetTensorFunction& getTensor)
+ template< class GC, class IV, class GetTensor >
+ void assembleWithGravity(Matrix& T, GC& g, Matrix& CA, IV& iv, const GetTensor& getTensor)
{
// assemble D into T & C into CA directly
assembleLocalMatrices_(iv.A(), iv.B(), CA, T, iv, getTensor);
@@ -190,6 +201,14 @@ public:
* 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 GetTensor 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
@@ -199,18 +218,15 @@ public:
* \param iv The mpfa-o interaction volume
* \param getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class GravityContainer,
- class OutsideGravityContainer,
- class OutsideTijContainer,
- class GetTensorFunction >
- void assembleWithGravity(OutsideTijContainer& outsideTij,
+ template< class GC, class GOut, class TOutside, class IV, class GetTensor >
+ void assembleWithGravity(TOutside& outsideTij,
Matrix& T,
- OutsideGravityContainer& outsideG,
- GravityContainer& g,
+ GOut& outsideG,
+ GC& g,
Matrix& CA,
Matrix& A,
- InteractionVolume& iv,
- const GetTensorFunction& getTensor)
+ IV& iv,
+ const GetTensor& getTensor)
{
// assemble D into T directly
assembleLocalMatrices_(iv.A(), iv.B(), iv.C(), T, iv, getTensor);
@@ -271,18 +287,17 @@ public:
* \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 u The vector to be filled with the cell unknowns
* \param iv The mpfa-o interaction volume
* \param getU Lambda to obtain the desired cell/Dirichlet value from grid index
*/
- template< class GetU >
- void assemble(Vector& u, const InteractionVolume& iv, const GetU& getU)
+ template< class IV, class GetU >
+ void assemble(Vector& u, const IV& iv, const GetU& getU)
{
- // resize given container
- u.resize(iv.numKnowns());
-
// put the cell pressures first
-
for (LocalIndexType i = 0; i < iv.numScvs(); ++i)
u[i] = getU( iv.localScv(i).globalScvIndex() );
@@ -300,17 +315,19 @@ public:
* 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 GetTensor 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 getTensor Lambda to evaluate scv-wise hydraulic conductivities
*/
- template< class GravityContainer, class GetTensorFunction >
- void assembleGravity(GravityContainer& g,
- const InteractionVolume& iv,
- const Matrix& CA,
- const GetTensorFunction& getTensor)
+ template< class GC, class IV, class GetTensor >
+ void assembleGravity(GC& g, const IV& iv, const Matrix& CA, const GetTensor& getTensor)
{
//! We require the gravity container to be a two-dimensional vector/array type, structured as follows:
//! - first index adresses the respective phases
@@ -329,8 +346,14 @@ public:
// reset gravity containers to zero
const auto numPhases = g.size();
- std::vector< Vector > sum_alphas(numPhases, Vector(iv.numUnknowns(), 0.0));
- std::for_each(g.begin(), g.end(), [&iv](auto& v) { v = 0.0; });
+ std::vector< Vector > 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)
{
@@ -412,6 +435,13 @@ public:
* 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 GetTensor 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
@@ -420,15 +450,13 @@ public:
* \param A Matrix needed for the "reconstruction" of face unknowns as a function of gravity
* \param getTensor Lambda to evaluate scv-wise hydraulic conductivities
*/
- template< class GravityContainer,
- class OutsideGravityContainer,
- class GetTensorFunction >
- void assembleGravity(GravityContainer& g,
- OutsideGravityContainer& outsideG,
- const InteractionVolume& iv,
+ template< class GC, class GOut, class IV, class GetTensor >
+ void assembleGravity(GC& g,
+ GOut& outsideG,
+ const IV& iv,
const Matrix& CA,
const Matrix& A,
- const GetTensorFunction& getTensor)
+ const GetTensor& getTensor)
{
//! We require the gravity container to be a two-dimensional vector/array type, structured as follows:
//! - first index adresses the respective phases
@@ -453,10 +481,12 @@ public:
// reset everything to zero
const auto numPhases = g.size();
- std::vector< Vector > sum_alphas(numPhases, Vector(iv.numUnknowns(), 0.0));
+ std::vector< Vector > 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; });
}
@@ -585,6 +615,10 @@ private:
*
* \note The matrices are expected to have been resized beforehand.
*
+ * \tparam IV The interaction volume type implementation
+ * \tparam GetTensor Lambda to obtain the tensor w.r.t.
+ * which the local system is to be solved
+ *
* \param A The A matrix of the iv-local equation system
* \param B The B matrix of the iv-local equation system
* \param C The C matrix of the iv-local flux expressions
@@ -592,10 +626,9 @@ private:
* \param iv The mpfa-o interaction volume
* \param getTensor Lambda to evaluate the scv-wise tensors
*/
- template< class GetTensorFunction >
+ template< class IV, class GetTensor >
void assembleLocalMatrices_(Matrix& A, Matrix& B, Matrix& C, Matrix& D,
- InteractionVolume& iv,
- const GetTensorFunction& getTensor)
+ IV& iv, const GetTensor& getTensor)
{
// Matrix D is assumed to have the right size already
assert(D.rows() == iv.numFaces() && D.cols() == iv.numKnowns() && "Matrix D does not have the correct size");
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
index 51ab083d0fa7e72d407da8d8925192a8edc490f5..7ddded18c3abba12cc578707e0358606da02597d 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
@@ -54,6 +54,9 @@ public:
static constexpr int myDimension = dim;
static constexpr int worldDimension = dimWorld;
+ //! The default constructor
+ CCMpfaOInteractionVolumeLocalScv() = default;
+
/*!
* \brief The constructor
*
@@ -69,7 +72,7 @@ public:
const SubControlVolume& scv,
const LocalIndexType localIndex,
const IvIndexSet& indexSet)
- : indexSet_(indexSet)
+ : indexSet_(&indexSet)
, globalScvIndex_(scv.dofIndex())
, localDofIndex_(localIndex)
{
@@ -103,7 +106,7 @@ public:
LocalIndexType scvfIdxLocal(unsigned int coordDir) const
{
assert(coordDir < myDimension);
- return indexSet_.scvfIdxLocal(localDofIndex_, coordDir);
+ return indexSet_->scvfIdxLocal(localDofIndex_, coordDir);
}
//! the nu vectors are needed for setting up the omegas of the iv
@@ -114,7 +117,7 @@ public:
}
private:
- const IvIndexSet& indexSet_;
+ const IvIndexSet* indexSet_;
GridIndexType globalScvIndex_;
LocalIndexType localDofIndex_;
LocalBasis nus_;
@@ -138,6 +141,9 @@ public:
using GridIndexType = typename IvIndexSet::GridIndexType;
using LocalIndexType = typename IvIndexSet::LocalIndexType;
+ //! The default constructor
+ CCMpfaOInteractionVolumeLocalScvf() = default;
+
/*!
* \brief The constructor
*
@@ -154,7 +160,7 @@ public:
: isDirichlet_(isDirichlet)
, scvfIdxGlobal_(scvf.index())
, localDofIndex_(localDofIdx)
- , neighborScvIndicesLocal_(localScvIndices)
+ , neighborScvIndicesLocal_(&localScvIndices)
{}
//! This is either the iv-local index of the intermediate unknown (interior/Neumann face)
@@ -165,7 +171,7 @@ public:
GridIndexType globalScvfIndex() const { return scvfIdxGlobal_; }
//! Returns the local indices of the scvs neighboring this scvf
- const LocalIndexContainer& neighboringLocalScvIndices() const { return neighborScvIndicesLocal_; }
+ const LocalIndexContainer& neighboringLocalScvIndices() const { return *neighborScvIndicesLocal_; }
//! states if this is scvf is on a Dirichlet boundary
bool isDirichlet() const { return isDirichlet_; }
@@ -174,7 +180,7 @@ private:
bool isDirichlet_;
GridIndexType scvfIdxGlobal_;
LocalIndexType localDofIndex_;
- const LocalIndexContainer& neighborScvIndicesLocal_;
+ const LocalIndexContainer* neighborScvIndicesLocal_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d4401e178678120e2d3c807aa683281c052f6a35
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
@@ -0,0 +1,317 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup CCMpfaDiscretization
+ * \brief Class for the interaction volume of the mpfa-o scheme to be used when the
+ * sizes are known at compile time, e.g. for non-boundary interaction volumes
+ * on structured grids.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_STATIC_INTERACTIONVOLUME_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_O_STATIC_INTERACTIONVOLUME_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/matrixvectorhelper.hh>
+
+#include <dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh>
+#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
+#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
+#include <dumux/discretization/cellcentered/mpfa/localfacedata.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+#include "localsubcontrolentities.hh"
+#include "interactionvolumeindexset.hh"
+
+namespace Dumux
+{
+//! Forward declaration of the o-method's static interaction volume
+template< class Traits, int localSize > class CCMpfaOStaticInteractionVolume;
+
+//! Specialization of the default interaction volume traits class for the mpfa-o method
+template< class TypeTag, int localSize >
+struct CCMpfaODefaultStaticInteractionVolumeTraits
+{
+private:
+ using GridIndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
+ static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
+ static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
+
+public:
+ //! export the problem type (needed for iv-local assembly)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ //! export the type of the local view on the finite volume grid geometry
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ //! export the type of the local view on the grid volume variables
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ //! export the type of grid view
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ //! export the type used for scalar values
+ using ScalarType = typename GET_PROP_TYPE(TypeTag, Scalar);
+ //! export the type of the mpfa helper class
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ //! export the type used for local indices
+ using LocalIndexType = std::uint8_t;
+ //! export the type for the interaction volume index set
+ using IndexSet = CCMpfaOInteractionVolumeIndexSet< DualGridNodalIndexSet<GridIndexType, LocalIndexType, dim> >;
+ //! export the type of interaction-volume local scvs
+ using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, ScalarType, dim, dimWorld >;
+ //! export the type of interaction-volume local scvfs
+ using LocalScvfType = CCMpfaOInteractionVolumeLocalScvf< IndexSet >;
+ //! export the type of used for the iv-local face data
+ using LocalFaceData = InteractionVolumeLocalFaceData<GridIndexType, LocalIndexType>;
+ //! export the type used for iv-local matrices
+ using Matrix = Dune::FieldMatrix< ScalarType, localSize, localSize >;
+ //! export the type used for iv-local vectors
+ using Vector = Dune::FieldVector< ScalarType, localSize >;
+ //! export the type used for the iv-stencils
+ using Stencil = std::vector< GridIndexType >;
+ //! export the data handle type for this iv
+ using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >;
+};
+
+/*!
+ * \ingroup CCMpfaDiscretization
+ * \brief Class for the interaction volume of the mpfa-o method.
+ * This implementation creates static objects of the local geometries
+ * and can only be used for interior interaction volumes with a static
+ * size known at compile time. This size has to match the sizes of the
+ * vectors/matrices defined in the traits class in case static types are used.
+ *
+ * \tparam Traits The type traits class to be used
+ * \tparam localSize The size of the local eq system
+ * This also is the number of local scvs/scvfs
+ */
+template< class Traits, int localSize >
+class CCMpfaOStaticInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >
+{
+ using ThisType = CCMpfaOInteractionVolume< Traits >;
+ using ParentType = CCMpfaInteractionVolumeBase< CCMpfaOInteractionVolume<Traits>, Traits >;
+
+ using Scalar = typename Traits::ScalarType;
+ using Helper = typename Traits::MpfaHelper;
+ using Problem = typename Traits::Problem;
+ using FVElementGeometry = typename Traits::FVElementGeometry;
+
+ using GridView = typename Traits::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ static constexpr int dim = GridView::dimension;
+ using DimVector = Dune::FieldVector<Scalar, dim>;
+
+ using Matrix = typename Traits::Matrix;
+ using LocalScvType = typename Traits::LocalScvType;
+ using LocalScvfType = typename Traits::LocalScvfType;
+ using LocalFaceData = typename Traits::LocalFaceData;
+
+ using IndexSet = typename Traits::IndexSet;
+ using GridIndexType = typename GridView::IndexSet::IndexType;
+ using LocalIndexType = typename Traits::LocalIndexType;
+ using Stencil = typename Traits::Stencil;
+
+ //! For the o method, the interaction volume stencil can be taken directly
+ //! from the nodal index set, which always uses dynamic types to be compatible
+ //! on the boundaries and unstructured grids. Thus, we have to make sure that
+ //! the type set for the stencils in the traits is castable.
+ static_assert( std::is_convertible<Stencil*, typename IndexSet::GridIndexContainer*>::value,
+ "The o-method uses the (dynamic) nodal index set's stencil as the interaction volume stencil. "
+ "Using a different type is not permissive here." );
+
+ //! Data attached to scvf touching Dirichlet boundaries.
+ //! For the default o-scheme, we only store the corresponding vol vars index.
+ class DirichletData
+ {
+ GridIndexType volVarIndex_;
+
+ public:
+ //! Constructor
+ DirichletData(const GridIndexType index) : volVarIndex_(index) {}
+
+ //! Return corresponding vol var index
+ GridIndexType volVarIndex() const { return volVarIndex_; }
+ };
+
+public:
+ //! publicly state the mpfa-scheme this interaction volume is associated with
+ static constexpr MpfaMethods MpfaMethod = MpfaMethods::oMethod;
+
+ //! Sets up the local scope for a given iv index set
+ void setUpLocalScope(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();
+
+ // set up local geometry containers
+ const auto& scvIndices = indexSet.globalScvIndices();
+ for (LocalIndexType localIdx = 0; localIdx < localSize; localIdx++)
+ {
+ // scv-related quantities
+ scvs_[localIdx] = LocalScvType(Helper(),
+ fvGeometry,
+ fvGeometry.scv( scvIndices[localIdx] ),
+ localIdx,
+ indexSet);
+ elements_[localIdx] = fvGeometry.fvGridGeometry().element( scvIndices[localIdx] );
+
+ // scvf-related quantities
+ const auto& scvf = fvGeometry.scvf(indexSet.scvfIdxGlobal(localIdx));
+
+ // this interaction volume implementation does not work on boundaries
+ assert(!scvf.boundary() && "static mpfa-o interaction volume cannot be used on boundaries");
+
+ const auto& neighborScvIndicesLocal = indexSet.neighboringLocalScvIndices(localIdx);
+ scvfs_[localIdx] = LocalScvfType(scvf, neighborScvIndicesLocal, localIdx, /*isDirichlet*/false);
+ localFaceData_[localIdx*2] = LocalFaceData(localIdx, neighborScvIndicesLocal[0], scvf.index());
+
+ const auto outsideLocalScvIdx = neighborScvIndicesLocal[1];
+ // loop over scvfs in outside scv until we find the one coinciding with current scvf
+ for (int coord = 0; coord < dim; ++coord)
+ {
+ if (indexSet.scvfIdxLocal(outsideLocalScvIdx, coord) == localIdx)
+ {
+ const auto globalScvfIdx = indexSet.nodalIndexSet().scvfIdxGlobal(outsideLocalScvIdx, coord);
+ const auto& flipScvf = fvGeometry.scvf(globalScvfIdx);
+ localFaceData_[localIdx*2 + 1] = LocalFaceData(localIdx, // iv-local scvf idx
+ outsideLocalScvIdx, // iv-local scv index
+ 0, // scvf-local index in outside faces
+ flipScvf.index()); // global scvf index
+ }
+ }
+ }
+
+ // Maybe resize local matrices if dynamic types are used
+ resizeMatrix(A_, localSize, localSize);
+ resizeMatrix(B_, localSize, localSize);
+ resizeMatrix(C_, localSize, localSize);
+ }
+
+ //! returns the number of primary scvfs of this interaction volume
+ static constexpr std::size_t numFaces() { return localSize; }
+
+ //! returns the number of intermediate unknowns within this interaction volume
+ static constexpr std::size_t numUnknowns() { return localSize; }
+
+ //! returns the number of (in this context) known solution values within this interaction volume
+ static constexpr std::size_t numKnowns() { return localSize; }
+
+ //! returns the number of scvs embedded in this interaction volume
+ static constexpr std::size_t numScvs() { return localSize; }
+
+ //! returns the cell-stencil of this interaction volume
+ const Stencil& stencil() const { return *stencil_; }
+
+ //! returns the grid element corresponding to a given iv-local scv idx
+ const Element& element(const LocalIndexType ivLocalScvIdx) const
+ {
+ assert(ivLocalScvIdx < localSize);
+ return elements_[ivLocalScvIdx];
+ }
+
+ //! returns the local scvf entity corresponding to a given iv-local scvf idx
+ const LocalScvfType& localScvf(const LocalIndexType ivLocalScvfIdx) const
+ {
+ assert(ivLocalScvfIdx < localSize);
+ return scvfs_[ivLocalScvfIdx];
+ }
+
+ //! returns the local scv entity corresponding to a given iv-local scv idx
+ const LocalScvType& localScv(const LocalIndexType ivLocalScvIdx) const
+ {
+ assert(ivLocalScvIdx < localSize);
+ return scvs_[ivLocalScvIdx];
+ }
+
+ //! returns a reference to the container with the local face data
+ const std::array<LocalFaceData, 2*localSize>& localFaceData() const { return localFaceData_; }
+
+ //! Returns a reference to the information container on Dirichlet BCs within this iv.
+ //! Here, we return an empty container as this implementation cannot be used on boundaries.
+ const std::array<DirichletData, 0>& dirichletData() const { return dirichletData_; }
+
+ //! returns the matrix associated with face unknowns in local equation system
+ const Matrix& A() const { return A_; }
+ Matrix& A() { return A_; }
+
+ //! returns the matrix associated with cell unknowns in local equation system
+ const Matrix& B() const { return B_; }
+ Matrix& B() { return B_; }
+
+ //! returns the matrix associated with face unknowns in flux expressions
+ const Matrix& C() const { return C_; }
+ Matrix& C() { return C_; }
+
+ //! returns container storing the transmissibilities for each face & coordinate
+ const std::array< std::array<DimVector, 2>, localSize >& omegas() const { return wijk_; }
+ std::array< std::array<DimVector, 2>, localSize >& omegas() { return wijk_; }
+
+ //! returns the number of interaction volumes living around a vertex
+ //! the mpfa-o scheme always constructs one iv per vertex
+ template< class NodalIndexSet >
+ static std::size_t numInteractionVolumesAtVertex(const NodalIndexSet& nodalIndexSet) { return 1; }
+
+ //! adds the iv index sets living around a vertex to a given container
+ //! and stores the the corresponding index in a map for each scvf
+ template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
+ static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
+ ScvfIndexMap& scvfIndexMap,
+ const NodalIndexSet& nodalIndexSet)
+ {
+ // 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);
+
+ // store the index mapping
+ for (const auto scvfIdx : nodalIndexSet.globalScvfIndices())
+ scvfIndexMap[scvfIdx] = curGlobalIndex;
+ }
+
+private:
+ // pointer to cell stencil (in iv index set)
+ const Stencil* stencil_;
+
+ // Variables defining the local scope
+ std::array<Element, localSize> elements_;
+ std::array<LocalScvType, localSize> scvs_;
+ std::array<LocalScvfType, localSize> scvfs_;
+ std::array<LocalFaceData, 2*localSize> localFaceData_;
+
+ // Empty Dirichlet container to be compatible with dynamic assembly
+ std::array<DirichletData, 0> dirichletData_;
+
+ // Matrices needed for computation of transmissibilities
+ Matrix A_;
+ Matrix B_;
+ Matrix C_;
+
+ // The omega factors are stored during assembly of local system
+ // we assume one "outside" face per scvf (does not work on surface grids)
+ std::array< std::array<DimVector, 2>, localSize > wijk_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh
index a2f6f6cba86b0d0300fc6a2dee44f4ed2409fa83..1a5e96d98a389f00c0976233008b749cc5f4e1d0 100644
--- a/dumux/discretization/cellcentered/mpfa/properties.hh
+++ b/dumux/discretization/cellcentered/mpfa/properties.hh
@@ -50,6 +50,7 @@
#include <dumux/discretization/cellcentered/mpfa/helper.hh>
#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
+#include <dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh>
namespace Dumux
{
@@ -73,7 +74,7 @@ SET_PROP(CCMpfaModel, MpfaMethod)
//! The mpfa helper class
SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>);
-//! Per default, we use the mpfa-o interaction volume
+//! Per default, we use the dynamic mpfa-o interaction volume
SET_PROP(CCMpfaModel, PrimaryInteractionVolume)
{
private:
@@ -83,10 +84,15 @@ public:
using type = CCMpfaOInteractionVolume< Traits >;
};
-//! Per default, we use the mpfa-o interaction volume everywhere (pure mpfa-o scheme)
-SET_TYPE_PROP(CCMpfaModel,
- SecondaryInteractionVolume,
- typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume));
+//! Per default, we use the dynamic mpfa-o interaction volume as secondary type
+SET_PROP(CCMpfaModel, SecondaryInteractionVolume)
+{
+private:
+ //! use the default traits
+ using Traits = CCMpfaODefaultInteractionVolumeTraits< TypeTag >;
+public:
+ using type = CCMpfaOInteractionVolume< Traits >;
+};
//! Set the default for the global finite volume geometry
SET_TYPE_PROP(CCMpfaModel,
diff --git a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
index 03c99d3a6bacdfc531e8fcf93ad6d78664be603e..08376c67351fd442482c051b2485d83af5ed76d1 100644
--- a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
+++ b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
@@ -123,21 +123,40 @@ public:
static_assert(FluidSystem::viscosityIsConstant(0),
"1p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!");
- const auto& fluxVarsCache = elemFluxVarsCache[scvf];
- const auto& volVarIndices = fluxVarsCache.advectionStencil();
- const auto& tij = fluxVarsCache.advectionTij();
-
// we know the "upwind factor" is constant, get inner one here and compute derivatives
static const Scalar up = curElemVolVars[scvf.insideScvIdx()].density()
/ curElemVolVars[scvf.insideScvIdx()].viscosity();
- // add partial derivatives to the respective given matrices
- for (unsigned int i = 0; i < volVarIndices.size();++i)
+ const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ if (fluxVarsCache.usesSecondaryIv())
{
- if (fluxVarsCache.advectionSwitchFluxSign())
- derivativeMatrices[volVarIndices[i]][conti0EqIdx][pressureIdx] -= tij[i]*up;
- else
- derivativeMatrices[volVarIndices[i]][conti0EqIdx][pressureIdx] += tij[i]*up;
+ const auto& stencil = fluxVarsCache.advectionStencilSecondaryIv();
+ const auto& tij = fluxVarsCache.advectionTijSecondaryIv();
+ 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;
+ }
+ }
+ else
+ {
+ const auto& stencil = fluxVarsCache.advectionStencilPrimaryIv();
+ const auto& tij = fluxVarsCache.advectionTijPrimaryIv();
+ 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;
+ }
}
}
@@ -210,22 +229,7 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
- const auto& fluxVarsCache = elemFluxVarsCache[scvf];
- const auto& volVarIndices = fluxVarsCache.advectionStencil();
- const auto& tij = fluxVarsCache.advectionTij();
-
- // we know the "upwind factor" is constant, get inner one here and compute derivatives
- static const Scalar up = curElemVolVars[scvf.insideScvIdx()].density()
- / curElemVolVars[scvf.insideScvIdx()].viscosity();
-
- // add partial derivatives to the respective given matrices
- for (unsigned int i = 0; i < volVarIndices.size();++i)
- {
- if (fluxVarsCache.advectionSwitchFluxSign())
- derivativeMatrices[volVarIndices[i]][conti0EqIdx][pressureIdx] -= tij[i]*up;
- else
- derivativeMatrices[volVarIndices[i]][conti0EqIdx][pressureIdx] += tij[i]*up;
- }
+ addFluxDerivatives(derivativeMatrices, problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
}
//! Robin-type flux derivatives
diff --git a/dumux/porousmediumflow/fluxvariablescache.hh b/dumux/porousmediumflow/fluxvariablescache.hh
index 4edfbaf182bc0229ca31e84e6d2395362a129239..51c5ae3bfac45e0b62969c96b6639984d8e28c5f 100644
--- a/dumux/porousmediumflow/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/fluxvariablescache.hh
@@ -151,13 +151,29 @@ class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethod
, public EnergyCacheChooser<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>
{
using GridIndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
+
+ using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ static constexpr bool considerSecondary = MpfaHelper::considerSecondaryIVs();
public:
//! Returns whether or not this cache has been updated
bool isUpdated() const { return isUpdated_; }
+ //! Returns whether or not this cache is associated with a secondary interaction volume
+ //! Specialization for deactivated secondary interaction volumes
+ template< bool doSecondary = considerSecondary, std::enable_if_t<!doSecondary, int> = 0>
+ constexpr bool usesSecondaryIv() const { return false; }
+
+ //! Returns whether or not this cache is associated with a secondary interaction volume
+ //! Specialization for activated secondary interaction volumes
+ template< bool doSecondary = considerSecondary, std::enable_if_t<doSecondary, int> = 0>
+ bool usesSecondaryIv() const { return usesSecondaryIv_; }
+
//! 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
+ 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; }
@@ -168,6 +184,9 @@ private:
//! indicates if cache has been fully updated
bool isUpdated_ = false;
+ //! indicates if cache is associated with secondary interaction volume
+ bool usesSecondaryIv_ = false;
+
//! the index of the iv (this scvf is embedded in) in its container
GridIndexType ivIndexInContainer_;
};