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

[mpfa] remove facetCoupling property from standard classes

we do not incorporate the handling of lower dimensional elements on the facets
in the base laws etc anymore. We want to overload functionalities required to
realize the coupling in the mixeddimension-specific classes. The upwind scheme
is still to be modified!
parent 63684f0d
......@@ -50,6 +50,8 @@ NEW_PROP_TAG(ProblemEnableGravity);
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
using Implementation = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
......@@ -69,7 +71,6 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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);
......@@ -91,10 +92,6 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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);
......@@ -116,11 +113,6 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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
......@@ -131,7 +123,6 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
Stencil advectionVolVarsStencil_;
PositionVector advectionVolVarsPositions_;
CoefficientVector advectionTij_;
CoefficientVector advectionCij_;
std::array<Scalar, numPhases> phaseNeumannFluxes_;
};
......@@ -174,7 +165,7 @@ public:
const unsigned int phaseIdx,
const ElementFluxVariablesCache& elemFluxVarsCache)
{
const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
const auto& volVarsStencil = fluxVarsCache.advectionVolVarsStencil();
......@@ -195,25 +186,7 @@ public:
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);
}
} ();
const auto rho = Implementation::interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
// calculate Tij*pj
Scalar flux(0.0);
......@@ -240,18 +213,65 @@ public:
if (!enableInteriorBoundaries)
return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
//////////////////////////////////////////////////////////////////
// Handle interior boundaries
//////////////////////////////////////////////////////////////////
flux += Implementation::computeInteriorBoundaryContribution(problem, fvGeometry, fluxVarsCache, phaseIdx, rho);
// For active facet coupling we will have to transform the interior flux vector
const auto& cij = fluxVarsCache.advectionCij();
// return overall resulting flux
return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
}
// The vector of interior neumann fluxes
CoefficientVector facetCouplingFluxes(cij.size(), 0.0);
static Scalar interpolateDensity(const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const FluxVariablesCache& fluxVarsCache,
const unsigned int phaseIdx,
const bool isInteriorBoundary)
{
static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
if (!gravity)
return Scalar(0.0);
else
{
// Treat interior Dirichlet boundaries differently
if (isInteriorBoundary)
{
const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
return data.facetVolVars(fvGeometry).density(phaseIdx);
}
// use arithmetic mean of the densities around the scvf
if (!scvf.boundary())
{
Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
for (auto outsideIdx : scvf.outsideScvIndices())
rho += elemVolVars[outsideIdx].density(phaseIdx);
return rho/(scvf.outsideScvIndices().size()+1);
}
else
return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
}
}
static Scalar computeInteriorBoundaryContribution(const Problem& problem,
const FVElementGeometry& fvGeometry,
const FluxVariablesCache& fluxVarsCache,
unsigned int phaseIdx, Scalar rho)
{
static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
// obtain the transmissibilites associated with all pressures
const auto& tij = fluxVarsCache.advectionTij();
// the interior dirichlet boundaries local indices start after
// the cell and the domain Dirichlet boundary pressures
const auto startIdx = fluxVarsCache.advectionVolVarsStencil().size();
// add interior Dirichlet boundary contributions
Scalar flux = 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);
......@@ -264,60 +284,11 @@ public:
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);
flux += tij[startIdx + data.localIndexInInteractionVolume()]*h;
}
}
// return overall resulting flux
const Scalar interiorNeumannFlux = facetCoupling ? cij*facetCouplingFluxes : 0.0;
return useTpfaBoundary ?
flux + interiorNeumannFlux :
flux + interiorNeumannFlux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
}
private:
static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const unsigned int phaseIdx)
{
// use arithmetic mean of the densities around the scvf
if (!scvf.boundary())
{
Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
for (auto outsideIdx : scvf.outsideScvIndices())
rho += elemVolVars[outsideIdx].density(phaseIdx);
return rho/(scvf.outsideScvIndices().size()+1);
}
else
return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
return flux;
}
};
......
......@@ -43,6 +43,8 @@ namespace Dumux
template <class TypeTag>
class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
using Implementation = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
......@@ -90,9 +92,6 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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.
......@@ -113,11 +112,6 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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
......@@ -130,7 +124,6 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
// 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_;
......@@ -201,63 +194,21 @@ public:
// get the scaling factor for the effective diffusive fluxes
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);
} ();
const auto effFactor = Implementation::computeEffectivityFactor(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
// if factor is zero, the flux will end up zero anyway
if (effFactor == 0.0)
return 0.0;
// lambda functions depending on if we use mole or mass fractions
auto getX = [useMoles, phaseIdx, compIdx] (const VolumeVariables& volVars)
auto getX = [useMoles, phaseIdx, compIdx] (const auto& volVars)
{ return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
auto getRho = [useMoles, phaseIdx] (const VolumeVariables& volVars)
auto getRho = [useMoles, phaseIdx] (const auto& 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);
} ();
const auto rho = Implementation::interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, getRho, isInteriorBoundary);
// calculate Tij*xj
Scalar flux(0.0);
......@@ -269,69 +220,30 @@ public:
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);
flux += Implementation::computeInteriorBoundaryContribution(fvGeometry, fluxVarsCache, getX, phaseIdx, compIdx);
// The vector of interior neumann fluxes
CoefficientVector 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;
// return overall resulting flux
return useTpfaBoundary ? flux*rho*effFactor : flux*rho*effFactor + fluxVarsCache.componentNeumannFlux(compIdx);
}
// add neumann fluxes for interior Neumann faces
if (facetCoupling && data.faceType() == MpfaFaceTypes::interiorNeumann)
template<typename GetRhoFunction>
static Scalar interpolateDensity(const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const FluxVariablesCache& fluxVarsCache,
const GetRhoFunction& getRho,
const bool isInteriorBoundary)
{
// 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 = [&] ()
// maybe use the density of the interior BC on the facet
if (isInteriorBoundary)
{
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);
const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
return getRho(data.facetVolVars(fvGeometry));
}
private:
template<typename GetRhoFunction>
static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const GetRhoFunction& getRho)
{
// use arithmetic mean of the densities around the scvf
if (!scvf.boundary())
{
......@@ -348,15 +260,40 @@ private:
//! scaled to get the effective diffusivity. For this we use the effective diffusivity with
//! a diffusion coefficient of 1.0 as input. Then we scale the transmissibilites during flux
//! calculation (above) with the harmonic average of the two factors
static Scalar calculateEffectiveDiffusivityFactor(const ElementVolumeVariables& elemVolVars,
static Scalar computeEffectivityFactor(const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const unsigned int phaseIdx)
const FluxVariablesCache& fluxVarsCache,
const unsigned int phaseIdx,
const bool isInteriorBoundary)
{
// Treat interior boundaries differently
if (isInteriorBoundary)
{
const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
// use harmonic mean between the interior and the facet volvars
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
{
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);
}
}
// use the harmonic mean between inside and outside
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(),
insideVolVars.saturation(phaseIdx),
/*Diffusion coefficient*/ 1.0);
if (!scvf.boundary())
{
// interpret outside factor as arithmetic mean
......@@ -376,6 +313,28 @@ private:
return factor;
}
template<typename GetXFunction>
static Scalar computeInteriorBoundaryContribution(const FVElementGeometry& fvGeometry,
const FluxVariablesCache& fluxVarsCache,
const GetXFunction& getX,
unsigned int phaseIdx, unsigned int compIdx)
{
// obtain the transmissibilites associated with all pressures
const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx);
// the interior dirichlet boundaries local indices start after
// the cell and the domain Dirichlet boundary pressures
const auto startIdx = fluxVarsCache.diffusionVolVarsStencil(phaseIdx, compIdx).size();
// add interior Dirichlet boundary contributions
Scalar flux = 0.0;
for (auto&& data : fluxVarsCache.interiorBoundaryData())
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
flux += tij[startIdx + data.localIndexInInteractionVolume()]*getX(data.facetVolVars(fvGeometry));
return flux;
}
};
} // end namespace
......
......@@ -41,6 +41,8 @@ namespace Dumux
template <class TypeTag>
class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
using Implementation = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
......@@ -59,7 +61,6 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
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);
......@@ -78,9 +79,6 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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
......@@ -94,11 +92,6 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
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
......@@ -108,7 +101,6 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
// Quantities associated with heat conduction
Stencil heatConductionVolVarsStencil_;
CoefficientVector heatConductionTij_;
CoefficientVector heatConductionCij_;
Scalar heatNeumannFlux_;
};
......@@ -177,68 +169,30 @@ public:
if (!enableInteriorBoundaries)
return useTpfaBoundary ? flux : flux + fluxVarsCache.heatNeumannFlux();
//////////////////////////////////////////////////////////////////
// Handle interior boundaries
//////////////////////////////////////////////////////////////////
flux += Implementation::computeInteriorBoundaryContribution(fvGeometry, fluxVarsCache);
// get coefficients to transform the vector of interior neumann boundary conditions
const auto& cij = fluxVarsCache.heatConductionCij();
// The Vector of interior neumann fluxes
CoefficientVector 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
flux += tij[localIdx + data.localIndexInInteractionVolume()]*data.facetVolVars(fvGeometry).temperature();
// return overall resulting flux
return useTpfaBoundary ? flux : flux + fluxVarsCache.heatNeumannFlux();
}
// add neumann fluxes for interior Neumann faces
if (data.faceType() == MpfaFaceTypes::interiorNeumann)
static Scalar computeInteriorBoundaryContribution(const FVElementGeometry& fvGeometry,
const FluxVariablesCache& fluxVarsCache)
{
if (facetCoupling)
{
// get the scvf corresponding to actual interior boundary face
const auto& curScvf = fvGeometry.scvf(data.scvfIndex());
// obtain the transmissibilites associated with all pressures
const auto& tij = fluxVarsCache.heatConductionTij();
// obtain the complete data on the facet element
const auto completeFacetData = data.completeCoupledFacetData();
// the interior dirichlet boundaries local indices start after
// the cell and the domain Dirichlet boundary pressures
const auto startIdx = fluxVarsCache.heatConductionVolVarsStencil().size();
// calculate "lekage factor"
const auto n = curScvf.unitOuterNormal();
const auto v = [&] ()
{
auto res = n;
res *= -0.5*completeFacetData.volVars().extrusionFactor();
res -= curScvf.ipGlobal();
res += curScvf.facetCorner();
res /= res.two_norm2();
return res;
} ();
// get the thermal conductivity in the facet element
const auto facetLambda = ThermalConductivityModel::effectiveThermalConductivity(completeFacetData.volVars(),
completeFacetData.spatialParams(),
completeFacetData.element(),
completeFacetData.fvGeometry(),
completeFacetData.scv());
// add value to vector of interior neumann fluxes