Skip to content

GitLab

Commit 1684060c authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[mpfa] change to new flux vars cache concept 

The flux variables cache now inherits from caches defined in the different laws.
These also define a respective filler class that is called when updating the flux
variables caches. That enables us to define different caches for different laws.
parent 968d1418
Expand all Hide whitespace changes
Inline Side-by-side
dumux/discretization/cellcentered/mpfa/darcyslaw.hh
View file @ 1684060c
......@@ -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
......
dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
View file @ 1684060c
......@@ -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
......
dumux/discretization/cellcentered/mpfa/fickslaw.hh
View file @ 1684060c
......@@ -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
......
dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh deleted 100644 → 0
View file @ 968d1418
// -*- 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,