From 1684060cf83584297e4833af7fc91b67c22fd2c1 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Sun, 29 Jan 2017 18:30:23 +0100
Subject: [PATCH] [mpfa] change to new flux vars cache concept
The flux variables cache now inherits from caches defined in the different laws.
These also define a respective filler class that is called when updating the flux
variables caches. That enables us to define different caches for different laws.
---
.../cellcentered/mpfa/darcyslaw.hh | 104 ++-
.../mpfa/elementfluxvariablescache.hh | 58 +-
.../cellcentered/mpfa/fickslaw.hh | 101 +-
.../mpfa/fluxvariablescachefiller.hh | 872 +++++++-----------
.../mpfa/fluxvariablescachefillerbase.hh | 247 -----
.../cellcentered/mpfa/fourierslaw.hh | 79 +-
.../mpfa/globalfluxvariablescache.hh | 20 +-
.../mpfa/interactionvolumebase.hh | 3 +
.../mpfa/omethod/interactionvolume.hh | 301 +++---
.../omethod/localsubcontrolentityseeds.hh | 2 +-
.../implicit/fluxvariablescache.hh | 258 +-----
11 files changed, 804 insertions(+), 1241 deletions(-)
delete mode 100644 dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh
diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 3aedae63d4..c26b0bca72 100644
--- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
@@ -53,38 +53,126 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
// Always use the dynamic type for vectors (compatibility with the boundary)
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
- using DynamicVector = typename BoundaryInteractionVolume::Vector;
-
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
- using Stencil = std::vector<IndexType>;
+ using CoefficientVector = typename BoundaryInteractionVolume::Vector;
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
static constexpr bool facetCoupling = GET_PROP_VALUE(TypeTag, MpfaFacetCoupling);
static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+ //! The cache used in conjunction with the mpfa Darcy's Law
+ class MpfaDarcysLawCache
+ {
+ // We always use the dynamic types here to be compatible on the boundary
+ using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
+ using PositionVector = typename BoundaryInteractionVolume::PositionVector;
+
+ public:
+ //! update cached objects
+ template<class InteractionVolume>
+ void updateAdvection(const InteractionVolume& iv, const SubControlVolumeFace &scvf)
+ {
+ const auto& localFaceData = iv.getLocalFaceData(scvf);
+ // update the quantities that are equal for all phases
+ advectionVolVarsStencil_ = iv.volVarsStencil();
+ advectionVolVarsPositions_ = iv.volVarsPositions();
+ advectionTij_ = iv.getTransmissibilities(localFaceData);
+
+ // we will need the neumann flux transformation only on interior boundaries with facet coupling
+ if (enableInteriorBoundaries && facetCoupling)
+ advectionCij_ = iv.getNeumannFluxTransformationCoefficients(localFaceData);
+
+ // The neumann fluxes always have to be set per phase
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ phaseNeumannFluxes_[phaseIdx] = iv.getNeumannFlux(localFaceData, phaseIdx);
+ }
+
+ //! Returns the volume variables indices necessary for flux computation
+ //! This includes all participating boundary volume variables. Since we
+ //! do not allow mixed BC for the mpfa this is the same for all phases.
+ const Stencil& advectionVolVarsStencil() const
+ { return advectionVolVarsStencil_; }
+
+ //! Returns the position on which the volume variables live. This is
+ //! necessary as we need to evaluate gravity also for the boundary volvars
+ const PositionVector& advectionVolVarsPositions() const
+ { return advectionVolVarsPositions_; }
+
+ //! Returns the transmissibilities associated with the volume variables
+ //! All phases flow through the same rock, thus, tij are equal for all phases
+ const CoefficientVector& advectionTij() const
+ { return advectionTij_; }
+
+ //! Returns the vector of coefficients with which the vector of neumann boundary conditions
+ //! has to be multiplied in order to transform them on the scvf this cache belongs to
+ const CoefficientVector& advectionCij() const
+ { return advectionCij_; }
+
+ //! If the useTpfaBoundary property is set to false, the boundary conditions
+ //! are put into the local systems leading to possible contributions on all faces
+ Scalar advectionNeumannFlux(unsigned int phaseIdx) const
+ { return phaseNeumannFluxes_[phaseIdx]; }
+
+ private:
+ // Quantities associated with advection
+ Stencil advectionVolVarsStencil_;
+ PositionVector advectionVolVarsPositions_;
+ CoefficientVector advectionTij_;
+ CoefficientVector advectionCij_;
+ std::array<Scalar, numPhases> phaseNeumannFluxes_;
+ };
+
+ //! Class that fills the cache corresponding to mpfa Darcy's Law
+ class MpfaDarcysLawCacheFiller
+ {
+ public:
+ //! Function to fill an MpfaDarcysLawCache of a given scvf
+ //! This interface has to be met by any advection-related cache filler class
+ template<class FluxVariablesCacheFiller>
+ static void fill(FluxVariablesCache& scvfFluxVarsCache,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const FluxVariablesCacheFiller& fluxVarsCacheFiller)
+ {
+ // get interaction volume from the flux vars cache filler & upate the cache
+ if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf))
+ scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.boundaryInteractionVolume(), scvf);
+ else
+ scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.interactionVolume(), scvf);
+ }
+ };
+
public:
// state the discretization method this implementation belongs to
static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;
+ // state the type for the corresponding cache and its filler
+ using Cache = MpfaDarcysLawCache;
+ using CacheFiller = MpfaDarcysLawCacheFiller;
+
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const unsigned int phaseIdx,
- const ElementFluxVarsCache& elemFluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache)
{
const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
@@ -160,7 +248,7 @@ public:
const auto& cij = fluxVarsCache.advectionCij();
// The vector of interior neumann fluxes
- DynamicVector facetCouplingFluxes(cij.size(), 0.0);
+ CoefficientVector facetCouplingFluxes(cij.size(), 0.0);
for (auto&& data : fluxVarsCache.interiorBoundaryData())
{
// Add additional Dirichlet fluxes for interior Dirichlet faces
diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
index 2bae60663e..954afca4c2 100644
--- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
@@ -140,7 +140,6 @@ public:
globalScvfIndices_.clear();
const auto& problem = globalFluxVarsCache().problem_();
- const auto& globalFvGeometry = problem.model().globalFvGeometry();
const auto globalI = problem.elementMapper().index(element);
const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI];
@@ -160,11 +159,17 @@ public:
for (auto scvfIdx : dataJ.scvfsJ)
globalScvfIndices_.push_back(scvfIdx);
- // prepare all the caches of the scvfs inside the corresponding interaction volumes using helper class
+ // helper class to fill flux variables caches
+ FluxVariablesCacheFiller filler(problem);
+
+ // prepare all the caches of the scvfs inside the corresponding interaction volumes
fluxVarsCache_.resize(globalScvfIndices_.size());
for (auto&& scvf : scvfs(fvGeometry))
- if (!(*this)[scvf].isUpdated())
- FluxVariablesCacheFiller::fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, *this);
+ {
+ auto& scvfCache = (*this)[scvf];
+ if (!scvfCache.isUpdated())
+ filler.fill(*this, scvfCache, element, fvGeometry, elemVolVars, scvf);
+ }
// prepare the caches in the remaining neighbors
for (auto&& dataJ : assemblyMapI)
@@ -172,10 +177,11 @@ public:
for (auto scvfIdx : dataJ.scvfsJ)
{
const auto& scvf = fvGeometry.scvf(scvfIdx);
- if (!(*this)[scvf].isUpdated())
+ auto& scvfCache = (*this)[scvf];
+ if (!scvfCache.isUpdated())
{
- auto elementJ = globalFvGeometry.element(dataJ.globalJ);
- FluxVariablesCacheFiller::fillFluxVarCache(problem, elementJ, fvGeometry, elemVolVars, scvf, *this);
+ auto elementJ = problem.model().globalFvGeometry().element(dataJ.globalJ);
+ filler.fill(*this, scvfCache, elementJ, fvGeometry, elemVolVars, scvf);
}
}
}
@@ -215,12 +221,40 @@ private:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
- for (auto&& scvf : scvfs(fvGeometry))
- (*this)[scvf].setUpdateStatus(false);
+ static const bool isSolIndependent = FluxVariablesCacheFiller::isSolutionIndependent();
- for (auto&& scvf : scvfs(fvGeometry))
- if (!(*this)[scvf].isUpdated())
- FluxVariablesCacheFiller::updateFluxVarCache(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf, *this);
+ if (!isSolIndependent)
+ {
+ const auto& problem = globalFluxVarsCache().problem_();
+
+ // helper class to fill flux variables caches
+ FluxVariablesCacheFiller filler(problem);
+
+ // set all the caches to "outdated"
+ for (auto& cache : fluxVarsCache_)
+ cache.setUpdateStatus(false);
+
+ // the global index of the element at hand
+ const auto globalI = problem.elementMapper().index(element);
+
+ // Let the filler do the update of the cache
+ for (unsigned int localScvfIdx = 0; localScvfIdx < globalScvfIndices_.size(); ++localScvfIdx)
+ {
+ const auto& scvf = fvGeometry.scvf(globalScvfIndices_[localScvfIdx]);
+
+ auto& scvfCache = fluxVarsCache_[localScvfIdx];
+ if (!scvfCache.isUpdated())
+ {
+ // obtain the corresponding element
+ const auto scvfInsideScvIdx = scvf.insideScvIdx();
+ const auto insideElement = scvfInsideScvIdx == globalI ?
+ element :
+ problem.model().globalFvGeometry().element(scvfInsideScvIdx);
+
+ filler.update(*this, scvfCache, insideElement, fvGeometry, elemVolVars, scvf);
+ }
+ }
+ }
}
// get index of an scvf in the local container
diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
index 6aaf11e188..34cdec1077 100644
--- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
@@ -54,10 +54,11 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
// Always use the dynamic type for vectors (compatibility with the boundary)
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
- using DynamicVector = typename BoundaryInteractionVolume::Vector;
+ using CoefficientVector = typename BoundaryInteractionVolume::Vector;
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
@@ -67,10 +68,106 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
static constexpr bool facetCoupling = GET_PROP_VALUE(TypeTag, MpfaFacetCoupling);
static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+ //! The cache used in conjunction with the mpfa Fick's Law
+ class MpfaFicksLawCache
+ {
+ static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+ // We always use the dynamic types here to be compatible on the boundary
+ using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
+ using PositionVector = typename BoundaryInteractionVolume::PositionVector;
+
+ public:
+ //! The constructor. Initializes the Neumann flux to zero
+ MpfaFicksLawCache() { componentNeumannFluxes_.fill(0.0); }
+
+ // update cached objects for the diffusive fluxes
+ template<typename InteractionVolume>
+ void updateDiffusion(const InteractionVolume& iv, const SubControlVolumeFace &scvf,
+ unsigned int phaseIdx, unsigned int compIdx)
+ {
+ const auto& localFaceData = iv.getLocalFaceData(scvf);
+ diffusionVolVarsStencils_[phaseIdx][compIdx] = iv.volVarsStencil();
+ diffusionTij_[phaseIdx][compIdx] = iv.getTransmissibilities(localFaceData);
+
+ if (enableInteriorBoundaries)
+ diffusionCij_[phaseIdx][compIdx] = iv.getNeumannFluxTransformationCoefficients(localFaceData);
+
+ //! For compositional models, we associate neumann fluxes with the phases (main components)
+ //! This is done in the AdvectionCache. However, in single-phasic models we solve the phase AND
+ //! the component mass balance equations. Thus, in this case we have diffusive neumann contributions.
+ //! we assume compIdx = eqIdx
+ if (numPhases == 1 && phaseIdx != compIdx)
+ componentNeumannFluxes_[compIdx] = iv.getNeumannFlux(localFaceData, compIdx);
+
+ }
+
+ //! Returns the volume variables indices necessary for diffusive flux
+ //! computation. This includes all participating boundary volume variables
+ //! and it can be different for the phases & components.
+ const Stencil& diffusionVolVarsStencil(unsigned int phaseIdx, unsigned int compIdx) const
+ { return diffusionVolVarsStencils_[phaseIdx][compIdx]; }
+
+ //! Returns the transmissibilities associated with the volume variables
+ //! This can be different for the phases & components.
+ const CoefficientVector& diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const
+ { return diffusionTij_[phaseIdx][compIdx]; }
+
+ //! Returns the vector of coefficients with which the vector of neumann boundary conditions
+ //! has to be multiplied in order to transform them on the scvf this cache belongs to
+ const CoefficientVector& diffusionCij(unsigned int phaseIdx, unsigned int compIdx) const
+ { return diffusionCij_[phaseIdx][compIdx]; }
+
+ //! If the useTpfaBoundary property is set to false, the boundary conditions
+ //! are put into the local systems leading to possible contributions on all faces
+ Scalar componentNeumannFlux(unsigned int compIdx) const
+ {
+ assert(numPhases == 1);
+ return componentNeumannFluxes_[compIdx];
+ }
+
+ private:
+ // Quantities associated with molecular diffusion
+ std::array< std::array<Stencil, numComponents>, numPhases> diffusionVolVarsStencils_;
+ std::array< std::array<CoefficientVector, numComponents>, numPhases> diffusionTij_;
+ std::array< std::array<CoefficientVector, numComponents>, numPhases> diffusionCij_;
+
+ // diffusive neumann flux for single-phasic models
+ std::array<Scalar, numComponents> componentNeumannFluxes_;
+ };
+
+ //! Class that fills the cache corresponding to mpfa Darcy's Law
+ class MpfaFicksLawCacheFiller
+ {
+ public:
+ //! Function to fill an MpfaFicksLawCache of a given scvf
+ //! This interface has to be met by any diffusion-related cache filler class
+ template<class FluxVariablesCacheFiller>
+ static void fill(FluxVariablesCache& scvfFluxVarsCache,
+ unsigned int phaseIdx, unsigned int compIdx,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const FluxVariablesCacheFiller& fluxVarsCacheFiller)
+ {
+ // get interaction volume from the flux vars cache filler & upate the cache
+ if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf))
+ scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.boundaryInteractionVolume(), scvf, phaseIdx, compIdx);
+ else
+ scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.interactionVolume(), scvf, phaseIdx, compIdx);
+ }
+ };
+
public:
// state the discretization method this implementation belongs to
static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;
+ // state the type for the corresponding cache and its filler
+ using Cache = MpfaFicksLawCache;
+ using CacheFiller = MpfaFicksLawCacheFiller;
+
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -180,7 +277,7 @@ public:
const auto& cij = fluxVarsCache.diffusionCij(phaseIdx, compIdx);
// The vector of interior neumann fluxes
- DynamicVector facetCouplingFluxes(cij.size(), 0.0);
+ CoefficientVector facetCouplingFluxes(cij.size(), 0.0);
for (auto&& data : fluxVarsCache.interiorBoundaryData())
{
// Add additional Dirichlet fluxes for interior Dirichlet faces
diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index 0ca84820d8..79504a910e 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -20,616 +20,422 @@
* \file
* \brief The flux variables cache filler class
*/
-#ifndef DUMUX_DISCRETIZATION_CCMPFA_GLOBAL_FLUXVARSCACHE_FILLER_HH
-#define DUMUX_DISCRETIZATION_CCMPFA_GLOBAL_FLUXVARSCACHE_FILLER_HH
+#ifndef DUMUX_DISCRETIZATION_CCMPFA_FLUXVARSCACHE_FILLER_HH
+#define DUMUX_DISCRETIZATION_CCMPFA_FLUXVARSCACHE_FILLER_HH
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/methods.hh>
-#include "fluxvariablescachefillerbase.hh"
namespace Dumux
{
-//! Forward declaration of the actual implementation
-template<class TypeTag, bool advection, bool diffusion, bool energy>
-class CCMpfaFluxVariablesCacheFillerImplementation;
+
+//! forward declaration of properties
+namespace Properties
+{
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(NumComponents);
+NEW_PROP_TAG(ThermalConductivityModel);
+};
/*!
* \ingroup ImplicitModel
* \brief Helper class to fill the flux var caches
*/
template<class TypeTag>
-using CCMpfaFluxVariablesCacheFiller = CCMpfaFluxVariablesCacheFillerImplementation<TypeTag,
- GET_PROP_VALUE(TypeTag, EnableAdvection),
- GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
- GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
-
-//! Implementation for purely advective problems
-template<class TypeTag>
-class CCMpfaFluxVariablesCacheFillerImplementation<TypeTag, true, false, false> : public CCMpfaAdvectionCacheFiller<TypeTag>
+class CCMpfaFluxVariablesCacheFiller
{
- using AdvectionFiller = CCMpfaAdvectionCacheFiller<TypeTag>;
-
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using Element = typename GridView::template Codim<0>::Entity;
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
- static constexpr bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection);
- static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+ static constexpr bool doAdvection = GET_PROP_VALUE(TypeTag, EnableAdvection);
+ static constexpr bool doDiffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion);
+ static constexpr bool doHeatConduction = GET_PROP_VALUE(TypeTag, EnableEnergyBalance);
+
+ static constexpr bool soldependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection);
+ static constexpr bool soldependentDiffusion = GET_PROP_VALUE(TypeTag, SolutionDependentMolecularDiffusion);
+ static constexpr bool soldependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction);
+
+ enum ProcessIndices : unsigned int
+ {
+ advectionIdx,
+ diffusionIdx,
+ heatConductionIdx
+ };
public:
- //! function to fill the flux var caches
- template<class FluxVarsCacheContainer>
- static void fillFluxVarCache(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- FluxVarsCacheContainer& fluxVarsCacheContainer)
+ //! The constructor. Sets the problem pointer
+ CCMpfaFluxVariablesCacheFiller(const Problem& problem) : problemPtr_(&problem) {}
+
+ /*!
+ * \brief function to fill the flux variables caches
+ *
+ * \param fluxVarsCacheContainer Either the element or global flux variables cache
+ * \param scvfFluxVarsCache The flux var cache to be updated corresponding to the given scvf
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The corresponding sub-control volume face
+ * \param doSubCaches Array of bools indicating which sub caches have to be updated
+ */
+ template<class FluxVariablesCacheContainer>
+ void fill(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const std::array<bool, 3>& doSubCaches = std::array<bool, 3>({true, true, true}))
{
- // Instantiate interaction volume depending on if scvf touches the boundary or not
- if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf))
+ // Set pointers
+ elementPtr_ = &element;
+ fvGeometryPtr_ = &fvGeometry;
+ elemVolVarsPtr_ = &elemVolVars;
+ scvfPtr_ = &scvf;
+
+ // prepare interaction volume and fill caches of all the scvfs connected to it
+ const auto& globalFvGeometry = problem().model().globalFvGeometry();
+ if (globalFvGeometry.touchesInteriorOrDomainBoundary(scvf))
{
- BoundaryInteractionVolume iv(problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf),
- problem,
- fvGeometry,
- elemVolVars);
-
- // forward to the filler for the advective quantities
- AdvectionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- //set update status and maybe update data on interior boundaries
- for (const auto scvfIdxJ : iv.globalScvfs())
- {
- if (enableInteriorBoundaries)
- fluxVarsCacheContainer[scvfIdxJ].updateInteriorBoundaryData(iv, fvGeometry.scvf(scvfIdxJ));
- fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true);
- }
+ bIv_ = std::make_unique<BoundaryInteractionVolume>(globalFvGeometry.boundaryInteractionVolumeSeed(scvf),
+ problem(),
+ fvGeometry,
+ elemVolVars);
+
+ // fill the caches for all the scvfs in the interaction volume
+ fillCachesInInteractionVolume_(fluxVarsCacheContainer, scvfFluxVarsCache, boundaryInteractionVolume(), doSubCaches);
}
else
{
- InteractionVolume iv(problem.model().globalFvGeometry().interactionVolumeSeed(scvf),
- problem,
- fvGeometry,
- elemVolVars);
-
- // forward to the filler for the advective quantities
- AdvectionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // set update status
- for (const auto scvfIdxJ : iv.globalScvfs())
- fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true);
+ iv_ = std::make_unique<InteractionVolume>(globalFvGeometry.interactionVolumeSeed(scvf),
+ problem(),
+ fvGeometry,
+ elemVolVars);
+
+ // fill the caches for all the scvfs in the interaction volume
+ fillCachesInInteractionVolume_(fluxVarsCacheContainer, scvfFluxVarsCache, interactionVolume(), doSubCaches);
}
}
- //! function to update the flux var caches during derivative calculation
- template<class FluxVarsCacheContainer>
- static void updateFluxVarCache(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- FluxVarsCacheContainer& fluxVarsCacheContainer)
+ /*!
+ * \brief function to update the flux variables caches during derivative calculation
+ *
+ * \copydoc fill
+ */
+ template<class FluxVariablesCacheContainer>
+ void update(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
{
- //! If advection is solution-independent, only update the caches for scvfs that touch
- //! a boundary as we have to update the boundary contributions (possibly solution-dependent)
- if (!solDependentAdvection)
+ // array of bool with which we indicate the sub-caches which have to be
+ // filled. During update, we only update solution-dependent quantities.
+ static const std::array<bool, 3> doSubCaches = []()
{
- const bool touchesBoundary = problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf);
-
- //! Do the whole update where a boundary is touched
- if (!useTpfaBoundary && touchesBoundary)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer);
- //! Simply tell the cache that it is up to date
- else
- fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true);
- }
- else
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer);
+ std::array<bool, 3> doCaches;
+ doCaches[ProcessIndices::advectionIdx] = doAdvection && soldependentAdvection;
+ doCaches[ProcessIndices::diffusionIdx] = doDiffusion && soldependentDiffusion;
+ doCaches[ProcessIndices::heatConductionIdx] = doHeatConduction && soldependentHeatConduction;
+ return doCaches;
+ } ();
+
+ // forward to fill routine
+ fill(fluxVarsCacheContainer, scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf, doSubCaches);
}
-};
-//! Implementation for problems considering advection & diffusion
-template<class TypeTag>
-class CCMpfaFluxVariablesCacheFillerImplementation<TypeTag, true, true, false> : public CCMpfaAdvectionCacheFiller<TypeTag>,
- public CCMpfaDiffusionCacheFiller<TypeTag>
-{
- using AdvectionFiller = CCMpfaAdvectionCacheFiller<TypeTag>;
- using DiffusionFiller = CCMpfaDiffusionCacheFiller<TypeTag>;
+ static bool isSolutionIndependent()
+ {
+ static const bool isSolDependent = (doAdvection && soldependentAdvection) ||
+ (doDiffusion && soldependentDiffusion) ||
+ (doHeatConduction && soldependentHeatConduction);
+ return !isSolDependent;
+ }
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
- using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
- using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+ const InteractionVolume& interactionVolume() const
+ { return *iv_.get(); }
- using Element = typename GridView::template Codim<0>::Entity;
+ const BoundaryInteractionVolume& boundaryInteractionVolume() const
+ { return *bIv_.get(); }
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+private:
- static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
- static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+ const Problem& problem() const
+ { return *problemPtr_; }
- static constexpr bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection);
- static constexpr bool advectionUsesMpfa = AdvectionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
+ const Element& element() const
+ { return *elementPtr_; }
- static constexpr bool solDependentDiffusion = GET_PROP_VALUE(TypeTag, SolutionDependentMolecularDiffusion);
- static constexpr bool diffusionUsesMpfa = MolecularDiffusionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
+ const FVElementGeometry& fvGeometry() const
+ { return *fvGeometryPtr_; }
- //! Whether or not the caches have to be updated after solution deflection depends on
- //! - the solution dependency of the parameters
- //! - if tpfa is used on the boundaries (advection has to be updated because of neumann boundaries)
- //! - if single-phasic compositional model is used
- //! (diffusion has to be updated because of neumann boundaries for the transported component)
- static constexpr bool doAdvectionUpdate = (solDependentAdvection || !useTpfaBoundary) && advectionUsesMpfa ;
- static constexpr bool doDiffusionUpdate = (solDependentDiffusion || (numPhases == 1 && !useTpfaBoundary)) && diffusionUsesMpfa;
+ const ElementVolumeVariables& elemVolVars() const
+ { return *elemVolVarsPtr_; }
- using BoolPair = std::pair<bool, bool>;
+ const SubControlVolumeFace& scvFace() const
+ { return *scvfPtr_; }
-public:
- //! function to fill the flux var caches
- template<class FluxVarsCacheContainer>
- static void fillFluxVarCache(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- FluxVarsCacheContainer& fluxVarsCacheContainer,
- const BoolPair doAdvectionOrDiffusion = BoolPair(true, true))
+ InteractionVolume& interactionVolume()
+ { return *iv_.get(); }
+
+ BoundaryInteractionVolume& boundaryInteractionVolume()
+ { return *bIv_.get(); }
+
+ //! Method to fill the flux var caches within an interaction volume
+ template<class FluxVariablesCacheContainer, class InteractionVolumeType>
+ void fillCachesInInteractionVolume_(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ InteractionVolumeType& iv,
+ const std::array<bool, 3>& doSubCaches)
{
- // Instantiate interaction volume depending on if scvf touches the boundary or not
- if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf))
+ // First we upate the interior boundary data and set the update status.
+ // We store pointers to the other flux var caches and the elements they are embedded in simultaneously.
+ // This way we have to obtain this data only once and can use it again in the sub-cache fillers.
+ std::vector<FluxVariablesCache*> otherFluxVarCaches(iv.globalLocalScvfPairedData().size()-1);
+ std::vector<Element> otherElements(iv.globalLocalScvfPairedData().size()-1);
+
+ scvfFluxVarsCache.updateInteriorBoundaryData(iv, scvFace());
+ scvfFluxVarsCache.setUpdateStatus(true);
+
+ const auto curScvfIdx = scvFace().index();
+ unsigned int otherScvfIdx = 0;
+ for (const auto& dataPair : iv.globalLocalScvfPairedData())
{
- BoundaryInteractionVolume iv(problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf),
- problem,
- fvGeometry,
- elemVolVars);
-
- // forward to the filler for the advective quantities
- if (doAdvectionOrDiffusion.first)
- AdvectionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // forward to the filler for the diffusive quantities
- if (doAdvectionOrDiffusion.second)
- DiffusionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- //set update status and maybe update data on interior boundaries
- for (const auto scvfIdxJ : iv.globalScvfs())
- {
- if (enableInteriorBoundaries)
- fluxVarsCacheContainer[scvfIdxJ].updateInteriorBoundaryData(iv, fvGeometry.scvf(scvfIdxJ));
- fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true);
- }
- }
- else
- {
- InteractionVolume iv(problem.model().globalFvGeometry().interactionVolumeSeed(scvf),
- problem,
- fvGeometry,
- elemVolVars);
-
- // forward to the filler for the advective quantities
- if (doAdvectionOrDiffusion.first)
- AdvectionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // forward to the filler for the diffusive quantities
- if (doAdvectionOrDiffusion.second)
- DiffusionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // set update status
- for (const auto scvfIdxJ : iv.globalScvfs())
- fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true);
+ const auto& scvfJ = *dataPair.first;
+ if (curScvfIdx == scvfJ.index())
+ continue;
+
+ // get the element scvfJ is embedded in
+ const auto scvfJInsideScvIndex = scvfJ.insideScvIdx();
+ otherElements[otherScvfIdx] = scvfJInsideScvIndex == scvFace().insideScvIdx() ?
+ element() :
+ problem().model().globalFvGeometry().element(scvfJInsideScvIndex);
+
+ // get the corresponding flux var cache
+ otherFluxVarCaches[otherScvfIdx] = &fluxVarsCacheContainer[scvfJ];
+ otherFluxVarCaches[otherScvfIdx]->updateInteriorBoundaryData(iv, scvfJ);
+ otherFluxVarCaches[otherScvfIdx]->setUpdateStatus(true);
+ otherScvfIdx++;
}
- }
- //! function to update the flux var caches during derivative calculation
- template<class FluxVarsCacheContainer>
- static void updateFluxVarCache(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- FluxVarsCacheContainer& fluxVarsCacheContainer)
- {
- const bool touchesDomainBoundary = problem.model().globalFvGeometry().touchesDomainBoundary(scvf);
-
- // maybe update Advection or diffusion or both
- const bool updateAdvection = doAdvectionUpdate && touchesDomainBoundary;
- const bool updateDiffusion = doDiffusionUpdate && touchesDomainBoundary;
+ //! Maybe update the advective quantities
+ if (doSubCaches[ProcessIndices::advectionIdx])
+ fillAdvection(fluxVarsCacheContainer, scvfFluxVarsCache, iv, otherFluxVarCaches, otherElements);
- if (updateAdvection && updateDiffusion)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolPair(true, true));
- else if (updateAdvection && !updateDiffusion)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolPair(true, false));
- else if (!updateAdvection && updateDiffusion)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolPair(false, true));
+ //! Maybe update the diffusive quantities
+ if (doSubCaches[ProcessIndices::diffusionIdx])
+ fillDiffusion(fluxVarsCacheContainer, scvfFluxVarsCache, iv, otherFluxVarCaches, otherElements);
- // tell the cache that it has been updated in any case
- fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true);
+ //! Maybe update quantities related to heat conduction
+ if (doSubCaches[ProcessIndices::heatConductionIdx])
+ fillHeatConduction(fluxVarsCacheContainer, scvfFluxVarsCache, iv, otherFluxVarCaches, otherElements);
}
-};
-
-//! Implementation for problems considering advection & heat conduction
-template<class TypeTag>
-class CCMpfaFluxVariablesCacheFillerImplementation<TypeTag, true, false, true> : public CCMpfaAdvectionCacheFiller<TypeTag>,
- public CCMpfaHeatConductionCacheFiller<TypeTag>
-{
- using AdvectionFiller = CCMpfaAdvectionCacheFiller<TypeTag>;
- using HeatConductionFiller = CCMpfaHeatConductionCacheFiller<TypeTag>;
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
- using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
+ //! method to fill the advective quantities
+ template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool advectionEnabled = doAdvection>
+ typename std::enable_if<advectionEnabled>::type
+ fillAdvection(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ InteractionVolumeType& iv,
+ const std::vector<FluxVariablesCache*> otherFluxVarCaches,
+ const std::vector<Element> otherElements)
+ {
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using AdvectionFiller = typename AdvectionType::CacheFiller;
- using Element = typename GridView::template Codim<0>::Entity;
+ static constexpr bool usesMpfa = AdvectionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
- static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
- static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+ // maybe solve the local system subject to K
+ if (usesMpfa)
+ {
+ auto K = [] (const Element& element, const auto& volVars, const SubControlVolume& scv)
+ { return volVars.permeability(); };
+ iv.solveLocalSystem(K);
+ }
- static constexpr bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection);
- static constexpr bool advectionUsesMpfa = AdvectionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
+ // fill the caches of all scvfs within this interaction volume
+ AdvectionFiller::fill(scvfFluxVarsCache, problem(), element(), fvGeometry(), elemVolVars(), scvFace(), *this);
- static constexpr bool solDependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction);
- static constexpr bool heatConductionUsesMpfa = HeatConductionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
+ unsigned int otherScvfIdx = 0;
+ const auto curScvfIdx = scvFace().index();
+ for (const auto& dataPair : iv.globalLocalScvfPairedData())
+ {
+ const auto& scvfJ = *dataPair.first;
+ if (curScvfIdx == scvfJ.index())
+ continue;
+
+ // fill corresponding cache
+ AdvectionFiller::fill(*otherFluxVarCaches[otherScvfIdx],
+ problem(),
+ otherElements[otherScvfIdx],
+ fvGeometry(),
+ elemVolVars(),
+ scvfJ,
+ *this);
+ otherScvfIdx++;
+ }
+ }
- //! Whether or not the caches have to be updated after solution deflection depends on
- //! - the solution dependency of the parameters
- //! - if tpfa is used on the boundaries (advection has to be updated because of neumann boundaries)
- static constexpr bool doAdvectionUpdate = (solDependentAdvection || !useTpfaBoundary) && advectionUsesMpfa ;
- static constexpr bool doHeatConductionUpdate = (solDependentHeatConduction || !useTpfaBoundary) && heatConductionUsesMpfa;
+ //! do nothing if advection is not enabled
+ template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool advectionEnabled = doAdvection>
+ typename std::enable_if<!advectionEnabled>::type
+ fillAdvection(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ InteractionVolumeType& iv,
+ const std::vector<FluxVariablesCache*> otherFluxVarCaches,
+ const std::vector<Element> otherElements)
+ {}
+
+ //! method to fill the diffusive quantities
+ template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool diffusionEnabled = doDiffusion>
+ typename std::enable_if<diffusionEnabled>::type
+ fillDiffusion(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ InteractionVolumeType& iv,
+ const std::vector<FluxVariablesCache*> otherFluxVarCaches,
+ const std::vector<Element> otherElements)
+ {
+ using DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+ using DiffusionFiller = typename DiffusionType::CacheFiller;
- using BoolPair = std::pair<bool, bool>;
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+ static constexpr bool usesMpfa = DiffusionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
-public:
- //! function to fill the flux var caches
- template<class FluxVarsCacheContainer>
- static void fillFluxVarCache(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- FluxVarsCacheContainer& fluxVarsCacheContainer,
- const BoolPair doAdvectionOrHeatConduction = BoolPair(true, true))
- {
- // Instantiate interaction volume depending on if scvf touches the boundary or not
- if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf))
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- BoundaryInteractionVolume iv(problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf),
- problem,
- fvGeometry,
- elemVolVars);
-
- // forward to the filler for the advective quantities
- if (doAdvectionOrHeatConduction.first)
- AdvectionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // forward to the filler for the heat conduction quantities
- if (doAdvectionOrHeatConduction.second)
- HeatConductionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- //set update status and maybe update data on interior boundaries
- for (const auto scvfIdxJ : iv.globalScvfs())
+ for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
{
- if (enableInteriorBoundaries)
- fluxVarsCacheContainer[scvfIdxJ].updateInteriorBoundaryData(iv, fvGeometry.scvf(scvfIdxJ));
- fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true);
+ if (phaseIdx == compIdx)
+ continue;
+
+ // solve the local system subject to the diffusion tensor (if uses mpfa)
+ if (usesMpfa)
+ {
+ auto D = [phaseIdx, compIdx](const Element& element, const auto& volVars, const SubControlVolume& scv)
+ { return volVars.diffusionCoefficient(phaseIdx, compIdx); };
+ iv.solveLocalSystem(D);
+ }
+
+ // fill the caches of all scvfs within this interaction volume
+ DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element(), fvGeometry(), elemVolVars(), scvFace(), *this);
+
+ unsigned int otherScvfIdx = 0;
+ const auto curScvfIdx = scvFace().index();
+ for (const auto& dataPair : iv.globalLocalScvfPairedData())
+ {
+ const auto& scvfJ = *dataPair.first;
+ if (curScvfIdx == scvfJ.index())
+ continue;
+
+ // fill corresponding cache
+ DiffusionFiller::fill(*otherFluxVarCaches[otherScvfIdx],
+ phaseIdx,
+ compIdx,
+ problem(),
+ otherElements[otherScvfIdx],
+ fvGeometry(),
+ elemVolVars(),
+ scvfJ,
+ *this);
+ otherScvfIdx++;
+ }
}
}
- else
- {
- InteractionVolume iv(problem.model().globalFvGeometry().interactionVolumeSeed(scvf),
- problem,
- fvGeometry,
- elemVolVars);
-
- // forward to the filler for the advective quantities
- if (doAdvectionOrHeatConduction.first)
- AdvectionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // forward to the filler for the heat conduction quantities
- if (doAdvectionOrHeatConduction.second)
- HeatConductionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // set update status
- for (const auto scvfIdxJ : iv.globalScvfs())
- fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true);
- }
}
- //! function to update the flux var caches during derivative calculation
- template<class FluxVarsCacheContainer>
- static void updateFluxVarCache(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- FluxVarsCacheContainer& fluxVarsCacheContainer)
+ //! do nothing if diffusion is not enabled
+ template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool diffusionEnabled = doDiffusion>
+ typename std::enable_if<!diffusionEnabled>::type
+ fillDiffusion(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ InteractionVolumeType& iv,
+ const std::vector<FluxVariablesCache*> otherFluxVarCaches,
+ const std::vector<Element> otherElements)
+ {}
+
+ //! method to fill the quantities related to heat conduction
+ template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool heatConductionEnabled = doHeatConduction>
+ typename std::enable_if<heatConductionEnabled>::type
+ fillHeatConduction(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ InteractionVolumeType& iv,
+ const std::vector<FluxVariablesCache*> otherFluxVarCaches,
+ const std::vector<Element> otherElements)
{
- const bool touchesDomainBoundary = problem.model().globalFvGeometry().touchesDomainBoundary(scvf);
-
- // maybe update Advection or diffusion or both
- const bool updateAdvection = doAdvectionUpdate && touchesDomainBoundary;
- const bool updateHeatConduction = doHeatConductionUpdate && touchesDomainBoundary;
-
- if (updateAdvection && updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolPair(true, true));
- else if (updateAdvection && !updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolPair(true, false));
- else if (!updateAdvection && updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolPair(false, true));
-
- // tell the cache that it has been updated in any case
- fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true);
- }
-};
+ using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
+ using HeatConductionFiller = typename HeatConductionType::CacheFiller;
+ using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
-//! Implementation for problems considering advection, diffusion & heat conduction
-template<class TypeTag>
-class CCMpfaFluxVariablesCacheFillerImplementation<TypeTag, true, true, true> : public CCMpfaAdvectionCacheFiller<TypeTag>,
- public CCMpfaDiffusionCacheFiller<TypeTag>,
- public CCMpfaHeatConductionCacheFiller<TypeTag>
-{
- using AdvectionFiller = CCMpfaAdvectionCacheFiller<TypeTag>;
- using DiffusionFiller = CCMpfaDiffusionCacheFiller<TypeTag>;
- using HeatConductionFiller = CCMpfaHeatConductionCacheFiller<TypeTag>;
+ static constexpr bool usesMpfa = HeatConductionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
- using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
- using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
- using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
-
- using Element = typename GridView::template Codim<0>::Entity;
-
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
-
- static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
- static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
-
- static constexpr bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection);
- static constexpr bool advectionUsesMpfa = AdvectionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
-
- static constexpr bool solDependentDiffusion = GET_PROP_VALUE(TypeTag, SolutionDependentMolecularDiffusion);
- static constexpr bool diffusionUsesMpfa = MolecularDiffusionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
-
- static constexpr bool solDependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction);
- static constexpr bool heatConductionUsesMpfa = HeatConductionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
-
- //! - if single-phasic compositional model is used
- //! (diffusion has to be updated because of neumann boundaries for the transported component)
- static constexpr bool doAdvectionUpdate = (solDependentAdvection || !useTpfaBoundary) && advectionUsesMpfa ;
- static constexpr bool doDiffusionUpdate = (solDependentDiffusion || (numPhases == 1 && !useTpfaBoundary)) && diffusionUsesMpfa;
- static constexpr bool doHeatConductionUpdate = (solDependentHeatConduction || !useTpfaBoundary) && heatConductionUsesMpfa;
-
- using BoolTriplet = std::array<bool, 3>;
-
-public:
- //! function to fill the flux var caches
- template<class FluxVarsCacheContainer>
- static void fillFluxVarCache(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- FluxVarsCacheContainer& fluxVarsCacheContainer,
- const BoolTriplet doAdvectionOrDiffusionOrHeatConduction = BoolTriplet({true, true, true}))
- {
- // Instantiate interaction volume depending on if scvf touches the boundary or not
- if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf))
+ // maybe solve the local system subject to K
+ if (usesMpfa)
{
- BoundaryInteractionVolume iv(problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf),
- problem,
- fvGeometry,
- elemVolVars);
-
- // forward to the filler for the advective quantities
- if (doAdvectionOrDiffusionOrHeatConduction[0])
- AdvectionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // forward to the filler for the diffusive quantities
- if (doAdvectionOrDiffusionOrHeatConduction[1])
- DiffusionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // forward to the filler for the diffusive quantities
- if (doAdvectionOrDiffusionOrHeatConduction[2])
- HeatConductionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- //set update status and maybe update data on interior boundaries
- for (const auto scvfIdxJ : iv.globalScvfs())
- {
- if (enableInteriorBoundaries)
- fluxVarsCacheContainer[scvfIdxJ].updateInteriorBoundaryData(iv, fvGeometry.scvf(scvfIdxJ));
- fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true);
- }
+ const auto& prob = problem();
+ const auto& fvGeom = fvGeometry();
+ auto lambda = [&prob, &fvGeom](const Element& element, const auto& volVars, const SubControlVolume& scv)
+ { return ThermalConductivityModel::effectiveThermalConductivity(volVars,
+ prob.spatialParams(),
+ element,
+ fvGeom,
+ scv); };
+ iv.solveLocalSystem(lambda);
}
- else
+
+ // fill the caches of all scvfs within this interaction volume
+ HeatConductionFiller::fill(scvfFluxVarsCache, problem(), element(), fvGeometry(), elemVolVars(), scvFace(), *this);
+
+ unsigned int otherScvfIdx = 0;
+ const auto curScvfIdx = scvFace().index();
+ for (const auto& dataPair : iv.globalLocalScvfPairedData())
{
- InteractionVolume iv(problem.model().globalFvGeometry().interactionVolumeSeed(scvf),
- problem,
- fvGeometry,
- elemVolVars);
-
- // forward to the filler for the advective quantities
- if (doAdvectionOrDiffusionOrHeatConduction[0])
- AdvectionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // forward to the filler for the diffusive quantities
- if (doAdvectionOrDiffusionOrHeatConduction[1])
- DiffusionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // forward to the filler for the diffusive quantities
- if (doAdvectionOrDiffusionOrHeatConduction[2])
- HeatConductionFiller::fillCaches(problem,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- iv,
- fluxVarsCacheContainer);
-
- // set update status
- for (const auto scvfIdxJ : iv.globalScvfs())
- fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true);
+ const auto& scvfJ = *dataPair.first;
+ if (curScvfIdx == scvfJ.index())
+ continue;
+
+ // fill corresponding cache
+ HeatConductionFiller::fill(*otherFluxVarCaches[otherScvfIdx],
+ problem(),
+ otherElements[otherScvfIdx],
+ fvGeometry(),
+ elemVolVars(),
+ scvfJ,
+ *this);
+ otherScvfIdx++;
}
}
- //! function to update the flux var caches during derivative calculation
- template<class FluxVarsCacheContainer>
- static void updateFluxVarCache(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- FluxVarsCacheContainer& fluxVarsCacheContainer)
- {
- const bool touchesDomainBoundary = problem.model().globalFvGeometry().touchesDomainBoundary(scvf);
-
- // maybe update Advection or diffusion or both
- const bool updateAdvection = doAdvectionUpdate && touchesDomainBoundary;
- const bool updateDiffusion = doDiffusionUpdate && touchesDomainBoundary;
- const bool updateHeatConduction = doHeatConductionUpdate && touchesDomainBoundary;
-
- if (updateAdvection && updateDiffusion && updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolTriplet({true, true, true}));
- else if (updateAdvection && updateDiffusion && !updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolTriplet({true, true, false}));
- else if (updateAdvection && !updateDiffusion && !updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolTriplet({true, false, false}));
- else if (updateAdvection && !updateDiffusion && updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolTriplet({true, false, true}));
- else if (!updateAdvection && updateDiffusion && updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolTriplet({false, true, true}));
- else if (!updateAdvection && updateDiffusion && !updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolTriplet({false, true, false}));
- else if (!updateAdvection && !updateDiffusion && !updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolTriplet({false, false, false}));
- else if (!updateAdvection && !updateDiffusion && updateHeatConduction)
- fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer, BoolTriplet({false, false, true}));
-
-
- // tell the cache that it has been updated in any case
- fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true);
- }
+ //! do nothing if heat conduction is disabled
+ template<class FluxVariablesCacheContainer, class InteractionVolumeType, bool heatConductionEnabled = doHeatConduction>
+ typename std::enable_if<!heatConductionEnabled>::type
+ fillHeatConduction(FluxVariablesCacheContainer& fluxVarsCacheContainer,
+ FluxVariablesCache& scvfFluxVarsCache,
+ InteractionVolumeType& iv,
+ const std::vector<FluxVariablesCache*> otherFluxVarCaches,
+ const std::vector<Element> otherElements)
+ {}
+
+ const Problem* problemPtr_;
+ const Element* elementPtr_;
+ const FVElementGeometry* fvGeometryPtr_;
+ const ElementVolumeVariables* elemVolVarsPtr_;
+ const SubControlVolumeFace* scvfPtr_;
+
+ // We store pointers to an inner and boundary interaction volume
+ // these are updated during the filling of the caches and the
+ // physics-related caches have access to them
+ std::unique_ptr<InteractionVolume> iv_;
+ std::unique_ptr<BoundaryInteractionVolume> bIv_;
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh
deleted file mode 100644
index 9ef0820f3a..0000000000
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh
+++ /dev/null
@@ -1,247 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Base classes for the flux variables cache filler
- */
-#ifndef DUMUX_DISCRETIZATION_CCMPFA_GLOBAL_FLUXVARSCACHE_FILLER_BASE_HH
-#define DUMUX_DISCRETIZATION_CCMPFA_GLOBAL_FLUXVARSCACHE_FILLER_BASE_HH
-
-#include <dumux/implicit/properties.hh>
-#include <dumux/discretization/methods.hh>
-
-namespace Dumux
-{
-
-//! Forward declaration of property tags
-namespace Properties
-{
-NEW_PROP_TAG(Indices);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(NumComponents);
-NEW_PROP_TAG(ThermalConductivityModel);
-NEW_PROP_TAG(EffectiveDiffusivityModel);
-}
-
-//! Fills the advective quantities in the caches
-template<class TypeTag>
-class CCMpfaAdvectionCacheFiller
-{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
-
- using Element = typename GridView::template Codim<0>::Entity;
-
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
- static constexpr bool enableDiffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion);
-
-public:
- //! function to fill the flux var caches
- template<class FluxVarsCacheContainer, class InteractionVolume>
- static void fillCaches(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- InteractionVolume& iv,
- FluxVarsCacheContainer& fluxVarsCacheContainer)
- {
- // lambda function to get the permeability tensor
- auto getK = [&problem](const Element& element,
- const VolumeVariables& volVars,
- const SubControlVolume& scv)
- { return volVars.permeability(); };
-
- // solve the local system subject to the permeability
- iv.solveLocalSystem(getK);
-
- // update the transmissibilities etc for each phase
- const auto scvfIdx = scvf.index();
- const auto scvfLocalFaceData = iv.getLocalFaceData(scvf);
- auto& scvfCache = fluxVarsCacheContainer[scvfIdx];
- scvfCache.updateAdvection(iv, scvf, scvfLocalFaceData);
-
- //! Update the transmissibilities in the other scvfs of the interaction volume
- for (auto scvfIdxJ : iv.globalScvfs())
- {
- if (scvfIdxJ != scvfIdx)
- {
- // update cache of scvfJ
- const auto& scvfJ = fvGeometry.scvf(scvfIdxJ);
- const auto scvfJLocalFaceData = iv.getLocalFaceData(scvfJ);
- fluxVarsCacheContainer[scvfIdxJ].updateAdvection(iv, scvfJ, scvfJLocalFaceData);
- }
- }
- }
-};
-
-//! Fills the diffusive quantities in the caches
-template<class TypeTag>
-class CCMpfaDiffusionCacheFiller
-{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-
- using Element = typename GridView::template Codim<0>::Entity;
-
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
- static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
-
-public:
- //! function to fill the flux var caches
- template<class FluxVarsCacheContainer, class InteractionVolume>
- static void fillCaches(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- InteractionVolume& iv,
- FluxVarsCacheContainer& fluxVarsCacheContainer)
- {
- //! update the cache for all components in all phases. We exclude the case
- //! phaseIdx = compIdx here, as diffusive fluxes of the major component in its
- //! own phase are not calculated explicitly during assembly (see compositional local residual)
- for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- if (phaseIdx == compIdx)
- continue;
-
- // lambda function to get the diffusion coefficient/tensor
- auto getD = [phaseIdx, compIdx](const Element& element,
- const VolumeVariables& volVars,
- const SubControlVolume& scv)
- { return volVars.diffusionCoefficient(phaseIdx, compIdx); };
-
- // solve for the transmissibilities
- iv.solveLocalSystem(getD);
-
- // update the caches
- const auto scvfIdx = scvf.index();
- const auto scvfLocalFaceData = iv.getLocalFaceData(scvf);
- auto& scvfCache = fluxVarsCacheContainer[scvfIdx];
- scvfCache.updateDiffusion(iv, scvf, scvfLocalFaceData, phaseIdx, compIdx);
-
- //! Update the caches in the other scvfs of the interaction volume
- for (auto scvfIdxJ : iv.globalScvfs())
- {
- if (scvfIdxJ != scvfIdx)
- {
- // store scvf pointer and local face data
- const auto& scvfJ = fvGeometry.scvf(scvfIdxJ);
- const auto scvfJLocalFaceData = iv.getLocalFaceData(scvfJ);
-
- // update cache
- fluxVarsCacheContainer[scvfIdxJ].updateDiffusion(iv, scvfJ, scvfJLocalFaceData, phaseIdx, compIdx);
- }
- }
- }
- }
- }
-};
-
-//! Fills the quantities for heat conduction in the caches
-template<class TypeTag>
-class CCMpfaHeatConductionCacheFiller
-{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
- using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
-
- using Element = typename GridView::template Codim<0>::Entity;
-
- static const int energyEqIdx = Indices::energyEqIdx;
- static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- static const bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
-
-public:
- //! function to fill the flux var caches
- template<class FluxVarsCacheContainer, class InteractionVolume>
- static void fillCaches(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- InteractionVolume& iv,
- FluxVarsCacheContainer& fluxVarsCacheContainer)
- {
- // lambda function to get the thermal conductivity
- auto getLambda = [&problem, &fvGeometry](const Element& element,
- const VolumeVariables& volVars,
- const SubControlVolume& scv)
- { return ThermalConductivityModel::effectiveThermalConductivity(volVars,
- problem.spatialParams(),
- element,
- fvGeometry,
- scv); };
-
- // solve for the transmissibilities
- iv.solveLocalSystem(getLambda);
-
- // update the caches
- const auto scvfIdx = scvf.index();
- const auto scvfLocalFaceData = iv.getLocalFaceData(scvf);
- auto& scvfCache = fluxVarsCacheContainer[scvfIdx];
- scvfCache.updateHeatConduction(iv, scvf, scvfLocalFaceData);
-
- //! Update the caches in the other scvfs of the interaction volume
- for (auto scvfIdxJ : iv.globalScvfs())
- {
- if (scvfIdxJ != scvfIdx)
- {
- // store scvf pointer and local face data
- const auto& scvfJ = fvGeometry.scvf(scvfIdxJ);
- const auto scvfJLocalFaceData = iv.getLocalFaceData(scvfJ);
-
- // update cache
- fluxVarsCacheContainer[scvfIdxJ].updateHeatConduction(iv, scvfJ, scvfJLocalFaceData);
- }
- }
- }
-};
-
-
-
-} // end namespace
-
-#endif
diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
index 03d00e0db4..76f15e411d 100644
--- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
@@ -48,11 +48,12 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
// Always use the dynamic type for vectors (compatibility with the boundary)
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
- using DynamicVector = typename BoundaryInteractionVolume::Vector;
+ using CoefficientVector = typename BoundaryInteractionVolume::Vector;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
@@ -64,10 +65,84 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
static constexpr int energyEqIdx = GET_PROP_TYPE(TypeTag, Indices)::energyEqIdx;
+ //! The cache used in conjunction with the mpfa Fourier's Law
+ class MpfaFouriersLawCache
+ {
+ using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
+ public:
+ // update cached objects for heat conduction
+ template<typename InteractionVolume>
+ void updateHeatConduction(const InteractionVolume& iv, const SubControlVolumeFace &scvf)
+ {
+ const auto& localFaceData = iv.getLocalFaceData(scvf);
+ heatConductionVolVarsStencil_ = iv.volVarsStencil();
+ heatConductionTij_ = iv.getTransmissibilities(localFaceData);
+ heatNeumannFlux_ = iv.getNeumannFlux(localFaceData, energyEqIdx);
+
+ if (enableInteriorBoundaries)
+ heatConductionCij_ = iv.getNeumannFluxTransformationCoefficients(localFaceData);
+ }
+
+ //! Returns the volume variables indices necessary for heat conduction flux
+ //! computation. This includes all participating boundary volume variables
+ //! and it can be different for the phases & components.
+ const Stencil& heatConductionVolVarsStencil() const
+ { return heatConductionVolVarsStencil_; }
+
+ //! Returns the transmissibilities associated with the volume variables
+ //! This can be different for the phases & components.
+ const CoefficientVector& heatConductionTij() const
+ { return heatConductionTij_; }
+
+ //! Returns the vector of coefficients with which the vector of neumann boundary conditions
+ //! has to be multiplied in order to transform them on the scvf this cache belongs to
+ const CoefficientVector& heatConductionCij() const
+ { return heatConductionCij_; }
+
+ //! If the useTpfaBoundary property is set to false, the boundary conditions
+ //! are put into the local systems leading to possible contributions on all faces
+ Scalar heatNeumannFlux() const
+ { return heatNeumannFlux_; }
+
+ private:
+ // Quantities associated with heat conduction
+ Stencil heatConductionVolVarsStencil_;
+ CoefficientVector heatConductionTij_;
+ CoefficientVector heatConductionCij_;
+ Scalar heatNeumannFlux_;
+ };
+
+ //! Class that fills the cache corresponding to mpfa Darcy's Law
+ class MpfaFouriersLawCacheFiller
+ {
+ public:
+ //! Function to fill an MpfaDarcysLawCache of a given scvf
+ //! This interface has to be met by any cache filler class for heat conduction quantities
+ template<class FluxVariablesCacheFiller>
+ static void fill(FluxVariablesCache& scvfFluxVarsCache,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const FluxVariablesCacheFiller& fluxVarsCacheFiller)
+ {
+ // get interaction volume from the flux vars cache filler & upate the cache
+ if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf))
+ scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.boundaryInteractionVolume(), scvf);
+ else
+ scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.interactionVolume(), scvf);
+ }
+ };
+
public:
// state the discretization method this implementation belongs to
static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;
+ // state the type for the corresponding cache and its filler
+ using Cache = MpfaFouriersLawCache;
+ using CacheFiller = MpfaFouriersLawCacheFiller;
+
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -110,7 +185,7 @@ public:
const auto& cij = fluxVarsCache.heatConductionCij();
// The Vector of interior neumann fluxes
- DynamicVector facetCouplingFluxes(cij.size(), 0.0);
+ CoefficientVector facetCouplingFluxes(cij.size(), 0.0);
for (auto&& data : fluxVarsCache.interiorBoundaryData())
{
// Add additional Dirichlet fluxes for interior Dirichlet faces
diff --git a/dumux/discretization/cellcentered/mpfa/globalfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/globalfluxvariablescache.hh
index 6bafe80a78..79ace006a3 100644
--- a/dumux/discretization/cellcentered/mpfa/globalfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/mpfa/globalfluxvariablescache.hh
@@ -24,8 +24,7 @@
#define DUMUX_DISCRETIZATION_CCMPFA_GLOBAL_FLUXVARSCACHE_HH
#include <dumux/implicit/properties.hh>
-
-#include "fluxvariablescachefiller.hh"
+#include <dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh>
namespace Dumux
{
@@ -45,8 +44,11 @@ class CCMpfaGlobalFluxVariablesCache;
template<class TypeTag>
class CCMpfaGlobalFluxVariablesCache<TypeTag, true>
{
- // the local class needs access to the problem
+ // the local class need access to the problem
friend typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ // the filler needs access to the operators
+ friend CCMpfaFluxVariablesCacheFiller<TypeTag>;
+
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
@@ -66,6 +68,10 @@ public:
void update(Problem& problem)
{
problemPtr_ = &problem;
+
+ // instantiate helper class to fill the caches
+ FluxVariablesCacheFiller filler(problem);
+
const auto& globalFvGeometry = problem.model().globalFvGeometry();
fluxVarsCache_.resize(globalFvGeometry.numScvf());
for (const auto& element : elements(problem.gridView()))
@@ -81,7 +87,7 @@ public:
for (auto&& scvf : scvfs(fvGeometry))
{
if (!fluxVarsCache_[scvf.index()].isUpdated())
- FluxVariablesCacheFiller::fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCache_);
+ filler.fill(*this, fluxVarsCache_[scvf.index()], element, fvGeometry, elemVolVars, scvf);
}
}
}
@@ -97,6 +103,12 @@ public:
private:
// access operators in the case of caching
+ const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const
+ { return fluxVarsCache_[scvf.index()]; }
+
+ FluxVariablesCache& operator [](const SubControlVolumeFace& scvf)
+ { return fluxVarsCache_[scvf.index()]; }
+
const FluxVariablesCache& operator [](IndexType scvfIdx) const
{ return fluxVarsCache_[scvfIdx]; }
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
index af8c2d793f..32af446ad2 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
@@ -79,6 +79,9 @@ public:
LocalIndexType localScvIndex;
bool isOutside;
+ //! default constructor
+ LocalFaceData() = default;
+
//! Constructor fully initializing the members
LocalFaceData(const LocalIndexType faceIndex,
const LocalIndexType scvIndex,
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 9450d4da51..7d3bb189bc 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -135,6 +135,9 @@ public:
using typename ParentType::PositionVector;
using typename ParentType::Seed;
+ // structure to store global and local face data
+ using GlobalLocalFaceDataPair = std::pair<const SubControlVolumeFace*, LocalFaceData>;
+
CCMpfaOInteractionVolume(const Seed& seed,
const Problem& problem,
const FVElementGeometry& fvGeometry,
@@ -142,69 +145,10 @@ public:
: problemPtr_(&problem),
fvGeometryPtr_(&fvGeometry),
elemVolVarsPtr_(&elemVolVars),
- onDomainOrInteriorBoundary_(seed.onDomainOrInteriorBoundary()),
- globalScvfIndices_(seed.globalScvfIndices())
+ onDomainOrInteriorBoundary_(seed.onDomainOrInteriorBoundary())
{
- // create local sub control entities from the seed
- createLocalEntities_(seed);
-
- // reserve memory
- auto numLocalScvs = localScvs_.size();
- auto numLocalScvfs = localScvfs_.size();
- auto maxNumVolVars = numLocalScvs + numLocalScvfs;
- volVarsStencil_ = seed.globalScvIndices(); // boundary vol vars are placed at the end
- volVarsStencil_.reserve(maxNumVolVars);
- volVarsPositions_.reserve(maxNumVolVars);
- dirichletFaceIndexSet_.reserve(numLocalScvfs);
- interiorDirichletFaceIndexSet_.reserve(numLocalScvfs);
- interiorBoundaryFaceIndexSet_.reserve(numLocalScvfs);
- interiorBoundaryData_.reserve(numLocalScvfs);
- fluxFaceIndexSet_.reserve(numLocalScvfs);
-
- // the positions where the vol vars are defined at (required for the gravitational acceleration)
- for (auto&& localScv : localScvs_)
- volVarsPositions_.push_back(localScv.center());
-
- // maybe add dirichlet vol var indices and set up local index sets of flux and dirichlet faces
- LocalIndexType localScvfIdx = 0;
- for (auto&& localScvf : localScvfs_)
- {
- auto faceType = localScvf.faceType();
-
- // eventually add vol var index and corresponding position
- if (faceType == MpfaFaceTypes::dirichlet)
- {
- volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex());
- volVarsPositions_.push_back(localScvf.ip());
- dirichletFaceIndexSet_.push_back(localScvfIdx++);
- }
- else if (faceType == MpfaFaceTypes::interior || faceType == MpfaFaceTypes::neumann)
- {
- fluxFaceIndexSet_.push_back(localScvfIdx++);
- }
- else if (faceType == MpfaFaceTypes::interiorNeumann)
- {
- interiorBoundaryData_.push_back(InteriorBoundaryData(problem,
- localScvf.insideGlobalScvIndex(),
- localScvf.insideGlobalScvfIndex(),
- fluxFaceIndexSet_.size(),
- faceType));
- fluxFaceIndexSet_.push_back(localScvfIdx);
- interiorBoundaryFaceIndexSet_.push_back(localScvfIdx++);
- }
- else if(faceType == MpfaFaceTypes::interiorDirichlet)
- {
- interiorBoundaryData_.push_back(InteriorBoundaryData(problem,
- localScvf.insideGlobalScvIndex(),
- localScvf.insideGlobalScvfIndex(),
- interiorDirichletFaceIndexSet_.size(),
- faceType));
- interiorDirichletFaceIndexSet_.push_back(localScvfIdx);
- interiorBoundaryFaceIndexSet_.push_back(localScvfIdx++);
- }
- else
- DUNE_THROW(Dune::InvalidStateException, "Unknown mpfa face type.");
- }
+ // set up the local scope of this interaction volume
+ bind(seed);
// initialize the vector containing the neumann fluxes
assembleNeumannFluxVector_();
@@ -239,123 +183,6 @@ public:
T_ += D;
}
- //! gets local data on a global scvf within the interaction volume
- //! specialization for dim = dimWorld
- template<int d = dim, int dw = dimWorld>
- typename std::enable_if< (d == dw), LocalFaceData >::type
- getLocalFaceData(const SubControlVolumeFace& scvf) const
- {
- auto scvfGlobalIdx = scvf.index();
-
- // if interior boundaries are disabled, do simplified search
- if (!enableInteriorBoundaries)
- {
- LocalIndexType localIdx = 0;
- for (auto&& localScvf : localScvfs_)
- {
- if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
- return LocalFaceData(localIdx, localScvf.insideLocalScvIndex(), false);
-
- if (!localScvf.boundary() && localScvf.outsideGlobalScvfIndex() == scvfGlobalIdx)
- return LocalFaceData(localIdx, localScvf.outsideLocalScvIndex(), true);
-
- localIdx++;
- }
- }
- else
- {
- // first check the inside data of the faces
- LocalIndexType localIdx = 0;
- for (auto&& localScvf : localScvfs_)
- {
- if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
- return LocalFaceData(localIdx, localScvf.insideLocalScvIndex(), false);
- localIdx++;
- }
-
- // then check the outside data of the faces
- localIdx = 0;
- for (auto&& localScvf : localScvfs_)
- {
- if (!localScvf.boundary() && localScvf.outsideGlobalScvfIndex() == scvfGlobalIdx)
- return LocalFaceData(localIdx, localScvf.outsideLocalScvIndex(), true);
- localIdx++;
- }
- }
-
- DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv face in the interaction volume for the given scvf with index: " << scvf.index());
- }
-
- //! gets local data on a global scvf within the interaction volume
- //! specialization for dim < dimWorld
- template<int d = dim, int dw = dimWorld>
- typename std::enable_if< (d < dw), LocalFaceData >::type
- getLocalFaceData(const SubControlVolumeFace& scvf) const
- {
- const auto scvfGlobalIdx = scvf.index();
-
- // if interior boundaries are disabled, do simplified search
- if (!enableInteriorBoundaries)
- {
- LocalIndexType localFaceIdx = 0;
- for (auto&& localScvf : localScvfs_)
- {
- if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
- return LocalFaceData(localFaceIdx, localScvf.insideLocalScvIndex(), false);
-
- if (!localScvf.boundary() && MpfaHelper::contains(localScvf.outsideGlobalScvfIndices(), scvfGlobalIdx))
- {
- if (localScvf.outsideLocalScvIndices().size() == 1)
- return LocalFaceData(localFaceIdx, localScvf.outsideLocalScvIndex(), true);
-
- // Now we have to find wich local scv is the one the global scvf is embedded in
- const auto scvGlobalIdx = scvf.insideScvIdx();
- for (auto localScvIdx : localScvf.outsideLocalScvIndices())
- if (localScv_(localScvIdx).globalIndex() == scvGlobalIdx)
- return LocalFaceData(localFaceIdx, localScvIdx, true);
-
- DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv in the interaction volume for the given scvf with index: " << scvfGlobalIdx);
- }
-
- localFaceIdx++;
- }
- }
- else
- {
- // first check the inside data of the faces
- LocalIndexType localFaceIdx = 0;
- for (auto&& localScvf : localScvfs_)
- {
- if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
- return LocalFaceData(localFaceIdx, localScvf.insideLocalScvIndex(), false);
- localFaceIdx++;
- }
-
- // then check the outside data
- localFaceIdx = 0;
- for (auto&& localScvf : localScvfs_)
- {
- if (!localScvf.boundary() && MpfaHelper::contains(localScvf.outsideGlobalScvfIndices(), scvfGlobalIdx))
- {
- if (localScvf.outsideLocalScvIndices().size() == 1)
- return LocalFaceData(localFaceIdx, localScvf.outsideLocalScvIndex(), true);
-
- // Now we have to find wich local scv is the one the global scvf is embedded in
- const auto scvGlobalIdx = scvf.insideScvIdx();
- for (auto localScvIdx : localScvf.outsideLocalScvIndices())
- if (localScv_(localScvIdx).globalIndex() == scvGlobalIdx)
- return LocalFaceData(localFaceIdx, localScvIdx, true);
-
- DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv in the interaction volume for the given scvf with index: " << scvfGlobalIdx);
- }
-
- localFaceIdx++;
- }
- }
-
- DUNE_THROW(Dune::InvalidStateException, "Could not find the local scv face in the interaction volume for the given scvf with index: " << scvfGlobalIdx);
- }
-
//! Gets the transmissibilities for a sub-control volume face within the interaction volume.
//! specialization for dim == dimWorld
template<int d = dim, int dw = dimWorld>
@@ -455,6 +282,9 @@ public:
return flux;
}
+ const LocalFaceData& getLocalFaceData(const SubControlVolumeFace& scvf) const
+ { return globalLocalScvfPairedData_[this->findIndexInVector(globalScvfIndices_, scvf.index())].second; }
+
bool onDomainOrInteriorBoundary() const
{ return onDomainOrInteriorBoundary_; }
@@ -464,8 +294,8 @@ public:
const PositionVector& volVarsPositions() const
{ return volVarsPositions_; }
- const GlobalIndexSet& globalScvfs() const
- { return globalScvfIndices_; }
+ const std::vector<GlobalLocalFaceDataPair>& globalLocalScvfPairedData() const
+ { return globalLocalScvfPairedData_; }
const std::vector<InteriorBoundaryData>& interiorBoundaryData() const
{ return interiorBoundaryData_; }
@@ -493,27 +323,107 @@ private:
const Element& localElement_(const LocalIndexType localScvIdx) const
{ return localElements_[localScvIdx]; }
- void createLocalEntities_(const Seed& seed)
+ void bind(const Seed& seed)
{
- const auto numScvs = seed.scvSeeds().size();
- const auto numScvfs = seed.scvfSeeds().size();
-
- localElements_.reserve(numScvs);
- localScvs_.reserve(numScvs);
- localScvfs_.reserve(numScvfs);
+ const auto numLocalScvs = seed.scvSeeds().size();
+ const auto numLocalScvfs = seed.scvfSeeds().size();
+ const auto numGlobalScvfs = seed.globalScvfIndices().size();
+ const auto maxNumVolVars = numLocalScvs + numLocalScvfs;
+
+ //! reserve memory for local entities
+ localElements_.reserve(numLocalScvs);
+ localScvs_.reserve(numLocalScvs);
+ localScvfs_.reserve(numLocalScvfs);
+ globalLocalScvfPairedData_.reserve(numGlobalScvfs);
+ globalScvfIndices_.reserve(numGlobalScvfs);
+
+ //! reserve memory for the index sets
+ volVarsStencil_ = seed.globalScvIndices(); // boundary vol vars are placed at the end
+ volVarsStencil_.reserve(maxNumVolVars);
+ volVarsPositions_.reserve(maxNumVolVars);
+ dirichletFaceIndexSet_.reserve(numLocalScvfs);
+ interiorDirichletFaceIndexSet_.reserve(numLocalScvfs);
+ interiorBoundaryFaceIndexSet_.reserve(numLocalScvfs);
+ interiorBoundaryData_.reserve(numLocalScvfs);
+ fluxFaceIndexSet_.reserve(numLocalScvfs);
+ // set up quantities related to sub-control volumes
for (auto&& scvSeed : seed.scvSeeds())
{
- auto element = problem_().model().globalFvGeometry().element(scvSeed.globalIndex());
- localScvs_.emplace_back(LocalScvType(problem_(), element, fvGeometry_(), scvSeed));
+ const auto element = problem_().model().globalFvGeometry().element(scvSeed.globalIndex());
+ localScvs_.emplace_back(problem_(), element, fvGeometry_(), scvSeed);
localElements_.emplace_back(std::move(element));
+ volVarsPositions_.push_back(localScvs_.back().center());
}
+ // set up quantitites related to sub-control volume faces
+ LocalIndexType localFaceIdx = 0;
for (auto&& scvfSeed : seed.scvfSeeds())
{
+ const auto faceType = scvfSeed.faceType();
+
// we have to use the "inside" scv face here
- auto&& scvf = fvGeometry_().scvf(scvfSeed.insideGlobalScvfIndex());
- localScvfs_.emplace_back(LocalScvfType(scvfSeed, scvf));
+ const auto& scvf = fvGeometry_().scvf(scvfSeed.insideGlobalScvfIndex());
+ localScvfs_.emplace_back(scvfSeed, scvf);
+
+ // create global/local face data for this face
+ // we simultaneously store the corresponding global scvf indices (allows global to local mapping later)
+ globalLocalScvfPairedData_.emplace_back(&scvf, LocalFaceData(localFaceIdx, scvfSeed.insideLocalScvIndex(), false));
+ globalScvfIndices_.push_back(scvf.index());
+
+ // set data depending on the face type
+ // obtain the local scvf entity just inserted
+ const auto& localScvf = localScvfs_.back();
+
+ // interior faces are flux faces
+ // also, set up outside global/local data for all the neighbors
+ if (faceType == MpfaFaceTypes::interior)
+ {
+ for (unsigned int i = 0; i < localScvf.outsideGlobalScvfIndices().size(); ++i)
+ {
+ const auto& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndex(i));
+ globalLocalScvfPairedData_.emplace_back( &outsideScvf,
+ LocalFaceData(localFaceIdx, scvfSeed.outsideLocalScvIndex(i), true) );
+ globalScvfIndices_.push_back(outsideScvf.index());
+ }
+ fluxFaceIndexSet_.push_back(localFaceIdx++);
+ }
+ // dirichlet faces are in the "stencil"
+ else if (faceType == MpfaFaceTypes::dirichlet)
+ {
+ volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex());
+ volVarsPositions_.push_back(localScvf.ip());
+ dirichletFaceIndexSet_.push_back(localFaceIdx++);
+ }
+ // neumann faces have an unknown associated with them
+ else if (faceType == MpfaFaceTypes::neumann)
+ {
+ fluxFaceIndexSet_.push_back(localFaceIdx++);
+ }
+ // interior neumann faces additionally produce interior boundary data
+ else if (faceType == MpfaFaceTypes::interiorNeumann)
+ {
+ interiorBoundaryData_.push_back(InteriorBoundaryData(problem_(),
+ localScvf.insideGlobalScvIndex(),
+ localScvf.insideGlobalScvfIndex(),
+ fluxFaceIndexSet_.size(),
+ faceType));
+ fluxFaceIndexSet_.push_back(localFaceIdx);
+ interiorBoundaryFaceIndexSet_.push_back(localFaceIdx++);
+ }
+ // as well as interior Dirichlet boundaries
+ else if (faceType == MpfaFaceTypes::interiorDirichlet)
+ {
+ interiorBoundaryData_.push_back(InteriorBoundaryData(problem_(),
+ localScvf.insideGlobalScvIndex(),
+ localScvf.insideGlobalScvfIndex(),
+ interiorDirichletFaceIndexSet_.size(),
+ faceType));
+ interiorDirichletFaceIndexSet_.push_back(localFaceIdx);
+ interiorBoundaryFaceIndexSet_.push_back(localFaceIdx++);
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Face type can not be handled by the mpfa o-method.");
}
}
@@ -867,11 +777,13 @@ private:
bool onDomainOrInteriorBoundary_;
+ // Variables defining the local scope
std::vector<Element> localElements_;
std::vector<LocalScvType> localScvs_;
std::vector<LocalScvfType> localScvfs_;
-
+ std::vector<GlobalLocalFaceDataPair> globalLocalScvfPairedData_;
GlobalIndexSet globalScvfIndices_;
+
GlobalIndexSet volVarsStencil_;
PositionVector volVarsPositions_;
@@ -881,12 +793,13 @@ private:
LocalIndexSet interiorBoundaryFaceIndexSet_;
std::vector<InteriorBoundaryData> interiorBoundaryData_;
+ // Quantities depending on the tensor the system is solved for
std::vector< std::vector< DimVector > > wijk_;
-
DynamicMatrix T_;
DynamicMatrix AinvB_;
DynamicMatrix CAinv_;
+ // stores all the neumann fluxes appearing in this interaction volume
std::vector<PrimaryVariables> neumannFluxes_;
};
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
index db5884a5ec..730017018c 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
@@ -165,7 +165,7 @@ public:
// for grids with dim < dimWorld, some outside indices might be doubled
// we want to make the outside indices unique, but, the i-th outside global scvf face
- // should correspond to the j-th outside local scv.Therefore we apply the same operations on both containers
+ // should correspond to the i-th outside local scv.Therefore we apply the same operations on both containers
void makeOutsideDataUnique()
{
auto scvIt = outsideLocalScvIndices_.begin();
diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
index 1e01be146b..181aa3240c 100644
--- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
@@ -162,240 +162,20 @@ private:
// forward declaration of the base class of the mpfa flux variables cache
template<class TypeTag, bool EnableAdvection, bool EnableMolecularDiffusion, bool EnableEnergyBalance>
-class MpfaPorousMediumFluxVariablesCache;
+class CCMpfaPorousMediumFluxVariablesCache;
-//! Base class for the advective cache in mpfa methods
-template<class TypeTag>
-class MpfaAdvectionCache
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
-
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
-
- static constexpr bool facetCoupling = GET_PROP_VALUE(TypeTag, MpfaFacetCoupling);
- static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
-
- // We always use the dynamic types here to be compatible on the boundary
- using IndexType = typename GridView::IndexSet::IndexType;
- using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
- using CoefficientVector = typename BoundaryInteractionVolume::Vector;
- using PositionVector = typename BoundaryInteractionVolume::PositionVector;
-
-public:
- //! update cached objects
- template<class InteractionVolume, class LocalFaceData>
- void updateAdvection(const InteractionVolume& interactionVolume,
- const SubControlVolumeFace &scvf,
- const LocalFaceData& scvfLocalFaceData)
- {
- // update the quantities that are equal for all phases
- advectionVolVarsStencil_ = interactionVolume.volVarsStencil();
- advectionVolVarsPositions_ = interactionVolume.volVarsPositions();
- advectionTij_ = interactionVolume.getTransmissibilities(scvfLocalFaceData);
-
- // we will need the neumann flux transformation only on interior boundaries with facet coupling
- if (enableInteriorBoundaries && facetCoupling)
- advectionCij_ = interactionVolume.getNeumannFluxTransformationCoefficients(scvfLocalFaceData);
-
- // The neumann fluxes always have to be set per phase
- for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- phaseNeumannFluxes_[phaseIdx] = interactionVolume.getNeumannFlux(scvfLocalFaceData, phaseIdx);
- }
-
- //! Returns the volume variables indices necessary for flux computation
- //! This includes all participating boundary volume variables. Since we
- //! do not allow mixed BC for the mpfa this is the same for all phases.
- const Stencil& advectionVolVarsStencil() const
- { return advectionVolVarsStencil_; }
-
- //! Returns the position on which the volume variables live. This is
- //! necessary as we need to evaluate gravity also for the boundary volvars
- const PositionVector& advectionVolVarsPositions() const
- { return advectionVolVarsPositions_; }
-
- //! Returns the transmissibilities associated with the volume variables
- //! All phases flow through the same rock, thus, tij are equal for all phases
- const CoefficientVector& advectionTij() const
- { return advectionTij_; }
-
- //! Returns the vector of coefficients with which the vector of neumann boundary conditions
- //! has to be multiplied in order to transform them on the scvf this cache belongs to
- const CoefficientVector& advectionCij() const
- { return advectionCij_; }
-
- //! If the useTpfaBoundary property is set to false, the boundary conditions
- //! are put into the local systems leading to possible contributions on all faces
- Scalar advectionNeumannFlux(unsigned int phaseIdx) const
- { return phaseNeumannFluxes_[phaseIdx]; }
-
-private:
- // Quantities associated with advection
- Stencil advectionVolVarsStencil_;
- PositionVector advectionVolVarsPositions_;
- CoefficientVector advectionTij_;
- CoefficientVector advectionCij_;
- std::array<Scalar, numPhases> phaseNeumannFluxes_;
-};
-
-//! Base class for the diffusive cache in mpfa methods
-template<class TypeTag>
-class MpfaDiffusionCache
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
-
- static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
- static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
-
- // We always use the dynamic types here to be compatible on the boundary
- using IndexType = typename GridView::IndexSet::IndexType;
- using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
- using CoefficientVector = typename BoundaryInteractionVolume::Vector;
- using PositionVector = typename BoundaryInteractionVolume::PositionVector;
-
-public:
- //! The constructor. Initializes the Neumann flux to zero
- MpfaDiffusionCache() { componentNeumannFluxes_.fill(0.0); }
-
- // update cached objects for the diffusive fluxes
- template<typename InteractionVolume, class LocalFaceData>
- void updateDiffusion(const InteractionVolume& interactionVolume,
- const SubControlVolumeFace &scvf,
- const LocalFaceData& scvfLocalFaceData,
- unsigned int phaseIdx,
- unsigned int compIdx)
- {
- diffusionVolVarsStencils_[phaseIdx][compIdx] = interactionVolume.volVarsStencil();
- diffusionTij_[phaseIdx][compIdx] = interactionVolume.getTransmissibilities(scvfLocalFaceData);
-
- if (enableInteriorBoundaries)
- diffusionCij_[phaseIdx][compIdx] = interactionVolume.getNeumannFluxTransformationCoefficients(scvfLocalFaceData);
-
- //! For compositional models, we associate neumann fluxes with the phases (main components)
- //! This is done in the AdvectionCache. However, in single-phasic models we solve the phase AND
- //! the component mass balance equations. Thus, in this case we have diffusive neumann contributions.
- //! we assume compIdx = eqIdx
- if (numPhases == 1 && phaseIdx != compIdx)
- componentNeumannFluxes_[compIdx] = interactionVolume.getNeumannFlux(scvfLocalFaceData, compIdx);
-
- }
-
- //! Returns the volume variables indices necessary for diffusive flux
- //! computation. This includes all participating boundary volume variables
- //! and it can be different for the phases & components.
- const Stencil& diffusionVolVarsStencil(unsigned int phaseIdx,
- unsigned int compIdx) const
- { return diffusionVolVarsStencils_[phaseIdx][compIdx]; }
-
- //! Returns the transmissibilities associated with the volume variables
- //! This can be different for the phases & components.
- const CoefficientVector& diffusionTij(unsigned int phaseIdx,
- unsigned int compIdx) const
- { return diffusionTij_[phaseIdx][compIdx]; }
-
- //! Returns the vector of coefficients with which the vector of neumann boundary conditions
- //! has to be multiplied in order to transform them on the scvf this cache belongs to
- const CoefficientVector& diffusionCij(unsigned int phaseIdx,
- unsigned int compIdx) const
- { return diffusionCij_[phaseIdx][compIdx]; }
-
- //! If the useTpfaBoundary property is set to false, the boundary conditions
- //! are put into the local systems leading to possible contributions on all faces
- Scalar componentNeumannFlux(unsigned int compIdx) const
- {
- assert(numPhases == 1);
- return componentNeumannFluxes_[compIdx];
- }
-
-private:
- // Quantities associated with molecular diffusion
- std::array< std::array<Stencil, numComponents>, numPhases> diffusionVolVarsStencils_;
- std::array< std::array<CoefficientVector, numComponents>, numPhases> diffusionTij_;
- std::array< std::array<CoefficientVector, numComponents>, numPhases> diffusionCij_;
-
- // diffusive neumann flux for single-phasic models
- std::array<Scalar, numComponents> componentNeumannFluxes_;
-};
-
-//! Base class for the heat conduction cache in mpfa methods
-template<class TypeTag>
-class MpfaHeatConductionCache
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
-
- // We always use the dynamic types here to be compatible on the boundary
- using IndexType = typename GridView::IndexSet::IndexType;
- using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
- using CoefficientVector = typename BoundaryInteractionVolume::Vector;
-
- static constexpr int energyEqIdx = GET_PROP_TYPE(TypeTag, Indices)::energyEqIdx;
- static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
-public:
- // update cached objects for heat conduction
- template<typename InteractionVolume, class LocalFaceData>
- void updateHeatConduction(const InteractionVolume& interactionVolume,
- const SubControlVolumeFace &scvf,
- const LocalFaceData& scvfLocalFaceData)
- {
- heatConductionVolVarsStencil_ = interactionVolume.volVarsStencil();
- heatConductionTij_ = interactionVolume.getTransmissibilities(scvfLocalFaceData);
- heatNeumannFlux_ = interactionVolume.getNeumannFlux(scvfLocalFaceData, energyEqIdx);
-
- if (enableInteriorBoundaries)
- heatConductionCij_ = interactionVolume.getNeumannFluxTransformationCoefficients(scvfLocalFaceData);
- }
-
- //! Returns the volume variables indices necessary for heat conduction flux
- //! computation. This includes all participating boundary volume variables
- //! and it can be different for the phases & components.
- const Stencil& heatConductionVolVarsStencil() const
- { return heatConductionVolVarsStencil_; }
-
- //! Returns the transmissibilities associated with the volume variables
- //! This can be different for the phases & components.
- const CoefficientVector& heatConductionTij() const
- { return heatConductionTij_; }
-
- //! Returns the vector of coefficients with which the vector of neumann boundary conditions
- //! has to be multiplied in order to transform them on the scvf this cache belongs to
- const CoefficientVector& heatConductionCij() const
- { return heatConductionCij_; }
-
- //! If the useTpfaBoundary property is set to false, the boundary conditions
- //! are put into the local systems leading to possible contributions on all faces
- Scalar heatNeumannFlux() const
- { return heatNeumannFlux_; }
-
-private:
- // Quantities associated with heat conduction
- Stencil heatConductionVolVarsStencil_;
- CoefficientVector heatConductionTij_;
- CoefficientVector heatConductionCij_;
- Scalar heatNeumannFlux_;
-};
-
-// specialization of the flux variables cache for cell centered mpfa methods
+// specialization of the flux variables cache for the cell centered finite volume mpfa scheme
template<class TypeTag>
class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::CCMpfa>
- : public MpfaPorousMediumFluxVariablesCache<TypeTag,
- GET_PROP_VALUE(TypeTag, EnableAdvection),
- GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
- GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>
+ : public CCMpfaPorousMediumFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableAdvection),
+ GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
+ GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>
{
using InteriorBoundaryData = typename GET_PROP_TYPE(TypeTag, InteriorBoundaryData);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using ParentType = MpfaPorousMediumFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableAdvection),
- GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
- GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
+ using ParentType = CCMpfaPorousMediumFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableAdvection),
+ GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
+ GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
@@ -421,12 +201,12 @@ public:
//! maybe update data on interior Dirichlet boundaries
template<class InteractionVolume>
- void updateInteriorBoundaryData(const InteractionVolume& interactionVolume,
+ void updateInteriorBoundaryData(const InteractionVolume& iv,
const SubControlVolumeFace &scvf)
{
if (enableInteriorBoundaries)
{
- interiorBoundaryData_ = interactionVolume.interiorBoundaryData();
+ interiorBoundaryData_ = iv.interiorBoundaryData();
// check if the actual scvf is on an interior Dirichlet boundary
const auto scvfIdx = scvf.index();
@@ -470,23 +250,25 @@ private:
// specialization for the case of pure advection
template<class TypeTag>
-class MpfaPorousMediumFluxVariablesCache<TypeTag, true, false, false> : public MpfaAdvectionCache<TypeTag> {};
+class CCMpfaPorousMediumFluxVariablesCache<TypeTag, true, false, false> : public GET_PROP_TYPE(TypeTag, AdvectionType)::Cache {};
// specialization for the case of advection & diffusion
template<class TypeTag>
-class MpfaPorousMediumFluxVariablesCache<TypeTag, true, true, false> : public MpfaAdvectionCache<TypeTag>,
- public MpfaDiffusionCache<TypeTag> {};
+class CCMpfaPorousMediumFluxVariablesCache<TypeTag, true, true, false> : public GET_PROP_TYPE(TypeTag, AdvectionType)::Cache,
+ public GET_PROP_TYPE(TypeTag, MolecularDiffusionType)::Cache {};
// specialization for the case of advection & heat conduction
template<class TypeTag>
-class MpfaPorousMediumFluxVariablesCache<TypeTag, true, false, true> : public MpfaAdvectionCache<TypeTag>,
- public MpfaHeatConductionCache<TypeTag> {};
+class CCMpfaPorousMediumFluxVariablesCache<TypeTag, true, false, true> : public GET_PROP_TYPE(TypeTag, AdvectionType)::Cache,
+ public GET_PROP_TYPE(TypeTag, HeatConductionType)::Cache {};
// specialization for the case of advection, diffusion & heat conduction
template<class TypeTag>
-class MpfaPorousMediumFluxVariablesCache<TypeTag, true, true, true> : public MpfaAdvectionCache<TypeTag>,
- public MpfaDiffusionCache<TypeTag>,
- public MpfaHeatConductionCache<TypeTag> {};
+class CCMpfaPorousMediumFluxVariablesCache<TypeTag, true, true, true> : public GET_PROP_TYPE(TypeTag, AdvectionType)::Cache,
+ public GET_PROP_TYPE(TypeTag, MolecularDiffusionType)::Cache,
+ public GET_PROP_TYPE(TypeTag, HeatConductionType)::Cache {};
+
+// TODO further specializations
} // end namespace
--
GitLab