diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 3aedae63d47b7268288c4c3d0439e20775228a23..c26b0bca72038035425e3a99e633f7b2ada6abb6 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 2bae60663e2c3bcc89f6e56db3790a9e3334032d..954afca4c268a26fe066b2466598c43a25fba82c 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 6aaf11e1881cc9b723bd65c3394dbc17843047ca..34cdec1077d514c4caa29d6ca5eafe8409c9ac3f 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 0ca84820d817e6e4201ff94ffbe873ad8649abf2..dfb3ad9b6f8bb79fa3ce5b220c8af48df4d4716e 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -20,616 +20,409 @@
* \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"
+#include <dumux/discretization/cellcentered/mpfa/tensorlambdafactory.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.
+ const auto numOtherScvfs = iv.globalLocalScvfPairedData().size()-1;
+ std::vector<FluxVariablesCache*> otherFluxVarCaches(numOtherScvfs);
+ std::vector<Element> otherElements(numOtherScvfs);
+
+ 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);
- }
+ 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++;
}
- 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);
- }
- }
-
- //! 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>;
+ //! 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 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);
+ static constexpr auto AdvectionMethod = AdvectionType::myDiscretizationMethod;
+ using LambdaFactory = TensorLambdaFactory<TypeTag, AdvectionMethod>;
- using Element = typename GridView::template Codim<0>::Entity;
+ // maybe solve the local system subject to K (if AdvectionType uses mpfa)
+ if (AdvectionMethod == DiscretizationMethods::CCMpfa)
+ iv.solveLocalSystem(LambdaFactory::getAdvectionLambda());
- static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
- static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+ // fill the caches of all scvfs within this interaction volume
+ AdvectionFiller::fill(scvfFluxVarsCache, problem(), element(), fvGeometry(), elemVolVars(), scvFace(), *this);
- static constexpr bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection);
- static constexpr bool advectionUsesMpfa = AdvectionType::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++;
+ }
+ }
- static constexpr bool solDependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction);
- static constexpr bool heatConductionUsesMpfa = HeatConductionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
+ //! 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;
- //! 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;
+ static constexpr auto DiffusionMethod = DiffusionType::myDiscretizationMethod;
+ using LambdaFactory = TensorLambdaFactory<TypeTag, DiffusionMethod>;
- using BoolPair = std::pair<bool, bool>;
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-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 (DiffusionMethod == DiscretizationMethods::CCMpfa)
+ iv.solveLocalSystem(LambdaFactory::getDiffusionLambda(phaseIdx, compIdx));
+
+ // 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);
+ using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
+ using HeatConductionFiller = typename HeatConductionType::CacheFiller;
- // maybe update Advection or diffusion or both
- const bool updateAdvection = doAdvectionUpdate && touchesDomainBoundary;
- const bool updateHeatConduction = doHeatConductionUpdate && touchesDomainBoundary;
+ static constexpr auto HeatConductionMethod = HeatConductionType::myDiscretizationMethod;
+ using LambdaFactory = TensorLambdaFactory<TypeTag, HeatConductionMethod>;
- 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);
- }
-};
+ // maybe solve the local system subject to fourier coefficient
+ if (HeatConductionMethod == DiscretizationMethods::CCMpfa)
+ iv.solveLocalSystem(LambdaFactory::getHeatConductionLambda());
-//! 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>;
+ // fill the caches of all scvfs within this interaction volume
+ HeatConductionFiller::fill(scvfFluxVarsCache, problem(), element(), fvGeometry(), elemVolVars(), scvFace(), *this);
- 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))
+ unsigned int otherScvfIdx = 0;
+ const auto curScvfIdx = scvFace().index();
+ 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 (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);
- }
- }
- else
- {
- 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 9ef0820f3a0142428ef6c9d16513b4ca0eb487f9..0000000000000000000000000000000000000000
--- 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 03d00e0db4d1a1370d0714c6edcade186a4a4b08..76f15e411da58c6bc0ee8e70e988e18d2c41653a 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 6bafe80a7834b2cc947d94662bce728e23aa6a5e..79ace006a3b5a534db1cd93b9de4f5d0884aee03 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 af8c2d793f151997880065581d45dc71462047d7..32af446ad2a96ed62aeadc991d5cd6a2222e5ba8 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/interiorboundarydata.hh b/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh
index c1cc25b6b088e56bd0b6191f481f54537418ace9..3244b8c74e3b2bf81d07001e2003fdaf3fe8b427 100644
--- a/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh
+++ b/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh
@@ -55,6 +55,11 @@ class InteriorBoundaryData
//! Note that the return types are also "wrong" (Here just to satisfy the compiler)
struct CompleteCoupledFacetData
{
+ const Problem& problem() const
+ {
+ DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided");
+ }
+
const SpatialParams& spatialParams() const
{
DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided");
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 9450d4da519ced549991478677a9e486d2277195..444fc6a823bf5866652f14127c7957da6f390c29 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.");
}
}
@@ -575,7 +485,7 @@ private:
const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
const auto& posVolVars = elemVolVars_()[posGlobalScv];
const auto& element = localElement_(posLocalScvIdx);
- const auto tensor = getTensor(element, posVolVars, posGlobalScv);
+ const auto tensor = getTensor(problem_(), element, posVolVars, fvGeometry_(), posGlobalScv);
// the omega factors of the "positive" sub volume
auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
@@ -613,23 +523,28 @@ private:
// get interior boundary data
const auto& data = interiorBoundaryData_[this->findIndexInVector(interiorBoundaryScvfIndexSet_(), curLocalScvfIdx)];
- // get the volvars on the actual interior boundary face
- const auto facetVolVars = data.facetVolVars(fvGeometry_());
+ // obtain the complete data on the facet element
+ const auto completeFacetData = data.completeCoupledFacetData();
// calculate "leakage factor"
const auto n = curLocalScvf.unitOuterNormal();
const auto v = [&] ()
{
auto res = n;
- res *= -0.5*facetVolVars.extrusionFactor();
+ res *= -0.5*completeFacetData.volVars().extrusionFactor();
res -= curLocalScvf.ip();
res += curLocalScvf.globalScvf().facetCorner();
res /= res.two_norm2();
return res;
} ();
- // substract from matrix entry
- A[idxInFluxFaces][curIdxInFluxFaces] -= MpfaHelper::nT_M_v(n, posVolVars.permeability(), v);
+ // substract n*T*v from diagonal matrix entry
+ const auto facetTensor = getTensor(completeFacetData.problem(),
+ completeFacetData.element(),
+ completeFacetData.volVars(),
+ completeFacetData.fvGeometry(),
+ completeFacetData.scv());
+ A[idxInFluxFaces][curIdxInFluxFaces] -= MpfaHelper::nT_M_v(n, facetTensor, v);
}
}
// this means we are on an interior face
@@ -683,7 +598,7 @@ private:
const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
const auto& negVolVars = elemVolVars_()[negGlobalScv];
const auto& negElement = localElement_(negLocalScvIdx);
- const auto negTensor = getTensor(negElement, negVolVars, negGlobalScv);
+ const auto negTensor = getTensor(problem_(), negElement, negVolVars, fvGeometry_(), negGlobalScv);
// the omega factors of the "negative" sub volume
DimVector negWijk;
@@ -771,7 +686,7 @@ private:
const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
const auto& posVolVars = elemVolVars_()[posGlobalScv];
const auto element = localElement_(posLocalScvIdx);
- const auto tensor = getTensor(element, posVolVars, posGlobalScv);
+ const auto tensor = getTensor(problem_(), element, posVolVars, fvGeometry_(), posGlobalScv);
// the omega factors of the "positive" sub volume
auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
@@ -867,11 +782,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 +798,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 db5884a5ec3595fcd713d55d2cd725ab5af228c8..730017018c7a344f78bc0332192a18bba40707a0 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/discretization/cellcentered/mpfa/tensorlambdafactory.hh b/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..548cd6598e072cd06ba515c80b4d60f6bff6d95a
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh
@@ -0,0 +1,149 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Helper class to be used to obtain lambda functions for the tensors
+ * involved in the laws that describe the different kind of fluxes that
+ * occur in DuMuX models (i.e. advective, diffusive and heat conduction fluxes).
+ * The local systems appearing in Mpfa methods have to be solved subject to
+ * the different tensors. This class returns lambda expressions to be used in the
+ * local systems. The specialization for other discretization methods allows
+ * compatibility with the TPFA scheme, which could be used for one or more of the tensors.
+ */
+#ifndef DUMUX_DISCRETIZATION_MPFA_TENSOR_LAMBDA_FACTORY_HH
+#define DUMUX_DISCRETIZATION_MPFA_TENSOR_LAMBDA_FACTORY_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh>
+
+namespace Dumux
+{
+
+//! forward declaration of properties
+namespace Properties
+{
+NEW_PROP_TAG(ThermalConductivityModel);
+};
+
+/*!
+ * \ingroup MpfaModel
+ * \brief Helper class to be used to obtain lambda functions for the tensors
+ * involved in the laws that describe the different kind of fluxes that
+ * occur in DuMuX models (i.e. advective, diffusive and heat conduction fluxes).
+ * The local systems appearing in Mpfa methods have to be solved subject to
+ * the different tensors. This class returns lambda expressions to be used in the
+ * local systems. The specialization for discretization methods other than mpfa allows
+ * compatibility with the TPFA scheme, which could be used for one or more of the tensors.
+ * The interfaces of the lambdas are chosen such that all involved tensors can be extracted
+ * with the given arguments.
+ */
+template<class TypeTag, DiscretizationMethods Method>
+class TensorLambdaFactory
+{
+public:
+
+ //! We return zero scalars here in the functions below.
+ //! We have to return something as the local systems expect a type
+ //! to perform actions on. We return 0.0 as this call should never happen
+ //! for a tensor which is not treated by an mpfa method anyway.
+
+ //! lambda for the law describing the advective term
+ static auto getAdvectionLambda()
+ {
+ return [] (const auto& problem,
+ const auto& element,
+ const auto& volVars,
+ const auto& fvGeometry,
+ const auto& scv)
+ { return 0.0; };
+ }
+
+ //! lambda for the diffusion coefficient/tensor
+ static auto getDiffusionLambda(unsigned int phaseIdx, unsigned int compIdx)
+ {
+ return [] (const auto& problem,
+ const auto& element,
+ const auto& volVars,
+ const auto& fvGeometry,
+ const auto& scv)
+ { return 0.0; };
+ }
+
+ //! lambda for the fourier coefficient
+ static auto getHeatConductionLambda()
+ {
+ return [] (const auto& problem,
+ const auto& element,
+ const auto& volVars,
+ const auto& fvGeometry,
+ const auto& scv)
+ { return 0.0; };
+ }
+};
+
+//! Specialization for mpfa schemes
+template<class TypeTag>
+class TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>
+{
+public:
+
+ //! return void lambda here
+ static auto getAdvectionLambda()
+ {
+ return [] (const auto& problem,
+ const auto& element,
+ const auto& volVars,
+ const auto& fvGeometry,
+ const auto& scv)
+ { return volVars.permeability(); };
+ }
+
+ //! return void lambda here
+ static auto getDiffusionLambda(unsigned int phaseIdx, unsigned int compIdx)
+ {
+ return [phaseIdx, compIdx] (const auto& problem,
+ const auto& element,
+ const auto& volVars,
+ const auto& fvGeometry,
+ const auto& scv)
+ { return volVars.diffusionCoefficient(phaseIdx, compIdx); };
+ }
+
+ //! return void lambda here
+ static auto getHeatConductionLambda()
+ {
+ using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
+
+ return [] (const auto& problem,
+ const auto& element,
+ const auto& volVars,
+ const auto& fvGeometry,
+ const auto& scv)
+ { return ThermalConductivityModel::effectiveThermalConductivity(volVars,
+ problem.spatialParams(),
+ element,
+ fvGeometry,
+ scv); };
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
index 3a3143739464163573a757f40d09403c5a765160..73e086080c224ac7cd550dbd4d3d52a979f89c90 100644
--- a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
@@ -51,6 +51,7 @@ NEW_PROP_TAG(ProblemEnableGravity);
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCTpfa>
{
+ using Implementation = DarcysLawImplementation<TypeTag, DiscretizationMethods::CCTpfa>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
@@ -60,6 +61,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCTpfa>
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
@@ -70,10 +72,52 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCTpfa>
using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ class TpfaDarcysLawCache
+ {
+ public:
+ void updateAdvection(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf)
+ {
+ tij_ = Implementation::calculateTransmissibilities(problem, element, fvGeometry, elemVolVars, scvf);
+ }
+
+ const Scalar& tij() const
+ { return tij_; }
+
+ private:
+ Scalar tij_;
+ };
+
+ //! Class that fills the cache corresponding to tpfa Darcy's Law
+ class TpfaDarcysLawCacheFiller
+ {
+ public:
+ //! Function to fill a TpfaDarcysLawCache 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)
+ {
+ scvfFluxVarsCache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf);
+ }
+ };
+
public:
// state the discretization method this implementation belongs to
static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCTpfa;
+ // state the type for the corresponding cache and its filler
+ using Cache = TpfaDarcysLawCache;
+ using CacheFiller = TpfaDarcysLawCacheFiller;
+
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
diff --git a/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh
index 343e7d13e7f24c10c755a36d30cfde0813cd7b41..f2360b74848afe0ac31467cdbe70723fdc4c8933 100644
--- a/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/tpfa/elementfluxvariablescache.hh
@@ -24,6 +24,7 @@
#define DUMUX_DISCRETIZATION_CCTPFA_ELEMENT_FLUXVARSCACHE_HH
#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/fluxvariablescachefiller.hh>
namespace Dumux
{
@@ -111,8 +112,7 @@ class CCTpfaElementFluxVariablesCache<TypeTag, false>
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using GlobalFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GlobalFluxVariablesCache);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-
- static const bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection);
+ using FluxVariablesCacheFiller = CCTpfaFluxVariablesCacheFiller<TypeTag>;
public:
CCTpfaElementFluxVariablesCache(const GlobalFluxVariablesCache& global)
@@ -128,12 +128,15 @@ public:
const auto numScvf = fvGeometry.numScvf();
fluxVarsCache_.resize(numScvf);
globalScvfIndices_.resize(numScvf);
- IndexType localScvfIdx = 0;
+ // instantiate helper class to fill the caches
+ FluxVariablesCacheFiller filler(globalFluxVarsCache().problem_());
+
+ IndexType localScvfIdx = 0;
// fill the containers
for (auto&& scvf : scvfs(fvGeometry))
{
- fluxVarsCache_[localScvfIdx].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf);
+ filler.fill(*this, fluxVarsCache_[localScvfIdx], element, fvGeometry, elemVolVars, scvf);
globalScvfIndices_[localScvfIdx] = scvf.index();
localScvfIdx++;
}
@@ -150,6 +153,9 @@ public:
const auto& assemblyMapI = problem.model().localJacobian().assemblyMap()[globalI];
const auto numNeighbors = assemblyMapI.size();
+ // instantiate helper class to fill the caches
+ FluxVariablesCacheFiller filler(problem);
+
// find the number of scv faces that need to be prepared
auto numScvf = fvGeometry.numScvf();
for (unsigned int localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
@@ -161,7 +167,7 @@ public:
unsigned int localScvfIdx = 0;
for (auto&& scvf : scvfs(fvGeometry))
{
- fluxVarsCache_[localScvfIdx].update(problem, element, fvGeometry, elemVolVars, scvf);
+ filler.fill(*this, fluxVarsCache_[localScvfIdx], element, fvGeometry, elemVolVars, scvf);
globalScvfIndices_[localScvfIdx] = scvf.index();
localScvfIdx++;
}
@@ -173,7 +179,7 @@ public:
for (auto scvfIdx : assemblyMapI[localIdxJ].scvfsJ)
{
auto&& scvfJ = fvGeometry.scvf(scvfIdx);
- fluxVarsCache_[localScvfIdx].update(problem, elementJ, fvGeometry, elemVolVars, scvfJ);
+ filler.fill(*this, fluxVarsCache_[localScvfIdx], elementJ, fvGeometry, elemVolVars, scvfJ);
globalScvfIndices_[localScvfIdx] = scvfJ.index();
localScvfIdx++;
}
@@ -188,7 +194,10 @@ public:
fluxVarsCache_.resize(1);
globalScvfIndices_.resize(1);
- fluxVarsCache_[0].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf);
+ // instantiate helper class to fill the caches
+ FluxVariablesCacheFiller filler(globalFluxVarsCache().problem_());
+
+ filler.fill(*this, fluxVarsCache_[0], element, fvGeometry, elemVolVars, scvf);
globalScvfIndices_[0] = scvf.index();
}
@@ -211,8 +220,29 @@ private:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
- if (solDependentAdvection)
- bind(element, fvGeometry, elemVolVars);
+ static const bool isSolIndependent = FluxVariablesCacheFiller::isSolutionIndependent();
+
+ if (!isSolIndependent)
+ {
+ const auto& problem = globalFluxVarsCache().problem_();
+ const auto globalI = problem.elementMapper().index(element);
+
+ // instantiate filler class
+ FluxVariablesCacheFiller filler(problem);
+
+ // let the filler class update the caches
+ for (unsigned int localScvfIdx = 0; localScvfIdx < fluxVarsCache_.size(); ++localScvfIdx)
+ {
+ const auto& scvf = fvGeometry.scvf(globalScvfIndices_[localScvfIdx]);
+
+ const auto scvfInsideScvIdx = scvf.insideScvIdx();
+ const auto& insideElement = scvfInsideScvIdx == globalI ?
+ element :
+ problem.model().globalFvGeometry().element(scvfInsideScvIdx);
+
+ filler.fill(*this, fluxVarsCache_[localScvfIdx], insideElement, fvGeometry, elemVolVars, scvf);
+ }
+ }
}
// get index of scvf in the local container
diff --git a/dumux/discretization/cellcentered/tpfa/fickslaw.hh b/dumux/discretization/cellcentered/tpfa/fickslaw.hh
index 7b05340bfad569bf72afc7fdff95d0d794c53d9d..a6edbd6a9c807f09a9958ee2c1e2d4d9577b5f9f 100644
--- a/dumux/discretization/cellcentered/tpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/fickslaw.hh
@@ -31,6 +31,7 @@
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/fluxvariablescaching.hh>
namespace Dumux
{
@@ -61,6 +62,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa >
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using Element = typename GridView::template Codim<0>::Entity;
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
@@ -73,6 +75,11 @@ public:
// state the discretization method this implementation belongs to
static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCTpfa;
+ //! state the type for the corresponding cache and its filler
+ //! We don't cache anything for this law
+ using Cache = FluxVariablesCaching::EmptyDiffusionCache;
+ using CacheFiller = FluxVariablesCaching::EmptyCacheFiller<TypeTag>;
+
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
diff --git a/dumux/discretization/cellcentered/tpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/tpfa/fluxvariablescachefiller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0fbddcc66f8bb22ed45e28e04e01173fd7cc6022
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/fluxvariablescachefiller.hh
@@ -0,0 +1,227 @@
+// -*- 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 The flux variables cache filler class for the cell-centered TPFA scheme
+ */
+#ifndef DUMUX_DISCRETIZATION_CCTPFA_FLUXVARSCACHE_FILLER_HH
+#define DUMUX_DISCRETIZATION_CCTPFA_FLUXVARSCACHE_FILLER_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Helper class to fill the flux var caches
+ */
+template<class TypeTag>
+class CCTpfaFluxVariablesCacheFiller
+{
+ 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 ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ 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:
+ //! The constructor. Sets the problem pointer
+ CCTpfaFluxVariablesCacheFiller(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}))
+ {
+ // fill the physics-related quantities of the caches
+ if (doSubCaches[ProcessIndices::advectionIdx])
+ fillAdvection(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
+ if (doSubCaches[ProcessIndices::diffusionIdx])
+ fillDiffusion(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
+ if (doSubCaches[ProcessIndices::heatConductionIdx])
+ fillHeatConduction(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
+ }
+
+ /*!
+ * \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)
+ {
+ // 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> updateSubCaches = []()
+ {
+ std::array<bool, 3> tmp;
+ tmp[ProcessIndices::advectionIdx] = doAdvection && soldependentAdvection;
+ tmp[ProcessIndices::diffusionIdx] = doDiffusion && soldependentDiffusion;
+ tmp[ProcessIndices::heatConductionIdx] = doHeatConduction && soldependentHeatConduction;
+ return tmp;
+ } ();
+
+ // forward to fill routine
+ fill(fluxVarsCacheContainer, scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf, updateSubCaches);
+ }
+
+ static bool isSolutionIndependent()
+ {
+ static const bool isSolDependent = (doAdvection && soldependentAdvection) ||
+ (doDiffusion && soldependentDiffusion) ||
+ (doHeatConduction && soldependentHeatConduction);
+ return !isSolDependent;
+ }
+
+private:
+
+ const Problem& problem() const
+ { return *problemPtr_; }
+
+ //! method to fill the advective quantities
+ template<bool advectionEnabled = doAdvection>
+ typename std::enable_if<advectionEnabled>::type
+ fillAdvection(FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+ using AdvectionFiller = typename AdvectionType::CacheFiller;
+
+ // forward to the filler for the advective quantities
+ AdvectionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this);
+ }
+
+ //! do nothing if advection is not enabled
+ template<bool advectionEnabled = doAdvection>
+ typename std::enable_if<!advectionEnabled>::type
+ fillAdvection(FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {}
+
+ //! method to fill the diffusive quantities
+ template<bool diffusionEnabled = doDiffusion>
+ typename std::enable_if<diffusionEnabled>::type
+ fillDiffusion(FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ using DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+ using DiffusionFiller = typename DiffusionType::CacheFiller;
+
+ static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+ // forward to the filler of the diffusive quantities
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
+ if (phaseIdx != compIdx)
+ DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this);
+ }
+
+ //! do nothing if diffusion is not enabled
+ template<bool diffusionEnabled = doDiffusion>
+ typename std::enable_if<!diffusionEnabled>::type
+ fillDiffusion(FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {}
+
+ //! method to fill the quantities related to heat conduction
+ template<bool heatConductionEnabled = doHeatConduction>
+ typename std::enable_if<heatConductionEnabled>::type
+ fillHeatConduction(FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {
+ using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
+ using HeatConductionFiller = typename HeatConductionType::CacheFiller;
+
+ // forward to the filler of the diffusive quantities
+ HeatConductionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this);
+ }
+
+ //! do nothing if heat conduction is disabled
+ template<bool heatConductionEnabled = doHeatConduction>
+ typename std::enable_if<!heatConductionEnabled>::type
+ fillHeatConduction(FluxVariablesCache& scvfFluxVarsCache,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
+ {}
+
+ const Problem* problemPtr_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/tpfa/fourierslaw.hh b/dumux/discretization/cellcentered/tpfa/fourierslaw.hh
index c46876547ae2ff4f9a4f60fb44861d615272fa46..016627964c2f17bf3d938afb2ac33ff42afa0bd7 100644
--- a/dumux/discretization/cellcentered/tpfa/fourierslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/fourierslaw.hh
@@ -30,6 +30,7 @@
#include <dumux/common/parameters.hh>
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/fluxvariablescaching.hh>
namespace Dumux
{
@@ -57,6 +58,7 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCTpfa>
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using Element = typename GridView::template Codim<0>::Entity;
using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
@@ -70,6 +72,11 @@ public:
// state the discretization method this implementation belongs to
static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCTpfa;
+ //! state the type for the corresponding cache and its filler
+ //! We don't cache anything for this law
+ using Cache = FluxVariablesCaching::EmptyHeatConductionCache;
+ using CacheFiller = FluxVariablesCaching::EmptyCacheFiller<TypeTag>;
+
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
diff --git a/dumux/discretization/cellcentered/tpfa/globalfluxvariablescache.hh b/dumux/discretization/cellcentered/tpfa/globalfluxvariablescache.hh
index 2718ffcf18bcf34b2b695aac548cb3c0e045cc8b..aba839d8eb379c9ce3568829792c3f6f324dc681 100644
--- a/dumux/discretization/cellcentered/tpfa/globalfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/tpfa/globalfluxvariablescache.hh
@@ -24,6 +24,7 @@
#define DUMUX_DISCRETIZATION_CCTPFA_GLOBAL_FLUXVARSCACHE_HH
#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/fluxvariablescachefiller.hh>
namespace Dumux
{
@@ -44,18 +45,26 @@ class CCTpfaGlobalFluxVariablesCache<TypeTag, true>
{
// the local class needs access to the problem
friend typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ // the filler class needs access to the access operators
+ friend CCTpfaFluxVariablesCacheFiller<TypeTag>;
+
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FluxVariablesCacheFiller = CCTpfaFluxVariablesCacheFiller<TypeTag>;
public:
// When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
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()))
@@ -69,7 +78,7 @@ public:
for (auto&& scvf : scvfs(fvGeometry))
{
- fluxVarsCache_[scvf.index()].update(problem, element, fvGeometry, elemVolVars, scvf);
+ filler.fill(*this, fluxVarsCache_[scvf.index()], element, fvGeometry, elemVolVars, scvf);
}
}
}
@@ -83,7 +92,13 @@ public:
{ return ElementFluxVariablesCache(global); }
private:
- // access operators in the case of caching
+ // // access operators in the case of caching
+ // const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const
+ // { return fluxVarsCache_[scvf.index()]; }
+
+ // FluxVariablesCache& operator [](const SubControlVolumeFace& scvf)
+ // { return fluxVarsCache_[scvf.index()]; }
+
const FluxVariablesCache& operator [](IndexType scvfIdx) const
{ return fluxVarsCache_[scvfIdx]; }
diff --git a/dumux/discretization/fluxvariablescaching.hh b/dumux/discretization/fluxvariablescaching.hh
new file mode 100644
index 0000000000000000000000000000000000000000..07e762fac061c7da115c774ba81e0f47bd200cd7
--- /dev/null
+++ b/dumux/discretization/fluxvariablescaching.hh
@@ -0,0 +1,83 @@
+// -*- 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 Classes related to flux variables caching
+ */
+#ifndef DUMUX_DISCRETIZATION_FLUXVAR_CACHING_HH
+#define DUMUX_DISCRETIZATION_FLUXVAR_CACHING_HH
+
+namespace Dumux
+{
+
+namespace FluxVariablesCaching
+{
+
+class _EmptyCache {};
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Empty caches to use in a law/process, e.g. Darcy's law
+ * \note Never use the _EmptyCache directly as it lead to ambiguous definitions
+ */
+class EmptyAdvectionCache : public _EmptyCache {};
+class EmptyDiffusionCache : public _EmptyCache {};
+class EmptyHeatConductionCache : public _EmptyCache {};
+
+//! The empty filler class corresponding to EmptyCache
+template<class TypeTag>
+class EmptyCacheFiller
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+public:
+ //! For advection filler
+ 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)
+ {}
+
+ //! For diffusion filler
+ 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)
+ {}
+};
+
+} // end namespace FluxVariablesCaching
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
index 7c4cfd92399e59e459fcf093253dd5ac42a5c551..5b4282d8b89740783ab9d549368b816683b06eb4 100644
--- a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
@@ -130,12 +130,7 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache)
{
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(this->problem(),
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- elemFluxVarsCache);
+ fluxVars.init(this->problem(), element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// get upwind weights into local scope
PrimaryVariables flux(0.0);
diff --git a/dumux/porousmediumflow/compositional/localresidual.hh b/dumux/porousmediumflow/compositional/localresidual.hh
index cc178e0cedb04f601e9ac7c7e743869fa8cc2ba9..01ff8dd50d0f9d3c9b4d9f5ee7abe717dd49320f 100644
--- a/dumux/porousmediumflow/compositional/localresidual.hh
+++ b/dumux/porousmediumflow/compositional/localresidual.hh
@@ -165,12 +165,7 @@ public:
bool useMoles = true)
{
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(this->problem(),
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- elemFluxVarsCache);
+ fluxVars.init(this->problem(), element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// get upwind weights into local scope
PrimaryVariables flux(0.0);
diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh
index d99d19b2987bd3fbbb54560b64a6fc45e82e8731..a7562d94db0f0d993803193b3e9536c5e8f7d90e 100644
--- a/dumux/porousmediumflow/immiscible/localresidual.hh
+++ b/dumux/porousmediumflow/immiscible/localresidual.hh
@@ -103,12 +103,7 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(this->problem(),
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- elemFluxVarsCache);
+ fluxVars.init(this->problem(), element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
PrimaryVariables flux;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
diff --git a/dumux/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh
index 87bd43686a0dadaaba111f73468990d28a286eb4..4e4f54a83b2c02112d17f5d06af8daa23c89d6ee 100644
--- a/dumux/porousmediumflow/implicit/fluxvariables.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariables.hh
@@ -77,16 +77,6 @@ class PorousMediumFluxVariablesImpl<TypeTag, true, false, false> : public FluxVa
public:
- void initAndComputeFluxes(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvFace,
- const ElementFluxVariablesCache& elemFluxVarsCache)
- {
- ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, elemFluxVarsCache);
- }
-
template<typename FunctionType>
Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindTerm)
{
@@ -124,18 +114,11 @@ class PorousMediumFluxVariablesImpl<TypeTag, true, true, false> : public FluxVar
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
public:
-
- void initAndComputeFluxes(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvFace,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ //! The constructor
+ PorousMediumFluxVariablesImpl()
{
advFluxCached_.reset();
advPreFlux_.fill(0.0);
-
- ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, elemFluxVarsCache);
}
template<typename FunctionType>
@@ -195,16 +178,11 @@ class PorousMediumFluxVariablesImpl<TypeTag, true, false, true> : public FluxVar
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
public:
-
- void initAndComputeFluxes(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvFace,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ //! The constructor
+ PorousMediumFluxVariablesImpl()
{
advFluxCached_.reset();
- ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, elemFluxVarsCache);
+ advPreFlux_.fill(0.0);
}
template<typename FunctionType>
@@ -264,16 +242,11 @@ class PorousMediumFluxVariablesImpl<TypeTag, true, true, true> : public FluxVari
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
public:
-
- void initAndComputeFluxes(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvFace,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ //! The constructor
+ PorousMediumFluxVariablesImpl()
{
advFluxCached_.reset();
- ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, elemFluxVarsCache);
+ advPreFlux_.fill(0.0);
}
template<typename FunctionType>
diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
index 1e01be146b4c4486c9bb63d930d105fa56e185da..558d8831545b83cc2f068abf3f317498fc0fa840 100644
--- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
@@ -40,15 +40,24 @@ namespace Properties
NEW_PROP_TAG(NumPhases);
}
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//! The cache is dependent on the active physical processes (advection, diffusion, heat conduction)
+//! For each type of process there is a base cache storing the data required to compute the respective fluxes
+//! Specializations of the overall cache are provided for combinations of processes
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
/*!
* \ingroup ImplicitModel
* \brief The flux variables cache classes for porous media.
- * Store flux stencils and data required for flux calculation
+ * Store data required for flux calculation. For each type of physical process (advection, diffusion, heat conduction)
+ * there is a base cache storing the data required to compute the respective fluxes. Specializations of the overall
+ * cache class are provided for different combinations of processes.
*/
template<class TypeTag>
using PorousMediumFluxVariablesCache = PorousMediumFluxVariablesCacheImplementation<TypeTag, GET_PROP_VALUE(TypeTag, DiscretizationMethod)>;
-// specialization for the Box Method
+//! We only store discretization-related quantities for the box method.
+//! Thus, we need no physics-dependent specialization.
template<class TypeTag>
class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::Box>
{
@@ -100,302 +109,61 @@ public:
const JacobianInverseTransposed& jacInvT() const
{ return jacInvT_; }
- // The stencil info is obsolete for the box method.
- // It is here for compatibility with cc methods
- const Stencil& stencil() const
- {
- return stencil_;
- }
-
private:
std::vector<ShapeJacobian> shapeJacobian_;
std::vector<ShapeValue> shapeValues_;
JacobianInverseTransposed jacInvT_;
-
- Stencil stencil_;
};
+// forward declaration of the base class of the tpfa flux variables cache
+template<class TypeTag, bool EnableAdvection, bool EnableMolecularDiffusion, bool EnableEnergyBalance>
+class CCTpfaPorousMediumFluxVariablesCache;
+
// specialization for the cell centered tpfa method
template<class TypeTag>
class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::CCTpfa>
-{
- 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 FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
- using Stencil = std::vector<IndexType>;
+ : public CCTpfaPorousMediumFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableAdvection),
+ GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
+ GET_PROP_VALUE(TypeTag, EnableEnergyBalance)> {};
-public:
- void update(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf)
- {
- stencil_ = FluxVariables::computeStencil(problem, element, fvGeometry, scvf);
- tij_ = AdvectionType::calculateTransmissibilities(problem, element, fvGeometry, elemVolVars, scvf);
- }
-
- const Stencil& stencil() const
- { return stencil_; }
-
- const Scalar& tij() const
- { return tij_; }
-
-private:
- Stencil stencil_;
- Scalar tij_;
-};
-
-///////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//! Classes building up the porous medium flux variables cache for mpfa methods
-//! The cache is dependent on the active physical processes (advection, diffusion, heat conduction)
-//! For each type of process there is a base cache storing the data required to compute the respective fluxes
-//! Specializations of the overall cache are provided for combinations of processes
-//////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-// forward declaration of the base class of the mpfa flux variables cache
-template<class TypeTag, bool EnableAdvection, bool EnableMolecularDiffusion, bool EnableEnergyBalance>
-class MpfaPorousMediumFluxVariablesCache;
-
-//! Base class for the advective cache in mpfa methods
+// specialization for the case of pure advection
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_; }
+class CCTpfaPorousMediumFluxVariablesCache<TypeTag, true, false, false> : public GET_PROP_TYPE(TypeTag, AdvectionType)::Cache {};
- //! 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
+// specialization for the case of advection & diffusion
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);
+class CCTpfaPorousMediumFluxVariablesCache<TypeTag, true, true, false> : public GET_PROP_TYPE(TypeTag, AdvectionType)::Cache,
+ public GET_PROP_TYPE(TypeTag, MolecularDiffusionType)::Cache {};
- 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
+// specialization for the case of advection & heat conduction
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_; }
+class CCTpfaPorousMediumFluxVariablesCache<TypeTag, true, false, true> : public GET_PROP_TYPE(TypeTag, AdvectionType)::Cache,
+ public GET_PROP_TYPE(TypeTag, HeatConductionType)::Cache {};
- //! 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_; }
+// specialization for the case of advection, diffusion & heat conduction
+template<class TypeTag>
+class CCTpfaPorousMediumFluxVariablesCache<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 {};
- //! 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_; }
+// TODO further specializations
-private:
- // Quantities associated with heat conduction
- Stencil heatConductionVolVarsStencil_;
- CoefficientVector heatConductionTij_;
- CoefficientVector heatConductionCij_;
- Scalar heatNeumannFlux_;
-};
+// forward declaration of the base class of the mpfa flux variables cache
+template<class TypeTag, bool EnableAdvection, bool EnableMolecularDiffusion, bool EnableEnergyBalance>
+class CCMpfaPorousMediumFluxVariablesCache;
-// 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 +189,11 @@ public:
//! maybe update data on interior Dirichlet boundaries
template<class InteractionVolume>
- void updateInteriorBoundaryData(const InteractionVolume& interactionVolume,
- const SubControlVolumeFace &scvf)
+ 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 +237,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
diff --git a/dumux/porousmediumflow/implicit/velocityoutput.hh b/dumux/porousmediumflow/implicit/velocityoutput.hh
index 5980213d44122d7061960440da626cd118ec850a..c2ade4f1eea2df3bc20529bbba6e1638df78234f 100644
--- a/dumux/porousmediumflow/implicit/velocityoutput.hh
+++ b/dumux/porousmediumflow/implicit/velocityoutput.hh
@@ -172,12 +172,7 @@ public:
// insantiate the flux variables
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(problem_,
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- elemFluxVarsCache);
+ fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// get the volume flux divided by the area of the
// subcontrolvolume face in the reference element
@@ -244,7 +239,7 @@ public:
if (!scvf.boundary())
{
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
scvfFluxes[scvfIndexInInside[localScvfIdx]] = fluxVars.advectiveFlux(phaseIdx, upwindTerm);
scvfFluxes[scvfIndexInInside[localScvfIdx]] /= problem_.extrusionFactor(element,
fvGeometry.scv(scvf.insideScvIdx()),
@@ -256,7 +251,7 @@ public:
if (bcTypes.hasOnlyDirichlet())
{
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
scvfFluxes[scvfIndexInInside[localScvfIdx]] = fluxVars.advectiveFlux(phaseIdx, upwindTerm);
scvfFluxes[scvfIndexInInside[localScvfIdx]] /= problem_.extrusionFactor(element,
fvGeometry.scv(scvf.insideScvIdx()),
diff --git a/dumux/porousmediumflow/richards/implicit/localresidual.hh b/dumux/porousmediumflow/richards/implicit/localresidual.hh
index 8e6e8a18e4686c41bec102627dc7bb4a01024ff4..5659a08f071c772ba1e286759ec405b2d4288d9b 100644
--- a/dumux/porousmediumflow/richards/implicit/localresidual.hh
+++ b/dumux/porousmediumflow/richards/implicit/localresidual.hh
@@ -102,12 +102,7 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(this->problem(),
- element,
- fvGeometry,
- elemVolVars,
- scvf,
- elemFluxVarsCache);
+ fluxVars.init(this->problem(), element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
PrimaryVariables flux;
// the physical quantities for which we perform upwinding