Skip to content

GitLab

Commit b5bf2f02 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[mpfa] enable interior boundary handling

parent e453c663
Expand all Hide whitespace changes
Inline Side-by-side
dumux/discretization/cellcentered/mpfa/darcyslaw.hh
View file @ b5bf2f02
......@@ -31,8 +31,8 @@
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
namespace Dumux
{
......@@ -50,21 +50,29 @@ NEW_PROP_TAG(ProblemEnableGravity);
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
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 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>;
static const bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
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);
public:
// state the discretization method this implementation belongs to
......@@ -81,12 +89,43 @@ public:
const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
const auto& volVarsStencil = fluxVarsCache.advectionVolVarsStencil(phaseIdx);
const auto& volVarsPositions = fluxVarsCache.advectionVolVarsPositions(phaseIdx);
const auto& tij = fluxVarsCache.advectionTij(phaseIdx);
// interface density as arithmetic mean of the neighbors (when gravity is on)
Scalar rho = gravity ? interpolateDensity(elemVolVars, scvf, phaseIdx) : 0.0;
const auto& volVarsStencil = fluxVarsCache.advectionVolVarsStencil();
const auto& volVarsPositions = fluxVarsCache.advectionVolVarsPositions();
const auto& tij = fluxVarsCache.advectionTij();
const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
// For interior Neumann boundaries when using Tpfa for Neumann boundary conditions, we simply
// return the user-specified flux. We assume phaseIdx = eqIdx here.
if (isInteriorBoundary
&& useTpfaBoundary
&& fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann)
return scvf.area()*
elemVolVars[scvf.insideScvIdx()].extrusionFactor()*
problem.neumann(element,
fvGeometry,
elemVolVars,
scvf)[phaseIdx];
// Calculate the interface density for gravity evaluation
const auto rho = [&] ()
{
if (!gravity)
return Scalar(0.0);
else
{
// Treat interior boundaries differently
if (enableInteriorBoundaries && isInteriorBoundary)
{
const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
if (facetCoupling || data.faceType() == MpfaFaceTypes::interiorDirichlet)
return data.facetVolVars(fvGeometry).density(phaseIdx);
else
return interpolateDensity(elemVolVars, scvf, phaseIdx);
}
else
return interpolateDensity(elemVolVars, scvf, phaseIdx);
}
} ();
// calculate Tij*pj
Scalar flux(0.0);
......@@ -105,13 +144,75 @@ public:
h -= rho*(g*x);
}
flux += tij[localIdx++]*h;
}
if (useTpfaBoundary)
return flux;
else
return flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
// if no interior boundaries are present, return the flux
if (!enableInteriorBoundaries)
return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
//////////////////////////////////////////////////////////////////
// Handle interior boundaries
//////////////////////////////////////////////////////////////////
// For active facet coupling we will have to transform the interior flux vector
const auto& cij = fluxVarsCache.advectionCij();
// The vector of interior neumann fluxes
DynamicVector facetCouplingFluxes(cij.size(), 0.0);
for (auto&& data : fluxVarsCache.interiorBoundaryData())
{
// Add additional Dirichlet fluxes for interior Dirichlet faces
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
{
Scalar h = data.facetVolVars(fvGeometry).pressure(phaseIdx);
if (gravity)
{
const auto x = fvGeometry.scvf(data.scvfIndex()).ipGlobal();
const auto g = problem.gravityAtPos(x);
h -= rho*(g*x);
}
// The transmissibilities of interior dirichlet boundaries are placed at the end
// So we simply keep incrementing the local index
flux += tij[localIdx + data.localIndexInInteractionVolume()]*h;
}
// add neumann fluxes for interior Neumann faces if facet coupling is active
if (facetCoupling && data.faceType() == MpfaFaceTypes::interiorNeumann)
{
// get the volvars of the actual interior neumann face
const auto facetVolVars = data.facetVolVars(fvGeometry);
// get the scvf corresponding to actual interior neumann face
const auto& curScvf = fvGeometry.scvf(data.scvfIndex());
// calculate "lekage factor"
const auto n = curScvf.unitOuterNormal();
const auto v = [&] ()
{
auto res = n;
res *= -0.5*facetVolVars.extrusionFactor();
res -= curScvf.ipGlobal();
res += curScvf.facetCorner();
res /= res.two_norm2();
return res;
} ();
// add value to vector of interior neumann fluxes
facetCouplingFluxes[data.localIndexInInteractionVolume()] -= facetVolVars.pressure(phaseIdx)*
MpfaHelper::nT_M_v(n, facetVolVars.permeability(), v);
}
}
// return overall resulting flux
const Scalar interiorNeumannFlux = facetCoupling ? cij*facetCouplingFluxes : 0.0;
return useTpfaBoundary ?
flux + interiorNeumannFlux :
flux + interiorNeumannFlux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
}
static Stencil stencil(const Problem& problem,
......@@ -122,7 +223,7 @@ public:
const auto& globalFvGeometry = problem.model().globalFvGeometry();
// return the scv (element) indices in the interaction region
if (globalFvGeometry.scvfTouchesBoundary(scvf))
if (globalFvGeometry.touchesInteriorOrDomainBoundary(scvf))
return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices();
else
return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices();
......
dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
View file @ b5bf2f02
......@@ -227,7 +227,7 @@ private:
int getLocalScvfIdx_(const int scvfIdx) const
{
auto it = std::find(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx);
assert(globalScvfIndices_[std::distance(globalScvfIndices_.begin(), it)] == scvfIdx && "Could not find the flux vars cache for scvfIdx");
assert(it != globalScvfIndices_.end() && "Could not find the flux vars cache for scvfIdx");
return std::distance(globalScvfIndices_.begin(), it);
}
......
dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
View file @ b5bf2f02
......@@ -194,8 +194,8 @@ public:
// Update boundary volume variables in the neighbors
for (auto&& scvf : scvfs(fvGeometry))
{
// skip the rest if the scvf does not touch a boundary
if (!globalFvGeometry.scvfTouchesBoundary(scvf))
// skip the rest if the scvf does not touch a domain boundary
if (!globalFvGeometry.touchesDomainBoundary(scvf))
continue;
// loop over all the scvfs in the interaction region
......@@ -304,7 +304,7 @@ private:
for (auto&& scvf : scvfs(fvGeometry))
{
bool boundary = scvf.boundary();
if (boundary || (!boundary && fvGeometry.globalFvGeometry().scvfTouchesBoundary(scvf)))
if (boundary || (!boundary && fvGeometry.globalFvGeometry().touchesDomainBoundary(scvf)))
bVolVarEstimate += dim-1;
}
......
dumux/discretization/cellcentered/mpfa/facetypes.hh
View file @ b5bf2f02
......@@ -29,7 +29,9 @@ namespace Dumux
{
interior,
neumann,
dirichlet
dirichlet,
interiorNeumann,
interiorDirichlet
};
} // end namespace
......
dumux/discretization/cellcentered/mpfa/fickslaw.hh
View file @ b5bf2f02
......@@ -46,19 +46,28 @@ class FicksLawImplementation<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 MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
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 InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
// 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>;
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
static constexpr bool facetCoupling = GET_PROP_VALUE(TypeTag, MpfaFacetCoupling);
static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
public:
// state the discretization method this implementation belongs to
static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;
......@@ -76,8 +85,55 @@ public:
const auto& volVarsStencil = fluxVarsCache.diffusionVolVarsStencil(phaseIdx, compIdx);
const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx);
const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
// For interior Neumann boundaries when using Tpfa for Neumann boundary conditions, we simply
// return the user-specified flux. Note that for compositional models we attribute the influxes
// to the major components, thus we do it per phase in Darcy's law. However, for single-phasic models
// wesolve the phase mass balance equation AND the transport equation, thus, in that case we incorporate
// the Neumann BCs here. We assume compIdx = eqIdx.
// Note that this way of including interior Neumann fluxes fails for mpnc models where n != m.
if (numPhases == 1
&& isInteriorBoundary
&& useTpfaBoundary
&& fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann)
return scvf.area()*
elemVolVars[scvf.insideScvIdx()].extrusionFactor()*
problem.neumann(element,
fvGeometry,
elemVolVars,
scvf)[compIdx];
// get the scaling factor for the effective diffusive fluxes
auto effFactor = calculateEffectiveDiffusivityFactor(elemVolVars, scvf, phaseIdx);
const auto effFactor = [&] ()
{
// Treat interior boundaries differently
if (isInteriorBoundary)
{
const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
// interpolate as usual for interior Neumann faces without facet coupling
if (data.faceType() == MpfaFaceTypes::interiorNeumann && !facetCoupling)
return calculateEffectiveDiffusivityFactor(elemVolVars, scvf, phaseIdx);
// use harmonic mean between the interior and the facet volvars
else
{
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(),
insideVolVars.saturation(phaseIdx),
/*Diffusion coefficient*/ 1.0);
const auto facetVolVars = data.facetVolVars(fvGeometry);
const auto outsideFactor = EffDiffModel::effectiveDiffusivity(facetVolVars.porosity(),
facetVolVars.saturation(phaseIdx),
/*Diffusion coefficient*/ 1.0);
return harmonicMean(factor, outsideFactor);
}
}
else
return calculateEffectiveDiffusivityFactor(elemVolVars, scvf, phaseIdx);
} ();
// if factor is zero, the flux will end up zero anyway
if (effFactor == 0.0)
......@@ -90,14 +146,88 @@ public:
auto getRho = [useMoles, phaseIdx] (const VolumeVariables& volVars)
{ return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
// calculate the density at the interface
const auto rho = [&] ()
{
// maybe use the density of the interior BC on the facet
if (isInteriorBoundary)
{
const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
if (facetCoupling || data.faceType() == MpfaFaceTypes::interiorDirichlet)
return useMoles ? data.facetVolVars(fvGeometry).molarDensity(phaseIdx) : data.facetVolVars(fvGeometry).density(phaseIdx);
else
return interpolateDensity(elemVolVars, scvf, getRho);
}
else
return interpolateDensity(elemVolVars, scvf, getRho);
} ();
// calculate Tij*xj
Scalar flux(0.0);
unsigned int localIdx = 0;
for (const auto volVarIdx : volVarsStencil)
flux += tij[localIdx++]*getX(elemVolVars[volVarIdx]);
// return effective mass flux
return flux*interpolateDensity(elemVolVars, scvf, getRho)*effFactor;
// if no interior boundaries are present, return effective mass flux
if (!enableInteriorBoundaries)
return useTpfaBoundary ? flux*rho*effFactor : flux*rho*effFactor + fluxVarsCache.componentNeumannFlux(compIdx);
//////////////////////////////////////////////////////////////////
// Handle interior boundaries
//////////////////////////////////////////////////////////////////
// get coefficients to transform the vector of interior neumann boundary conditions
const auto& cij = fluxVarsCache.diffusionCij(phaseIdx, compIdx);
// The vector of interior neumann fluxes
DynamicVector facetCouplingFluxes(cij.size(), 0.0);
for (auto&& data : fluxVarsCache.interiorBoundaryData())
{
// Add additional Dirichlet fluxes for interior Dirichlet faces
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
{
// The transmissibilities of interior dirichlet boundaries are placed at the end
// So we simply keep incrementing the local index
const auto x = useMoles ?
data.facetVolVars(fvGeometry).moleFraction(phaseIdx, compIdx) :
data.facetVolVars(fvGeometry).massFraction(phaseIdx, compIdx);
flux += tij[localIdx + data.localIndexInInteractionVolume()]*x;
}
// add neumann fluxes for interior Neumann faces
if (facetCoupling && data.faceType() == MpfaFaceTypes::interiorNeumann)
{
// get the scvf corresponding to actual interior boundary face
const auto& curScvf = fvGeometry.scvf(data.scvfIndex());
// get the volvars of the actual interior neumann face
const auto facetVolVars = data.facetVolVars(fvGeometry);
// calculate "lekage factor"
const auto n = curScvf.unitOuterNormal();
const auto v = [&] ()
{
auto res = n;
res *= -0.5*facetVolVars.extrusionFactor();
res -= curScvf.ipGlobal();
res += curScvf.facetCorner();
res /= res.two_norm2();
return res;
} ();
// add value to vector of interior neumann fluxes
facetCouplingFluxes[data.localIndexInInteractionVolume()] += MpfaHelper::nT_M_v(n,
facetVolVars.diffusionCoefficient(phaseIdx, compIdx),
v);
}
}
// return overall resulting flux
const Scalar interiorNeumannFlux = facetCoupling ? cij*facetCouplingFluxes : 0.0;
return useTpfaBoundary ?
flux + interiorNeumannFlux :
flux + interiorNeumannFlux + fluxVarsCache.componentNeumannFlux(compIdx);
}
static Stencil stencil(const Problem& problem,
......@@ -108,7 +238,7 @@ public:
const auto& globalFvGeometry = problem.model().globalFvGeometry();
// return the scv (element) indices in the interaction region
if (globalFvGeometry.scvfTouchesBoundary(scvf))
if (globalFvGeometry.touchesInteriorOrDomainBoundary(scvf))
return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices();
else
return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices();
......@@ -141,19 +271,25 @@ private:
const unsigned int phaseIdx)
{
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(),
const auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(),
insideVolVars.saturation(phaseIdx),
/*Diffusion coefficient*/ 1.0);
if (!scvf.boundary())
{
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
auto outsideFactor = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(),
// interpret outside factor as arithmetic mean
Scalar outsideFactor = 0.0;
for (auto outsideIdx : scvf.outsideScvIndices())
{
const auto& outsideVolVars = elemVolVars[outsideIdx];
outsideFactor += EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(),
outsideVolVars.saturation(phaseIdx),
/*Diffusion coefficient*/ 1.0);
}
outsideFactor /= scvf.outsideScvIndices().size();
// use the harmonic mean of the two
factor = harmonicMean(factor, outsideFactor);
return harmonicMean(factor, outsideFactor);
}
return factor;
......
dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh
View file @ b5bf2f02
......@@ -47,8 +47,6 @@ class CCMpfaAdvectionCacheFiller
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 InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
......@@ -56,18 +54,19 @@ class CCMpfaAdvectionCacheFiller
using Element = typename GridView::template Codim<0>::Entity;
static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
static const bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
static const bool enableDiffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion);
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>
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
......@@ -76,68 +75,24 @@ public:
const SubControlVolume& scv)
{ return volVars.permeability(); };
// update flux var caches
if (problem.model().globalFvGeometry().scvfTouchesBoundary(scvf))
{
const auto& seed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf);
BoundaryInteractionVolume iv(seed, problem, fvGeometry, elemVolVars);
iv.solveLocalSystem(getK);
// update the transmissibilities etc for each phase
const auto scvfIdx = scvf.index();
auto& cache = fluxVarsCacheContainer[scvfIdx];
cache.updateAdvection(scvf, iv);
cache.setUpdateStatus(true);
// solve the local system subject to the permeability
iv.solveLocalSystem(getK);
//! set the neumann boundary conditions in case we do not use tpfa on the
//! boundary and diffusion is not enabled (then we assume neumann BCs to be diffusive)
if (!useTpfaBoundary && !enableDiffusion)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
iv.assembleNeumannFluxes(phaseIdx);
cache.updatePhaseNeumannFlux(scvf, iv, phaseIdx);
}
// 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);
for (const auto scvfIdxJ : iv.globalScvfs())
{
if (scvfIdxJ != scvfIdx)
{
const auto& scvfJ = fvGeometry.scvf(scvfIdxJ);
auto& cacheJ = fluxVarsCacheContainer[scvfIdxJ];
cacheJ.updateAdvection(scvfJ, iv);
cacheJ.setUpdateStatus(true);
if (!useTpfaBoundary && !enableDiffusion)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
iv.assembleNeumannFluxes(phaseIdx);
cacheJ.updatePhaseNeumannFlux(scvfJ, iv, phaseIdx);
}
}
}
}
else
//! Update the transmissibilities in the other scvfs of the interaction volume
for (auto scvfIdxJ : iv.globalScvfs())
{
const auto& seed = problem.model().globalFvGeometry().interactionVolumeSeed(scvf);
InteractionVolume iv(seed, problem, fvGeometry, elemVolVars);
iv.solveLocalSystem(getK);
// update flux variables cache
const auto scvfIdx = scvf.index();
auto& cache = fluxVarsCacheContainer[scvfIdx];
cache.updateAdvection(scvf, iv);
cache.setUpdateStatus(true);
// update flux variable caches of the other scvfs of the interaction volume
for (const auto scvfIdxJ : iv.globalScvfs())
if (scvfIdxJ != scvfIdx)
{
if (scvfIdxJ != scvfIdx)
{
const auto& scvfJ = fvGeometry.scvf(scvfIdxJ);
auto& cacheJ = fluxVarsCacheContainer[scvfIdxJ];
cacheJ.updateAdvection(scvfJ, iv);
cacheJ.setUpdateStatus(true);
}
// update cache of scvfJ
const auto& scvfJ = fvGeometry.scvf(scvfIdxJ);
const auto scvfJLocalFaceData = iv.getLocalFaceData(scvfJ);
fluxVarsCacheContainer[scvfIdxJ].updateAdvection(iv, scvfJ, scvfJLocalFaceData);
}
}
}
......@@ -152,8 +107,6 @@ class CCMpfaDiffusionCacheFiller
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 InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
......@@ -161,107 +114,57 @@ class CCMpfaDiffusionCacheFiller
using Element = typename GridView::template Codim<0>::Entity;
static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
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);