From 0dd5db6664e96a848a46a1a4876ee3bc68fd0bd9 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 Jan 2017 12:00:11 +0100 Subject: [PATCH 01/10] [mpfa] remove obsolete stencil class --- .../cellcentered/mpfa/stencils.hh | 41 ------------------- 1 file changed, 41 deletions(-) delete mode 100644 dumux/discretization/cellcentered/mpfa/stencils.hh diff --git a/dumux/discretization/cellcentered/mpfa/stencils.hh b/dumux/discretization/cellcentered/mpfa/stencils.hh deleted file mode 100644 index 734951864e..0000000000 --- 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 -- GitLab From f76ff2592bbdf8f2398778e8ef7232761b668276 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 Jan 2017 12:00:57 +0100 Subject: [PATCH 02/10] [mpfa][elemfluxvarcache] change throw statement --- .../cellcentered/mpfa/elementfluxvariablescache.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh index 01d8be6c03..098703639f 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 -- GitLab From 381105ebd29caf3ffc7ebde4190fa3dfb37ab69e Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 Jan 2017 12:25:39 +0100 Subject: [PATCH 03/10] [mpfa][iavol] remove brackets from forward declaration --- dumux/discretization/cellcentered/mpfa/interactionvolume.hh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/interactionvolume.hh index 8d8f233b95..729fe9ed94 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 -- GitLab From 3007a8e517aad648503bdf38286534ec004e2a7c Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 Jan 2017 12:25:57 +0100 Subject: [PATCH 04/10] [mpfa-o][localscvf] store pointer to scvf instead of copies The interaction volumes have to be instantiated with a valid fvGeometry etc and during assembly they only exist temporarily for the transmissibility calculations. Interaction volumes can now not be instantiated and stored globally anymore, but, why would one like to do that. If global flux var caching is enabled, then the transmissibilities are stored over the whole simulation time. This is the main thing one wants from them. --- .../mpfa/omethod/localsubcontrolentities.hh | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh index c3d8b4c98e..f3aa1a5ba7 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,13 +170,13 @@ public: { return scvfSeed_().faceType(); } GlobalPosition ip() const - { return ip_; } + { return scvf_().ipGlobal(); } GlobalPosition unitOuterNormal() const - { return normal_; } + { return scvf_().unitOuterNormal(); } Scalar area() const - { return area_; } + { return scvf_().area(); } bool boundary() const { return scvfSeed_().boundary(); } @@ -187,10 +185,11 @@ private: const LocalScvfSeed& scvfSeed_() const { return *seedPtr_; } + const SubControlVolumeFace& scvf_() const + { return *scvfPtr_; } + const LocalScvfSeed* seedPtr_; - GlobalPosition ip_; - GlobalPosition normal_; - Scalar area_; + const SubControlVolumeFace* scvfPtr_; }; } // end namespace -- GitLab From b9d97eb41889bc392d08301f37e4a561864f8d23 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 23 Jan 2017 19:02:06 +0100 Subject: [PATCH 05/10] [mpfa] enable interior boundary handling --- .../cellcentered/mpfa/darcyslaw.hh | 135 ++++- .../mpfa/elementfluxvariablescache.hh | 2 +- .../mpfa/elementvolumevariables.hh | 6 +- .../cellcentered/mpfa/facetypes.hh | 4 +- .../cellcentered/mpfa/fickslaw.hh | 154 ++++- .../mpfa/fluxvariablescachefiller.hh | 533 ++++++++++++++---- .../mpfa/fluxvariablescachefillerbase.hh | 290 +++------- .../cellcentered/mpfa/fourierslaw.hh | 95 +++- .../cellcentered/mpfa/fvelementgeometry.hh | 14 +- .../cellcentered/mpfa/globalfvgeometrybase.hh | 343 +++++++---- .../mpfa/globalinteractionvolumeseedsbase.hh | 4 +- .../mpfa/globalvolumevariables.hh | 4 +- .../cellcentered/mpfa/helper.hh | 89 ++- .../mpfa/interactionvolumebase.hh | 5 +- .../cellcentered/mpfa/interiorboundarydata.hh | 166 ++++++ .../lmethod/globalinteractionvolumeseeds.hh | 24 +- .../cellcentered/mpfa/lmethod/helper.hh | 67 ++- .../mpfa/lmethod/interactionregions.hh | 43 +- .../mpfa/lmethod/interactionvolume.hh | 19 +- .../omethod/globalinteractionvolumeseeds.hh | 10 +- .../cellcentered/mpfa/omethod/helper.hh | 239 +++++--- .../mpfa/omethod/interactionvolume.hh | 498 +++++++++++----- .../mpfa/omethod/interactionvolumeseed.hh | 18 +- .../mpfa/omethod/localsubcontrolentities.hh | 12 +- .../mpfa/omethodfps/interactionvolume.hh | 41 +- .../implicit/cellcentered/mpfa/properties.hh | 4 + .../cellcentered/mpfa/propertydefaults.hh | 13 + dumux/implicit/problem.hh | 16 + .../implicit/fluxvariablescache.hh | 363 ++++++------ 29 files changed, 2237 insertions(+), 974 deletions(-) create mode 100644 dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh index c024ebbf6a..d82738c41d 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,13 +144,75 @@ public: h -= rho*(g*x); } + flux += tij[localIdx++]*h; } - if (useTpfaBoundary) - return flux; - else - return flux + fluxVarsCache.advectionNeumannFlux(phaseIdx); + // if no interior boundaries are present, return the flux + if (!enableInteriorBoundaries) + return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx); + + ////////////////////////////////////////////////////////////////// + // Handle interior boundaries + ////////////////////////////////////////////////////////////////// + + // For active facet coupling we will have to transform the interior flux vector + const auto& cij = fluxVarsCache.advectionCij(); + + // The vector of interior neumann fluxes + DynamicVector facetCouplingFluxes(cij.size(), 0.0); + for (auto&& data : fluxVarsCache.interiorBoundaryData()) + { + // Add additional Dirichlet fluxes for interior Dirichlet faces + if (data.faceType() == MpfaFaceTypes::interiorDirichlet) + { + Scalar h = data.facetVolVars(fvGeometry).pressure(phaseIdx); + + if (gravity) + { + const auto x = fvGeometry.scvf(data.scvfIndex()).ipGlobal(); + const auto g = problem.gravityAtPos(x); + + h -= rho*(g*x); + } + + // The transmissibilities of interior dirichlet boundaries are placed at the end + // So we simply keep incrementing the local index + flux += tij[localIdx + data.localIndexInInteractionVolume()]*h; + } + + // add neumann fluxes for interior Neumann faces if facet coupling is active + if (facetCoupling && data.faceType() == MpfaFaceTypes::interiorNeumann) + { + // get the volvars of the actual interior neumann face + const auto facetVolVars = data.facetVolVars(fvGeometry); + + // get the scvf corresponding to actual interior neumann face + const auto& curScvf = fvGeometry.scvf(data.scvfIndex()); + + // calculate "lekage factor" + const auto n = curScvf.unitOuterNormal(); + const auto v = [&] () + { + auto res = n; + res *= -0.5*facetVolVars.extrusionFactor(); + res -= curScvf.ipGlobal(); + res += curScvf.facetCorner(); + res /= res.two_norm2(); + return res; + } (); + + // add value to vector of interior neumann fluxes + facetCouplingFluxes[data.localIndexInInteractionVolume()] -= facetVolVars.pressure(phaseIdx)* + MpfaHelper::nT_M_v(n, facetVolVars.permeability(), v); + } + } + + // return overall resulting flux + const Scalar interiorNeumannFlux = facetCoupling ? cij*facetCouplingFluxes : 0.0; + return useTpfaBoundary ? + flux + interiorNeumannFlux : + flux + interiorNeumannFlux + fluxVarsCache.advectionNeumannFlux(phaseIdx); } static Stencil stencil(const Problem& problem, @@ -122,7 +223,7 @@ public: const auto& globalFvGeometry = problem.model().globalFvGeometry(); // return the scv (element) indices in the interaction region - if (globalFvGeometry.scvfTouchesBoundary(scvf)) + if (globalFvGeometry.touchesInteriorOrDomainBoundary(scvf)) return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices(); else return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices(); diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh index 098703639f..2bae60663e 100644 --- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh +++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh @@ -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 6620dee3b0..e0470acfc4 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 87d526e167..c278de9e39 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 fb27d52806..09ad7b2801 100644 --- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh @@ -46,19 +46,28 @@ 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 static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa; @@ -76,8 +85,55 @@ public: const auto& volVarsStencil = fluxVarsCache.diffusionVolVarsStencil(phaseIdx, compIdx); const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx); + const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary(); + // For interior Neumann boundaries when using Tpfa for Neumann boundary conditions, we simply + // return the user-specified flux. Note that for compositional models we attribute the influxes + // to the major components, thus we do it per phase in Darcy's law. However, for single-phasic models + // wesolve the phase mass balance equation AND the transport equation, thus, in that case we incorporate + // the Neumann BCs here. We assume compIdx = eqIdx. + // Note that this way of including interior Neumann fluxes fails for mpnc models where n != m. + if (numPhases == 1 + && isInteriorBoundary + && useTpfaBoundary + && fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann) + return scvf.area()* + elemVolVars[scvf.insideScvIdx()].extrusionFactor()* + problem.neumann(element, + fvGeometry, + elemVolVars, + scvf)[compIdx]; + + // get the scaling factor for the effective diffusive fluxes - auto effFactor = calculateEffectiveDiffusivityFactor(elemVolVars, scvf, phaseIdx); + const auto effFactor = [&] () + { + // Treat interior boundaries differently + if (isInteriorBoundary) + { + const auto& data = fluxVarsCache.interiorBoundaryDataSelf(); + + // interpolate as usual for interior Neumann faces without facet coupling + if (data.faceType() == MpfaFaceTypes::interiorNeumann && !facetCoupling) + return calculateEffectiveDiffusivityFactor(elemVolVars, scvf, phaseIdx); + // use harmonic mean between the interior and the facet volvars + else + { + const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; + const auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), + insideVolVars.saturation(phaseIdx), + /*Diffusion coefficient*/ 1.0); + + const auto facetVolVars = data.facetVolVars(fvGeometry); + const auto outsideFactor = EffDiffModel::effectiveDiffusivity(facetVolVars.porosity(), + facetVolVars.saturation(phaseIdx), + /*Diffusion coefficient*/ 1.0); + + return harmonicMean(factor, outsideFactor); + } + } + else + return calculateEffectiveDiffusivityFactor(elemVolVars, scvf, phaseIdx); + } (); // if factor is zero, the flux will end up zero anyway if (effFactor == 0.0) @@ -90,14 +146,88 @@ public: auto getRho = [useMoles, phaseIdx] (const VolumeVariables& volVars) { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); }; + // calculate the density at the interface + const auto rho = [&] () + { + // maybe use the density of the interior BC on the facet + if (isInteriorBoundary) + { + const auto& data = fluxVarsCache.interiorBoundaryDataSelf(); + + if (facetCoupling || data.faceType() == MpfaFaceTypes::interiorDirichlet) + return useMoles ? data.facetVolVars(fvGeometry).molarDensity(phaseIdx) : data.facetVolVars(fvGeometry).density(phaseIdx); + else + return interpolateDensity(elemVolVars, scvf, getRho); + } + else + return interpolateDensity(elemVolVars, scvf, getRho); + } (); + // calculate Tij*xj Scalar flux(0.0); unsigned int localIdx = 0; for (const auto volVarIdx : volVarsStencil) flux += tij[localIdx++]*getX(elemVolVars[volVarIdx]); - // return effective mass flux - return flux*interpolateDensity(elemVolVars, scvf, getRho)*effFactor; + // if no interior boundaries are present, return effective mass flux + if (!enableInteriorBoundaries) + return useTpfaBoundary ? flux*rho*effFactor : flux*rho*effFactor + fluxVarsCache.componentNeumannFlux(compIdx); + + ////////////////////////////////////////////////////////////////// + // Handle interior boundaries + ////////////////////////////////////////////////////////////////// + + // get coefficients to transform the vector of interior neumann boundary conditions + const auto& cij = fluxVarsCache.diffusionCij(phaseIdx, compIdx); + + // The vector of interior neumann fluxes + DynamicVector facetCouplingFluxes(cij.size(), 0.0); + for (auto&& data : fluxVarsCache.interiorBoundaryData()) + { + // Add additional Dirichlet fluxes for interior Dirichlet faces + if (data.faceType() == MpfaFaceTypes::interiorDirichlet) + { + // The transmissibilities of interior dirichlet boundaries are placed at the end + // So we simply keep incrementing the local index + const auto x = useMoles ? + data.facetVolVars(fvGeometry).moleFraction(phaseIdx, compIdx) : + data.facetVolVars(fvGeometry).massFraction(phaseIdx, compIdx); + flux += tij[localIdx + data.localIndexInInteractionVolume()]*x; + } + + // add neumann fluxes for interior Neumann faces + if (facetCoupling && data.faceType() == MpfaFaceTypes::interiorNeumann) + { + // get the scvf corresponding to actual interior boundary face + const auto& curScvf = fvGeometry.scvf(data.scvfIndex()); + + // get the volvars of the actual interior neumann face + const auto facetVolVars = data.facetVolVars(fvGeometry); + + // calculate "lekage factor" + const auto n = curScvf.unitOuterNormal(); + const auto v = [&] () + { + auto res = n; + res *= -0.5*facetVolVars.extrusionFactor(); + res -= curScvf.ipGlobal(); + res += curScvf.facetCorner(); + res /= res.two_norm2(); + return res; + } (); + + // add value to vector of interior neumann fluxes + facetCouplingFluxes[data.localIndexInInteractionVolume()] += MpfaHelper::nT_M_v(n, + facetVolVars.diffusionCoefficient(phaseIdx, compIdx), + v); + } + } + + // return overall resulting flux + const Scalar interiorNeumannFlux = facetCoupling ? cij*facetCouplingFluxes : 0.0; + return useTpfaBoundary ? + flux + interiorNeumannFlux : + flux + interiorNeumannFlux + fluxVarsCache.componentNeumannFlux(compIdx); } static Stencil stencil(const Problem& problem, @@ -108,7 +238,7 @@ public: const auto& globalFvGeometry = problem.model().globalFvGeometry(); // return the scv (element) indices in the interaction region - if (globalFvGeometry.scvfTouchesBoundary(scvf)) + if (globalFvGeometry.touchesInteriorOrDomainBoundary(scvf)) return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices(); else return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices(); @@ -141,19 +271,25 @@ private: const unsigned int phaseIdx) { const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; - auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), + const auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), /*Diffusion coefficient*/ 1.0); if (!scvf.boundary()) { - const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; - auto outsideFactor = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), + // interpret outside factor as arithmetic mean + Scalar outsideFactor = 0.0; + for (auto outsideIdx : scvf.outsideScvIndices()) + { + const auto& outsideVolVars = elemVolVars[outsideIdx]; + outsideFactor += EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), /*Diffusion coefficient*/ 1.0); + } + outsideFactor /= scvf.outsideScvIndices().size(); // use the harmonic mean of the two - factor = harmonicMean(factor, outsideFactor); + return harmonicMean(factor, outsideFactor); } return factor; diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh index 6b05dd47b2..0ca84820d8 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 9790e04e1d..9ef0820f3a 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 fcca346bcc..5c770db6b7 100644 --- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh @@ -43,16 +43,28 @@ 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 static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa; @@ -68,14 +80,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(); + + ////////////////////////////////////////////////////////////////// + // 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 overall resulting flux + const Scalar interiorNeumannFlux = facetCoupling ? cij*facetCouplingFluxes : 0.0; + return useTpfaBoundary ? + flux + interiorNeumannFlux : + flux + interiorNeumannFlux + fluxVarsCache.heatNeumannFlux(); } static Stencil stencil(const Problem& problem, @@ -86,7 +175,7 @@ public: const auto& globalFvGeometry = problem.model().globalFvGeometry(); // return the scv (element) indices in the interaction region - if (globalFvGeometry.scvfTouchesBoundary(scvf)) + if (globalFvGeometry.touchesInteriorOrDomainBoundary(scvf)) return globalFvGeometry.boundaryInteractionVolumeSeed(scvf).globalScvIndices(); else return globalFvGeometry.interactionVolumeSeed(scvf).globalScvIndices(); diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh index 05e95ae2ee..e79a41db5e 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 b275a73873..fe3f89a5f4 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 cadbd60826..06b4056d92 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 3869317a56..54ce225e4d 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 23fe038d1a..5f27a14dbd 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/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh index ea3ced18cf..af8c2d793f 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 0000000000..c1cc25b6b0 --- /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 7b36d027cd..8a3da8af4f 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 f11cb30f2e..f5d88a6897 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 2a9a4287fc..fc4c039ea1 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 5cf4968104..51e7929705 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 f009866d9c..22effd4f00 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 12ae0c9781..9c7aca1b01 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 e49ba03208..9450d4da51 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 3a4c5ccf35..a3b11e1add 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 f3aa1a5ba7..959706e166 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh @@ -170,24 +170,24 @@ public: { return scvfSeed_().faceType(); } GlobalPosition ip() const - { return scvf_().ipGlobal(); } + { return globalScvf().ipGlobal(); } GlobalPosition unitOuterNormal() const - { return scvf_().unitOuterNormal(); } + { return globalScvf().unitOuterNormal(); } Scalar area() const - { return scvf_().area(); } + { return globalScvf().area(); } bool boundary() const { return scvfSeed_().boundary(); } + const SubControlVolumeFace& globalScvf() const + { return *scvfPtr_; } + private: const LocalScvfSeed& scvfSeed_() const { return *seedPtr_; } - const SubControlVolumeFace& scvf_() const - { return *scvfPtr_; } - const LocalScvfSeed* seedPtr_; const SubControlVolumeFace* scvfPtr_; }; diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh index 2262c42f33..6a905987d6 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/implicit/cellcentered/mpfa/properties.hh b/dumux/implicit/cellcentered/mpfa/properties.hh index 274783410e..0118bbd4d8 100644 --- a/dumux/implicit/cellcentered/mpfa/properties.hh +++ b/dumux/implicit/cellcentered/mpfa/properties.hh @@ -49,6 +49,10 @@ 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(EnableInteriorBoundaries); //! Enables or disables the use of internal boundaries on interior scvfs +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 3278dafaef..8276f91ef2 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,18 @@ public: // By default, we use tpfa on the boundaries SET_BOOL_PROP(CCMpfaModel, UseTpfaBoundary, true); +// By default, we disable interior boundaries +SET_BOOL_PROP(CCMpfaModel, EnableInteriorBoundaries, false); + +// 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/problem.hh b/dumux/implicit/problem.hh index 9bdc259051..1008db47d6 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -1018,6 +1018,22 @@ public: std::vector getAdditionalDofDependencies(IndexType globalIdx) const { return std::vector(); } + /*! + * \brief Function to mark intersections as interior boundaries. This functionality is only + * available for models using a cell-centered Mpfa scheme. 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 + typename std::enable_if::type + 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/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh index 0b5a82bd1e..dbc18f3e6b 100644 --- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh +++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh @@ -174,64 +174,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 +252,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 +336,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 +363,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 +379,8 @@ public: private: // Quantities associated with heat conduction Stencil heatConductionVolVarsStencil_; - TransmissibilityVector heatConductionTij_; + CoefficientVector heatConductionTij_; + CoefficientVector heatConductionCij_; Scalar heatNeumannFlux_; }; @@ -349,168 +390,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_; -}; - -// specialization for the case of advection & heat conduction -template -class MpfaPorousMediumFluxVariablesCache - : public MpfaAdvectionCache, - public MpfaHeatConductionCache -{ - using AdvectionCache = MpfaAdvectionCache; - using HeatConductionCache = MpfaHeatConductionCache; - -public: - // the constructor - MpfaPorousMediumFluxVariablesCache() - : AdvectionCache(), - HeatConductionCache(), - isUpdated_(false) - {} + bool isInteriorBoundary() const + { return interiorBoundaryInfoSelf_.first; } - //! Returns whether or not this cache has been updated - bool isUpdated() const - { return isUpdated_; } + const std::vector& interiorBoundaryData() const + { return interiorBoundaryData_; } - //! 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 -- GitLab From 11ae6bb119165ff4e92fb1912e81bbbd7194011d Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 26 Jan 2017 19:53:13 +0100 Subject: [PATCH 06/10] [fluxvars] specialize upwind scheme to different discretization schemes --- dumux/discretization/fluxvariablesbase.hh | 133 ++---- dumux/discretization/upwindscheme.hh | 401 ++++++++++++++++++ .../implicit/fluxvariables.hh | 28 +- 3 files changed, 437 insertions(+), 125 deletions(-) create mode 100644 dumux/discretization/upwindscheme.hh diff --git a/dumux/discretization/fluxvariablesbase.hh b/dumux/discretization/fluxvariablesbase.hh index 216edff503..c6d0012011 100644 --- a/dumux/discretization/fluxvariablesbase.hh +++ b/dumux/discretization/fluxvariablesbase.hh @@ -24,6 +24,9 @@ #define DUMUX_DISCRETIZATION_FLUXVARIABLESBASE_HH #include +#include +#include +#include namespace Dumux { @@ -31,16 +34,16 @@ namespace Dumux namespace Properties { // forward declaration -NEW_PROP_TAG(ImplicitUpwindWeight); } - /*! * \ingroup Discretization - * \brief Base class for the flux variables - * Actual flux variables inherit from this class + * \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 +51,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 +71,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 +91,34 @@ 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) - { - 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)); - } + Stencil computeStencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + { DUNE_THROW(Dune::InvalidStateException, "computeStencil() routine is not provided by the implementation."); } 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_; }; +/*! + * \ingroup Discretization + * \brief The flux variables base class class + * The upwind scheme is chosen depending on the discretization method + */ +template +using FluxVariablesBase = FluxVariablesBaseImplementation>; + } // end namespace #endif diff --git a/dumux/discretization/upwindscheme.hh b/dumux/discretization/upwindscheme.hh new file mode 100644 index 0000000000..94d1cdc60d --- /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/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh index d2840c174a..66ffbdb8ab 100644 --- a/dumux/porousmediumflow/implicit/fluxvariables.hh +++ b/dumux/porousmediumflow/implicit/fluxvariables.hh @@ -57,10 +57,9 @@ 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; @@ -100,7 +99,7 @@ public: phaseIdx, this->elemFluxVarsCache()); - return this->upwindScheme(flux, phaseIdx, upwindTerm); + return this->applyUpwindScheme(upwindTerm, flux, phaseIdx); } Stencil computeStencil(const Problem& problem, @@ -113,10 +112,9 @@ public: // 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; @@ -164,7 +162,7 @@ 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) @@ -209,10 +207,9 @@ 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; @@ -258,7 +255,7 @@ public: advFluxCached_.set(phaseIdx, true); } - return this->upwindScheme(advPreFlux_[phaseIdx], phaseIdx, upwindTerm); + return this->applyUpwindScheme(upwindTerm, advPreFlux_[phaseIdx], phaseIdx); } Scalar heatConductionFlux() @@ -303,10 +300,9 @@ 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; @@ -353,7 +349,7 @@ 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) -- GitLab From fda8a3e5cd0d0ff95bae5d93c15d218d0329ffba Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 26 Jan 2017 19:59:55 +0100 Subject: [PATCH 07/10] [tpfa] enable interior boundaries --- .../cellcentered/tpfa/fvelementgeometry.hh | 19 +++++++++++---- .../cellcentered/tpfa/globalfvgeometry.hh | 24 ++++++++++++------- .../cellcentered/tpfa/subcontrolvolumeface.hh | 15 +++++++++--- .../implicit/cellcentered/mpfa/properties.hh | 1 - .../cellcentered/mpfa/propertydefaults.hh | 3 --- dumux/implicit/cellcentered/properties.hh | 1 + .../implicit/cellcentered/propertydefaults.hh | 3 +++ dumux/implicit/model.hh | 13 ++++++---- dumux/implicit/problem.hh | 7 +++--- 9 files changed, 58 insertions(+), 28 deletions(-) diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh index 6d9efed3fd..64d46973be 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 59a5758639..7e48089aea 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 d77fcb98c6..30a287bfca 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/implicit/cellcentered/mpfa/properties.hh b/dumux/implicit/cellcentered/mpfa/properties.hh index 0118bbd4d8..f78c330e81 100644 --- a/dumux/implicit/cellcentered/mpfa/properties.hh +++ b/dumux/implicit/cellcentered/mpfa/properties.hh @@ -49,7 +49,6 @@ 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(EnableInteriorBoundaries); //! Enables or disables the use of internal boundaries on interior scvfs 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) diff --git a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh index 8276f91ef2..483f831e83 100644 --- a/dumux/implicit/cellcentered/mpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/mpfa/propertydefaults.hh @@ -109,9 +109,6 @@ public: // By default, we use tpfa on the boundaries SET_BOOL_PROP(CCMpfaModel, UseTpfaBoundary, true); -// By default, we disable interior boundaries -SET_BOOL_PROP(CCMpfaModel, EnableInteriorBoundaries, false); - // By default, interior boundaries are static SET_BOOL_PROP(CCMpfaModel, MpfaFacetCoupling, false); diff --git a/dumux/implicit/cellcentered/properties.hh b/dumux/implicit/cellcentered/properties.hh index bea7da122b..c5757fcbc3 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 20cadb848c..702259c31e 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 a067d55a69..6124493f88 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 1008db47d6..7b862f6dde 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -1019,8 +1019,8 @@ public: { return std::vector(); } /*! - * \brief Function to mark intersections as interior boundaries. This functionality is only - * available for models using a cell-centered Mpfa scheme. The corresponding boundary + * \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 @@ -1030,8 +1030,7 @@ public: * Per default we don't have interior boundaries */ template - typename std::enable_if::type - isInteriorBoundary(const Element& element, const Intersection& intersection) const + bool isInteriorBoundary(const Element& element, const Intersection& intersection) const { return false; } /*! -- GitLab From b2f7e854a6df123720fc6808194008e0c67af2e7 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 26 Jan 2017 21:38:02 +0100 Subject: [PATCH 08/10] [stencil] Improve flux stencils Fix the flux stencil for every method independent of the stencils used in the different flux laws. The matrix entries have to be reserved for all possible occuring stencils if the method has parameter or solution dependent stencils, otherwise the matrix needs to be rebuild every time. This implements this unified flux stencil for box, cc, ccmpfa. Remove a lot of duplicate code. Separate the physical laws from the dicretization stencil. --- dumux/discretization/box/darcyslaw.hh | 9 -- dumux/discretization/box/fickslaw.hh | 8 -- dumux/discretization/box/fourierslaw.hh | 8 -- .../cellcentered/mpfa/darcyslaw.hh | 14 -- .../cellcentered/mpfa/fickslaw.hh | 15 -- .../cellcentered/mpfa/fourierslaw.hh | 15 -- .../cellcentered/tpfa/darcyslaw.hh | 19 --- .../cellcentered/tpfa/fickslaw.hh | 12 -- .../cellcentered/tpfa/fourierslaw.hh | 12 -- dumux/discretization/fluxstencil.hh | 133 ++++++++++++++++++ dumux/discretization/fluxvariablesbase.hh | 39 ++--- .../implicit/fluxvariables.hh | 80 ----------- 12 files changed, 156 insertions(+), 208 deletions(-) create mode 100644 dumux/discretization/fluxstencil.hh diff --git a/dumux/discretization/box/darcyslaw.hh b/dumux/discretization/box/darcyslaw.hh index d432223e5f..f4b1bc04b9 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/fickslaw.hh b/dumux/discretization/box/fickslaw.hh index fc2d13eedb..881072dd3d 100644 --- a/dumux/discretization/box/fickslaw.hh +++ b/dumux/discretization/box/fickslaw.hh @@ -65,7 +65,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 +133,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 c7288999b3..d77e89a689 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 d82738c41d..3aedae63d4 100644 --- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh @@ -215,20 +215,6 @@ public: flux + interiorNeumannFlux + fluxVarsCache.advectionNeumannFlux(phaseIdx); } - 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(); - } - private: static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf, diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh index 09ad7b2801..6aaf11e188 100644 --- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh @@ -61,7 +61,6 @@ class FicksLawImplementation 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); @@ -230,20 +229,6 @@ public: flux + interiorNeumannFlux + fluxVarsCache.componentNeumannFlux(compIdx); } - 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(); - } - private: template static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars, diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh index 5c770db6b7..03d00e0db4 100644 --- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh +++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh @@ -57,7 +57,6 @@ class FouriersLawImplementation 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); @@ -166,20 +165,6 @@ public: flux + interiorNeumannFlux : flux + interiorNeumannFlux + fluxVarsCache.heatNeumannFlux(); } - - 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 diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh index 63cedea3ba..3a31437394 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 8fcd9ce4a0..7b05340bfa 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 1806ae63da..c46876547a 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/fluxstencil.hh b/dumux/discretization/fluxstencil.hh new file mode 100644 index 0000000000..72a49236c3 --- /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 c6d0012011..b20dd3f15e 100644 --- a/dumux/discretization/fluxvariablesbase.hh +++ b/dumux/discretization/fluxvariablesbase.hh @@ -26,15 +26,23 @@ #include #include #include +#include #include namespace Dumux { -namespace Properties -{ -// forward declaration -} +template +class FluxVariablesBaseImplementation; + +/*! + * \ingroup Discretization + * \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 @@ -42,7 +50,7 @@ namespace Properties * \param TypeTag The type tag * \param UpwindScheme The type used for the upwinding of the advective fluxes */ -template +template class FluxVariablesBaseImplementation { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); @@ -99,8 +107,15 @@ public: return UpwindScheme::apply(*this, upwindTerm, flux, phaseIdx); } - Stencil computeStencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) - { DUNE_THROW(Dune::InvalidStateException, "computeStencil() routine is not provided by the implementation."); } + static Stencil computeStencil(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf) + { + //! Give the upwind scheme access to the cached variables + //! Forward to the discretization specific implementation + return FluxStencil::stencil(problem, element, fvGeometry, scvf); + } private: const Problem* problemPtr_; //! Pointer to the problem @@ -111,14 +126,6 @@ private: const ElementFluxVariablesCache* elemFluxVarsCachePtr_; }; -/*! - * \ingroup Discretization - * \brief The flux variables base class class - * The upwind scheme is chosen depending on the discretization method - */ -template -using FluxVariablesBase = FluxVariablesBaseImplementation>; - -} // end namespace +} // end namespace Dumux #endif diff --git a/dumux/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh index 66ffbdb8ab..a60ad9901c 100644 --- a/dumux/porousmediumflow/implicit/fluxvariables.hh +++ b/dumux/porousmediumflow/implicit/fluxvariables.hh @@ -65,7 +65,6 @@ class PorousMediumFluxVariablesImpl : public FluxVa 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); @@ -101,12 +100,6 @@ public: 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); } }; @@ -120,7 +113,6 @@ class PorousMediumFluxVariablesImpl : public FluxVar 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); @@ -177,28 +169,6 @@ public: 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 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_; @@ -215,7 +185,6 @@ class PorousMediumFluxVariablesImpl : public FluxVar 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); @@ -269,29 +238,6 @@ public: 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; - } - private: //! simple caching if advection flux is used twice with different upwind function std::bitset advFluxCached_; @@ -308,7 +254,6 @@ class PorousMediumFluxVariablesImpl : public FluxVari 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); @@ -375,31 +320,6 @@ public: 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; - } - private: //! simple caching if advection flux is used twice with different upwind function std::bitset advFluxCached_; -- GitLab From cc139828193d216ffce5cdd62f6ae0eaf8a73977 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 26 Jan 2017 22:28:15 +0100 Subject: [PATCH 09/10] [box] Fix forgotten change to new interface for flux cache --- dumux/discretization/box/elementfluxvariablescache.hh | 4 ++-- dumux/discretization/box/fickslaw.hh | 2 -- dumux/porousmediumflow/implicit/fluxvariablescache.hh | 11 +++++------ 3 files changed, 7 insertions(+), 10 deletions(-) diff --git a/dumux/discretization/box/elementfluxvariablescache.hh b/dumux/discretization/box/elementfluxvariablescache.hh index e0b699c9f9..97fb569f4b 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 881072dd3d..0ae115f46c 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); diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh index dbc18f3e6b..1e01be146b 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); } -- GitLab From c11ede3d329413ad903f7a2e77b14da958eeb5f3 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 26 Jan 2017 22:28:50 +0100 Subject: [PATCH 10/10] [pmfluxvariables][cleanup] directly return --- .../implicit/fluxvariables.hh | 56 +++++++++---------- 1 file changed, 26 insertions(+), 30 deletions(-) diff --git a/dumux/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh index a60ad9901c..87bd43686a 100644 --- a/dumux/porousmediumflow/implicit/fluxvariables.hh +++ b/dumux/porousmediumflow/implicit/fluxvariables.hh @@ -159,14 +159,13 @@ public: 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()); } private: @@ -229,13 +228,12 @@ public: Scalar heatConductionFlux() { - Scalar flux = HeatConductionType::flux(this->problem(), - this->element(), - this->fvGeometry(), - this->elemVolVars(), - this->scvFace(), - this->elemFluxVarsCache()); - return flux; + return HeatConductionType::flux(this->problem(), + this->element(), + this->fvGeometry(), + this->elemVolVars(), + this->scvFace(), + this->elemFluxVarsCache()); } private: @@ -299,25 +297,23 @@ public: 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; + return HeatConductionType::flux(this->problem(), + this->element(), + this->fvGeometry(), + this->elemVolVars(), + this->scvFace(), + this->elemFluxVarsCache()); } private: -- GitLab