diff --git a/dumux/discretization/box/darcyslaw.hh b/dumux/discretization/box/darcyslaw.hh index d432223e5fe4a6adfb7cd1ba20a5dcf06c09ce9f..f4b1bc04b99f9ecc2e7eb18234902dced45e6354 100644 --- a/dumux/discretization/box/darcyslaw.hh +++ b/dumux/discretization/box/darcyslaw.hh @@ -64,7 +64,6 @@ class DarcysLawImplementation using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; using CoordScalar = typename GridView::ctype; - using Stencil = std::vector; enum { dim = GridView::dimension}; enum { dimWorld = GridView::dimensionworld}; @@ -137,14 +136,6 @@ public: return -1.0*(KGradP*scvf.unitOuterNormal())*scvf.area(); } - // This is for compatibility with the cc methods. The flux stencil info is obsolete for the box method. - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { - return Stencil(0); - } private: static DimVector applyPermeability_(const DimWorldMatrix& K, const DimVector& gradI) { diff --git a/dumux/discretization/box/elementfluxvariablescache.hh b/dumux/discretization/box/elementfluxvariablescache.hh index e0b699c9f9fdf63d26e1fcf6619adcdb16b2d4f4..97fb569f4b4fb451e2909c90a64a7970984ee67a 100644 --- a/dumux/discretization/box/elementfluxvariablescache.hh +++ b/dumux/discretization/box/elementfluxvariablescache.hh @@ -132,7 +132,7 @@ public: // temporary resizing of the cache fluxVarsCache_.resize(fvGeometry.numScvf()); for (auto&& scvf : scvfs(fvGeometry)) - (*this)[scvf].update(globalFluxVarsCache().problem_(), element, fvGeometry, scvf); + (*this)[scvf].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf); } void bindScvf(const Element& element, @@ -141,7 +141,7 @@ public: const SubControlVolumeFace& scvf) { fluxVarsCache_.resize(fvGeometry.numScvf()); - (*this)[scvf].update(globalFluxVarsCache().problem_(), element, fvGeometry, scvf); + (*this)[scvf].update(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf); } // access operator diff --git a/dumux/discretization/box/fickslaw.hh b/dumux/discretization/box/fickslaw.hh index fc2d13eedb07e8cd102d16f55829e7835d979b2d..0ae115f46c3c862313e63923addb6d894ed33a63 100644 --- a/dumux/discretization/box/fickslaw.hh +++ b/dumux/discretization/box/fickslaw.hh @@ -39,7 +39,6 @@ namespace Properties { // forward declaration of properties NEW_PROP_TAG(NumPhases); -NEW_PROP_TAG(FluidState); NEW_PROP_TAG(FluidSystem); NEW_PROP_TAG(EffectiveDiffusivityModel); } @@ -53,7 +52,6 @@ class FicksLawImplementation { using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState); using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); @@ -65,7 +63,6 @@ class FicksLawImplementation using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using IndexType = typename GridView::IndexSet::IndexType; - using Stencil = typename std::vector; using Element = typename GridView::template Codim<0>::Entity; @@ -134,13 +131,6 @@ public: return -1.0*(DGradX*scvf.unitOuterNormal())*scvf.area(); } - // This is for compatibility with the cc methods. The flux stencil info is obsolete for the box method. - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { return Stencil(0); } - private: static GlobalPosition applyDiffusionTensor_(const DimWorldMatrix& D, const GlobalPosition& gradI) { diff --git a/dumux/discretization/box/fourierslaw.hh b/dumux/discretization/box/fourierslaw.hh index c7288999b3cb16a703de35eca14f528032d6da86..d77e89a689cd78e22c0b4dfb388fe9793876716d 100644 --- a/dumux/discretization/box/fourierslaw.hh +++ b/dumux/discretization/box/fourierslaw.hh @@ -65,7 +65,6 @@ class FouriersLawImplementation using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using IndexType = typename GridView::IndexSet::IndexType; - using Stencil = typename std::vector; using Element = typename GridView::template Codim<0>::Entity; @@ -123,13 +122,6 @@ public: return -1.0*(lambdaGradT*scvf.unitOuterNormal())*scvf.area(); } - // This is for compatibility with the cc methods. The flux stencil info is obsolete for the box method. - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { return Stencil(0); } - private: static GlobalPosition applyHeatConductivityTensor_(const DimWorldMatrix& Lambda, const GlobalPosition& gradI) { diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh index c024ebbf6ae14ef534f1ea1cf6c8c1610d576bbd..3aedae63d47b7268288c4c3d0439e20775228a23 100644 --- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh @@ -31,8 +31,8 @@ #include #include - #include +#include namespace Dumux { @@ -50,21 +50,29 @@ NEW_PROP_TAG(ProblemEnableGravity); template class DarcysLawImplementation { + 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; - 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,27 +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); - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { - const auto& globalFvGeometry = problem.model().globalFvGeometry(); + ////////////////////////////////////////////////////////////////// + // Handle interior boundaries + ////////////////////////////////////////////////////////////////// - // return the scv (element) indices in the interaction region - if (globalFvGeometry.scvfTouchesBoundary(scvf)) - return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices(); - else - return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices(); + // 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); } private: diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh index 01d8be6c03fbb1763b8549f722d3797dac7d7b74..2bae60663e2c3bcc89f6e56db3790a9e3334032d 100644 --- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh +++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh @@ -127,7 +127,7 @@ public: const ElementVolumeVariables& elemVolVars) { // TODO - DUNE_THROW(Dune::InvalidStateException, "Does local flux var cache binding make sense in general?"); + DUNE_THROW(Dune::NotImplemented, "Local element binding of the flux variables cache in mpfa schemes"); } // This function is called by the CCLocalResidual before flux calculations during assembly. @@ -187,7 +187,7 @@ public: const SubControlVolumeFace& scvf) { // TODO - DUNE_THROW(Dune::InvalidStateException, "Does binding of one scvf make sense in general?"); + DUNE_THROW(Dune::NotImplemented, "Local element binding of the flux variables cache in mpfa schemes"); } // access operators in the case of no caching @@ -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); } diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index 6620dee3b07cec0cd8df9f3d2c667576fe699ae1..e0470acfc44d4a3044250574227aab05ccb9264d 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -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; } diff --git a/dumux/discretization/cellcentered/mpfa/facetypes.hh b/dumux/discretization/cellcentered/mpfa/facetypes.hh index 87d526e167db11fbf695e0a603483bd62b4e1e74..c278de9e3913b3938248805445b65d5cdc696280 100644 --- a/dumux/discretization/cellcentered/mpfa/facetypes.hh +++ b/dumux/discretization/cellcentered/mpfa/facetypes.hh @@ -29,7 +29,9 @@ namespace Dumux { interior, neumann, - dirichlet + dirichlet, + interiorNeumann, + interiorDirichlet }; } // end namespace diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh index fb27d52806dec57ffc6cfc30be9e20a92d33bbf1..6aaf11e1881cc9b723bd65c3394dbc17843047ca 100644 --- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh @@ -46,18 +46,26 @@ class FicksLawImplementation 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; + + 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 @@ -76,8 +84,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,28 +145,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); - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { - const auto& globalFvGeometry = problem.model().globalFvGeometry(); + ////////////////////////////////////////////////////////////////// + // Handle interior boundaries + ////////////////////////////////////////////////////////////////// - // return the scv (element) indices in the interaction region - if (globalFvGeometry.scvfTouchesBoundary(scvf)) - return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices(); - else - return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices(); + // 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); } private: @@ -141,19 +256,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; diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh index 6b05dd47b2d4b1a9fe30834d375298a1c5745d06..0ca84820d817e6e4201ff94ffbe873ad8649abf2 100644 --- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh +++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh @@ -24,6 +24,7 @@ #define DUMUX_DISCRETIZATION_CCMPFA_GLOBAL_FLUXVARSCACHE_FILLER_HH #include +#include #include "fluxvariablescachefillerbase.hh" namespace Dumux @@ -44,8 +45,7 @@ using CCMpfaFluxVariablesCacheFiller = CCMpfaFluxVariablesCacheFillerImplementat //! Implementation for purely advective problems template -class CCMpfaFluxVariablesCacheFillerImplementation - : public CCMpfaAdvectionCacheFiller +class CCMpfaFluxVariablesCacheFillerImplementation : public CCMpfaAdvectionCacheFiller { using AdvectionFiller = CCMpfaAdvectionCacheFiller; @@ -56,12 +56,15 @@ class CCMpfaFluxVariablesCacheFillerImplementation 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 InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume); + using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume); 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 solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection); + 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); public: //! function to fill the flux var caches @@ -73,13 +76,51 @@ public: const SubControlVolumeFace& scvf, FluxVarsCacheContainer& fluxVarsCacheContainer) { - // forward to the filler for the advective quantities - AdvectionFiller::fillCaches(problem, - element, - fvGeometry, - elemVolVars, - scvf, - fluxVarsCacheContainer); + // Instantiate interaction volume depending on if scvf touches the boundary or not + if (problem.model().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); + } + } + 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); + } } //! function to update the flux var caches during derivative calculation @@ -91,13 +132,19 @@ public: const SubControlVolumeFace& scvf, FluxVarsCacheContainer& fluxVarsCacheContainer) { - // Do basically nothing if advection is solution-independent. Although we have - // to set the update status to true here as it has been set to false before. - // This is for compatibility reasons with compositional models. - if (!solDependentAdvection && useTpfaBoundary) - fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true); - // If we don't use tpfa on the boundaries we have to update the whole thing anyway, - // as we need the local matrices to assemble the neumann fluxes (which could be solution-dependent) + //! 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) + { + 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); } @@ -105,9 +152,8 @@ public: //! Implementation for problems considering advection & diffusion template -class CCMpfaFluxVariablesCacheFillerImplementation - : public CCMpfaAdvectionCacheFiller, - public CCMpfaDiffusionCacheFiller +class CCMpfaFluxVariablesCacheFillerImplementation : public CCMpfaAdvectionCacheFiller, + public CCMpfaDiffusionCacheFiller { using AdvectionFiller = CCMpfaAdvectionCacheFiller; using DiffusionFiller = CCMpfaDiffusionCacheFiller; @@ -119,12 +165,33 @@ class CCMpfaFluxVariablesCacheFillerImplementation 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); using Element = typename GridView::template Codim<0>::Entity; - static const bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary); - static const bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection); - static const bool solDependentDiffusion = GET_PROP_VALUE(TypeTag, SolutionDependentMolecularDiffusion); + 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; + + //! 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; + + using BoolPair = std::pair; public: //! function to fill the flux var caches @@ -134,10 +201,76 @@ public: const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf, - FluxVarsCacheContainer& fluxVarsCacheContainer) + FluxVarsCacheContainer& fluxVarsCacheContainer, + const BoolPair doAdvectionOrDiffusion = BoolPair(true, true)) { - AdvectionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); - DiffusionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + // Instantiate interaction volume depending on if scvf touches the boundary or not + if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf)) + { + BoundaryInteractionVolume iv(problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf), + problem, + fvGeometry, + elemVolVars); + + // forward to the filler for the advective quantities + if (doAdvectionOrDiffusion.first) + AdvectionFiller::fillCaches(problem, + element, + fvGeometry, + elemVolVars, + scvf, + iv, + fluxVarsCacheContainer); + + // forward to the filler for the diffusive quantities + if (doAdvectionOrDiffusion.second) + DiffusionFiller::fillCaches(problem, + element, + fvGeometry, + elemVolVars, + scvf, + iv, + fluxVarsCacheContainer); + + //set update status and maybe update data on interior boundaries + for (const auto scvfIdxJ : iv.globalScvfs()) + { + if (enableInteriorBoundaries) + fluxVarsCacheContainer[scvfIdxJ].updateInteriorBoundaryData(iv, fvGeometry.scvf(scvfIdxJ)); + fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true); + } + } + else + { + InteractionVolume iv(problem.model().globalFvGeometry().interactionVolumeSeed(scvf), + problem, + fvGeometry, + elemVolVars); + + // forward to the filler for the advective quantities + if (doAdvectionOrDiffusion.first) + AdvectionFiller::fillCaches(problem, + element, + fvGeometry, + elemVolVars, + scvf, + iv, + fluxVarsCacheContainer); + + // forward to the filler for the diffusive quantities + if (doAdvectionOrDiffusion.second) + DiffusionFiller::fillCaches(problem, + element, + fvGeometry, + elemVolVars, + scvf, + iv, + fluxVarsCacheContainer); + + // set update status + for (const auto scvfIdxJ : iv.globalScvfs()) + fluxVarsCacheContainer[scvfIdxJ].setUpdateStatus(true); + } } //! function to update the flux var caches during derivative calculation @@ -149,35 +282,28 @@ public: const SubControlVolumeFace& scvf, FluxVarsCacheContainer& fluxVarsCacheContainer) { - // Do basically nothing if the parameters are solution-independent. Although we have - // to set the update status to true here as it has been set to false before. - // This is for compatibility reasons with compositional models. - - // TODO: How to treat !useTpfaBoundary??? - if (!solDependentAdvection && !solDependentDiffusion) - { - if (useTpfaBoundary) - fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true); - else - DiffusionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + const bool touchesDomainBoundary = problem.model().globalFvGeometry().touchesDomainBoundary(scvf); - // we're done here - return; - } + // maybe update Advection or diffusion or both + const bool updateAdvection = doAdvectionUpdate && touchesDomainBoundary; + const bool updateDiffusion = doDiffusionUpdate && touchesDomainBoundary; - if (solDependentAdvection) - AdvectionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + 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)); - if (solDependentDiffusion || !useTpfaBoundary) - DiffusionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + // tell the cache that it has been updated in any case + fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true); } }; //! Implementation for problems considering advection & heat conduction template -class CCMpfaFluxVariablesCacheFillerImplementation - : public CCMpfaAdvectionCacheFiller, - public CCMpfaHeatConductionCacheFiller +class CCMpfaFluxVariablesCacheFillerImplementation : public CCMpfaAdvectionCacheFiller, + public CCMpfaHeatConductionCacheFiller { using AdvectionFiller = CCMpfaAdvectionCacheFiller; using HeatConductionFiller = CCMpfaHeatConductionCacheFiller; @@ -188,13 +314,30 @@ class CCMpfaFluxVariablesCacheFillerImplementation 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); using Element = typename GridView::template Codim<0>::Entity; - static const bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary); - static const bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection); - static const bool solDependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction); + 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 solDependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction); + static constexpr bool heatConductionUsesMpfa = HeatConductionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa; + + //! 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; + + using BoolPair = std::pair; public: //! function to fill the flux var caches @@ -204,10 +347,76 @@ public: const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf, - FluxVarsCacheContainer& fluxVarsCacheContainer) + FluxVarsCacheContainer& fluxVarsCacheContainer, + const BoolPair doAdvectionOrHeatConduction = BoolPair(true, true)) { - AdvectionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); - HeatConductionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + // Instantiate interaction volume depending on if scvf touches the boundary or not + if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf)) + { + 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()) + { + 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 (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 @@ -219,36 +428,29 @@ public: const SubControlVolumeFace& scvf, FluxVarsCacheContainer& fluxVarsCacheContainer) { - // Do basically nothing if the parameters are solution-independent. Although we have - // to set the update status to true here as it has been set to false before. - // This is for compatibility reasons with compositional models. + const bool touchesDomainBoundary = problem.model().globalFvGeometry().touchesDomainBoundary(scvf); - // TODO: How to treat !useTpfaBoundary??? - if (!solDependentAdvection && !solDependentHeatConduction) - { - if (useTpfaBoundary) - fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true); - else - fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + // maybe update Advection or diffusion or both + const bool updateAdvection = doAdvectionUpdate && touchesDomainBoundary; + const bool updateHeatConduction = doHeatConductionUpdate && touchesDomainBoundary; - // we're done here - return; - } - - if (solDependentAdvection || !useTpfaBoundary) - AdvectionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + 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)); - if (solDependentHeatConduction || !useTpfaBoundary) - HeatConductionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + // tell the cache that it has been updated in any case + fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true); } }; //! Implementation for problems considering advection, diffusion & heat conduction template -class CCMpfaFluxVariablesCacheFillerImplementation - : public CCMpfaAdvectionCacheFiller, - public CCMpfaDiffusionCacheFiller, - public CCMpfaHeatConductionCacheFiller +class CCMpfaFluxVariablesCacheFillerImplementation : public CCMpfaAdvectionCacheFiller, + public CCMpfaDiffusionCacheFiller, + public CCMpfaHeatConductionCacheFiller { using AdvectionFiller = CCMpfaAdvectionCacheFiller; using DiffusionFiller = CCMpfaDiffusionCacheFiller; @@ -261,13 +463,35 @@ class CCMpfaFluxVariablesCacheFillerImplementation 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 const bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary); - static const bool solDependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection); - static const bool solDependentDiffusion = GET_PROP_VALUE(TypeTag, SolutionDependentMolecularDiffusion); - static const bool solDependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction); + 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; public: //! function to fill the flux var caches @@ -277,11 +501,96 @@ public: const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf, - FluxVarsCacheContainer& fluxVarsCacheContainer) + FluxVarsCacheContainer& fluxVarsCacheContainer, + const BoolTriplet doAdvectionOrDiffusionOrHeatConduction = BoolTriplet({true, true, true})) { - AdvectionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); - DiffusionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); - HeatConductionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); + // Instantiate interaction volume depending on if scvf touches the boundary or not + if (problem.model().globalFvGeometry().touchesInteriorOrDomainBoundary(scvf)) + { + 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); + } } //! function to update the flux var caches during derivative calculation @@ -293,33 +602,33 @@ public: const SubControlVolumeFace& scvf, FluxVarsCacheContainer& fluxVarsCacheContainer) { - // Do basically nothing if the parameters are solution-independent. Although we have - // to set the update status to true here as it has been set to false before. - // This is for compatibility reasons with compositional models. - - // TODO: How to treat !useTpfaBoundary??? - if (!solDependentAdvection && !solDependentDiffusion && !solDependentHeatConduction) - { - if (useTpfaBoundary) - fluxVarsCacheContainer[scvf.index()].setUpdateStatus(true); - else - { - DiffusionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); - HeatConductionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); - } - - // we're done here - return; - } - - if (solDependentAdvection) - AdvectionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); - - if (solDependentDiffusion || !useTpfaBoundary) - DiffusionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCacheContainer); - - if (solDependentHeatConduction || !useTpfaBoundary) - HeatConductionFiller::fillCaches(problem, element, fvGeometry, elemVolVars, scvf, 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); } }; diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh index 9790e04e1daa551bfcafadff70c4170694185fe6..9ef0820f3a0142428ef6c9d16513b4ca0eb487f9 100644 --- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh +++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefillerbase.hh @@ -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 + template 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); public: //! function to fill the flux var caches - template + template static void fillCaches(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf, + InteractionVolume& iv, FluxVarsCacheContainer& fluxVarsCacheContainer) { - //! If the problem does not use an mpfa method for diffusion, do nothing - //! This is known at compile time so it gets optimized away - if (MolecularDiffusionType::myDiscretizationMethod != DiscretizationMethods::CCMpfa) - return; - - if (problem.model().globalFvGeometry().scvfTouchesBoundary(scvf)) + //! 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) { - const auto& seed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf); - BoundaryInteractionVolume iv(seed, problem, fvGeometry, elemVolVars); - - //! 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) { - 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 (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(); - auto& cache = fluxVarsCacheContainer[scvfIdx]; - cache.updateDiffusion(scvf, iv, phaseIdx, compIdx); - cache.setUpdateStatus(true); - - for (const auto scvfIdxJ : iv.globalScvfs()) + if (scvfIdxJ != scvfIdx) { - if (scvfIdxJ != scvfIdx) - { - const auto& scvfJ = fvGeometry.scvf(scvfIdxJ); - auto& cacheJ = fluxVarsCacheContainer[scvfIdxJ]; - cacheJ.updateDiffusion(scvfJ, iv, phaseIdx, compIdx); - cacheJ.setUpdateStatus(true); - } - } - } - } - } - else - { - const auto& seed = problem.model().globalFvGeometry().interactionVolumeSeed(scvf); - InteractionVolume iv(seed, problem, fvGeometry, elemVolVars); - - //! 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 - // TODO: How to include effective diffusivity? - 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(); - auto& cache = fluxVarsCacheContainer[scvfIdx]; - cache.updateDiffusion(scvf, iv, phaseIdx, compIdx); - cache.setUpdateStatus(true); + // store scvf pointer and local face data + const auto& scvfJ = fvGeometry.scvf(scvfIdxJ); + const auto scvfJLocalFaceData = iv.getLocalFaceData(scvfJ); - for (const auto scvfIdxJ : iv.globalScvfs()) - { - if (scvfIdxJ != scvfIdx) - { - const auto& scvfJ = fvGeometry.scvf(scvfIdxJ); - auto& cacheJ = fluxVarsCacheContainer[scvfIdxJ]; - cacheJ.updateDiffusion(scvfJ, iv, phaseIdx, compIdx); - cacheJ.setUpdateStatus(true); - } + // update cache + fluxVarsCacheContainer[scvfIdxJ].updateDiffusion(iv, scvfJ, scvfJLocalFaceData, phaseIdx, compIdx); } } } @@ -280,7 +183,6 @@ class CCMpfaHeatConductionCacheFiller 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 ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel); @@ -294,19 +196,15 @@ class CCMpfaHeatConductionCacheFiller public: //! function to fill the flux var caches - template + template static void fillCaches(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf, + InteractionVolume& iv, FluxVarsCacheContainer& fluxVarsCacheContainer) { - //! If the problem does not use an mpfa method for diffusion, do nothing - //! This is known at compile time so it gets optimized away - if (HeatConductionType::myDiscretizationMethod != DiscretizationMethods::CCMpfa) - return; - // lambda function to get the thermal conductivity auto getLambda = [&problem, &fvGeometry](const Element& element, const VolumeVariables& volVars, @@ -317,68 +215,26 @@ public: fvGeometry, scv); }; - if (problem.model().globalFvGeometry().scvfTouchesBoundary(scvf)) - { - const auto& seed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf); - BoundaryInteractionVolume iv(seed, problem, fvGeometry, elemVolVars); - - // solve for the transmissibilities - iv.solveLocalSystem(getLambda); - - // update the caches - const auto scvfIdx = scvf.index(); - auto& cache = fluxVarsCacheContainer[scvfIdx]; - cache.updateHeatConduction(scvf, iv); - cache.setUpdateStatus(true); - - //! set the neumann boundary conditions in case we do not use tpfa on the boundary - if (!useTpfaBoundary) - { - iv.assembleNeumannFluxes(energyEqIdx); - cache.updateHeatNeumannFlux(scvf, iv); - } + // solve for the transmissibilities + iv.solveLocalSystem(getLambda); - for (const auto scvfIdxJ : iv.globalScvfs()) - { - if (scvfIdxJ != scvfIdx) - { - const auto& scvfJ = fvGeometry.scvf(scvfIdxJ); - auto& cacheJ = fluxVarsCacheContainer[scvfIdxJ]; - cacheJ.updateHeatConduction(scvfJ, iv); - cacheJ.setUpdateStatus(true); + // update the caches + const auto scvfIdx = scvf.index(); + const auto scvfLocalFaceData = iv.getLocalFaceData(scvf); + auto& scvfCache = fluxVarsCacheContainer[scvfIdx]; + scvfCache.updateHeatConduction(iv, scvf, scvfLocalFaceData); - //! set the neumann boundary conditions in case we do not use tpfa on the boundary - if (!useTpfaBoundary) - { - iv.assembleNeumannFluxes(energyEqIdx); - cache.updateHeatNeumannFlux(scvf, iv); - } - } - } - } - else + //! Update the caches 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); - - // solve for the transmissibilities - iv.solveLocalSystem(getLambda); - - // update the caches - const auto scvfIdx = scvf.index(); - auto& cache = fluxVarsCacheContainer[scvfIdx]; - cache.updateHeatConduction(scvf, iv); - cache.setUpdateStatus(true); - - for (const auto scvfIdxJ : iv.globalScvfs()) + if (scvfIdxJ != scvfIdx) { - if (scvfIdxJ != scvfIdx) - { - const auto& scvfJ = fvGeometry.scvf(scvfIdxJ); - auto& cacheJ = fluxVarsCacheContainer[scvfIdxJ]; - cacheJ.updateHeatConduction(scvfJ, iv); - cacheJ.setUpdateStatus(true); - } + // 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); } } } diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh index fcca346bcc23cf14e7aad66705c0f4f5cce965aa..03d00e0db4d1a1370d0714c6edcade186a4a4b08 100644 --- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh @@ -43,15 +43,26 @@ class FouriersLawImplementation { using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); 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 ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache); + 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 GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; - using Stencil = typename std::vector; + + 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); + + static constexpr int energyEqIdx = GET_PROP_TYPE(TypeTag, Indices)::energyEqIdx; public: // state the discretization method this implementation belongs to @@ -68,28 +79,91 @@ public: const auto& volVarsStencil = fluxVarsCache.heatConductionVolVarsStencil(); const auto& tij = fluxVarsCache.heatConductionTij(); + const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary(); + // For interior Neumann boundaries when using Tpfa for Neumann boundary conditions, we simply + // return the user-specified flux + if (isInteriorBoundary + && useTpfaBoundary + && fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann) + return scvf.area()* + elemVolVars[scvf.insideScvIdx()].extrusionFactor()* + problem.neumann(element, + fvGeometry, + elemVolVars, + scvf)[energyEqIdx]; + // calculate Tij*tj Scalar flux(0.0); unsigned int localIdx = 0; for (const auto volVarIdx : volVarsStencil) flux += tij[localIdx++]*elemVolVars[volVarIdx].temperature(); - // return heat conduction flux - return flux + fluxVarsCache.heatNeumannFlux(); - } + // if no interior boundaries are present, return heat conduction flux + if (!enableInteriorBoundaries) + return useTpfaBoundary ? flux : flux + fluxVarsCache.heatNeumannFlux(); - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { - const auto& globalFvGeometry = problem.model().globalFvGeometry(); + ////////////////////////////////////////////////////////////////// + // Handle interior boundaries + ////////////////////////////////////////////////////////////////// + + // get coefficients to transform the vector of interior neumann boundary conditions + const auto& cij = fluxVarsCache.heatConductionCij(); + + // 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 + flux += tij[localIdx + data.localIndexInInteractionVolume()]*data.facetVolVars(fvGeometry).temperature(); + } + + // add neumann fluxes for interior Neumann faces + if (data.faceType() == MpfaFaceTypes::interiorNeumann) + { + if (facetCoupling) + { + // get the scvf corresponding to actual interior boundary face + const auto& curScvf = fvGeometry.scvf(data.scvfIndex()); + + // obtain the complete data on the facet element + const auto completeFacetData = data.completeCoupledFacetData(); + + // 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 + facetCouplingFluxes[data.localIndexInInteractionVolume()] += MpfaHelper::nT_M_v(n, + facetLambda, + v); + } + } + } - // return the scv (element) indices in the interaction region - if (globalFvGeometry.scvfTouchesBoundary(scvf)) - return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices(); - else - return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices(); + // return overall resulting flux + const Scalar interiorNeumannFlux = facetCoupling ? cij*facetCouplingFluxes : 0.0; + return useTpfaBoundary ? + flux + interiorNeumannFlux : + flux + interiorNeumannFlux + fluxVarsCache.heatNeumannFlux(); } }; diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh index 05e95ae2ee9a5f97fe39457041d0dfc24ff239c2..e79a41db5eca80335665f9585cb48af9f002a811 100644 --- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh +++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh @@ -351,9 +351,6 @@ private: const auto& scvFaceIndices = globalFvGeometry().scvfIndicesOfScv(eIdx); const auto& neighborVolVarIndices = globalFvGeometry().neighborVolVarIndices(eIdx); - // Instantiate helper class to pass it to scvf constructors - MpfaHelper helper; - // the quadrature point to be used on the scvf const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q); @@ -407,8 +404,8 @@ private: if (globalFvGeometry().isGhostVertex(vIdxGlobal)) continue; - scvfs_.emplace_back(helper, - helper.getScvfCorners(isCorners, c), + scvfs_.emplace_back(MpfaHelper(), + MpfaHelper::getScvfCorners(isCorners, c), is.centerUnitOuterNormal(), vIdxGlobal, vIdxLocal, @@ -437,9 +434,6 @@ private: const auto& scvFaceIndices = globalFvGeometry().scvfIndicesOfScv(eIdxGlobal); const auto& neighborVolVarIndices = globalFvGeometry().neighborVolVarIndices(eIdxGlobal); - // Instantiate a helper class to pass it to scvf constructors - MpfaHelper helper; - // the quadrature point to be used on the scvf const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q); @@ -491,8 +485,8 @@ private: // only build the scvf if it is in the list of necessary indices if (MpfaHelper::contains(scvfIndices, scvFaceIndices[scvfCounter])) { - neighborScvfs_.emplace_back(helper, - helper.getScvfCorners(isCorners, c), + neighborScvfs_.emplace_back(MpfaHelper(), + MpfaHelper::getScvfCorners(isCorners, c), is.centerUnitOuterNormal(), vIdxGlobal, vIdxLocal, diff --git a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh index b275a738732c055a5b2106873064d050632b2a4a..fe3f89a5f4bea2e9a336267d331b1e103c9edc19 100644 --- a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh +++ b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh @@ -72,12 +72,13 @@ class CCMpfaGlobalFVGeometryBase using IndexType = typename GridView::IndexSet::IndexType; using LocalIndexType = typename InteractionVolume::LocalIndexType; - static const int dim = GridView::dimension; - static const int dimWorld = GridView::dimensionworld; - + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; using GlobalPosition = Dune::FieldVector; using ReferenceElements = typename Dune::ReferenceElements; + static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries); + public: //! Constructor CCMpfaGlobalFVGeometryBase(const GridView gridView) @@ -87,17 +88,29 @@ public: std::size_t numScv() const { return scvs_.size(); } - //! The total number of sun control volume faces + //! The total number of sub control volume faces std::size_t numScvf() const { return scvfs_.size(); } - //! The total number of boundary sub control volume faces - std::size_t numBoundaryScvf() const - { return numBoundaryScvf_; } + //! The total number of scvfs on the domain boundaries + std::size_t numDomainBoundaryScvf() const + { return numDomainBoundaryScvf_; } - //! The total number of vertices on the boundary - std::size_t numBoundaryVertices() const - { return numBoundaryVertices_; } + //! The total number of scvfs on the interior boundaries + std::size_t numInteriorBoundaryScvf() const + { return numInteriorBoundaryScvf_; } + + //! The total number of scvfs connected to branching points + std::size_t numBranchingPointScvf() const + { return numBranchingPointScvf_; } + + //! The total number of vertices connected to branching points + std::size_t numBranchingPointVertices() const + { return numBranchingPointVertices_; } + + //! The total number of vertices on interior and domain boundaries + std::size_t numInteriorOrDomainBoundaryVertices() const + { return numInteriorOrDomainBoundaryVertices_; } // Get an element from a sub control volume contained in it Element element(const SubControlVolume& scv) const @@ -115,60 +128,70 @@ public: const BoundaryInteractionVolumeSeed& boundaryInteractionVolumeSeed(const SubControlVolumeFace& scvf) const { return globalInteractionVolumeSeeds_.boundarySeed(scvf); } - //! Returns whether or not an scvf touches the boundary (has to be called before getting an interaction volume) - bool scvfTouchesBoundary(const SubControlVolumeFace& scvf) const - { return boundaryVertices_[scvf.vertexIndex()]; } + //! Returns whether or not an scvf is on an interior boundary + bool isOnInteriorBoundary(const SubControlVolumeFace& scvf) const + { return enableInteriorBoundaries ? interiorBoundaryScvfs_[scvf.index()] : false; } + + //! Returns whether or not an scvf touches the domain boundary + bool touchesDomainBoundary(const SubControlVolumeFace& scvf) const + { return domainBoundaryVertices_[scvf.vertexIndex()]; } + + //! Returns whether or not an scvf touches the domain boundary + bool touchesInteriorBoundary(const SubControlVolumeFace& scvf) const + { return enableInteriorBoundaries ? interiorBoundaryVertices_[scvf.vertexIndex()] : false; } + + //! Returns whether or not an scvf touches an interior or a domain boundary (has to be called before getting an interaction volume) + bool touchesInteriorOrDomainBoundary(const SubControlVolumeFace& scvf) const + { return touchesDomainBoundary(scvf) || touchesInteriorBoundary(scvf); } //! Returns whether or not an scvf touches a branching point (for dim < dimWorld) - bool scvfTouchesBranchingPoint(const SubControlVolumeFace& scvf) const - { - if (dim == dimWorld) - return false; - else - return branchingVertices_[scvf.vertexIndex()]; - } + bool touchesBranchingPoint(const SubControlVolumeFace& scvf) const + { return dim == dimWorld ? false : branchingVertices_[scvf.vertexIndex()]; } //! update all fvElementGeometries (do this again after grid adaption) void update(const Problem& problem) { problemPtr_ = &problem; - // clear containers (necessary after grid refinement) - scvs_.clear(); - scvfs_.clear(); - scvfIndicesOfScv_.clear(); - boundaryVertices_.clear(); - elementMap_.clear(); + // resize containers + const auto numVert = gridView_.size(dim); + const auto numScvs = gridView_.size(0); + const auto numScvfs = MpfaHelper::getGlobalNumScvf(gridView_); - // reserve memory - std::size_t numScvs = gridView_.size(0); - std::size_t numScvfs = MpfaHelper::getGlobalNumScvf(gridView_); scvs_.resize(numScvs); scvfs_.reserve(numScvfs); scvfIndicesOfScv_.resize(numScvs); - boundaryVertices_.resize(gridView_.size(dim), false); elementMap_.resize(numScvs); - // keep track of branching points + // Keep track of domain (and mabybe interior) boundaries + domainBoundaryVertices_.resize(numVert, false); + std::vector interiorOrDomainBoundaryVertices(numVert, false); + if (enableInteriorBoundaries) + { + interiorBoundaryVertices_.resize(numVert, false); + interiorBoundaryScvfs_.resize(numScvfs, false); + } + + // Maybe keep track of branching points if (dim < dimWorld) branchingVertices_.resize(gridView_.size(dim), false); // find vertices on processor boundaries - auto isGhostVertex = MpfaHelper::findGhostVertices(problem, gridView_); - - // the quadrature point to be used on the scvf - const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q); + const auto isGhostVertex = MpfaHelper::findGhostVertices(problem, gridView_); - // Instantiate a helper class to pass it to the scvf constructors - MpfaHelper helper; + // the quadrature point to be used on the scvfs + static const Scalar q = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Q); // Build the SCVs and SCV faces IndexType scvfIdx = 0; - numBoundaryScvf_ = 0; - numBoundaryVertices_ = 0; + numDomainBoundaryScvf_ = 0; + numInteriorBoundaryScvf_ = 0; + numBranchingPointScvf_ = 0; + numBranchingPointVertices_ = 0; + numInteriorOrDomainBoundaryVertices_ = 0; for (const auto& element : elements(gridView_)) { - auto eIdx = problem.elementMapper().index(element); + const auto eIdx = problem.elementMapper().index(element); // the element geometry auto elementGeometry = element.geometry(); @@ -200,25 +223,26 @@ public: // construct the sub control volume faces for (const auto& is : intersections(gridView_, element)) { - auto indexInInside = is.indexInInside(); - bool boundary = is.boundary(); - bool neighbor = is.neighbor(); + const auto indexInInside = is.indexInInside(); + const bool boundary = is.boundary(); + const bool neighbor = is.neighbor(); + const bool interiorBoundary = enableInteriorBoundaries ? problem.isInteriorBoundary(element, is) : false; // for network grids, skip the rest if handled already if (dim < dimWorld && neighbor && outsideIndices[indexInInside].empty()) continue; // determine the outside volvar indices - std::vector nIndices; - if (neighbor && dim == dimWorld) - nIndices = std::vector( {problem.elementMapper().index(is.outside())} ); + const std::vector nIndices = neighbor && dim == dimWorld ? + std::vector({problem.elementMapper().index(is.outside())}) : + std::vector(); // get the intersection corners according to generic numbering - auto numCorners = is.geometry().corners(); + const auto numCorners = is.geometry().corners(); std::vector isCorners(numCorners); // if outside level > inside level, use the outside element in the following - bool useNeighbor = neighbor && is.outside().level() > element.level(); + const bool useNeighbor = neighbor && is.outside().level() > element.level(); const auto& e = useNeighbor ? is.outside() : element; const auto indexInElement = useNeighbor ? is.indexInOutside() : indexInInside; const auto eg = e.geometry(); @@ -231,31 +255,29 @@ public: for (unsigned int c = 0; c < numCorners; ++c) { // get the global vertex index the scv face is connected to - auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim); - auto vIdxGlobal = problem.vertexMapper().subIndex(e, vIdxLocal, dim); + const auto vIdxLocal = refElement.subEntity(indexInElement, 1, c, dim); + const auto vIdxGlobal = problem.vertexMapper().subIndex(e, vIdxLocal, dim); // do not build scvfs connected to a processor boundary if (isGhostVertex[vIdxGlobal]) continue; - // store info on which vertices are on the domain boundary - if (boundary && !boundaryVertices_[vIdxGlobal]) - { - boundaryVertices_[vIdxGlobal] = true; - numBoundaryVertices_++; - } - // is vertex on a branching point? if (dim < dimWorld && outsideIndices[indexInInside].size() > 1) + { + if (!branchingVertices_[vIdxGlobal]) + numBranchingPointVertices_++; branchingVertices_[vIdxGlobal] = true; + numBranchingPointScvf_++; + } - // make the scv face (for inside scvfs on network grids, use precalculated outside indices) + // make the scv face (for non-boundary scvfs on network grids, use precalculated outside indices) if (!boundary) { const auto& outsideScvIndices = dim == dimWorld ? nIndices : outsideIndices[indexInInside]; scvfIndexSet.push_back(scvfIdx); - scvfs_.emplace_back(helper, - helper.getScvfCorners(isCorners, c), + scvfs_.emplace_back(MpfaHelper(), + MpfaHelper::getScvfCorners(isCorners, c), is.centerUnitOuterNormal(), vIdxGlobal, vIdxLocal, @@ -268,22 +290,42 @@ public: } else { - nIndices.resize(1); - nIndices[0] = numScvs + numBoundaryScvf_++; + const std::vector boundaryIdx = [&] () + { + IndexType bIdx = numScvs + numDomainBoundaryScvf_++; + return std::vector({bIdx}); + } (); scvfIndexSet.push_back(scvfIdx); - scvfs_.emplace_back(helper, - helper.getScvfCorners(isCorners, c), + scvfs_.emplace_back(MpfaHelper(), + MpfaHelper::getScvfCorners(isCorners, c), is.centerUnitOuterNormal(), vIdxGlobal, vIdxLocal, scvfIdx, eIdx, - nIndices, + boundaryIdx, q, boundary ); } + // if a new interior or domain boundary has been found, increase counter + if ((boundary || interiorBoundary) && !interiorOrDomainBoundaryVertices[vIdxGlobal]) + { + interiorOrDomainBoundaryVertices[vIdxGlobal] = true; + numInteriorOrDomainBoundaryVertices_++; + } + + // store info on which vertices are on the domain/interior boundary + if (boundary) + domainBoundaryVertices_[vIdxGlobal] = true; + if (enableInteriorBoundaries && interiorBoundary) + { + interiorBoundaryVertices_[vIdxGlobal] = true; + interiorBoundaryScvfs_[scvfIdx] = true; + numInteriorBoundaryScvf_++; + } + // increment scvf counter scvfIdx++; } @@ -333,8 +375,11 @@ public: } } + // make sure we found as many scvfs as previously estimated + assert(scvfIdx == numScvfs); + // Initialize the interaction volume seeds - globalInteractionVolumeSeeds_.update(problem, boundaryVertices_); + globalInteractionVolumeSeeds_.update(problem, interiorOrDomainBoundaryVertices); } /*! @@ -379,10 +424,15 @@ private: // vectors that store the global data std::vector> scvfIndicesOfScv_; - std::vector boundaryVertices_; + std::vector domainBoundaryVertices_; + std::vector interiorBoundaryScvfs_; + std::vector interiorBoundaryVertices_; std::vector branchingVertices_; - std::size_t numBoundaryScvf_; - std::size_t numBoundaryVertices_; + std::size_t numInteriorOrDomainBoundaryVertices_; + std::size_t numDomainBoundaryScvf_; + std::size_t numInteriorBoundaryScvf_; + std::size_t numBranchingPointScvf_; + std::size_t numBranchingPointVertices_; // needed for embedded surface and network grids (dim < dimWorld) std::vector> flipScvfIndices_; @@ -416,12 +466,14 @@ class CCMpfaGlobalFVGeometryBase using IndexType = typename GridView::IndexSet::IndexType; using LocalIndexType = typename InteractionVolume::LocalIndexType; - static const int dim = GridView::dimension; - static const int dimWorld = GridView::dimensionworld; + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; using GlobalPosition = Dune::FieldVector; using ReferenceElements = typename Dune::ReferenceElements; + static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries); + public: //! Constructor CCMpfaGlobalFVGeometryBase(const GridView gridView) @@ -435,13 +487,25 @@ public: std::size_t numScvf() const { return numScvf_; } - //! The total number of boundary sub control volume faces - std::size_t numBoundaryScvf() const - { return numBoundaryScvf_; } + //! The total number of scvfs on the interior boundaries + std::size_t numDomainBoundaryScvf() const + { return numDomainBoundaryScvf_; } + + //! The total number of scvfs on the interior boundaries + std::size_t numInteriorBoundaryScvf() const + { return numInteriorBoundaryScvf_; } + + //! The total number of scvfs connected to branching points + std::size_t numBranchingPointScvf() const + { return numBranchingPointScvf_; } + + //! The total number of vertices connected to branching points + std::size_t numBranchingPointVertices() const + { return numBranchingPointVertices_; } //! The total number of vertices on the boundary - std::size_t numBoundaryVertices() const - { return numBoundaryVertices_; } + std::size_t numInteriorOrDomainBoundaryVertices() const + { return numInteriorOrDomainBoundaryVertices_; } // Get an element from a sub control volume contained in it Element element(const SubControlVolume& scv) const @@ -463,38 +527,52 @@ public: const BoundaryInteractionVolumeSeed& boundaryInteractionVolumeSeed(const SubControlVolumeFace& scvf) const { return globalInteractionVolumeSeeds_.boundarySeed(scvf); } - //! Returns whether or not an scvf touches the boundary (has to be called before getting an interaction volume) - bool scvfTouchesBoundary(const SubControlVolumeFace& scvf) const - { return boundaryVertices_[scvf.vertexIndex()]; } + //! Returns whether or not an scvf is on an interior boundary + bool isOnInteriorBoundary(const SubControlVolumeFace& scvf) const + { return enableInteriorBoundaries ? interiorBoundaryScvfs_[scvf.index()] : false; } + + //! Returns whether or not an scvf touches the domain boundary + bool touchesDomainBoundary(const SubControlVolumeFace& scvf) const + { return domainBoundaryVertices_[scvf.vertexIndex()]; } + + //! Returns whether or not an scvf touches the domain boundary + bool touchesInteriorBoundary(const SubControlVolumeFace& scvf) const + { return enableInteriorBoundaries ? interiorBoundaryVertices_[scvf.vertexIndex()] : false; } + + //! Returns whether or not an scvf touches an interior or a domain boundary (has to be called before getting an interaction volume) + bool touchesInteriorOrDomainBoundary(const SubControlVolumeFace& scvf) const + { return touchesDomainBoundary(scvf) || touchesInteriorBoundary(scvf); } //! Returns whether or not a vertex is on a processor boundary bool isGhostVertex(const IndexType vIdxGlobal) const { return ghostVertices_[vIdxGlobal]; } //! Returns whether or not an scvf touches a branching point (for dim < dimWorld) - bool scvfTouchesBranchingPoint(const SubControlVolumeFace& scvf) const - { return branchingVertices_[scvf.vertexIndex()]; } + bool touchesBranchingPoint(const SubControlVolumeFace& scvf) const + { return dim < dimWorld ? branchingVertices_[scvf.vertexIndex()] : false; } //! update all fvElementGeometries (do this again after grid adaption) void update(const Problem& problem) { problemPtr_ = &problem; - // clear containers (necessary after grid refinement) - scvfIndicesOfScv_.clear(); - neighborVolVarIndices_.clear(); - boundaryVertices_.clear(); - elementMap_.clear(); - - // reserve memory or resize the containers + // resize the containers + const auto numScvfs = MpfaHelper::getGlobalNumScvf(gridView_); + const auto numVert = gridView_.size(dim); numScvs_ = gridView_.size(0); - numScvf_ = 0; - numBoundaryScvf_ = 0; - numBoundaryVertices_ = 0; + elementMap_.resize(numScvs_); scvfIndicesOfScv_.resize(numScvs_); neighborVolVarIndices_.resize(numScvs_); - boundaryVertices_.resize(gridView_.size(dim), false); + + // keep track of the domain (and maybe interior) boundaries + domainBoundaryVertices_.resize(numVert, false); + std::vector interiorOrDomainBoundaryVertices(numVert, false); + if (enableInteriorBoundaries) + { + interiorBoundaryScvfs_.resize(numScvfs, false); + interiorBoundaryVertices_.resize(numVert, false); + } // keep track of branching points if (dim < dimWorld) @@ -504,9 +582,16 @@ public: ghostVertices_ = MpfaHelper::findGhostVertices(problem, gridView_); // Store necessary info on SCVs and SCV faces + // reset counters for the tracking of the indices + numScvf_ = 0; + numDomainBoundaryScvf_ = 0; + numInteriorBoundaryScvf_ = 0; + numBranchingPointScvf_ = 0; + numBranchingPointVertices_ = 0; + numInteriorOrDomainBoundaryVertices_ = 0; for (const auto& element : elements(gridView_)) { - auto eIdx = problem.elementMapper().index(element); + const auto eIdx = problem.elementMapper().index(element); // fill the element map with seeds elementMap_[eIdx] = element.seed(); @@ -515,7 +600,7 @@ public: auto eg = element.geometry(); // the element-wise index sets for finite volume geometry - auto numLocalFaces = MpfaHelper::getNumLocalScvfs(eg.type()); + const auto numLocalFaces = MpfaHelper::getNumLocalScvfs(eg.type()); std::vector scvfsIndexSet; std::vector< std::vector > neighborVolVarIndexSet; scvfsIndexSet.reserve(numLocalFaces); @@ -542,20 +627,22 @@ public: // construct the sub control volume faces for (const auto& is : intersections(gridView_, element)) { - auto indexInInside = is.indexInInside(); - bool boundary = is.boundary(); - bool neighbor = is.neighbor(); + const auto indexInInside = is.indexInInside(); + const bool boundary = is.boundary(); + const bool neighbor = is.neighbor(); + const bool interiorBoundary = enableInteriorBoundaries ? problem.isInteriorBoundary(element, is) : false; + // for network grids, skip the rest if handled already if (dim < dimWorld && neighbor && outsideIndices[indexInInside].empty()) continue; // determine the outside volvar indices - std::vector nIndices; - if (neighbor && dim == dimWorld) - nIndices = std::vector( {problem.elementMapper().index(is.outside())} ); + const std::vector nIndices = neighbor && dim == dimWorld ? + std::vector({problem.elementMapper().index(is.outside())}) : + std::vector(); // if outside level > inside level, use the outside element in the following - bool useNeighbor = neighbor && is.outside().level() > element.level(); + const bool useNeighbor = neighbor && is.outside().level() > element.level(); const auto& e = useNeighbor ? is.outside() : element; const auto indexInElement = useNeighbor ? is.indexInOutside() : indexInInside; const auto eg = e.geometry(); @@ -572,16 +659,14 @@ public: if (ghostVertices_[vIdxGlobal]) continue; - // store info on which vertices are on the domain boundary - if (boundary && !boundaryVertices_[vIdxGlobal]) - { - boundaryVertices_[vIdxGlobal] = true; - numBoundaryVertices_++; - } - // is vertex on a branching point? if (dim < dimWorld && outsideIndices[indexInInside].size() > 1) + { + if (!branchingVertices_[vIdxGlobal]) + numBranchingPointVertices_++; branchingVertices_[vIdxGlobal] = true; + numBranchingPointScvf_++; + } // store information on the scv face (for inner scvfs on network grids use precalculated outside indices) if (!boundary) @@ -592,10 +677,30 @@ public: } else { - nIndices.resize(1); - nIndices[0] = numScvs_ + numBoundaryScvf_++; + const std::vector boundaryIdx = [&] () + { + IndexType bIdx = numScvs_ + numDomainBoundaryScvf_++; + return std::vector({bIdx}); + } (); scvfsIndexSet.push_back(numScvf_++); - neighborVolVarIndexSet.push_back(nIndices); + neighborVolVarIndexSet.push_back(boundaryIdx); + } + + // if a new interior or domain boundary has been found, increase counter + if ((boundary || interiorBoundary) && !interiorOrDomainBoundaryVertices[vIdxGlobal]) + { + interiorOrDomainBoundaryVertices[vIdxGlobal] = true; + numInteriorOrDomainBoundaryVertices_++; + } + + // store info on which vertices are on the domain/interior boundary + if (boundary) + domainBoundaryVertices_[vIdxGlobal] = true; + if (enableInteriorBoundaries && interiorBoundary) + { + interiorBoundaryVertices_[vIdxGlobal] = true; + interiorBoundaryScvfs_[numScvf_-1] = true; + numInteriorBoundaryScvf_++; } // increment counter @@ -612,8 +717,11 @@ public: neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet; } + // make sure we found as many scvfs as previously estimated + assert(numScvf_ == numScvfs); + // Initialize the interaction volume seeds - globalInteractionVolumeSeeds_.update(problem, boundaryVertices_); + globalInteractionVolumeSeeds_.update(problem, interiorOrDomainBoundaryVertices); } const std::vector& scvfIndicesOfScv(IndexType scvIdx) const @@ -642,14 +750,19 @@ private: // Information on the global number of geometries std::size_t numScvs_; std::size_t numScvf_; - std::size_t numBoundaryScvf_; - std::size_t numBoundaryVertices_; + std::size_t numDomainBoundaryScvf_; + std::size_t numInteriorBoundaryScvf_; + std::size_t numBranchingPointScvf_; + std::size_t numBranchingPointVertices_; + std::size_t numInteriorOrDomainBoundaryVertices_; // vectors that store the global data Dumux::ElementMap elementMap_; std::vector> scvfIndicesOfScv_; std::vector< std::vector< std::vector > > neighborVolVarIndices_; - std::vector boundaryVertices_; + std::vector domainBoundaryVertices_; + std::vector interiorBoundaryScvfs_; + std::vector interiorBoundaryVertices_; std::vector ghostVertices_; std::vector branchingVertices_; diff --git a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh index cadbd608266eda58b1d1f3acafdcf1d5480d21c0..06b4056d9214b291fb8e5ed619e790cdd29bda88 100644 --- a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh +++ b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh @@ -52,12 +52,12 @@ public: CCMpfaGlobalInteractionVolumeSeedsBase(const GridView& gridView) : gridView_(gridView) {} // initializes the interaction volumes or the seeds - void update(const Problem& p, const std::vector& boundaryVertices) + void update(const Problem& p, const std::vector& interiorOrDomainBoundaryVertices) { problemPtr_ = &p; // initialize the seeds according to the mpfa method - asImp_().initializeSeeds(boundaryVertices, + asImp_().initializeSeeds(interiorOrDomainBoundaryVertices, scvfIndexMap_, seeds_, boundarySeeds_); diff --git a/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh index 3869317a56738415b409db51fb470d891142b6bd..54ce225e4df8b0d9de67b144fc988800300b4fe5 100644 --- a/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh @@ -68,9 +68,9 @@ public: problemPtr_ = &problem; auto numScv = problem.model().globalFvGeometry().numScv(); - auto numBoundaryScvf = problem.model().globalFvGeometry().numBoundaryScvf(); + auto numDomainBoundaryScvf = problem.model().globalFvGeometry().numDomainBoundaryScvf(); - volumeVariables_.resize(numScv + numBoundaryScvf); + volumeVariables_.resize(numScv + numDomainBoundaryScvf); for (const auto& element : elements(problem.gridView())) { auto fvGeometry = localView(problem.model().globalFvGeometry()); diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh index 23fe038d1acbd21db62ddd97bf80e91a38c45050..5f27a14dbd4466d2c5de4c93825e4d22f1f9f13c 100644 --- a/dumux/discretization/cellcentered/mpfa/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/helper.hh @@ -445,6 +445,7 @@ class CCMpfaHelperImplementation : public MpfaDimensionHelper; + using DimWorldMatrix = Dune::FieldMatrix; using ScvfVector = std::array; using LocalBasis = std::array; using CoordScalar = typename GridView::ctype; using ReferenceElements = typename Dune::ReferenceElements; + static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries); + public: // returns shared pointers to the scv faces that share a vertex in the order of a right hand system static ScvfVector getScvFacesAtVertex(const GlobalIndexType vIdxGlobal, @@ -522,22 +526,51 @@ public: const Element& element, const SubControlVolumeFace& scvf) { - if (!scvf.boundary()) - return MpfaFaceTypes::interior; + // We do simplified checks if interior boundaries are disabled + if (!enableInteriorBoundaries) + { + if (!scvf.boundary()) + return MpfaFaceTypes::interior; - auto bcTypes = problem.boundaryTypes(element, scvf); - if (bcTypes.hasOnlyNeumann()) - return MpfaFaceTypes::neumann; - if (bcTypes.hasOnlyDirichlet()) - return MpfaFaceTypes::dirichlet; + const auto bcTypes = problem.boundaryTypes(element, scvf); + if (bcTypes.hasOnlyNeumann()) + return MpfaFaceTypes::neumann; + if (bcTypes.hasOnlyDirichlet()) + return MpfaFaceTypes::dirichlet; - // throw for outflow or mixed boundary conditions - if (bcTypes.hasOutflow()) - DUNE_THROW(Dune::NotImplemented, "outflow BC for mpfa schemes"); - if (bcTypes.hasDirichlet() && bcTypes.hasNeumann()) - DUNE_THROW(Dune::InvalidStateException, "Mixed BC are not allowed for cellcentered schemes"); + // throw exception + return throwBoundaryExceptionMessage(bcTypes); + } + else + { + const auto bcTypes = problem.boundaryTypes(element, scvf); - DUNE_THROW(Dune::InvalidStateException, "unknown boundary condition type"); + // if we are on an interior boundary return interior types + if (problem.model().globalFvGeometry().isOnInteriorBoundary(scvf)) + { + if (bcTypes.hasOnlyNeumann()) + return MpfaFaceTypes::interiorNeumann; + if (bcTypes.hasOnlyDirichlet()) + return MpfaFaceTypes::interiorDirichlet; + + // throw exception + return throwBoundaryExceptionMessage(bcTypes); + } + + if (scvf.boundary()) + { + if (bcTypes.hasOnlyNeumann()) + return MpfaFaceTypes::neumann; + if (bcTypes.hasOnlyDirichlet()) + return MpfaFaceTypes::dirichlet; + + // throw exception + return throwBoundaryExceptionMessage(bcTypes); + } + + // This is an interior scvf + return MpfaFaceTypes::interior; + } } // returns a vector, which maps a bool (true if ghost vertex) to each vertex index @@ -582,6 +615,36 @@ public: template static bool contains(const std::vector& vector, const V2 value) { return std::find(vector.begin(), vector.end(), value) != vector.end(); } + + // calculates the product of a transposed vector n, a Matrix M and another vector v (n^T M v) + static Scalar nT_M_v(const GlobalPosition& n, + const DimWorldMatrix& M, + const GlobalPosition& v) + { + GlobalPosition tmp; + M.mv(v, tmp); + return n*tmp; + } + + // calculates the product of a transposed vector n, a Scalar M and another vector v (n^T M v) + static Scalar nT_M_v(const GlobalPosition& n, + const Scalar M, + const GlobalPosition& v) + { + return M*(n*v); + } + +private: + static MpfaFaceTypes throwBoundaryExceptionMessage(const BoundaryTypes& bcTypes) + { + // throw for outflow or mixed boundary conditions + if (bcTypes.hasOutflow()) + DUNE_THROW(Dune::NotImplemented, "outflow BC for mpfa schemes"); + if (bcTypes.hasDirichlet() && bcTypes.hasNeumann()) + DUNE_THROW(Dune::InvalidStateException, "Mixed BC are not allowed for cellcentered schemes"); + + DUNE_THROW(Dune::InvalidStateException, "unknown boundary condition type"); + } }; /*! diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/interactionvolume.hh index 8d8f233b95d10045126e4b55e0c4542659f85c19..729fe9ed949c14f15a2f3729836ef262f7b4e05c 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolume.hh @@ -29,8 +29,7 @@ namespace Dumux { // forward declaration of the base class template -class CCMpfaInteractionVolumeImplementation -{}; +class CCMpfaInteractionVolumeImplementation; /*! * \ingroup Mpfa diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh index ea3ced18cfc90177fd6445e5f77732c45f3806dd..af8c2d793f151997880065581d45dc71462047d7 100644 --- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh +++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh @@ -79,6 +79,7 @@ public: LocalIndexType localScvIndex; bool isOutside; + //! Constructor fully initializing the members LocalFaceData(const LocalIndexType faceIndex, const LocalIndexType scvIndex, bool isOut) @@ -113,12 +114,12 @@ public: { DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getTransmissibilities() method."); } //! returns the neumann flux corresponding to a local scvf - Scalar getNeumannFlux(const LocalFaceData& localFaceData) const + Scalar getNeumannFlux(const LocalFaceData& localFaceData, unsigned int eqIdx) const { DUNE_THROW(Dune::NotImplemented, "Actual interaction volume implementation does not provide a getNeumannFlux() method."); } //! returns the local index in a vector for a given global index template - LocalIndexType findLocalIndex(const std::vector& vector, const IdxType2 globalIdx) const + LocalIndexType findIndexInVector(const std::vector& vector, const IdxType2 globalIdx) const { auto it = std::find(vector.begin(), vector.end(), globalIdx); assert(it != vector.end() && "could not find local index in the vector for the given global index!"); diff --git a/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh b/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh new file mode 100644 index 0000000000000000000000000000000000000000..c1cc25b6b088e56bd0b6191f481f54537418ace9 --- /dev/null +++ b/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh @@ -0,0 +1,166 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * \brief A class to store info on interior boundaries + */ +#ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERIORBOUNDARYDATA_HH +#define DUMUX_DISCRETIZATION_CC_MPFA_INTERIORBOUNDARYDATA_HH + +#include + +namespace Dumux +{ + +//! forward declaration of property tags +namespace Properties +{ + NEW_PROP_TAG(SpatialParams); + NEW_PROP_TAG(FacetCoupling); +} + +template +class InteriorBoundaryData +{ + 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 SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams); + using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); + using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume); + + using IndexType = typename GridView::IndexSet::IndexType; + using LocalIndexType = typename BoundaryInteractionVolume::LocalIndexType; + + //! Dummy type for the CompleteCoupledFacetData struct. + //! Implementations need to have at least the provided interfaces. + //! Note that the return types are also "wrong" (Here just to satisfy the compiler) + struct CompleteCoupledFacetData + { + const SpatialParams& spatialParams() const + { + DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided"); + } + + const VolumeVariables volVars() const + { + DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided"); + } + + const Element element() const + { + DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided"); + } + + const FVElementGeometry fvGeometry() const + { + DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided"); + } + + const FVElementGeometry scv() const + { + DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided"); + } + }; + +public: + //! the constructor + InteriorBoundaryData(const Problem& problem, + IndexType elementIndex, + IndexType scvfIndex, + LocalIndexType localIndex, + MpfaFaceTypes faceType) + : problemPtr_(&problem), + elementIndex_(elementIndex), + scvfIndex_(scvfIndex), + localIndex_(localIndex), + faceType_(faceType) + {} + + //! returns the global index of the element/scv connected to the interior boundary + IndexType elementIndex() const + { return elementIndex_; } + + //! returns the global index of the scvf connected to the interior boundary + IndexType scvfIndex() const + { return scvfIndex_; } + + //! returns the local index i of the scvf within the interaction volume. + //! This is either: + //! - the i-th flux face index (interior neumann boundaries) + //! - the i-th interior dirichlet face (interior dirichlet boundaries) + LocalIndexType localIndexInInteractionVolume() const + { return localIndex_; } + + //! returns the face type of this scvf + MpfaFaceTypes faceType() const + { return faceType_; } + + //! returns the volume variables for interior dirichlet boundaries + VolumeVariables facetVolVars(const FVElementGeometry& fvGeometry) const + { + //! This cannot be called when FacetCoupling is active + assert(!GET_PROP_VALUE(TypeTag, MpfaFacetCoupling) && "For models with a coupled problem on the element facets you have to" + "provide a suitable implementation of the InteriorBoundaryData class"); + + //! This can only be called for interior Dirichlet boundaries + assert(faceType_ == MpfaFaceType::interiorDirichlet && "requesting Dirichlet vol vars for a face which is" + "not marked as interior Dirichlet face."); + + auto element = problem_().model().globalFvGeometry().element(elementIndex()); + auto priVars = problem_().dirichlet(element, fvGeometry.scvf(scvfIndex())); + + VolumeVariables volVars; + volVars.update(ElementSolutionVector({priVars}), + problem_(), + element, + fvGeometry.scv(elementIndex())); + + return volVars; + } + + //! The following interface is to be overloaded for problems using facet coupling. + //! prepares all the necessary variables of the other domain. + //! Note that also an implementation of the CompleteFacetData structure has to be provided. + CompleteCoupledFacetData completeCoupledFacetData() const + { + if (GET_PROP_VALUE(TypeTag, MpfaFacetCoupling)) + DUNE_THROW(Dune::InvalidStateException, "You have to use an InteriorBoundaryData class designed " + "for handling a coupled problem on the element facets."); + else + DUNE_THROW(Dune::InvalidStateException, "FacetCoupling is not active. " + "Calling completeCoupledFacetData() for uncoupled problems is invalid."); + } + +private: + const Problem& problem_() const + { return *problemPtr_; } + + const Problem* problemPtr_; + + IndexType elementIndex_; + IndexType scvfIndex_; + LocalIndexType localIndex_; + MpfaFaceTypes faceType_; +}; +} // end namespace + +#endif diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh index 7b36d027cd6cf5e2f25d40b30bc6ed2449f1bb9e..8a3da8af4f4bdc17d2e75722b525f3dc270592e9 100644 --- a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh +++ b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh @@ -52,7 +52,7 @@ public: CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView& gridView) : ParentType(gridView) {} template - void initializeSeeds(const std::vector& boundaryVertices, + void initializeSeeds(const std::vector& interiorOrDomainBoundaryVertices, std::vector& scvfIndexMap, SeedVector& seeds, BoundarySeedVector& boundarySeeds) @@ -63,12 +63,17 @@ public: // reserve memory const auto numScvf = this->problem().model().globalFvGeometry().numScvf(); - const auto numBoundaryScvf = this->problem().model().globalFvGeometry().numBoundaryScvf(); - const int numInteriorScvf = (numScvf-numBoundaryScvf)/2; + const auto numIBScvf = this->problem().model().globalFvGeometry().numInteriorBoundaryScvf(); + const auto numDBScvf = this->problem().model().globalFvGeometry().numDomainBoundaryScvf(); + const auto numBPScvf = this->problem().model().globalFvGeometry().numBranchingPointScvf(); - if (numInteriorScvf > 0) - seeds.reserve( numInteriorScvf ); - boundarySeeds.reserve(numBoundaryScvf); + const int numLMethodIVs = (numScvf-numIBScvf-numDBScvf-numBPScvf)/2; + const auto numOMethodIVs = this->problem().model().globalFvGeometry().numInteriorOrDomainBoundaryVertices() + + this->problem().model().globalFvGeometry().numBranchingPointVertices(); + + if (numLMethodIVs > 0) + seeds.reserve( numLMethodIVs ); + boundarySeeds.reserve(numOMethodIVs); scvfIndexMap.resize(numScvf); // Keep track of which faces have been handled already @@ -87,7 +92,12 @@ public: if (isFaceHandled[scvf.index()]) continue; - if (boundaryVertices[scvf.vertexIndex()]) + // make boundary interaction volume seeds (o-method seeds) on: + // - interior boundaries + // - domain boundaries + // - branching points (l-method can not handle this) + if (interiorOrDomainBoundaryVertices[scvf.vertexIndex()] + || this->problem().model().globalFvGeometry().touchesBranchingPoint(scvf)) { // make the boundary interaction volume seed boundarySeeds.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem(), diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh index f11cb30f2e2313d3415a4c7fdf39015c8d955bc3..f5d88a68972c485f9b8111d617c90e6b3f744612 100644 --- a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh @@ -45,7 +45,7 @@ class MpfaMethodHelper static const int dimWorld = 2; // The mpfa-o helper class used to construct the boundary interaction volume seeds - using oMethodHelper = CCMpfaHelperImplementation; + using oMethodHelper = CCMpfaHelperImplementation; using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); @@ -64,7 +64,6 @@ class MpfaMethodHelper using GlobalIndexSet = typename InteractionVolume::GlobalIndexSet; using GlobalIndexType = typename InteractionVolume::GlobalIndexType; - using LocalIndexSet = typename InteractionVolume::LocalIndexSet; using LocalIndexType = typename InteractionVolume::LocalIndexType; using Matrix = typename InteractionVolume::Matrix; public: @@ -177,6 +176,70 @@ private: } }; +/*! + * \ingroup Mpfa + * \brief Helper class to get the required information on an interaction volume. + * Specialization for the Mpfa-L method for 2d embedded in 3d space. + */ +template +class MpfaMethodHelper +{ + using Implementation = typename GET_PROP_TYPE(TypeTag, MpfaHelper); + + static const int dim = 2; + static const int dimWorld = 3; + + // The mpfa-o helper class used to construct the boundary interaction volume seeds + using oMethodHelper = CCMpfaHelperImplementation; + + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + 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 Element = typename GridView::template Codim<0>::Entity; + + using InteractionVolumeSeed = typename InteractionVolume::Seed; + using BoundaryInteractionVolumeSeed = typename BoundaryInteractionVolume::Seed; + + using LocalIndexType = typename InteractionVolume::LocalIndexType; + using Matrix = typename InteractionVolume::Matrix; +public: + + //! We know that this is only called for scvfs NOT touching: + //! - branching points + //! - domain boundaries + //! - interior boundaries + //! Thus, here we can use the same algorithm as in the 2d case. + static InteractionVolumeSeed makeInnerInteractionVolumeSeed(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) + { + return MpfaMethodHelper::makeInnerInteractionVolumeSeed(problem, + element, + fvGeometry, + scvf); + } + + //! Here we simply forward to the o-method helper (we use o-interaction volumes on the boundary) + static BoundaryInteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) + { return oMethodHelper::makeBoundaryInteractionVolumeSeed(problem, element, fvGeometry, scvf); } + + //! The selection criterion is the same as in the dim = 2 = dimWorld case + template + static LocalIndexType selectionCriterion(const InteractionRegion& I1, + const InteractionRegion& I2, + const Matrix& M1, + const Matrix& M2) + { return MpfaMethodHelper::selectionCriterion(I1, I2, M1, M2); } +}; + } // end namespace #endif diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh index 2a9a4287fcee3d73ca46574cabfda23419ffd1f4..fc4c039ea1d955f95c7aa93de31bbc978e0217dd 100644 --- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh +++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh @@ -36,6 +36,7 @@ private: 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 FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume); @@ -46,6 +47,7 @@ private: using LocalIndexType = typename InteractionVolume::LocalIndexType; using GlobalIndexType = typename InteractionVolume::GlobalIndexType; using GlobalPosition = Dune::FieldVector; + using LocalBasis = std::array; public: LocalIndexType contiFaceLocalIdx; @@ -58,7 +60,7 @@ public: std::array globalScvfs; - // Constructor for dim == 2 + // Constructor signature for dim == 2 template InteractionRegion(const Problem& problem, const FVElementGeometry& fvGeometry, @@ -107,21 +109,36 @@ public: scvIndices[1] = scv2.index(); scvIndices[2] = scv3.index(); - // calculate nus and detXs - static const Dune::FieldMatrix R = {{0.0, 1.0}, {-1.0, 0.0}}; + // set up the local basis of the elements + LocalBasis basis1, basis2, basis3, basis4; + basis1[0] = f1 - scvCenters[0]; + basis1[1] = f2 - scvCenters[0]; + basis2[0] = v - scvCenters[1]; + basis2[1] = f1 - scvCenters[1]; + basis3[0] = f2 - scvCenters[2]; + basis3[1] = v - scvCenters[2]; + basis4[0] = basis1[0]; + basis4[1] = v - scvCenters[0]; + + // calculate nus + const auto nus1 = MpfaHelper::calculateInnerNormals(basis1); + const auto nus2 = MpfaHelper::calculateInnerNormals(basis2); + const auto nus3 = MpfaHelper::calculateInnerNormals(basis3); + const auto nus4 = MpfaHelper::calculateInnerNormals(basis1); + nu.resize(7); - R.mv(f2-scvCenters[0], nu[0]); - R.mv(scvCenters[0]-f1, nu[1]); - R.mv(f1-scvCenters[1], nu[2]); - R.mv(scvCenters[1]-v, nu[3]); - R.mv(v-scvCenters[2], nu[4]); - R.mv(scvCenters[2]-f2, nu[5]); - R.mv(v-scvCenters[0], nu[6]); + nu[0] = nus1[0]; + nu[1] = nus1[1]; + nu[2] = nus2[0]; + nu[3] = nus2[1]; + nu[4] = nus3[0]; + nu[5] = nus3[1]; + nu[6] = nus4[0]; detX.resize(3); - detX[0] = (f1-scvCenters[0])*nu[0]; - detX[1] = (v-scvCenters[1])*nu[2]; - detX[2] = (f2-scvCenters[2])*nu[4]; + detX[0] = MpfaHelper::calculateDetX(basis1); + detX[1] = MpfaHelper::calculateDetX(basis2); + detX[2] = MpfaHelper::calculateDetX(basis3); } }; }; //end namespace diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh index 5cf49681045d604dcfd5ef06777e22918b01304b..51e7929705249eb21b8656195fdd125552499394 100644 --- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh @@ -85,9 +85,10 @@ class CCMpfaInteractionVolumeImplementation : pub static const int dimWorld = GridView::dimensionworld; using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = Dune::FieldVector; + using LocalBasis = std::array; using InteractionRegion = Dumux::InteractionRegion; - using Vector = typename Traits::Vector; + using CoefficientVector = typename Traits::Vector; using Tensor = typename Traits::Tensor; public: @@ -151,7 +152,7 @@ public: } //! returns the transmissibilities corresponding to the bound scvf - Vector getTransmissibilities(const LocalFaceData& localFaceData) const + CoefficientVector getTransmissibilities(const LocalFaceData& localFaceData) const { assert(systemSolved_ && "Transmissibilities not calculated yet. You have to call solveLocalSystem() beforehand."); @@ -179,6 +180,15 @@ public: return globalScvfIndices_; } + //! Boundaries will be treated by a different mpfa method (e.g. o method). Thus, on + //! faces in l-method interaction volumes there will never be a Neumann flux contribution. + Scalar getNeumannFlux(const LocalFaceData& localFaceData, unsigned int eqIdx) const + { return 0.0; } + + //! See comment of getNeumannFlux() + CoefficientVector getNeumannFluxTransformationCoefficients(const LocalFaceData& localFaceData) const + { return CoefficientVector(0.0); } + const Matrix& matrix() const { return T_; } @@ -331,7 +341,10 @@ private: calculateXi_(const GlobalPosition& nu1, const GlobalPosition& nu2, const Scalar detX) const - { return crossProduct(nu1, nu2)/detX; } + { + LocalBasis basis({nu1, nu2}); + return MpfaHelper::calculateDetX(basis)/detX; + } const Seed& seed_() const { return seedPtr_; } diff --git a/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh index f009866d9c5c925a0addff22f00d8f5ca32f38bf..22effd4f005ab8e45a6f2aeeb1c598453cd1d7c8 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh @@ -55,7 +55,7 @@ public: CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView& gridView) : ParentType(gridView) {} template - void initializeSeeds(const std::vector& boundaryVertices, + void initializeSeeds(const std::vector& interiorOrDomainBoundaryVertices, std::vector& scvfIndexMap, SeedVector& seeds, BoundarySeedVector& boundarySeeds) @@ -66,12 +66,12 @@ public: // reserve memory const auto numScvf = this->problem().model().globalFvGeometry().numScvf(); - const auto numBoundaryVertices = this->problem().model().globalFvGeometry().numBoundaryVertices(); - const int numInteriorVertices = this->gridView().size(dim) - numBoundaryVertices; + const auto numInteriorOrDomainBoundaryVertices = this->problem().model().globalFvGeometry().numInteriorOrDomainBoundaryVertices(); + const int numInteriorVertices = this->gridView().size(dim) - numInteriorOrDomainBoundaryVertices; if (numInteriorVertices > 0) seeds.reserve(numInteriorVertices); - boundarySeeds.reserve(numBoundaryVertices); + boundarySeeds.reserve(numInteriorOrDomainBoundaryVertices); scvfIndexMap.resize(numScvf); // Keep track of which faces have been handled already @@ -90,7 +90,7 @@ public: if (isFaceHandled[scvf.index()]) continue; - if (boundaryVertices[scvf.vertexIndex()]) + if (interiorOrDomainBoundaryVertices[scvf.vertexIndex()]) { // make the boundary interaction volume seed boundarySeeds.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem(), diff --git a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh index 12ae0c978118b4980f81a7b997219716f2a1378b..9c7aca1b0175be1a7f358f6cdce2ab046d6bba6f 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh @@ -42,8 +42,8 @@ class MpfaMethodHelper { using Implementation = typename GET_PROP_TYPE(TypeTag, MpfaHelper); - static const int dim = 2; - static const int dimWorld = 2; + static constexpr int dim = 2; + static constexpr int dimWorld = 2; using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); @@ -67,6 +67,9 @@ class MpfaMethodHelper using GlobalIndexSet = typename InteractionVolume::GlobalIndexSet; using LocalIndexSet = typename InteractionVolume::LocalIndexSet; + + static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries); + public: static InteractionVolumeSeed makeInnerInteractionVolumeSeed(const Problem& problem, const Element& element, @@ -107,7 +110,9 @@ public: scvSeeds.shrink_to_fit(); scvfSeeds.shrink_to_fit(); - return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true); + // check if touches domain boundary (not only interior boundary) + const bool boundary = problem.model().globalFvGeometry().touchesDomainBoundary(scvf); + return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), boundary); } private: @@ -120,13 +125,13 @@ private: const SubControlVolumeFace& scvf) { // Check whether or not we are touching the boundary here - bool onBoundary = problem.model().globalFvGeometry().scvfTouchesBoundary(scvf); + const bool onBoundary = problem.model().globalFvGeometry().touchesDomainBoundary(scvf); // Get the two scv faces in the first scv - auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry); + const auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry); // The global index of the first scv of the interaction region - auto scvIdx0 = scvf.insideScvIdx(); + const auto scvIdx0 = scvf.insideScvIdx(); // rotate counter clockwise and create the entities performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0); @@ -134,7 +139,7 @@ private: if (onBoundary) { // the local scvf index of the second local face of the first local scv - LocalIndexType storeIdx = scvfSeeds.size(); + const LocalIndexType storeIdx = scvfSeeds.size(); // clockwise rotation until hitting the boundary again performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0, /*clockwise*/true); @@ -170,7 +175,7 @@ private: auto outsideFvGeometry = localView(problem.model().globalFvGeometry()); // Start/continue interaction region construction from the given scv face - LocalIndexType startScvfIdx = clockWise ? 1 : 0; + const LocalIndexType startScvfIdx = clockWise ? 1 : 0; auto curScvf = *scvfVector[startScvfIdx]; bool firstIteration = true; bool finished = false; @@ -178,17 +183,16 @@ private: while (!finished) { // Get some indices beforehand - GlobalIndexType outsideGlobalScvIdx = curScvf.outsideScvIdx(); - LocalIndexType insideLocalScvIdx = firstIteration ? 0 : localScvIdx; + const GlobalIndexType outsideGlobalScvIdx = curScvf.outsideScvIdx(); + const LocalIndexType insideLocalScvIdx = firstIteration ? 0 : localScvIdx; // the current element inside of the scv face - auto insideElement = problem.model().globalFvGeometry().element(curScvf.insideScvIdx()); - auto faceType = Implementation::getMpfaFaceType(problem, insideElement, curScvf); + const auto insideElement = problem.model().globalFvGeometry().element(curScvf.insideScvIdx()); + const auto faceType = Implementation::getMpfaFaceType(problem, insideElement, curScvf); // if the face touches the boundary, create a boundary scvf entity if (curScvf.boundary()) { - assert(faceType == MpfaFaceTypes::neumann || faceType == MpfaFaceTypes::dirichlet); scvfSeeds.emplace_back( curScvf, insideLocalScvIdx, LocalIndexSet(), @@ -208,23 +212,31 @@ private: GlobalIndexSet({scvfVector[1]->index()}), faceType ); + // create duplicate face if it is an interior boundary, + if (faceType != MpfaFaceTypes::interior) + scvfSeeds.emplace_back( *scvfVector[1], + 0, + LocalIndexSet({insideLocalScvIdx}), + GlobalIndexSet({curScvf.index()}), + faceType ); + // rotation loop is finished finished = true; return; } // If we get here, there are outside entities - auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx); + const auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx); outsideFvGeometry.bindElement(outsideElement); // get the two scv faces in the outside element that share the vertex - auto outsideScvfVector = Implementation::getCommonAndNextScvFace(curScvf, outsideFvGeometry, clockWise); - GlobalIndexType commonFaceCoordIdx = clockWise ? 0 : 1; - GlobalIndexType nextFaceCoordIdx = clockWise ? 1 : 0; - auto&& commonScvf = *outsideScvfVector[commonFaceCoordIdx]; - auto&& nextScvf = *outsideScvfVector[nextFaceCoordIdx]; + const auto outsideScvfVector = Implementation::getCommonAndNextScvFace(curScvf, outsideFvGeometry, clockWise); + const GlobalIndexType commonFaceCoordIdx = clockWise ? 0 : 1; + const GlobalIndexType nextFaceCoordIdx = clockWise ? 1 : 0; + const auto& commonScvf = *outsideScvfVector[commonFaceCoordIdx]; + const auto& nextScvf = *outsideScvfVector[nextFaceCoordIdx]; // create local scv face entity of the current scvf - LocalIndexType outsideLocalScvIdx = localScvIdx+1; + const LocalIndexType outsideLocalScvIdx = localScvIdx+1; scvfSeeds.emplace_back( curScvf, insideLocalScvIdx, LocalIndexSet({outsideLocalScvIdx}), @@ -232,6 +244,18 @@ private: faceType ); localScvfIdx++; + // create duplicate face if it is an interior boundary, + if (faceType != MpfaFaceTypes::interior) + { + // face from "outside" to "inside", index sets change + scvfSeeds.emplace_back( commonScvf, + outsideLocalScvIdx, + LocalIndexSet({insideLocalScvIdx}), + GlobalIndexSet({curScvf.index()}), + faceType ); + localScvfIdx++; + } + // create index set storing the two local scvf indices LocalIndexSet localScvfs(2); localScvfs[commonFaceCoordIdx] = localScvfIdx-1; @@ -239,7 +263,9 @@ private: // create "outside" scv GlobalIndexSet globalScvfIndices({outsideScvfVector[0]->index(), outsideScvfVector[1]->index()}); - scvSeeds.emplace_back(std::move(globalScvfIndices), std::move(localScvfs), outsideGlobalScvIdx); + scvSeeds.emplace_back(std::move(globalScvfIndices), + std::move(localScvfs), + outsideGlobalScvIdx); localScvIdx++; // create the next scvf in the following iteration @@ -291,21 +317,24 @@ public: const SubControlVolumeFace& scvf) { // if the scvf does not touch a branching point, use simplified algorithm to create interaction volume seed - if (!problem.model().globalFvGeometry().scvfTouchesBranchingPoint(scvf)) - return MpfaMethodHelper::makeInnerInteractionVolumeSeed(problem, element, fvGeometry, scvf); + if (!problem.model().globalFvGeometry().touchesBranchingPoint(scvf)) + return MpfaMethodHelper::makeInnerInteractionVolumeSeed(problem, + element, + fvGeometry, + scvf); std::vector scvSeeds; std::vector scvfSeeds; // reserve sufficient memory - scvSeeds.reserve(20); - scvfSeeds.reserve(20); + scvSeeds.reserve(50); + scvfSeeds.reserve(50); // The vertex index around which we construct the interaction volume - auto vIdxGlobal = scvf.vertexIndex(); + const auto vIdxGlobal = scvf.vertexIndex(); // Get the two scv faces in the scv - auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry); + const auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry); // fill the entity seed data fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvfVector, vIdxGlobal); @@ -326,21 +355,24 @@ public: const SubControlVolumeFace& scvf) { // if the scvf does not touch a branching point, use simplified algorithm to create interaction volume seed - if (!problem.model().globalFvGeometry().scvfTouchesBranchingPoint(scvf)) - return MpfaMethodHelper::makeBoundaryInteractionVolumeSeed(problem, element, fvGeometry, scvf); + if (!problem.model().globalFvGeometry().touchesBranchingPoint(scvf)) + return MpfaMethodHelper::makeBoundaryInteractionVolumeSeed(problem, + element, + fvGeometry, + scvf); std::vector scvSeeds; std::vector scvfSeeds; // reserve sufficient memory - scvSeeds.reserve(20); - scvfSeeds.reserve(20); + scvSeeds.reserve(50); + scvfSeeds.reserve(50); // The vertex index around which we construct the interaction volume - auto vIdxGlobal = scvf.vertexIndex(); + const auto vIdxGlobal = scvf.vertexIndex(); // Get the two scv faces in the scv - auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry); + const auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry); // fill the entity seed data fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvfVector, vIdxGlobal); @@ -373,11 +405,11 @@ private: // make the scvf seeds for the two scvfs connected to the scv auto& actualScvSeed = scvSeeds.back(); - LocalIndexType actualLocalScvIdx = scvSeeds.size()-1; + const LocalIndexType actualLocalScvIdx = scvSeeds.size()-1; for (int coordDir = 0; coordDir < dim; ++coordDir) { - auto&& actualScvf = *scvfVector[coordDir]; + const auto& actualScvf = *scvfVector[coordDir]; // if scvf is on a boundary, we create the scvfSeed and make no neighbor if (actualScvf.boundary()) @@ -398,15 +430,15 @@ private: for (auto outsideGlobalScvIdx : actualScvf.outsideScvIndices()) { // get outside element, fvgeometry etc. - auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx); + const auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx); auto outsideFvGeometry = localView(problem.model().globalFvGeometry()); outsideFvGeometry.bindElement(outsideElement); // find scvf in outside corresponding to the actual scvf - auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry); - auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector); - auto&& outsideScvf = *outsideScvfVector[commonFaceLocalIdx]; - auto outsideScvfIdx = outsideScvf.index(); + const auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry); + const auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector); + const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx]; + const auto outsideScvfIdx = outsideScvf.index(); // check if the outside scv already exists and get its local index bool outsideScvExists = false; @@ -415,7 +447,8 @@ private: { if (scvSeed.globalIndex() == outsideGlobalScvIdx) { - outsideScvExists = true; break; + outsideScvExists = true; + break; } // keep track of local index outsideLocalScvIdx++; @@ -426,31 +459,36 @@ private: bool insideScvfExists = false; LocalIndexType outsideLocalScvfIdx = 0; LocalIndexType insideLocalScvfIdx = 0; + + // If the outside scvf exists, we need to find out its face type + // We initialize it here to please the compiler + MpfaFaceTypes outsideScvfFaceType = MpfaFaceTypes::interior; + for (auto&& scvfSeed : scvfSeeds) { if (Implementation::contains(actualScvf.outsideScvIndices(), scvfSeed.insideGlobalScvIndex()) && Implementation::contains(scvfSeed.outsideGlobalScvIndices(), actualScvf.insideScvIdx())) + { + outsideScvfFaceType = scvfSeed.faceType(); outsideScvfExists = true; - else - if (!outsideScvfExists) - outsideLocalScvfIdx++; + } + else if (!outsideScvfExists) + outsideLocalScvfIdx++; if (scvfSeed.insideGlobalScvIndex() == actualScvf.insideScvIdx() && scvfSeed.outsideGlobalScvIndices().size() == actualScvf.outsideScvIndices().size() && std::equal(scvfSeed.outsideGlobalScvIndices().begin(), scvfSeed.outsideGlobalScvIndices().end(), actualScvf.outsideScvIndices().begin())) + { insideScvfExists = true; - else - if (!insideScvfExists) - insideLocalScvfIdx++; + } + else if (!insideScvfExists) + insideLocalScvfIdx++; } - // we should never have two local scv faces existing inside and in an outside element - assert(!(outsideScvfExists && insideScvfExists) && "The scv face seems to exist twice!"); - // outside scv has to be created if (!outsideScvExists) { - // The face does not yet exist + // No face exists yet if (!insideScvfExists && !outsideScvfExists) { // set the local scvfIndex of the face that is about to created @@ -461,7 +499,7 @@ private: actualLocalScvIdx, Implementation::getMpfaFaceType(problem, element, actualScvf)); - // pass the actual outside indices to the new scvf seed + // pass the outside index that is about to be created to the new scvf seed scvfSeeds.back().addOutsideData(outsideScvfIdx, static_cast(scvSeeds.size())); } else if (insideScvfExists && !outsideScvfExists) @@ -471,8 +509,24 @@ private: } else if (!insideScvfExists && outsideScvfExists) { - // set the local scvfIndex of the outside face - actualScvSeed.setLocalScvfIndex(coordDir, outsideLocalScvfIdx); + // if outside is an interior boundary, we create a duplicate local scvf here for inside + if (outsideScvfFaceType == MpfaFaceTypes::interiorDirichlet || + outsideScvfFaceType == MpfaFaceTypes::interiorNeumann) + { + // set the local scvfIndex of the face that is about to created + actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size()); + + // create scvf seed + scvfSeeds.emplace_back(actualScvf, + actualLocalScvIdx, + Implementation::getMpfaFaceType(problem, element, actualScvf)); + + // pass the outside index that is about to be created to the new scvf seed + scvfSeeds.back().addOutsideData(outsideScvfIdx, static_cast(scvSeeds.size())); + } + else + // set the local scvfIndex of the outside face + actualScvSeed.setLocalScvfIndex(coordDir, outsideLocalScvfIdx); // pass info on inside to the outside scvf seed scvfSeeds[outsideLocalScvfIdx].addOutsideData(actualScvf.index(), actualLocalScvIdx); @@ -499,8 +553,24 @@ private: } else if (!insideScvfExists && outsideScvfExists) { - // set the local scvfIndex of the found outside local scv face seed - actualScvSeed.setLocalScvfIndex(coordDir, outsideLocalScvfIdx); + // if outside is an interior boundary, we create a duplicate local scvf here for inside + if (outsideScvfFaceType == MpfaFaceTypes::interiorDirichlet || + outsideScvfFaceType == MpfaFaceTypes::interiorNeumann) + { + // set the local scvfIndex of the face that is about to created + actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size()); + + // create scvf seed + scvfSeeds.emplace_back(actualScvf, + actualLocalScvIdx, + Implementation::getMpfaFaceType(problem, element, actualScvf)); + + // pass the actual outside indices to the new scvf seed + scvfSeeds.back().addOutsideData(outsideScvfIdx, outsideLocalScvIdx); + } + else + // set the local scvfIndex of the found outside local scv face seed + actualScvSeed.setLocalScvfIndex(coordDir, outsideLocalScvfIdx); // pass info on inside to the outside found local scvf seed scvfSeeds[outsideLocalScvfIdx].addOutsideData(actualScvf.index(), actualLocalScvIdx); @@ -562,10 +632,10 @@ public: scvfSeeds.reserve(50); // The vertex index around which we construct the interaction volume - auto vIdxGlobal = scvf.vertexIndex(); + const auto vIdxGlobal = scvf.vertexIndex(); // Get the three scv faces in the scv - auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry); + const auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry); // create the scv entity seeds fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvfVector, vIdxGlobal); @@ -590,10 +660,10 @@ public: scvfSeeds.reserve(50); // The vertex index around which we construct the interaction volume - auto vIdxGlobal = scvf.vertexIndex(); + const auto vIdxGlobal = scvf.vertexIndex(); // Get the three scv faces in the scv - auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry); + const auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry); // create the scv entity seeds fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvfVector, vIdxGlobal); @@ -624,7 +694,7 @@ private: // make the scvf seeds for the three scvfs connected to the scv auto& actualScvSeed = scvSeeds.back(); - LocalIndexType actualLocalScvIdx = scvSeeds.size()-1; + const LocalIndexType actualLocalScvIdx = scvSeeds.size()-1; for (int coordDir = 0; coordDir < dim; ++coordDir) { @@ -637,17 +707,16 @@ private: actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size()); // create the scvf seed - auto faceType = Implementation::getMpfaFaceType(problem, element, actualScvf); scvfSeeds.emplace_back( actualScvf, actualLocalScvIdx, LocalIndexSet(), GlobalIndexSet(), - faceType ); + Implementation::getMpfaFaceType(problem, element, actualScvf) ); } else { - auto outsideGlobalScvIdx = actualScvf.outsideScvIdx(); - auto globalScvfIndex = actualScvf.index(); + const auto outsideGlobalScvIdx = actualScvf.outsideScvIdx(); + const auto globalScvfIndex = actualScvf.index(); // check if the outside scv already exists and get its local index bool outsideExists = false; @@ -656,7 +725,8 @@ private: { if (scvSeed.globalIndex() == outsideGlobalScvIdx) { - outsideExists = true; break; + outsideExists = true; + break; } // keep track of local index outsideLocalScvIdx++; @@ -669,14 +739,14 @@ private: actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size()); // get outside element, fvgeometry etc. - auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx); + const auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx); auto outsideFvGeometry = localView(problem.model().globalFvGeometry()); outsideFvGeometry.bindElement(outsideElement); // find scvf in outside corresponding to the actual scvf - auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry); - auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector); - auto&& outsideScvf = *outsideScvfVector[commonFaceLocalIdx]; + const auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry); + const auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector); + const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx]; // create scvf seed scvfSeeds.emplace_back(actualScvf, @@ -699,8 +769,23 @@ private: // boundary scvf seeds have no outside scvf if (!scvfSeed.boundary() && scvfSeed.outsideGlobalScvfIndex() == globalScvfIndex) { + const auto faceType = scvfSeed.faceType(); + // on interior boundaries we have to create a new face + if (faceType != MpfaFaceTypes::interior) + { + // set the local scvfIndex of the face that is about to created + actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size()); + + // create local scvf + scvfSeeds.emplace_back(actualScvf, + actualLocalScvIdx, + LocalIndexSet({scvfSeed.insideLocalScvIndex()}), + GlobalIndexSet({scvfSeed.insideGlobalScvfIndex()}), + faceType); + } // pass local scvf index to local scv - actualScvSeed.setLocalScvfIndex(coordDir, localScvfIdx); + else + actualScvSeed.setLocalScvfIndex(coordDir, localScvfIdx); // we found the corresponding face found = true; break; @@ -716,20 +801,20 @@ private: actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size()); // get outside element, fvgeometry etc. - auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx); + const auto outsideElement = problem.model().globalFvGeometry().element(outsideGlobalScvIdx); auto outsideFvGeometry = localView(problem.model().globalFvGeometry()); outsideFvGeometry.bindElement(outsideElement); // find scvf in outside corresponding to the actual scvf - auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry); - auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector); - auto&& outsideScvf = *outsideScvfVector[commonFaceLocalIdx]; + const auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry); + const auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector); + const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx]; // make scv face seed scvfSeeds.emplace_back(actualScvf, actualLocalScvIdx, - LocalIndexSet( {outsideLocalScvIdx} ), - GlobalIndexSet( {outsideScvf.index()} ), + LocalIndexSet({outsideLocalScvIdx}), + GlobalIndexSet({outsideScvf.index()}), Implementation::getMpfaFaceType(problem, element, actualScvf)); } } diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh index e49ba0320806727645946de5dd354041b4d47752..9450d4da519ced549991478677a9e486d2277195 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh @@ -97,15 +97,22 @@ class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase; using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper); + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); 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 InteriorBoundaryData = typename GET_PROP_TYPE(TypeTag, InteriorBoundaryData); + + 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); static const int dim = GridView::dimension; static const int dimWorld = GridView::dimensionworld; @@ -135,7 +142,7 @@ public: : problemPtr_(&problem), fvGeometryPtr_(&fvGeometry), elemVolVarsPtr_(&elemVolVars), - onBoundary_(seed.onBoundary()), + onDomainOrInteriorBoundary_(seed.onDomainOrInteriorBoundary()), globalScvfIndices_(seed.globalScvfIndices()) { // create local sub control entities from the seed @@ -149,42 +156,71 @@ public: 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()); - // eventually add dirichlet vol var indices and set up local index sets of flux and dirichlet faces + // 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 (localScvf.faceType() == MpfaFaceTypes::dirichlet) + if (faceType == MpfaFaceTypes::dirichlet) { volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex()); volVarsPositions_.push_back(localScvf.ip()); dirichletFaceIndexSet_.push_back(localScvfIdx++); } - else + 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."); } - // initialize the neumann fluxes vector to zero - neumannFluxes_ = DynamicVector(fluxFaceIndexSet_.size(), 0.0); + // initialize the vector containing the neumann fluxes + assembleNeumannFluxVector_(); } template void solveLocalSystem(const GetTensorFunction& getTensor) { - const std::size_t numFluxFaces = fluxScvfIndexSet_().size(); + const auto numFluxFaces = fluxScvfIndexSet_().size(); // if only dirichlet faces are present, assemble T_ directly if (numFluxFaces == 0) return assemblePureDirichletSystem_(getTensor); - const std::size_t numFaces = localScvfs_.size(); - const std::size_t numPotentials = volVarsStencil().size(); + const auto numFaces = localScvfs_.size(); + const auto numPotentials = volVarsStencil().size() + interiorDirichletScvfIndexSet_().size(); // the local matrices DynamicMatrix A(numFluxFaces, numFluxFaces, 0.0); @@ -203,31 +239,6 @@ public: T_ += D; } - void assembleNeumannFluxes(unsigned int eqIdx) - { - if (!onBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary)) - return; - - LocalIndexType fluxFaceIdx = 0; - for (auto localFluxFaceIdx : fluxFaceIndexSet_) - { - auto&& localScvf = localScvf_(localFluxFaceIdx); - const auto faceType = localScvf.faceType(); - if (faceType == MpfaFaceTypes::neumann) - { - auto&& element = localElement_(localScvf.insideLocalScvIndex()); - auto&& globalScvf = fvGeometry_().scvf(localScvf.insideGlobalScvfIndex()); - auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx]; - neumannFlux *= globalScvf.area(); - - // The flux is assumed to be prescribed in the form of -D*gradU - neumannFluxes_[fluxFaceIdx] = neumannFlux; - } - - fluxFaceIdx++; - } - } - //! gets local data on a global scvf within the interaction volume //! specialization for dim = dimWorld template @@ -236,16 +247,40 @@ public: { auto scvfGlobalIdx = scvf.index(); - LocalIndexType localIdx = 0; - for (auto&& localScvf : localScvfs_) + // if interior boundaries are disabled, do simplified search + if (!enableInteriorBoundaries) { - if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx) - return LocalFaceData(localIdx, localScvf.insideLocalScvIndex(), false); + 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); - 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++; + } - 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()); @@ -257,85 +292,131 @@ public: typename std::enable_if< (d < dw), LocalFaceData >::type getLocalFaceData(const SubControlVolumeFace& scvf) const { - auto scvfGlobalIdx = scvf.index(); + const auto scvfGlobalIdx = scvf.index(); - LocalIndexType localFaceIdx = 0; - for (auto&& localScvf : localScvfs_) + // if interior boundaries are disabled, do simplified search + if (!enableInteriorBoundaries) { - if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx) - return LocalFaceData(localFaceIdx, localScvf.insideLocalScvIndex(), false); - - if (!localScvf.boundary() && MpfaHelper::contains(localScvf.outsideGlobalScvfIndices(), scvfGlobalIdx)) + LocalIndexType localFaceIdx = 0; + for (auto&& localScvf : localScvfs_) { - if (localScvf.outsideLocalScvIndices().size() == 1) - return LocalFaceData(localFaceIdx, localScvf.outsideLocalScvIndex(), true); + if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx) + return LocalFaceData(localFaceIdx, localScvf.insideLocalScvIndex(), false); - // Now we have to find wich local scv is the one the global scvf is embedded in - auto scvGlobalIdx = scvf.insideScvIdx(); - for (auto localScvIdx : localScvf.outsideLocalScvIndices()) - if (localScv_(localScvIdx).globalIndex() == scvGlobalIdx) - return LocalFaceData(localFaceIdx, localScvIdx, true); + 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); + } - 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); - localFaceIdx++; + 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 + //! Gets the transmissibilities for a sub-control volume face within the interaction volume. //! specialization for dim == dimWorld template typename std::enable_if< (d == dw), DynamicVector >::type getTransmissibilities(const LocalFaceData& localFaceData) const { - auto tij = T_[localFaceData.localScvfIndex]; if (localFaceData.isOutside) + { + auto tij = T_[localFaceData.localScvfIndex]; tij *= -1.0; - return tij; + return tij; + } + else + return T_[localFaceData.localScvfIndex]; } - //! gets the transmissibilities for a sub-control volume face within the interaction volume - //! specialization for dim < dimWorld + //! Gets the transmissibilities for a sub-control volume face within the interaction volume. + //! specialization for dim < dimWorld. template typename std::enable_if< (d < dw), DynamicVector >::type getTransmissibilities(const LocalFaceData& localFaceData) const { - // This means we are on a branching point. If we come from the inside, simply return tij + // If we come from the inside, simply return tij if (!localFaceData.isOutside) return T_[localFaceData.localScvfIndex]; - // otherwise compute the outside transmissibilities - DynamicVector tij(volVarsStencil().size(), 0.0); + // compute the outside transmissibilities + DynamicVector tij(volVarsStencil().size() + interiorDirichletScvfIndexSet_().size(), 0.0); // get the local scv and iterate over local coordinates - std::size_t numLocalScvs = localScvs_.size(); - auto&& localScv = localScv_(localFaceData.localScvIndex); - auto&& localScvf = localScvf_(localFaceData.localScvfIndex); + const auto numLocalScvs = localScvs_.size(); + const auto numDirichletScvfs = dirichletScvfIndexSet_().size(); + const auto& localScv = localScv_(localFaceData.localScvIndex); + const auto& localScvf = localScvf_(localFaceData.localScvfIndex); - auto idxInOutside = this->findLocalIndex(localScvf.outsideLocalScvIndices(), localFaceData.localScvIndex); - auto&& wijk = wijk_[localFaceData.localScvfIndex][idxInOutside+1]; + const auto idxInOutside = this->findIndexInVector(localScvf.outsideLocalScvIndices(), localFaceData.localScvIndex); + const auto& wijk = wijk_[localFaceData.localScvfIndex][idxInOutside+1]; for (int localDir = 0; localDir < dim; localDir++) { - auto localScvfIdx = localScv.localScvfIndex(localDir); - auto&& localScvf = localScvf_(localScvfIdx); + const auto localScvfIdx = localScv.localScvfIndex(localDir); + const auto& localScvf = localScvf_(localScvfIdx); - // add entries from the face unknowns - if (localScvf.faceType() != MpfaFaceTypes::dirichlet) + const auto faceType = localScvf.faceType(); + if (faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet) { - auto fluxFaceIndex = this->findLocalIndex(fluxScvfIndexSet_(), localScvfIdx); + const auto fluxFaceIndex = this->findIndexInVector(fluxScvfIndexSet_(), localScvfIdx); auto tmp = AinvB_[fluxFaceIndex]; tmp *= wijk[localDir]; tij += tmp; } - else + else if (faceType == MpfaFaceTypes::dirichlet) { - auto idxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), localScvfIdx); + const auto idxInDiriFaces = this->findIndexInVector(dirichletScvfIndexSet_(), localScvfIdx); tij[numLocalScvs + idxInDiriFaces] += wijk[localDir]; } + else if (faceType == MpfaFaceTypes::interiorDirichlet) + { + const auto idxInInteriorDiriFaces = this->findIndexInVector(interiorDirichletScvfIndexSet_(), localScvfIdx); + tij[numLocalScvs + numDirichletScvfs + idxInInteriorDiriFaces] += wijk[localDir]; + } // add entry from the scv unknown tij[localFaceData.localScvIndex] -= wijk[localDir]; @@ -344,19 +425,38 @@ public: return tij; } - Scalar getNeumannFlux(const LocalFaceData& localFaceData) const + //! 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 + DynamicVector getNeumannFluxTransformationCoefficients(const LocalFaceData& localFaceData) const + { + auto cij = CAinv_[localFaceData.localScvfIndex]; + if (localFaceData.isOutside) + cij *= -1.0; + return cij; + } + + Scalar getNeumannFlux(const LocalFaceData& localFaceData, unsigned int eqIdx) const { - if (!onBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary) || fluxScvfIndexSet_().size() == 0 ) + if (!onDomainOrInteriorBoundary() || useTpfaBoundary || fluxScvfIndexSet_().size() == 0 ) return 0.0; - auto flux = CAinv_[localFaceData.localScvfIndex] * neumannFluxes_; + // Do the scalar product CAinv_*neumannFluxes[eqIdx] + assert(CAinv_[localFaceData.localScvfIndex].size() == neumannFluxes_.size() && + "Number of columns of matrix does not correspond to number entries in vector!"); + + Scalar flux(0.0); + for (unsigned int i = 0; i < neumannFluxes_.size(); ++i) + flux += CAinv_[localFaceData.localScvfIndex][i] * neumannFluxes_[i][eqIdx]; + + // flip sign if we are coming from the outside if (localFaceData.isOutside) return -1.0*flux; + return flux; } - bool onBoundary() const - { return onBoundary_; } + bool onDomainOrInteriorBoundary() const + { return onDomainOrInteriorBoundary_; } const GlobalIndexSet& volVarsStencil() const { return volVarsStencil_; } @@ -367,6 +467,9 @@ public: const GlobalIndexSet& globalScvfs() const { return globalScvfIndices_; } + const std::vector& interiorBoundaryData() const + { return interiorBoundaryData_; } + private: const LocalScvfType& localScvf_(const LocalIndexType localScvfIdx) const @@ -381,13 +484,19 @@ private: const LocalIndexSet& dirichletScvfIndexSet_() const { return dirichletFaceIndexSet_; } + const LocalIndexSet& interiorDirichletScvfIndexSet_() const + { return interiorDirichletFaceIndexSet_; } + + const LocalIndexSet& interiorBoundaryScvfIndexSet_() const + { return interiorBoundaryFaceIndexSet_; } + const Element& localElement_(const LocalIndexType localScvIdx) const { return localElements_[localScvIdx]; } void createLocalEntities_(const Seed& seed) { - auto numScvs = seed.scvSeeds().size(); - auto numScvfs = seed.scvfSeeds().size(); + const auto numScvs = seed.scvSeeds().size(); + const auto numScvfs = seed.scvfSeeds().size(); localElements_.reserve(numScvs); localScvs_.reserve(numScvs); @@ -408,6 +517,36 @@ private: } } + void assembleNeumannFluxVector_() + { + // initialize the neumann fluxes vector to zero + neumannFluxes_.resize(fluxFaceIndexSet_.size(), PrimaryVariables(0.0)); + + if (!onDomainOrInteriorBoundary() || useTpfaBoundary) + return; + + LocalIndexType fluxFaceIdx = 0; + for (auto localFluxFaceIdx : fluxFaceIndexSet_) + { + const auto& localScvf = localScvf_(localFluxFaceIdx); + const auto faceType = localScvf.faceType(); + + if (faceType == MpfaFaceTypes::neumann || (faceType == MpfaFaceTypes::interiorNeumann && !facetCoupling)) + { + const auto& element = localElement_(localScvf.insideLocalScvIndex()); + const auto& globalScvf = fvGeometry_().scvf(localScvf.insideGlobalScvfIndex()); + auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf); + neumannFlux *= globalScvf.area(); + neumannFlux *= elemVolVars_()[globalScvf.insideScvIdx()].extrusionFactor(); + + // The flux is assumed to be prescribed in the form of -D*gradU + neumannFluxes_[fluxFaceIdx] = neumannFlux; + } + + fluxFaceIdx++; + } + } + template void assembleLocalMatrices_(const GetTensorFunction& getTensor, DynamicMatrix& A, @@ -415,26 +554,28 @@ private: DynamicMatrix& C, DynamicMatrix& D) { - const std::size_t numLocalScvs = localScvs_.size(); + static const auto xi = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, Xi); + const auto numLocalScvs = localScvs_.size(); + const auto numDirichletScvfs = dirichletScvfIndexSet_().size(); // reserve space for the omegas wijk_.resize(localScvfs_.size()); // loop over the local faces - LocalIndexType rowIdx = 0; + unsigned int rowIdx = 0; for (auto&& localScvf : localScvfs_) { - auto faceType = localScvf.faceType(); - bool hasUnknown = faceType != MpfaFaceTypes::dirichlet; - LocalIndexType idxInFluxFaces = hasUnknown ? this->findLocalIndex(fluxScvfIndexSet_(), rowIdx) : -1; + const auto faceType = localScvf.faceType(); + const bool hasUnknown = faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet; + const LocalIndexType idxInFluxFaces = hasUnknown ? this->findIndexInVector(fluxScvfIndexSet_(), rowIdx) : -1; // get diffusion tensor in "positive" sub volume - auto posLocalScvIdx = localScvf.insideLocalScvIndex(); - auto&& posLocalScv = localScv_(posLocalScvIdx); - auto&& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex()); - auto&& posVolVars = elemVolVars_()[posGlobalScv]; - auto element = localElement_(posLocalScvIdx); - auto tensor = getTensor(element, posVolVars, posGlobalScv); + const auto posLocalScvIdx = localScvf.insideLocalScvIndex(); + const auto& posLocalScv = localScv_(posLocalScvIdx); + const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex()); + const auto& posVolVars = elemVolVars_()[posGlobalScv]; + const auto& element = localElement_(posLocalScvIdx); + const auto tensor = getTensor(element, posVolVars, posGlobalScv); // the omega factors of the "positive" sub volume auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor); @@ -443,51 +584,106 @@ private: // Check the local directions of the positive sub volume for (int localDir = 0; localDir < dim; localDir++) { - auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir); - auto&& curLocalScvf = localScvf_(curLocalScvfIdx); + const auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir); + const auto& curLocalScvf = localScvf_(curLocalScvfIdx); + const auto curFaceType = curLocalScvf.faceType(); - // First, add the entries associated with face pressures (unkown or dirichlet) - if (curLocalScvf.faceType() != MpfaFaceTypes::dirichlet) + // First, add the entries associated with unknown face pressures + if (curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet) { // we need the index of the current local scvf in the flux face indices - auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx); + auto curIdxInFluxFaces = this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx); + // this creates an entry in matrix C C[rowIdx][curIdxInFluxFaces] += posWijk[localDir]; + + // proceed depending on if the current face has an unknown if (hasUnknown) - A[idxInFluxFaces][curIdxInFluxFaces] += posWijk[localDir]; + { + if (faceType == MpfaFaceTypes::interiorNeumann) + { + // on interior neumann faces, apply xi factor + // However, if this is a boundary face at the same time, don't! + A[idxInFluxFaces][curIdxInFluxFaces] += localScvf.globalScvf().boundary() ? + posWijk[localDir] : xi*posWijk[localDir]; + + // add values from other domain in case of facet coupling + if (facetCoupling && curIdxInFluxFaces == idxInFluxFaces) + { + // 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_()); + + // calculate "leakage factor" + const auto n = curLocalScvf.unitOuterNormal(); + const auto v = [&] () + { + auto res = n; + res *= -0.5*facetVolVars.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); + } + } + // this means we are on an interior face + else + A[idxInFluxFaces][curIdxInFluxFaces] += posWijk[localDir]; + } } - else + else if (curFaceType == MpfaFaceTypes::dirichlet) { // the current face is a Dirichlet face and creates entries in D & eventually B - auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx); + auto curIdxInDiriFaces = this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx); D[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir]; if (hasUnknown) B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] -= posWijk[localDir]; } + else if (curFaceType == MpfaFaceTypes::interiorDirichlet) + { + // the current face is an interior Dirichlet face and creates entries in D & eventually B + auto curIdxInInteriorDiriFaces = this->findIndexInVector(interiorDirichletScvfIndexSet_(), curLocalScvfIdx); + + D[rowIdx][numLocalScvs + numDirichletScvfs + curIdxInInteriorDiriFaces] += posWijk[localDir]; + if (hasUnknown) + B[idxInFluxFaces][numLocalScvs + numDirichletScvfs + curIdxInInteriorDiriFaces] -= posWijk[localDir]; + } // add entries related to pressures at the scv centers (dofs) D[rowIdx][posLocalScvIdx] -= posWijk[localDir]; if (hasUnknown) - B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir]; + { + if (faceType == MpfaFaceTypes::interiorNeumann && !localScvf.globalScvf().boundary()) + B[idxInFluxFaces][posLocalScvIdx] += xi*posWijk[localDir]; + else + B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir]; + } } // store the omegas wijk_[rowIdx].emplace_back(std::move(posWijk)); - // If face is not on a boundary or an interior dirichlet face, add entries for the "negative" scv - if (faceType == MpfaFaceTypes::interior) + // If we are on an interior or interior neumann face, add values from negative sub volume + if (!localScvf.globalScvf().boundary() && + (faceType == MpfaFaceTypes::interior || faceType == MpfaFaceTypes::interiorNeumann)) { // loop over all the outside neighbors of this face and add entries unsigned int indexInOutsideData = 0; for (auto negLocalScvIdx : localScvf.outsideLocalScvIndices()) { - auto&& negLocalScv = localScv_(negLocalScvIdx); - auto&& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex()); - auto&& negVolVars = elemVolVars_()[negGlobalScv]; - auto negElement = localElement_(negLocalScvIdx); - auto negTensor = getTensor(negElement, negVolVars, negGlobalScv); + const auto& negLocalScv = localScv_(negLocalScvIdx); + const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex()); + const auto& negVolVars = elemVolVars_()[negGlobalScv]; + const auto& negElement = localElement_(negLocalScvIdx); + const auto negTensor = getTensor(negElement, negVolVars, negGlobalScv); // the omega factors of the "negative" sub volume DimVector negWijk; @@ -496,7 +692,7 @@ private: if (dim < dimWorld) { // outside scvf - auto&& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndex(indexInOutsideData)); + const auto& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndex(indexInOutsideData)); auto negNormal = outsideScvf.unitOuterNormal(); negNormal *= -1.0; negWijk = calculateOmegas_(negLocalScv, negNormal, localScvf.area(), negTensor); @@ -510,24 +706,30 @@ private: // Check local directions of negative sub volume for (int localDir = 0; localDir < dim; localDir++) { - auto curLocalScvfIdx = negLocalScv.localScvfIndex(localDir); - auto&& curLocalScvf = localScvf_(curLocalScvfIdx); + const auto curLocalScvfIdx = negLocalScv.localScvfIndex(localDir); + const auto& curLocalScvf = localScvf_(curLocalScvfIdx); + const auto curFaceType = curLocalScvf.faceType(); - if (curLocalScvf.faceType() != MpfaFaceTypes::dirichlet) + if (curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet) { - // we need the index of the current local scvf in the flux face indices - auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx); - A[idxInFluxFaces][curIdxInFluxFaces] -= negWijk[localDir]; + if (faceType == MpfaFaceTypes::interiorNeumann) + A[idxInFluxFaces][this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx)] -= (1-xi)*negWijk[localDir]; + else + A[idxInFluxFaces][this->findIndexInVector(fluxScvfIndexSet_(), curLocalScvfIdx)] -= negWijk[localDir]; } - else + else if (curFaceType == MpfaFaceTypes::interiorDirichlet) { - // the current face is a Dirichlet face and creates entries in B - auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx); - B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] += negWijk[localDir]; + const auto idxInInteriorDiriFaces = this->findIndexInVector(interiorDirichletScvfIndexSet_(), curLocalScvfIdx); + B[idxInFluxFaces][numLocalScvs + numDirichletScvfs + idxInInteriorDiriFaces] += negWijk[localDir]; } + else // Dirichlet face + B[idxInFluxFaces][numLocalScvs + this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx)] += negWijk[localDir]; // add entries to matrix B - B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir]; + if (faceType == MpfaFaceTypes::interiorNeumann) + B[idxInFluxFaces][negLocalScvIdx] -= (1-xi)*negWijk[localDir]; + else + B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir]; } // store the omegas (negative, because of normal vector sign switch) @@ -546,9 +748,10 @@ private: template void assemblePureDirichletSystem_(const GetTensorFunction& getTensor) { - const std::size_t numLocalScvs = localScvs_.size(); - const std::size_t numFaces = localScvfs_.size(); - const std::size_t numPotentials = volVarsStencil().size(); + const auto numLocalScvs = localScvs_.size(); + const auto numFaces = localScvfs_.size(); + const auto numInteriorDirichletFaces = interiorDirichletScvfIndexSet_().size(); + const auto numPotentials = volVarsStencil().size() + numInteriorDirichletFaces; // resize matrices, only T_ will have entries T_.resize(numFaces, numPotentials, 0.0); @@ -563,12 +766,12 @@ private: for (auto&& localScvf : localScvfs_) { // get diffusion tensor in "positive" sub volume - auto posLocalScvIdx = localScvf.insideLocalScvIndex(); - auto&& posLocalScv = localScv_(posLocalScvIdx); - auto&& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex()); - auto&& posVolVars = elemVolVars_()[posGlobalScv]; - auto element = localElement_(posLocalScvIdx); - auto tensor = getTensor(element, posVolVars, posGlobalScv); + const auto posLocalScvIdx = localScvf.insideLocalScvIndex(); + const auto& posLocalScv = localScv_(posLocalScvIdx); + const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex()); + const auto& posVolVars = elemVolVars_()[posGlobalScv]; + const auto element = localElement_(posLocalScvIdx); + const auto tensor = getTensor(element, posVolVars, posGlobalScv); // the omega factors of the "positive" sub volume auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor); @@ -576,11 +779,27 @@ private: for (int localDir = 0; localDir < dim; localDir++) { - auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir); - auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx); + // When interior boundaries are disabled, all faces will be of dirichlet type + if (!enableInteriorBoundaries) + { + const auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir); + const auto curIdxInDiriFaces = this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx); + T_[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir]; + T_[rowIdx][posLocalScvIdx] -= posWijk[localDir]; + } + else + { + const auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir); + const auto& curLocalScvf = localScvf_(curLocalScvfIdx); + const auto curFaceType = curLocalScvf.faceType(); + + const auto curIdxInDiriFaces = curFaceType == MpfaFaceTypes::dirichlet ? + this->findIndexInVector(dirichletScvfIndexSet_(), curLocalScvfIdx) : + numInteriorDirichletFaces + this->findIndexInVector(interiorDirichletScvfIndexSet_(), curLocalScvfIdx); - T_[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir]; - T_[rowIdx][posLocalScvIdx] -= posWijk[localDir]; + T_[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir]; + T_[rowIdx][posLocalScvIdx] -= posWijk[localDir]; + } } // store the omegas @@ -646,7 +865,7 @@ private: const FVElementGeometry* fvGeometryPtr_; const ElementVolumeVariables* elemVolVarsPtr_; - bool onBoundary_; + bool onDomainOrInteriorBoundary_; std::vector localElements_; std::vector localScvs_; @@ -658,6 +877,9 @@ private: LocalIndexSet fluxFaceIndexSet_; LocalIndexSet dirichletFaceIndexSet_; + LocalIndexSet interiorDirichletFaceIndexSet_; + LocalIndexSet interiorBoundaryFaceIndexSet_; + std::vector interiorBoundaryData_; std::vector< std::vector< DimVector > > wijk_; @@ -665,7 +887,7 @@ private: DynamicMatrix AinvB_; DynamicMatrix CAinv_; - DynamicVector neumannFluxes_; + std::vector neumannFluxes_; }; } // end namespace diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh index 3a4c5ccf35092ed127a68546801ad554267d6902..a3b11e1add00ca1e89888715af4f59428ce9c8a5 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh @@ -24,6 +24,7 @@ #define DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUMESEED_HH #include +#include #include "localsubcontrolentityseeds.hh" namespace Dumux @@ -45,13 +46,13 @@ public: CCMpfaOInteractionVolumeSeed(std::vector&& scvSeeds, std::vector&& scvfSeeds, - bool onBoundary) - : onBoundary_(onBoundary), + bool onDomainOrInteriorBoundary) + : onDomainOrInteriorBoundary_(onDomainOrInteriorBoundary), scvSeeds_(std::move(scvSeeds)), scvfSeeds_(std::move(scvfSeeds)) {} - bool onBoundary() const - { return onBoundary_; } + bool onDomainOrInteriorBoundary() const + { return onDomainOrInteriorBoundary_; } const std::vector& scvSeeds() const { return scvSeeds_; } @@ -78,15 +79,18 @@ public: for (auto&& localScvfSeed : scvfSeeds()) { globalIndices.push_back(localScvfSeed.insideGlobalScvfIndex()); - for (auto scvfIdxGlobal : localScvfSeed.outsideGlobalScvfIndices()) - globalIndices.push_back(scvfIdxGlobal); + + // add outside indices for interior faces + if (localScvfSeed.faceType() == MpfaFaceTypes::interior) + for (auto scvfIdxGlobal : localScvfSeed.outsideGlobalScvfIndices()) + globalIndices.push_back(scvfIdxGlobal); } return globalIndices; } private: - bool onBoundary_; + bool onDomainOrInteriorBoundary_; std::vector scvSeeds_; std::vector scvfSeeds_; }; diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh index c3d8b4c98e5785857375389fe67f20b9d1249c96..959706e166409d6e5fffd1a9fb9beac26f299fb5 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh @@ -136,9 +136,7 @@ public: CCMpfaOLocalScvf(const LocalScvfSeed& scvfSeed, const SubControlVolumeFace& scvf) : seedPtr_(&scvfSeed), - ip_(scvf.ipGlobal()), - normal_(scvf.unitOuterNormal()), - area_(scvf.area()) + scvfPtr_(&scvf) {} GlobalIndexType insideGlobalScvfIndex() const @@ -172,25 +170,26 @@ public: { return scvfSeed_().faceType(); } GlobalPosition ip() const - { return ip_; } + { return globalScvf().ipGlobal(); } GlobalPosition unitOuterNormal() const - { return normal_; } + { return globalScvf().unitOuterNormal(); } Scalar area() const - { return area_; } + { return globalScvf().area(); } bool boundary() const { return scvfSeed_().boundary(); } + const SubControlVolumeFace& globalScvf() const + { return *scvfPtr_; } + private: const LocalScvfSeed& scvfSeed_() const { return *seedPtr_; } const LocalScvfSeed* seedPtr_; - GlobalPosition ip_; - GlobalPosition normal_; - Scalar area_; + const SubControlVolumeFace* scvfPtr_; }; } // end namespace diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh index 2262c42f330e1fcc66647d51a19fa8f226c0f977..6a905987d6bbff25c1c9560da92e0b3e266d2ac9 100644 --- a/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh @@ -41,8 +41,14 @@ template class CCMpfaOFpsInteractionVolumeTraits : public CCMpfaOInteractionVolumeTraits { public: - // the fps method uses its own interaction volumes at the boundary - using BoundaryInteractionVolume = CCMpfaInteractionVolumeImplementation; + // Interior boundaries can not yet be handled by the currend o-method fps implementation + // In that case we use the o-interactionvolume, otherwise we use its own interaction volumes at the boundary + // TODO Fix the std::conditional + using BoundaryInteractionVolume = typename CCMpfaOInteractionVolumeTraits::BoundaryInteractionVolume; + // using BoundaryInteractionVolume = typename std::conditional::BoundaryInteractionVolume, + // typename CCMpfaInteractionVolumeImplementation + // >::type; // The local sub-control volume type differs from the standard mpfa-o method using LocalScvType = CCMpfaOFpsLocalScv; @@ -61,6 +67,7 @@ class CCMpfaInteractionVolumeImplementation : using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); using LocalScvType = typename Traits::LocalScvType; @@ -116,6 +123,9 @@ public: { if (dim == 3) DUNE_THROW(Dune::NotImplemented, "Fps scheme in 3d"); + + //! add entry to the vector of neumann fluxes + addAuxiliaryCellNeumannFlux_(); } template @@ -147,20 +157,19 @@ public: this->T_ += mc.BF; } - void assembleNeumannFluxes(unsigned int eqIdx) - { - ParentType::assembleNeumannFluxes(eqIdx); +private: - if (!this->onBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary)) + void addAuxiliaryCellNeumannFlux_() + { + if (!this->onDomainOrInteriorBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary)) return; + // add one entry to the vector of neumann fluxes auto& neumannFluxes = this->neumannFluxes_; neumannFluxes.resize(this->fluxScvfIndexSet_().size() + 1); - neumannFluxes[neumannFluxes.size()-1] = std::accumulate(neumannFluxes.begin(), neumannFluxes.end()-1, 0.0); + neumannFluxes[neumannFluxes.size()-1] = std::accumulate(neumannFluxes.begin(), neumannFluxes.end()-1, PrimaryVariables(0.0)); } -private: - template void assembleLocalMatrices_(const GetTensorFunction& getTensor, LocalMatrixContainer& mc) { @@ -277,7 +286,7 @@ private: bool isRHS = false) { // In case we're on a flux continuity face, get local index - LocalIndexType eqSystemIdx = isFluxEq ? this->findLocalIndex(this->fluxScvfIndexSet_(), rowIdx) : -1; + LocalIndexType eqSystemIdx = isFluxEq ? this->findIndexInVector(this->fluxScvfIndexSet_(), rowIdx) : -1; // Fluxes stemming from the RHS have to have the opposite sign Scalar factor = isRHS ? 1.0 : -1.0; @@ -298,7 +307,7 @@ private: localD *= localScv.geometry().integrationElement(ipLocal); // add matrix entries for the pressure in the cell center - auto cellPressureIdx = this->findLocalIndex(this->volVarsStencil(), localScv.globalIndex()); + auto cellPressureIdx = this->findIndexInVector(this->volVarsStencil(), localScv.globalIndex()); Scalar bi0 = factor*(localD[normalDir]*shapeJacobian[0][0]); if (!isRHS) mc.BF[rowIdx][cellPressureIdx] += bi0; if (isFluxEq) B[eqSystemIdx][cellPressureIdx] += bi0; @@ -310,14 +319,14 @@ private: if (this->localScvf_(localScvfIdx).faceType() != MpfaFaceTypes::dirichlet) { Scalar aij = factor*(localD[normalDir]*shapeJacobian[localDir+1][0]); - auto colIdx = this->findLocalIndex(this->fluxScvfIndexSet_(), localScvfIdx); + auto colIdx = this->findIndexInVector(this->fluxScvfIndexSet_(), localScvfIdx); if (!isRHS) mc.AF[rowIdx][colIdx] += aij; if (isFluxEq) A[eqSystemIdx][colIdx] += aij; } else { Scalar bij = factor*(localD[normalDir]*shapeJacobian[localDir+1][0]); - auto colIdx = this->localScvs_.size() + this->findLocalIndex(this->dirichletScvfIndexSet_(), localScvfIdx); + auto colIdx = this->localScvs_.size() + this->findIndexInVector(this->dirichletScvfIndexSet_(), localScvfIdx); if (!isRHS) mc.BF[rowIdx][colIdx] += bij; if (isFluxEq) B[eqSystemIdx][colIdx] += bij; } @@ -355,7 +364,7 @@ private: localD *= localScv.geometry().integrationElement(ipLocal); // add matrix entries for the pressure in the cell center - auto cellPressureIdx = this->findLocalIndex(this->volVarsStencil(), localScv.globalIndex()); + auto cellPressureIdx = this->findIndexInVector(this->volVarsStencil(), localScv.globalIndex()); mc.BL[divEqIdx_][cellPressureIdx] += factor*(localD[normalDir]*shapeJacobian[0][0]); // Add entries from the local scv faces @@ -365,12 +374,12 @@ private: if (this->localScvf_(localScvfIdx).faceType() != MpfaFaceTypes::dirichlet) { - auto colIdx = this->findLocalIndex(this->fluxScvfIndexSet_(), localScvfIdx); + auto colIdx = this->findIndexInVector(this->fluxScvfIndexSet_(), localScvfIdx); mc.AL[divEqIdx_][colIdx] += factor*(localD[normalDir]*shapeJacobian[localDir+1][0]); } else { - auto colIdx = this->localScvs_.size() + this->findLocalIndex(this->dirichletScvfIndexSet_(), localScvfIdx); + auto colIdx = this->localScvs_.size() + this->findIndexInVector(this->dirichletScvfIndexSet_(), localScvfIdx); mc.BL[divEqIdx_][colIdx] += factor*(localD[normalDir]*shapeJacobian[localDir+1][0]); } } diff --git a/dumux/discretization/cellcentered/mpfa/stencils.hh b/dumux/discretization/cellcentered/mpfa/stencils.hh deleted file mode 100644 index 734951864ea814c656c4bfad1d40b752c0e4e6a7..0000000000000000000000000000000000000000 --- a/dumux/discretization/cellcentered/mpfa/stencils.hh +++ /dev/null @@ -1,41 +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 . * - *****************************************************************************/ -/*! - * \file - * \brief Implements the notion of stencils for cell-centered models - */ -#ifndef DUMUX_DISCRETIZATION_CC_MPFA_STENCILS_HH -#define DUMUX_DISCRETIZATION_CC_MPFA_STENCILS_HH - -#include - -namespace Dumux -{ -//! Method-specific implementation of the stencils vector -//! By default, we inherit from the symmetric stencils vector class -//! For non-symmetric mpfa methods a corresponding specialization has to be provided below -template -class CCMpfaStencilsVectorImplementation : public CCSymmetricStencilsVector {}; - -template -using CCMpfaStencilsVector = CCMpfaStencilsVectorImplementation; - -} // end namespace - -#endif diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh index 63cedea3ba391784a6380b80445ef72e89b6a478..3a3143739464163573a757f40d09403c5a765160 100644 --- a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh +++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh @@ -63,7 +63,6 @@ class DarcysLawImplementation using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; - using Stencil = std::vector; static const int dim = GridView::dimension; static const int dimWorld = GridView::dimensionworld; @@ -198,24 +197,6 @@ public: } } - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { - if (scvf.boundary()) - return Stencil({scvf.insideScvIdx()}); - else if (scvf.numOutsideScvs() > 1) - { - Stencil stencil({scvf.insideScvIdx()}); - for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) - stencil.push_back(scvf.outsideScvIdx(i)); - return stencil; - } - else - return Stencil({scvf.insideScvIdx(), scvf.outsideScvIdx()}); - } - // The flux variables cache has to be bound to an element prior to flux calculations // During the binding, the transmissibilities will be computed and stored using the method below. static Scalar calculateTransmissibilities(const Problem& problem, diff --git a/dumux/discretization/cellcentered/tpfa/fickslaw.hh b/dumux/discretization/cellcentered/tpfa/fickslaw.hh index 8fcd9ce4a014a2d3ad31f9501ed918c3ab8e9918..7b05340bfad569bf72afc7fdff95d0d794c53d9d 100644 --- a/dumux/discretization/cellcentered/tpfa/fickslaw.hh +++ b/dumux/discretization/cellcentered/tpfa/fickslaw.hh @@ -56,7 +56,6 @@ class FicksLawImplementation using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using IndexType = typename GridView::IndexSet::IndexType; - using Stencil = typename std::vector; using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); @@ -109,17 +108,6 @@ public: return rho*tij*(xInside - xOutside); } - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { - if (!scvf.boundary()) - return Stencil({scvf.insideScvIdx(), scvf.outsideScvIdx()}); - else - return Stencil({scvf.insideScvIdx()}); - } - private: //! compute the mole/mass fraction at branching facets for network grids diff --git a/dumux/discretization/cellcentered/tpfa/fourierslaw.hh b/dumux/discretization/cellcentered/tpfa/fourierslaw.hh index 1806ae63dac1eb10aae445adda614003c4f53680..c46876547ae2ff4f9a4f60fb44861d615272fa46 100644 --- a/dumux/discretization/cellcentered/tpfa/fourierslaw.hh +++ b/dumux/discretization/cellcentered/tpfa/fourierslaw.hh @@ -53,7 +53,6 @@ class FouriersLawImplementation using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using IndexType = typename GridView::IndexSet::IndexType; - using Stencil = typename std::vector; using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); using Element = typename GridView::template Codim<0>::Entity; @@ -89,17 +88,6 @@ public: return tij*(tInside - tOutside); } - static Stencil stencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf) - { - if (!scvf.boundary()) - return Stencil({scvf.insideScvIdx(), scvf.outsideScvIdx()}); - else - return Stencil({scvf.insideScvIdx()}); - } - private: //! compute the temperature at branching facets for network grids diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh index 6d9efed3fd510502d1730de2f74454fd9fa85c50..64d46973bef124dd676e489c6727cffcb92863f8 100644 --- a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh +++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh @@ -375,6 +375,9 @@ private: int scvfCounter = 0; for (const auto& intersection : intersections(globalFvGeometry().gridView(), element)) { + // check if intersection is on interior boundary + const auto isInteriorBoundary = globalFvGeometry().problem_().isInteriorBoundary(element, intersection); + if (dim < dimWorld) if (handledScvf[intersection.indexInInside()]) continue; @@ -386,7 +389,8 @@ private: scvfs_.emplace_back(intersection, intersection.geometry(), scvFaceIndices[scvfCounter], - scvIndices); + scvIndices, + intersection.boundary() || isInteriorBoundary); scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]); scvfCounter++; @@ -418,11 +422,14 @@ private: int scvfCounter = 0; for (const auto& intersection : intersections(globalFvGeometry().gridView(), element)) { + // check if intersection is on interior boundary + const auto isInteriorBoundary = globalFvGeometry().problem_().isInteriorBoundary(element, intersection); + if (dim < dimWorld) if (handledScvf[intersection.indexInInside()]) continue; - if (intersection.neighbor()) + if (intersection.neighbor() && !isInteriorBoundary) { // only create subcontrol faces where the outside element is the bound element if (dim == dimWorld) @@ -434,7 +441,8 @@ private: neighborScvfs_.emplace_back(intersection, intersection.geometry(), scvFaceIndices[scvfCounter], - scvIndices); + scvIndices, + false); neighborScvfIndices_.push_back(scvFaceIndices[scvfCounter]); } @@ -454,7 +462,8 @@ private: neighborScvfs_.emplace_back(intersection, intersection.geometry(), scvFaceIndices[scvfCounter], - scvIndices); + scvIndices, + false); neighborScvfIndices_.push_back(scvFaceIndices[scvfCounter]); break; @@ -468,7 +477,7 @@ private: handledScvf[intersection.indexInInside()] = true; scvfCounter++; } - else if (intersection.boundary()) + else if (intersection.boundary() || isInteriorBoundary) { // for surface and network grids mark that we handled this face if (dim < dimWorld) diff --git a/dumux/discretization/cellcentered/tpfa/globalfvgeometry.hh b/dumux/discretization/cellcentered/tpfa/globalfvgeometry.hh index 59a5758639fbadc8252139765a9d28c6d70eb9fb..7e48089aeae5ee7c4cfebaaa1de2392d33640e23 100644 --- a/dumux/discretization/cellcentered/tpfa/globalfvgeometry.hh +++ b/dumux/discretization/cellcentered/tpfa/globalfvgeometry.hh @@ -156,8 +156,11 @@ public: for (const auto& intersection : intersections(gridView_, element)) { + // check if intersection is on interior boundary + const auto isInteriorBoundary = problem.isInteriorBoundary(element, intersection); + // inner sub control volume faces - if (intersection.neighbor()) + if (intersection.neighbor() && !isInteriorBoundary) { if (dim == dimWorld) { @@ -165,7 +168,8 @@ public: scvfs_.emplace_back(intersection, intersection.geometry(), scvfIdx, - std::vector({eIdx, nIdx})); + std::vector({eIdx, nIdx}), + false); scvfsIndexSet.push_back(scvfIdx++); } // this is for network grids @@ -183,20 +187,21 @@ public: scvfs_.emplace_back(intersection, intersection.geometry(), scvfIdx, - scvIndices); + scvIndices, + false); scvfsIndexSet.push_back(scvfIdx++); outsideIndices[indexInInside].clear(); } } } // boundary sub control volume faces - else if (intersection.boundary()) + else if (intersection.boundary() || isInteriorBoundary) { scvfs_.emplace_back(intersection, intersection.geometry(), scvfIdx, - std::vector({eIdx, gridView_.size(0) + numBoundaryScvf_++}) - ); + std::vector({eIdx, gridView_.size(0) + numBoundaryScvf_++}), + true); scvfsIndexSet.push_back(scvfIdx++); } } @@ -391,8 +396,11 @@ public: for (const auto& intersection : intersections(gridView_, element)) { + // check if intersection is on interior boundary + const auto isInteriorBoundary = problem.isInteriorBoundary(element, intersection); + // inner sub control volume faces - if (intersection.neighbor()) + if (intersection.neighbor() && !isInteriorBoundary) { if (dim == dimWorld) { @@ -417,7 +425,7 @@ public: } } // boundary sub control volume faces - else if (intersection.boundary()) + else if (intersection.boundary() || isInteriorBoundary) { scvfsIndexSet.push_back(numScvf_++); neighborVolVarIndexSet.push_back({numScvs_ + numBoundaryScvf_++}); diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh index d77fcb98c6c1cacdc3b10c8882141d3eeec84619..30a287bfcaf020eb32c6debe0ae69c4b5f08ee12 100644 --- a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh +++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh @@ -54,12 +54,21 @@ public: // the default constructor CCTpfaSubControlVolumeFace() = default; - //! Constructor with intersection + /*! + * \brief Constructor with intersection + * + * \param is The intersection + * \param isGeometry The geometry of the intersection + * \param scvfIndex The global index of this scv face + * \param scvIndices The inside/outside scv indices connected to this face + * \param isBoundary Bool to specify whether or not the scvf is on an interior or the domain boundary + */ template CCTpfaSubControlVolumeFace(const Intersection& is, const typename Intersection::Geometry& isGeometry, IndexType scvfIndex, - const std::vector& scvIndices) + const std::vector& scvIndices, + bool isBoundary) : ParentType(), geomType_(isGeometry.type()), area_(isGeometry.volume()), @@ -67,7 +76,7 @@ public: unitOuterNormal_(is.centerUnitOuterNormal()), scvfIndex_(scvfIndex), scvIndices_(scvIndices), - boundary_(is.boundary()) + boundary_(isBoundary) { corners_.resize(isGeometry.corners()); for (int i = 0; i < isGeometry.corners(); ++i) diff --git a/dumux/discretization/fluxstencil.hh b/dumux/discretization/fluxstencil.hh new file mode 100644 index 0000000000000000000000000000000000000000..72a49236c3c4064ef831bd66d9b459b6310d1e70 --- /dev/null +++ b/dumux/discretization/fluxstencil.hh @@ -0,0 +1,133 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * \brief Base class the flux stencil + */ +#ifndef DUMUX_DISCRETIZATION_FLUXSTENCIL_HH +#define DUMUX_DISCRETIZATION_FLUXSTENCIL_HH + +#include +#include + +namespace Dumux +{ + +//! Forward declaration of the upwind scheme implementation +template +class FluxStencilImplementation; + +/*! + * \ingroup Discretization + * \brief The flux stencil specialized for each discretization method + * \note There might be different stencils used for e.g. advection and diffusion for schemes + * where the stencil depends on variables. Also schemes might even have solution dependent + * stencil. However, we always reserve the stencil or all DOFs that are possibly involved + * since we use the flux stencil for matrix and assembly. This might lead to some zeros stored + * in the matrix. + */ +template +using FluxStencil = FluxStencilImplementation; + +//! Flux stencil for the box method +template +class FluxStencilImplementation +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using Element = typename GridView::template Codim<0>::Entity; + using IndexType = typename GridView::IndexSet::IndexType; + using Stencil = std::vector; + +public: + // This is for compatibility with the cc methods. The flux stencil info is obsolete for the box method. + static Stencil stencil(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) + { + return Stencil(); + } +}; + +//! Flux stencil for the cell-centered TPFA scheme +template +class FluxStencilImplementation +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using Element = typename GridView::template Codim<0>::Entity; + using IndexType = typename GridView::IndexSet::IndexType; + using Stencil = std::vector; + +public: + static Stencil stencil(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) + { + if (scvf.boundary()) + return Stencil({scvf.insideScvIdx()}); + else if (scvf.numOutsideScvs() > 1) + { + Stencil stencil({scvf.insideScvIdx()}); + for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) + stencil.push_back(scvf.outsideScvIdx(i)); + return stencil; + } + else + return Stencil({scvf.insideScvIdx(), scvf.outsideScvIdx()}); + } +}; + +//! Specialization for cell-centered MPFA schemes +template +class FluxStencilImplementation +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using Element = typename GridView::template Codim<0>::Entity; + using IndexType = typename GridView::IndexSet::IndexType; + using Stencil = std::vector; + +public: + static Stencil stencil(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) + { + const auto& globalFvGeometry = problem.model().globalFvGeometry(); + + // return the scv (element) indices in the interaction region + if (globalFvGeometry.touchesInteriorOrDomainBoundary(scvf)) + return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices(); + else + return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices(); + } +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/discretization/fluxvariablesbase.hh b/dumux/discretization/fluxvariablesbase.hh index 216edff503ad123ec3d8c3e20dc9cc43f757c716..b20dd3f15e97211e66fc547fa802c2a6e41adfd6 100644 --- a/dumux/discretization/fluxvariablesbase.hh +++ b/dumux/discretization/fluxvariablesbase.hh @@ -24,23 +24,34 @@ #define DUMUX_DISCRETIZATION_FLUXVARIABLESBASE_HH #include +#include +#include +#include +#include namespace Dumux { -namespace Properties -{ -// forward declaration -NEW_PROP_TAG(ImplicitUpwindWeight); -} +template +class FluxVariablesBaseImplementation; /*! * \ingroup Discretization - * \brief Base class for the flux variables - * Actual flux variables inherit from this class + * \brief The flux variables base class class + * The upwind scheme is chosen depending on the discretization method + */ +template +using FluxVariablesBase = FluxVariablesBaseImplementation, FluxStencil>; + +/*! + * \ingroup Discretization + * \brief Implementation of the base class of the flux variables + * + * \param TypeTag The type tag + * \param UpwindScheme The type used for the upwinding of the advective fluxes */ -template -class FluxVariablesBase +template +class FluxVariablesBaseImplementation { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); @@ -48,7 +59,6 @@ class FluxVariablesBase using IndexType = typename GridView::IndexSet::IndexType; using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using Stencil = std::vector; - using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); @@ -69,10 +79,6 @@ public: fvGeometryPtr_ = &fvGeometry; elemVolVarsPtr_ = &elemVolVars; elemFluxVarsCachePtr_ = &elemFluxVarsCache; - // retrieve the upwind weight for the mass conservation equations. Use the value - // specified via the property system as default, and overwrite - // it by the run-time parameter from the Dune::ParameterTree - upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, UpwindWeight); } const Problem& problem() const @@ -93,117 +99,33 @@ public: const ElementFluxVariablesCache& elemFluxVarsCache() const { return *elemFluxVarsCachePtr_; } - Stencil computeStencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) - { DUNE_THROW(Dune::InvalidStateException, "computeStencil() routine is not provided by the implementation."); } - - // For cell-centered surface and network grids (dim < dimWorld) we have to do a special upwind scheme - template - typename std::enable_if::type - upwindScheme(Scalar flux, int phaseIdx, const FunctionType& upwindTerm) + //! Applies the upwind scheme to precalculated fluxes + template + Scalar applyUpwindScheme(const UpwindTermFunction& upwindTerm, Scalar flux, int phaseIdx) { - const auto& insideScv = this->fvGeometry().scv(this->scvFace().insideScvIdx()); - const auto& insideVolVars = this->elemVolVars()[insideScv]; - - // check if this is a branching point - if (this->scvFace().numOutsideScvs() > 1) - { - // more complicated upwind scheme - // we compute a flux-weighted average of all inflowing branches - Scalar branchingPointUpwindTerm = 0.0; - Scalar sumUpwindFluxes = 0.0; - - // if the inside flux is positive (outflow) do fully upwind and return flux - if (!std::signbit(flux)) - return upwindTerm(insideVolVars)*flux; - else - sumUpwindFluxes += flux; - - for (unsigned int i = 0; i < this->scvFace().numOutsideScvs(); ++i) - { - // compute the outside flux - const auto outsideScvIdx = this->scvFace().outsideScvIdx(i); - const auto outsideElement = this->fvGeometry().globalFvGeometry().element(outsideScvIdx); - const auto& flippedScvf = this->fvGeometry().flipScvf(this->scvFace().index(), i); - - const auto outsideFlux = AdvectionType::flux(this->problem(), - outsideElement, - this->fvGeometry(), - this->elemVolVars(), - flippedScvf, - phaseIdx, - this->elemFluxVarsCache()); - - if (!std::signbit(outsideFlux)) - branchingPointUpwindTerm += upwindTerm(this->elemVolVars()[outsideScvIdx])*outsideFlux; - else - sumUpwindFluxes += outsideFlux; - } - - // the flux might be zero - if (sumUpwindFluxes != 0.0) - branchingPointUpwindTerm /= -sumUpwindFluxes; - else - branchingPointUpwindTerm = 0.0; - - // upwind scheme (always do fully upwind at branching points) - // a weighting here would lead to an error since the derivation is based on a fully upwind scheme - // TODO How to implement a weight of e.g. 0.5 - if (std::signbit(flux)) - return flux*branchingPointUpwindTerm; - else - return flux*upwindTerm(insideVolVars); - } - // non-branching points and boundaries - else - { - // upwind scheme - const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()]; - if (std::signbit(flux)) - return flux*(upwindWeight_*upwindTerm(outsideVolVars) - + (1.0 - upwindWeight_)*upwindTerm(insideVolVars)); - else - return flux*(upwindWeight_*upwindTerm(insideVolVars) - + (1.0 - upwindWeight_)*upwindTerm(outsideVolVars)); - } + //! Give the upwind scheme access to the cached variables + return UpwindScheme::apply(*this, upwindTerm, flux, phaseIdx); } - // For grids with dim == dimWorld or the box-method we use a simple upwinding scheme - template - typename std::enable_if::type - upwindScheme(Scalar flux, int phaseIdx, const FunctionType& upwindTerm) + static Stencil computeStencil(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) { - const auto& insideScv = this->fvGeometry().scv(this->scvFace().insideScvIdx()); - const auto& insideVolVars = this->elemVolVars()[insideScv]; - - // upwind scheme - const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()]; - if (std::signbit(flux)) - return flux*(upwindWeight_*upwindTerm(outsideVolVars) - + (1.0 - upwindWeight_)*upwindTerm(insideVolVars)); - else - return flux*(upwindWeight_*upwindTerm(insideVolVars) - + (1.0 - upwindWeight_)*upwindTerm(outsideVolVars)); + //! Give the upwind scheme access to the cached variables + //! Forward to the discretization specific implementation + return FluxStencil::stencil(problem, element, fvGeometry, scvf); } private: - - Implementation &asImp_() - { return *static_cast(this); } - - const Implementation &asImp_() const - { return *static_cast(this); } - const Problem* problemPtr_; //! Pointer to the problem const Element* elementPtr_; //! Pointer to the element at hand const FVElementGeometry* fvGeometryPtr_; const SubControlVolumeFace* scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created const ElementVolumeVariables* elemVolVarsPtr_; const ElementFluxVariablesCache* elemFluxVarsCachePtr_; - Scalar upwindWeight_; }; -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/discretization/upwindscheme.hh b/dumux/discretization/upwindscheme.hh new file mode 100644 index 0000000000000000000000000000000000000000..94d1cdc60d745338a881f657ef2cafc862a0eeab --- /dev/null +++ b/dumux/discretization/upwindscheme.hh @@ -0,0 +1,401 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * \brief Base class for the upwind scheme + */ +#ifndef DUMUX_DISCRETIZATION_UPWINDSCHEME_HH +#define DUMUX_DISCRETIZATION_UPWINDSCHEME_HH + +#include +#include +#include + +namespace Dumux +{ + +namespace Properties +{ +// forward declaration +NEW_PROP_TAG(ImplicitUpwindWeight); +NEW_PROP_TAG(EnableInteriorBoundaries); +NEW_PROP_TAG(MpfaFacetCoupling); +NEW_PROP_TAG(UseTpfaBoundary); +} + +//! Forward declaration of the upwind scheme implementation +template +class UpwindSchemeImplementation; + + +//! Upwind scheme for the box method +template +class UpwindSchemeImplementation +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + +public: + // applies a simple upwind scheme to the precalculated advective flux + template + static Scalar apply(const FluxVariables& fluxVars, + const UpwindTermFunction& upwindTerm, + Scalar flux, int phaseIdx) + { + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + static const Scalar upwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, UpwindWeight); + + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()]; + if (std::signbit(flux)) + return flux*(upwindWeight*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); + else + return flux*(upwindWeight*upwindTerm(insideVolVars) + + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); + } +}; + +//! Upwind scheme for the cell-centered TPFA scheme +template +class UpwindSchemeImplementation +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); + + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + +public: + // For surface and network grids (dim < dimWorld) we have to do a special upwind scheme + template + static typename std::enable_if<(d < dw), Scalar>::type + apply(const FluxVariables& fluxVars, + const UpwindTermFunction& upwindTerm, + Scalar flux, int phaseIdx) + { + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + static const Scalar upwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, UpwindWeight); + + // the volume variables of the inside sub-control volume + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + + // check if this is a branching point + if (fluxVars.scvFace().numOutsideScvs() > 1) + { + // more complicated upwind scheme + // we compute a flux-weighted average of all inflowing branches + Scalar branchingPointUpwindTerm = 0.0; + Scalar sumUpwindFluxes = 0.0; + + // if the inside flux is positive (outflow) do fully upwind and return flux + if (!std::signbit(flux)) + return upwindTerm(insideVolVars)*flux; + else + sumUpwindFluxes += flux; + + for (unsigned int i = 0; i < fluxVars.scvFace().numOutsideScvs(); ++i) + { + // compute the outside flux + const auto outsideScvIdx = fluxVars.scvFace().outsideScvIdx(i); + const auto outsideElement = fluxVars.fvGeometry().globalFvGeometry().element(outsideScvIdx); + const auto& flippedScvf = fluxVars.fvGeometry().flipScvf(fluxVars.scvFace().index(), i); + + const auto outsideFlux = AdvectionType::flux(fluxVars.problem(), + outsideElement, + fluxVars.fvGeometry(), + fluxVars.elemVolVars(), + flippedScvf, + phaseIdx, + fluxVars.elemFluxVarsCache()); + + if (!std::signbit(outsideFlux)) + branchingPointUpwindTerm += upwindTerm(fluxVars.elemVolVars()[outsideScvIdx])*outsideFlux; + else + sumUpwindFluxes += outsideFlux; + } + + // the flux might be zero + if (sumUpwindFluxes != 0.0) + branchingPointUpwindTerm /= -sumUpwindFluxes; + else + branchingPointUpwindTerm = 0.0; + + // upwind scheme (always do fully upwind at branching points) + // a weighting here would lead to an error since the derivation is based on a fully upwind scheme + // TODO How to implement a weight of e.g. 0.5 + if (std::signbit(flux)) + return flux*branchingPointUpwindTerm; + else + return flux*upwindTerm(insideVolVars); + } + // non-branching points and boundaries + else + { + // upwind scheme + const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()]; + if (std::signbit(flux)) + return flux*(upwindWeight*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); + else + return flux*(upwindWeight*upwindTerm(insideVolVars) + + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); + } + } + + // For grids with dim == dimWorld we use a simple upwinding scheme + template + static typename std::enable_if<(d == dw), Scalar>::type + apply(const FluxVariables& fluxVars, + const UpwindTermFunction& upwindTerm, + Scalar flux, int phaseIdx) + { + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + static const Scalar upwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, UpwindWeight); + + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()]; + if (std::signbit(flux)) + return flux*(upwindWeight*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); + else + return flux*(upwindWeight*upwindTerm(insideVolVars) + + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); + } +}; + +//! Specialization for cell-centered MPFA schemes +template +class UpwindSchemeImplementation +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); + + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + + 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: + + // For surface and network grids (dim < dimWorld) we have to do a special upwind scheme + template + static typename std::enable_if<(d < dw), Scalar>::type + apply(const FluxVariables& fluxVars, + const UpwindTermFunction& upwindTerm, + Scalar flux, int phaseIdx) + { + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + static const Scalar upwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, UpwindWeight); + + // check if we have to handle a branching point + const auto isInteriorBoundary = enableInteriorBoundaries ? + fluxVars.elemFluxVarsCache()[fluxVars.scvFace()].isInteriorBoundary() : + false; + + // on branching points (which are not interior boundaries) we use a more complicated upwind scheme + if (fluxVars.scvFace().numOutsideScvs() > 1 && !isInteriorBoundary) + { + // the volume variables of the inside sub-control volume + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + + // we compute a flux-weighted average of all inflowing branches + Scalar branchingPointUpwindTerm = 0.0; + Scalar sumUpwindFluxes = 0.0; + + // if the inside flux is positive (outflow) do fully upwind and return flux + if (!std::signbit(flux)) + return upwindTerm(insideVolVars)*flux; + else + sumUpwindFluxes += flux; + + for (unsigned int i = 0; i < fluxVars.scvFace().numOutsideScvs(); ++i) + { + // compute the outside flux + const auto outsideScvIdx = fluxVars.scvFace().outsideScvIdx(i); + const auto outsideElement = fluxVars.fvGeometry().globalFvGeometry().element(outsideScvIdx); + const auto& flippedScvf = fluxVars.fvGeometry().flipScvf(fluxVars.scvFace().index(), i); + + const auto outsideFlux = AdvectionType::flux(fluxVars.problem(), + outsideElement, + fluxVars.fvGeometry(), + fluxVars.elemVolVars(), + flippedScvf, + phaseIdx, + fluxVars.elemFluxVarsCache()); + + if (!std::signbit(outsideFlux)) + branchingPointUpwindTerm += upwindTerm(fluxVars.elemVolVars()[outsideScvIdx])*outsideFlux; + else + sumUpwindFluxes += outsideFlux; + } + + // the flux might be zero + if (sumUpwindFluxes != 0.0) + branchingPointUpwindTerm /= -sumUpwindFluxes; + else + branchingPointUpwindTerm = 0.0; + + // upwind scheme (always do fully upwind at branching points) + // a weighting here would lead to an error since the derivation is based on a fully upwind scheme + // TODO How to implement a weight of e.g. 0.5 + if (std::signbit(flux)) + return flux*branchingPointUpwindTerm; + else + return flux*upwindTerm(insideVolVars); + } + // non-branching points and domain boundaries + else if (!isInteriorBoundary) + { + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()]; + if (std::signbit(flux)) + return flux*(upwindWeight*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); + else + return flux*(upwindWeight*upwindTerm(insideVolVars) + + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); + } + // interior boundaries + else + { + const auto& scvfFluxVarsCache = fluxVars.elemFluxVarsCache()[fluxVars.scvFace()]; + + // on interior Dirichlet Boundaries or for active facetCoupling we use the facet vol vars + if (facetCoupling || scvfFluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorDirichlet) + { + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + const auto& outsideVolVars = scvfFluxVarsCache.interiorBoundaryDataSelf().facetVolVars(fluxVars.fvGeometry()); + if (std::signbit(flux)) + return flux*(upwindWeight*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); + else + return flux*(upwindWeight*upwindTerm(insideVolVars) + + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); + } + + // This is an interior neumann boundary. Thus, the flux at hand carries the user specified + // boundary conditions and we do not have to do any upwinding here. + // Note: + // Mpfa on the boundary can essentially only be used for flux terms of the form F = -DgradU, + // without an appearing upwind term. Thus, only for rather mathematical than realistic problems. + // Assuming that -DgradU has been specified by the user, we simply return the precalculated flux, + // as it has the boundary fluxes incorporated. Note that for Neumann boundaries above we did do + // the upwinding, thus, when using Mpfa on Neumann boundaries an upwind term of 1.0 has to be set + // for it to work. We usually achieve this by using the Unit fluid Component. Here, we can't apply + // the upwinding though because the outside vol vars are not defined. + // TODO: Implement Neumann BCs using additional Dofs for !useTpfaBoundary property + return flux; + } + } + + // For grids with dim == dimWorld we use a simple upwinding scheme + template + static typename std::enable_if<(d == dw), Scalar>::type + apply(const FluxVariables& fluxVars, + const UpwindTermFunction& upwindTerm, + Scalar flux, int phaseIdx) + { + // retrieve the upwind weight for the mass conservation equations. Use the value + // specified via the property system as default, and overwrite + // it by the run-time parameter from the Dune::ParameterTree + static const Scalar upwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, UpwindWeight); + + if (!enableInteriorBoundaries) + { + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()]; + if (std::signbit(flux)) + return flux*(upwindWeight*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); + else + return flux*(upwindWeight*upwindTerm(insideVolVars) + + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); + } + else + { + const auto& scvfFluxVarsCache = fluxVars.elemFluxVarsCache()[fluxVars.scvFace()]; + + if (!scvfFluxVarsCache.isInteriorBoundary()) + { + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()]; + if (std::signbit(flux)) + return flux*(upwindWeight*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); + else + return flux*(upwindWeight*upwindTerm(insideVolVars) + + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); + } + // on interior boundaries, the outside volume variables might be facet volume variables + else + { + // on interior Dirichlet Boundaries or for active facetCoupling we use the facet vol vars + if (facetCoupling || scvfFluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorDirichlet) + { + const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()]; + const auto& outsideVolVars = scvfFluxVarsCache.interiorBoundaryDataSelf().facetVolVars(fluxVars.fvGeometry()); + if (std::signbit(flux)) + return flux*(upwindWeight*upwindTerm(outsideVolVars) + + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); + else + return flux*(upwindWeight*upwindTerm(insideVolVars) + + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); + } + + // This is an interior neumann boundary. Thus, the flux at hand carries the user specified + // boundary conditions and we do not have to do any upwinding here. + // Note: + // Mpfa on the boundary can essentially only be used for flux terms of the form F = -DgradU, + // without an appearing upwind term. Thus, only for rather mathematical than realistic problems. + // Assuming that -DgradU has been specified by the user, we simply return the precalculated flux, + // as it has the boundary fluxes incorporated. Note that for Neumann boundaries above we did do + // the upwinding, thus, when using Mpfa on Neumann boundaries an upwind term of 1.0 has to be set + // for it to work. We usually achieve this by using the Unit fluid Component. Here, we can't apply + // the upwinding though because the outside vol vars are not defined. + // TODO: Implement Neumann BCs using additional Dofs for !useTpfaBoundary property + return flux; + } + } + } +}; + +/*! + * \ingroup Discretization + * \brief The upwind scheme used for the advective fluxes. + * This depends on the chosen discretization method. + */ +template +using UpwindScheme = UpwindSchemeImplementation; + + +} // end namespace + +#endif diff --git a/dumux/implicit/cellcentered/mpfa/properties.hh b/dumux/implicit/cellcentered/mpfa/properties.hh index 274783410ea5828db3a85df22b34435cb398f20d..f78c330e81eb1f0fea0c3d883be694784f7ab48e 100644 --- a/dumux/implicit/cellcentered/mpfa/properties.hh +++ b/dumux/implicit/cellcentered/mpfa/properties.hh @@ -49,6 +49,9 @@ NEW_PROP_TAG(InteractionVolume); //! The inner interaction volume type NEW_PROP_TAG(BoundaryInteractionVolume); //! The interaction volume type used on the boundaries NEW_PROP_TAG(GlobalInteractionVolumeSeeds); //! Class storing and managing the interaction volume seeds NEW_PROP_TAG(UseTpfaBoundary); //! This property specifies whether or not tpfa is to be used to handle the boundary fluxes +NEW_PROP_TAG(InteriorBoundaryData); //! Stores and obtains additional data on interior boundaries +NEW_PROP_TAG(MpfaFacetCoupling); //! Specifies if the interior boundaries are static or coupled to another domain +NEW_PROP_TAG(MpfaXi); //! Parameter for the coupling conditions when using mpfa in a context with coupling on facets (0 <= xi <= 1) NEW_PROP_TAG(MpfaQ); //! The quadrature point on the sub control volume faces (0.0 <= q <= 1.0) NEW_PROP_TAG(MpfaC); //! Parameterisation of the size of the auxiliary volume in the FPS scheme NEW_PROP_TAG(MpfaP); //! The quadrature point on the auxiliary sub faces in the FPS scheme diff --git a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh index 3278dafaef82cb8bda520f5522d7cd82cacf0ff0..483f831e8367bbe97157c1a43e74944574674f9d 100644 --- a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh @@ -41,6 +41,7 @@ #include #include #include +#include #include #include @@ -108,6 +109,15 @@ public: // By default, we use tpfa on the boundaries SET_BOOL_PROP(CCMpfaModel, UseTpfaBoundary, true); +// By default, interior boundaries are static +SET_BOOL_PROP(CCMpfaModel, MpfaFacetCoupling, false); + +// The default interior Dirichlet boundary data +SET_TYPE_PROP(CCMpfaModel, InteriorBoundaryData, InteriorBoundaryData); + +// By default, we use simple coupling conditions (Xi = 1) +SET_SCALAR_PROP(CCMpfaModel, MpfaXi, 1.0); + // By default, we set the quadrature point to the mid point of the element facets SET_SCALAR_PROP(CCMpfaModel, MpfaQ, 0.0); diff --git a/dumux/implicit/cellcentered/properties.hh b/dumux/implicit/cellcentered/properties.hh index bea7da122be7702eefd8d15cad7991fe9f63aafb..c5757fcbc3db4f28f83c9d6d633e339c4e60c65f 100644 --- a/dumux/implicit/cellcentered/properties.hh +++ b/dumux/implicit/cellcentered/properties.hh @@ -44,6 +44,7 @@ namespace Properties NEW_TYPE_TAG(CCModel, INHERITS_FROM(ImplicitBase)); NEW_PROP_TAG(AssemblyMap); +NEW_PROP_TAG(EnableInteriorBoundaries); //! Enables or disables the use of internal boundaries on interior scvfs } } diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh index 20cadb848c32b332e8e8a0f1914571b823d2f078..702259c31e8b32641980917887ad55f61c822ec8 100644 --- a/dumux/implicit/cellcentered/propertydefaults.hh +++ b/dumux/implicit/cellcentered/propertydefaults.hh @@ -79,6 +79,9 @@ SET_TYPE_PROP(CCModel, BaseLocalResidual, CCLocalResidual); //! Set the AssemblyMap to the CCAssemblyMap (default: symmetric occupation pattern) SET_TYPE_PROP(CCModel, AssemblyMap, Dumux::CCSymmetricAssemblyMap); +// By default, we disable interior boundaries +SET_BOOL_PROP(CCModel, EnableInteriorBoundaries, false); + //! indicate that this is no box discretization SET_BOOL_PROP(CCModel, ImplicitIsBox, false); diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index a067d55a698951c230053405bb4eccfca9c7f7be..6124493f88251dd3a750c93ed34007a104ab9281 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -786,9 +786,11 @@ protected: Dune::GeometryType geomType = element.geometry().type(); const ReferenceElement &refElement = ReferenceElements::general(geomType); - for (const auto& intersection : intersections(gridView_(), element)) { - if (intersection.boundary()) { - if (isBox) + for (const auto& intersection : intersections(gridView_(), element)) + { + if (isBox) + { + if (intersection.boundary()) { // add all vertices on the intersection to the set of // boundary vertices @@ -806,7 +808,10 @@ protected: boundaryIndices_[vIdxGlobal] = true; } } - else + } + else + { + if (intersection.boundary() || problem_().isInteriorBoundary(element, intersection)) { int eIdxGlobal = elementMapper().index(element); boundaryIndices_[eIdxGlobal] = true; diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh index 9bdc2590519b1c83fce42a971e5153db87cfc961..7b862f6ddec42a85a357ae9dfb4edc078fb66e94 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -1018,6 +1018,21 @@ public: std::vector getAdditionalDofDependencies(IndexType globalIdx) const { return std::vector(); } + /*! + * \brief Function to set intersections as interior boundaries. This functionality is only + * available for models using cell-centered schemes. The corresponding boundary + * types and conditions are obtained from the standard methods. + * + * \param element The finite element + * \param intersection The intersection within the element + * \return boolean to mark an intersection as an interior boundary + * + * Per default we don't have interior boundaries + */ + template + bool isInteriorBoundary(const Element& element, const Intersection& intersection) const + { return false; } + /*! * \brief Capability to introduce problem-specific routines at the * beginning of the grid adaptation diff --git a/dumux/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh index d2840c174ad24982fcdda768a40a3361209ca1dd..87bd43686a0dadaaba111f73468990d28a286eb4 100644 --- a/dumux/porousmediumflow/implicit/fluxvariables.hh +++ b/dumux/porousmediumflow/implicit/fluxvariables.hh @@ -57,16 +57,14 @@ using PorousMediumFluxVariables = PorousMediumFluxVariablesImpl -class PorousMediumFluxVariablesImpl -: public FluxVariablesBase> +class PorousMediumFluxVariablesImpl : public FluxVariablesBase { - using ParentType = FluxVariablesBase>; + using ParentType = FluxVariablesBase; 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 IndexType = typename GridView::IndexSet::IndexType; using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using Stencil = std::vector; using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); @@ -100,29 +98,21 @@ public: phaseIdx, this->elemFluxVarsCache()); - return this->upwindScheme(flux, phaseIdx, upwindTerm); + return this->applyUpwindScheme(upwindTerm, flux, phaseIdx); } - - Stencil computeStencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvFace) - { return AdvectionType::stencil(problem, element, fvGeometry, scvFace); } }; // specialization for isothermal advection molecularDiffusion equations template -class PorousMediumFluxVariablesImpl -: public FluxVariablesBase> +class PorousMediumFluxVariablesImpl : public FluxVariablesBase { - using ParentType = FluxVariablesBase>; + using ParentType = FluxVariablesBase; 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 IndexType = typename GridView::IndexSet::IndexType; using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using Stencil = std::vector; using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); @@ -164,43 +154,20 @@ public: advFluxCached_.set(phaseIdx, true); } - return this->upwindScheme(advPreFlux_[phaseIdx], phaseIdx, upwindTerm); + return this->applyUpwindScheme(upwindTerm, advPreFlux_[phaseIdx], phaseIdx); } Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) { - Scalar flux = MolecularDiffusionType::flux(this->problem(), - this->element(), - this->fvGeometry(), - this->elemVolVars(), - this->scvFace(), - phaseIdx, compIdx, - this->elemFluxVarsCache()); - return flux; + return MolecularDiffusionType::flux(this->problem(), + this->element(), + this->fvGeometry(), + this->elemVolVars(), + this->scvFace(), + phaseIdx, compIdx, + this->elemFluxVarsCache()); } - Stencil computeStencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvFace) - { - // In the case of cctpfa or box the stencils for all laws are the same... - if (GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::CCTpfa - || GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box) - { - return AdvectionType::stencil(problem, element, fvGeometry, scvFace); - } - - // ...in general: unifiy advective and diffusive stencil - Stencil stencil = AdvectionType::stencil(problem, element, fvGeometry, scvFace); - Stencil diffusionStencil = MolecularDiffusionType::stencil(problem, element, fvGeometry, scvFace); - - stencil.insert(stencil.end(), diffusionStencil.begin(), diffusionStencil.end()); - std::sort(stencil.begin(), stencil.end()); - stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end()); - - return stencil; - } private: //! simple caching if advection flux is used twice with different upwind function std::bitset advFluxCached_; @@ -209,16 +176,14 @@ private: // specialization for non-isothermal advective flow (e.g. non-isothermal one-phase darcy equation) template -class PorousMediumFluxVariablesImpl -: public FluxVariablesBase> +class PorousMediumFluxVariablesImpl : public FluxVariablesBase { - using ParentType = FluxVariablesBase>; + using ParentType = FluxVariablesBase; 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 IndexType = typename GridView::IndexSet::IndexType; using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using Stencil = std::vector; using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); @@ -258,41 +223,17 @@ public: advFluxCached_.set(phaseIdx, true); } - return this->upwindScheme(advPreFlux_[phaseIdx], phaseIdx, upwindTerm); + return this->applyUpwindScheme(upwindTerm, advPreFlux_[phaseIdx], phaseIdx); } Scalar heatConductionFlux() { - Scalar flux = HeatConductionType::flux(this->problem(), - this->element(), - this->fvGeometry(), - this->elemVolVars(), - this->scvFace(), - this->elemFluxVarsCache()); - return flux; - } - - Stencil computeStencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvFace) - { - // In the case of cctpfa or box the stencils for all laws are the same... - if (GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::CCTpfa - || GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box) - { - return AdvectionType::stencil(problem, element, fvGeometry, scvFace); - } - - // ...in general: unifiy advective and heat conduction stencil - Stencil stencil = AdvectionType::stencil(problem, element, fvGeometry, scvFace); - Stencil energyStencil = HeatConductionType::stencil(problem, element, fvGeometry, scvFace); - - stencil.insert(stencil.end(), energyStencil.begin(), energyStencil.end()); - std::sort(stencil.begin(), stencil.end()); - stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end()); - - return stencil; + return HeatConductionType::flux(this->problem(), + this->element(), + this->fvGeometry(), + this->elemVolVars(), + this->scvFace(), + this->elemFluxVarsCache()); } private: @@ -303,16 +244,14 @@ private: // specialization for non-isothermal advection and difussion equations (e.g. non-isothermal three-phase three-component flow) template -class PorousMediumFluxVariablesImpl -: public FluxVariablesBase> +class PorousMediumFluxVariablesImpl : public FluxVariablesBase { - using ParentType = FluxVariablesBase>; + using ParentType = FluxVariablesBase; 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 IndexType = typename GridView::IndexSet::IndexType; using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using Stencil = std::vector; using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); @@ -353,55 +292,28 @@ public: advFluxCached_.set(phaseIdx, true); } - return this->upwindScheme(advPreFlux_[phaseIdx], phaseIdx, upwindTerm); + return this->applyUpwindScheme(upwindTerm, advPreFlux_[phaseIdx], phaseIdx); } Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) { - Scalar flux = MolecularDiffusionType::flux(this->problem(), - this->element(), - this->fvGeometry(), - this->elemVolVars(), - this->scvFace(), - phaseIdx, compIdx, - this->elemFluxVarsCache()); - return flux; + return MolecularDiffusionType::flux(this->problem(), + this->element(), + this->fvGeometry(), + this->elemVolVars(), + this->scvFace(), + phaseIdx, compIdx, + this->elemFluxVarsCache()); } Scalar heatConductionFlux() { - Scalar flux = HeatConductionType::flux(this->problem(), - this->element(), - this->fvGeometry(), - this->elemVolVars(), - this->scvFace(), - this->elemFluxVarsCache()); - return flux; - } - - Stencil computeStencil(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvFace) - { - // In the case of cctpfa or box the stencils for all laws are the same... - if (GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::CCTpfa - || GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box) - { - return AdvectionType::stencil(problem, element, fvGeometry, scvFace); - } - - // ...in general: unifiy advective, diffusive and heat conduction stencil - Stencil stencil = AdvectionType::stencil(problem, element, fvGeometry, scvFace); - Stencil diffusionStencil = MolecularDiffusionType::stencil(problem, element, fvGeometry, scvFace); - Stencil energyStencil = HeatConductionType::stencil(problem, element, fvGeometry, scvFace); - - stencil.insert(stencil.end(), diffusionStencil.begin(), diffusionStencil.end()); - stencil.insert(stencil.end(), energyStencil.begin(), energyStencil.end()); - std::sort(stencil.begin(), stencil.end()); - stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end()); - - return stencil; + return HeatConductionType::flux(this->problem(), + this->element(), + this->fvGeometry(), + this->elemVolVars(), + this->scvFace(), + this->elemFluxVarsCache()); } private: diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh index 0b5a82bd1efba7b822371e96efea5638a337515e..1e01be146b4c4486c9bb63d930d105fa56e185da 100644 --- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh +++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh @@ -57,6 +57,7 @@ class PorousMediumFluxVariablesCacheImplementation::Entity; using IndexType = typename GridView::IndexSet::IndexType; @@ -77,6 +78,7 @@ public: void update(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace &scvf) { const auto geometry = element.geometry(); @@ -87,10 +89,6 @@ public: jacInvT_ = geometry.jacobianInverseTransposed(ipLocal); localBasis.evaluateJacobian(ipLocal, shapeJacobian_); localBasis.evaluateFunction(ipLocal, shapeValues_); // do we need the shapeValues for the flux? - - // The stencil info is obsolete for the box method. - // It is here for compatibility with cc methods - stencil_ = Stencil(0); } const std::vector& shapeJacobian() const @@ -102,6 +100,8 @@ 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_; @@ -138,8 +138,7 @@ public: const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace &scvf) { - FluxVariables fluxVars; - stencil_ = fluxVars.computeStencil(problem, element, fvGeometry, scvf); + stencil_ = FluxVariables::computeStencil(problem, element, fvGeometry, scvf); tij_ = AdvectionType::calculateTransmissibilities(problem, element, fvGeometry, elemVolVars, scvf); } @@ -174,64 +173,70 @@ class MpfaAdvectionCache using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume); - static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); + static 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 TransmissibilityVector = typename BoundaryInteractionVolume::Vector; + using CoefficientVector = typename BoundaryInteractionVolume::Vector; using PositionVector = typename BoundaryInteractionVolume::PositionVector; public: - MpfaAdvectionCache() { phaseNeumannFluxes_.fill(0.0); } - //! update cached objects - template - void updateAdvection(const SubControlVolumeFace &scvf, - const InteractionVolume& interactionVolume) + template + void updateAdvection(const InteractionVolume& interactionVolume, + const SubControlVolumeFace &scvf, + const LocalFaceData& scvfLocalFaceData) { - const auto& localFaceData = interactionVolume.getLocalFaceData(scvf); - volVarsStencil_ = interactionVolume.volVarsStencil(); - volVarsPositions_ = interactionVolume.volVarsPositions(); - tij_ = interactionVolume.getTransmissibilities(localFaceData); - } - - //! update cached neumann boundary flux - template - void updatePhaseNeumannFlux(const SubControlVolumeFace &scvf, - const InteractionVolume& interactionVolume, - const unsigned int phaseIdx) - { - const auto& localFaceData = interactionVolume.getLocalFaceData(scvf); - phaseNeumannFluxes_[phaseIdx] = interactionVolume.getNeumannFlux(localFaceData); + // 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 unsigned int phaseIdx) const - { return volVarsStencil_; } + 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 unsigned int phaseIdx) const - { return volVarsPositions_; } + 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 TransmissibilityVector& advectionTij(const unsigned int phaseIdx) const - { return tij_; } + 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(const unsigned int phaseIdx) const + Scalar advectionNeumannFlux(unsigned int phaseIdx) const { return phaseNeumannFluxes_[phaseIdx]; } private: // Quantities associated with advection - Stencil volVarsStencil_; - PositionVector volVarsPositions_; - TransmissibilityVector tij_; + Stencil advectionVolVarsStencil_; + PositionVector advectionVolVarsPositions_; + CoefficientVector advectionTij_; + CoefficientVector advectionCij_; std::array phaseNeumannFluxes_; }; @@ -246,43 +251,76 @@ class MpfaDiffusionCache 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 TransmissibilityVector = typename BoundaryInteractionVolume::Vector; + 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 - void updateDiffusion(const SubControlVolumeFace &scvf, - const InteractionVolume& interactionVolume, - const unsigned int phaseIdx, - const unsigned int compIdx) + template + void updateDiffusion(const InteractionVolume& interactionVolume, + const SubControlVolumeFace &scvf, + const LocalFaceData& scvfLocalFaceData, + unsigned int phaseIdx, + unsigned int compIdx) { - const auto& localFaceData = interactionVolume.getLocalFaceData(scvf); diffusionVolVarsStencils_[phaseIdx][compIdx] = interactionVolume.volVarsStencil(); - diffusionTij_[phaseIdx][compIdx] = interactionVolume.getTransmissibilities(localFaceData); + 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(const unsigned int phaseIdx, - const unsigned int compIdx) const + 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 TransmissibilityVector& diffusionTij(const unsigned int phaseIdx, - const unsigned int compIdx) const + 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, numPhases> diffusionVolVarsStencils_; - std::array< std::array, numPhases> diffusionTij_; + std::array< std::array, numPhases> diffusionTij_; + std::array< std::array, numPhases> diffusionCij_; + + // diffusive neumann flux for single-phasic models + std::array componentNeumannFluxes_; }; //! Base class for the heat conduction cache in mpfa methods @@ -297,27 +335,23 @@ class MpfaHeatConductionCache // 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 TransmissibilityVector = typename BoundaryInteractionVolume::Vector; -public: - MpfaHeatConductionCache() : heatNeumannFlux_(0.0) {} + 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 - void updateHeatConduction(const SubControlVolumeFace &scvf, - const InteractionVolume& interactionVolume) + template + void updateHeatConduction(const InteractionVolume& interactionVolume, + const SubControlVolumeFace &scvf, + const LocalFaceData& scvfLocalFaceData) { - const auto& localFaceData = interactionVolume.getLocalFaceData(scvf); heatConductionVolVarsStencil_ = interactionVolume.volVarsStencil(); - heatConductionTij_ = interactionVolume.getTransmissibilities(localFaceData); - } + heatConductionTij_ = interactionVolume.getTransmissibilities(scvfLocalFaceData); + heatNeumannFlux_ = interactionVolume.getNeumannFlux(scvfLocalFaceData, energyEqIdx); - // update cached objects for heat conduction - template - void updateHeatNeumannFlux(const SubControlVolumeFace &scvf, - const InteractionVolume& interactionVolume) - { - const auto& localFaceData = interactionVolume.getLocalFaceData(scvf); - heatNeumannFlux_ = interactionVolume.getNeumannFlux(localFaceData); + if (enableInteriorBoundaries) + heatConductionCij_ = interactionVolume.getNeumannFluxTransformationCoefficients(scvfLocalFaceData); } //! Returns the volume variables indices necessary for heat conduction flux @@ -328,9 +362,14 @@ public: //! Returns the transmissibilities associated with the volume variables //! This can be different for the phases & components. - const TransmissibilityVector& heatConductionTij() const + 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 @@ -339,7 +378,8 @@ public: private: // Quantities associated with heat conduction Stencil heatConductionVolVarsStencil_; - TransmissibilityVector heatConductionTij_; + CoefficientVector heatConductionTij_; + CoefficientVector heatConductionCij_; Scalar heatNeumannFlux_; }; @@ -349,168 +389,104 @@ class PorousMediumFluxVariablesCacheImplementation {}; - -// specialization for the case of pure advection -template -class MpfaPorousMediumFluxVariablesCache - : public MpfaAdvectionCache + GET_PROP_VALUE(TypeTag, EnableEnergyBalance)> { - using AdvectionCache = MpfaAdvectionCache; + using InteriorBoundaryData = typename GET_PROP_TYPE(TypeTag, InteriorBoundaryData); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using ParentType = MpfaPorousMediumFluxVariablesCache; + + static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries); public: //! the constructor - MpfaPorousMediumFluxVariablesCache() - : AdvectionCache(), - isUpdated_(false) + PorousMediumFluxVariablesCacheImplementation() + : ParentType(), + isUpdated_(false), + interiorBoundaryInfoSelf_(false, -1) {} //! Returns whether or not this cache has been updated bool isUpdated() const { return isUpdated_; } - //! Sets the update status from outside. Allows an update of the cache specific - //! to processes that have solution dependent parameters, e.g. only updating - //! the diffusion transmissibilities leaving the advective ones untouched + //! Sets the update status from outside. This is used to only update + //! the cache once after solving the local system. When visiting an scvf + //! of the same interaction region again, the update is skipped. void setUpdateStatus(const bool status) { isUpdated_ = status; } -private: - bool isUpdated_; -}; - -// specialization for the case of advection & diffusion -template -class MpfaPorousMediumFluxVariablesCache - : public MpfaAdvectionCache, - public MpfaDiffusionCache -{ - using AdvectionCache = MpfaAdvectionCache; - using DiffusionCache = MpfaDiffusionCache; - - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - -public: - // the constructor - MpfaPorousMediumFluxVariablesCache() - : AdvectionCache(), - DiffusionCache(), - isUpdated_(false) - {} - - //! For compositional problems, neumann fluxes are not associated with a phase anymore - //! TODO: How to implement neumann fluxes for !useTpfa - template - void updatePhaseNeumannFlux(const SubControlVolumeFace &scvf, - const InteractionVolume& interactionVolume, - const unsigned int phaseIdx) {} - - //! TODO: How to implement neumann fluxes for !useTpfa - Scalar advectionNeumannFlux(const unsigned int phaseIdx) const - { return 0.0; } - - //! Returns whether or not this cache has been updated - bool isUpdated() const - { return isUpdated_; } - - //! Sets the update status from outside. Allows an update of the cache specific - //! to processes that have solution dependent parameters, e.g. only updating - //! the diffusion transmissibilities leaving the advective ones untouched - void setUpdateStatus(const bool status) + //! maybe update data on interior Dirichlet boundaries + template + void updateInteriorBoundaryData(const InteractionVolume& interactionVolume, + const SubControlVolumeFace &scvf) { - isUpdated_ = status; + if (enableInteriorBoundaries) + { + interiorBoundaryData_ = interactionVolume.interiorBoundaryData(); + + // check if the actual scvf is on an interior Dirichlet boundary + const auto scvfIdx = scvf.index(); + unsigned int indexInData = 0; + for (auto&& data : interiorBoundaryData_) + { + if (data.scvfIndex() == scvfIdx) + { + interiorBoundaryInfoSelf_ = std::make_pair(true, indexInData); + break; + } + + indexInData++; + } + } } -private: - bool isUpdated_; -}; + bool isInteriorBoundary() const + { return interiorBoundaryInfoSelf_.first; } -// specialization for the case of advection & heat conduction -template -class MpfaPorousMediumFluxVariablesCache - : public MpfaAdvectionCache, - public MpfaHeatConductionCache -{ - using AdvectionCache = MpfaAdvectionCache; - using HeatConductionCache = MpfaHeatConductionCache; + const std::vector& interiorBoundaryData() const + { return interiorBoundaryData_; } -public: - // the constructor - MpfaPorousMediumFluxVariablesCache() - : AdvectionCache(), - HeatConductionCache(), - isUpdated_(false) - {} - - //! Returns whether or not this cache has been updated - bool isUpdated() const - { return isUpdated_; } - - //! Sets the update status from outside. Allows an update of the cache specific - //! to processes that have solution dependent parameters, e.g. only updating - //! the diffusion transmissibilities leaving the advective ones untouched - void setUpdateStatus(const bool status) + const InteriorBoundaryData& interiorBoundaryDataSelf() const { - isUpdated_ = status; + assert(interiorBoundaryInfoSelf_.first && "Trying to obtain interior boundary data on a face that is not marked as such"); + assert(interiorBoundaryInfoSelf_.second != -1 && "The index to the interior boundary data of this face has not been set"); + return interiorBoundaryData_[interiorBoundaryInfoSelf_.second]; } private: + // indicates whether or not this cache has been fully updated bool isUpdated_; -}; - -// specialization for the case of advection, diffusion & heat conduction -template -class MpfaPorousMediumFluxVariablesCache - : public MpfaAdvectionCache, - public MpfaDiffusionCache, - public MpfaHeatConductionCache -{ - using AdvectionCache = MpfaAdvectionCache; - using DiffusionCache = MpfaDiffusionCache; - using HeatConductionCache = MpfaHeatConductionCache; - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - -public: - // the constructor - MpfaPorousMediumFluxVariablesCache() - : AdvectionCache(), - DiffusionCache(), - HeatConductionCache(), - isUpdated_(false) - {} - - //! TODO: How to implement neumann fluxes for !useTpfa when diffusion/heat conduction is active? - template - void updatePhaseNeumannFlux(const SubControlVolumeFace &scvf, - const InteractionVolume& interactionVolume, - const unsigned int phaseIdx) {} + // if this face is an interior Dirichlet boundary itself, store additional data + std::pair interiorBoundaryInfoSelf_; + // contains all the interior Dirichlet boundary data within the stencil of this face + std::vector interiorBoundaryData_; +}; - //! This method is needed for compatibility reasons - //! TODO: How to implement neumann fluxes for !useTpfa when diffusion is active? - Scalar advectionNeumannFlux(const unsigned int phaseIdx) const - { return 0.0; } +// specialization for the case of pure advection +template +class MpfaPorousMediumFluxVariablesCache : public MpfaAdvectionCache {}; - //! Returns whether or not this cache has been updated - bool isUpdated() const - { return isUpdated_; } +// specialization for the case of advection & diffusion +template +class MpfaPorousMediumFluxVariablesCache : public MpfaAdvectionCache, + public MpfaDiffusionCache {}; - //! Sets the update status from outside. Allows an update of the cache specific - //! to processes that have solution dependent parameters, e.g. only updating - //! the diffusion transmissibilities leaving the advective ones untouched - void setUpdateStatus(const bool status) - { - isUpdated_ = status; - } +// specialization for the case of advection & heat conduction +template +class MpfaPorousMediumFluxVariablesCache : public MpfaAdvectionCache, + public MpfaHeatConductionCache {}; -private: - bool isUpdated_; -}; +// specialization for the case of advection, diffusion & heat conduction +template +class MpfaPorousMediumFluxVariablesCache : public MpfaAdvectionCache, + public MpfaDiffusionCache, + public MpfaHeatConductionCache {}; } // end namespace