diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh
index c3b3278d1eb8ba25e566ad0bc1d1d24f2ee5d633..2b48bbb80e88fdcb21d90dfac0f6a0a0184b744c 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh
@@ -7,448 +7,38 @@
/*!
* \file
* \ingroup FreeFlowPorousMediumCoupling
- * \copydoc Dumux::FreeFlowPorousMediumCouplingConditions
+ * \brief Coupling conditions specialized for different discretization schemes
*/
#ifndef DUMUX_MD_FREEFLOW_POROUSMEDIUM_COUPLINGCONDITIONS_HH
#define DUMUX_MD_FREEFLOW_POROUSMEDIUM_COUPLINGCONDITIONS_HH
-#include <numeric>
+#include <dumux/discretization/method.hh>
-#include <dune/common/exceptions.hh>
-
-#include <dumux/common/properties.hh>
-#include <dumux/common/math.hh>
-
-#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
-
-#include <dumux/flux/darcyslaw_fwd.hh>
-#include <dumux/flux/fickslaw_fwd.hh>
-#include <dumux/flux/forchheimerslaw_fwd.hh>
-
-#include <dumux/multidomain/boundary/freeflowporousmedium/traits.hh>
-#include <dumux/multidomain/boundary/freeflowporousmedium/indexhelper.hh>
+#include "couplingconditions_staggered_cctpfa.hh"
namespace Dumux {
-/*!
- * \ingroup FreeFlowPorousMediumCoupling
- * \brief This structs holds a set of options which allow to modify the Stokes-Darcy
- * coupling mechanism during runtime.
- */
-struct FreeFlowPorousMediumCouplingOptions
-{
- /*!
- * \brief Defines which kind of averanging of diffusion coefficients
- * (moleculat diffusion or thermal conductance) at the interface
- * between free flow and porous medium shall be used.
- */
- enum class DiffusionCoefficientAveragingType
- {
- harmonic, arithmetic, ffOnly, pmOnly
- };
-
- /*!
- * \brief Convenience function to convert user input given as std::string to the corresponding enum class used for choosing the type
- * of averaging of the diffusion/conduction parameter at the interface between the two domains.
- */
- static DiffusionCoefficientAveragingType stringToEnum(DiffusionCoefficientAveragingType, const std::string& diffusionCoefficientAveragingType)
- {
- if (diffusionCoefficientAveragingType == "Harmonic")
- return DiffusionCoefficientAveragingType::harmonic;
- else if (diffusionCoefficientAveragingType == "Arithmetic")
- return DiffusionCoefficientAveragingType::arithmetic;
- else if (diffusionCoefficientAveragingType == "FreeFlowOnly")
- return DiffusionCoefficientAveragingType::ffOnly;
- else if (diffusionCoefficientAveragingType == "PorousMediumOnly")
- return DiffusionCoefficientAveragingType::pmOnly;
- else
- DUNE_THROW(Dune::IOError, "Unknown DiffusionCoefficientAveragingType");
- }
-
-};
+#ifndef DOXYGEN
+namespace FreeFlowPorousMediumDetail {
-template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional>
-class FreeFlowPorousMediumCouplingConditionsImplementation;
+// declaration (specialize for different discretization types)
+template<class MDTraits, class CouplingManager,
+ class DiscFFMomentum = typename MDTraits::template SubDomain<CouplingManager::freeFlowMomentumIndex>::GridGeometry::DiscretizationMethod,
+ class DiscFFMass = typename MDTraits::template SubDomain<CouplingManager::freeFlowMassIndex>::GridGeometry::DiscretizationMethod,
+ class DiscPM = typename MDTraits::template SubDomain<CouplingManager::porousMediumIndex>::GridGeometry::DiscretizationMethod
+ >
+struct FreeFlowPorousMediumCouplingConditionsSelector;
-/*!
-* \ingroup FreeFlowPorousMediumCoupling
-* \brief Data for the coupling of a Darcy model (cell-centered finite volume)
-* with a (Navier-)Stokes model (staggerd grid).
-*/
template<class MDTraits, class CouplingManager>
-using FreeFlowPorousMediumCouplingConditions
- = FreeFlowPorousMediumCouplingConditionsImplementation<
- MDTraits, CouplingManager,
- GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(),
- (GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::numFluidComponents() > 1)
- >;
-
-/*!
- * \ingroup FreeFlowPorousMediumCoupling
- * \brief A base class which provides some common methods used for Stokes-Darcy coupling.
- */
-template<class MDTraits, class CouplingManager>
-class FreeFlowPorousMediumCouplingConditionsImplementationBase
-{
- using Scalar = typename MDTraits::Scalar;
-
- template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
- template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
- template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
- template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
- template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
- template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
- template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
- template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
- template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
- template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
- template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>;
- template<std::size_t id> using ModelTraits = GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>;
- template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate;
- template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector;
-
-public:
- static constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex;
- static constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex;
- static constexpr auto porousMediumIndex = CouplingManager::porousMediumIndex;
-private:
-
- using AdvectionType = GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::AdvectionType>;
- using DarcysLaw = Dumux::DarcysLaw<SubDomainTypeTag<porousMediumIndex>>;
- using ForchheimersLaw = Dumux::ForchheimersLaw<SubDomainTypeTag<porousMediumIndex>>;
-
- static constexpr bool adapterUsed = ModelTraits<porousMediumIndex>::numFluidPhases() > 1;
- using IndexHelper = FreeFlowPorousMediumCoupling::IndexHelper<freeFlowMassIndex, porousMediumIndex, FluidSystem<freeFlowMassIndex>, adapterUsed>;
-
- static constexpr int enableEnergyBalance = GetPropType<SubDomainTypeTag<freeFlowMassIndex>, Properties::ModelTraits>::enableEnergyBalance();
- static_assert(GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::ModelTraits>::enableEnergyBalance() == enableEnergyBalance,
- "All submodels must both be either isothermal or non-isothermal");
-
- static_assert(FreeFlowPorousMediumCoupling::IsSameFluidSystem<FluidSystem<freeFlowMassIndex>,
- FluidSystem<porousMediumIndex>>::value,
- "All submodels must use the same fluid system");
-
- using DiffusionCoefficientAveragingType = typename FreeFlowPorousMediumCouplingOptions::DiffusionCoefficientAveragingType;
-
-public:
- /*!
- * \brief Returns the corresponding phase index needed for coupling.
- */
- template<std::size_t i>
- static constexpr auto couplingPhaseIdx(Dune::index_constant<i> id, int coupledPhaseIdx = 0)
- { return IndexHelper::couplingPhaseIdx(id, coupledPhaseIdx); }
-
- /*!
- * \brief Returns the corresponding component index needed for coupling.
- */
- template<std::size_t i>
- static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompIdx)
- { return IndexHelper::couplingCompIdx(id, coupledCompIdx); }
-
- /*!
- * \brief Returns the intrinsic permeability of the coupled Darcy element.
- */
- template<class Context>
- static auto darcyPermeability(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
- const Context& context)
- { return context.volVars.permeability(); }
-
- /*!
- * \brief Returns the momentum flux across the coupling boundary.
- *
- * For the normal momentum coupling, the porous medium side of the coupling condition
- * is evaluated, i.e. -[p n]^pm.
- *
- */
- template<class Context>
- static NumEqVector<freeFlowMomentumIndex> momentumCouplingCondition(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
- const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars,
- const Context& context)
- {
- static constexpr auto numPhasesDarcy = GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::ModelTraits>::numFluidPhases();
- NumEqVector<freeFlowMomentumIndex> momentumFlux(0.0);
-
- if (!scvf.isFrontal())
- return momentumFlux;
-
- const auto pmPhaseIdx = couplingPhaseIdx(porousMediumIndex);
-
- if(numPhasesDarcy > 1)
- momentumFlux[scvf.normalAxis()] = context.volVars.pressure(pmPhaseIdx);
- else // use pressure reconstruction for single phase models
- momentumFlux[scvf.normalAxis()] = pressureAtInterface_(fvGeometry, scvf, elemVolVars, context);
-
- // TODO: generalize for permeability tensors
-
- // normalize pressure
- momentumFlux[scvf.normalAxis()] -= elemVolVars.gridVolVars().problem().referencePressure(fvGeometry.element(), fvGeometry, scvf);
-
- momentumFlux[scvf.normalAxis()] *= scvf.directionSign();
-
- return momentumFlux;
- }
-
- /*!
- * \brief Evaluate an advective flux across the interface and consider upwinding.
- */
- static Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream)
- {
- const Scalar upwindWeight = 1.0; //TODO use Flux.UpwindWeight or something like Coupling.UpwindWeight?
-
- if(insideIsUpstream)
- return (upwindWeight * insideQuantity + (1.0 - upwindWeight) * outsideQuantity) * volumeFlow;
- else
- return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow;
- }
+struct FreeFlowPorousMediumCouplingConditionsSelector<MDTraits, CouplingManager, DiscretizationMethods::FCStaggered, DiscretizationMethods::CCTpfa, DiscretizationMethods::CCTpfa>
+{ using type = FFPMCouplingConditionsStaggeredCCTpfa<MDTraits, CouplingManager>; };
-protected:
+} // end namespace FreeFlowPorousMediumDetail
+#endif // DOXYGEN
- /*!
- * \brief Returns the transmissibility used for either molecular diffusion or thermal conductivity.
- */
- template<std::size_t i, std::size_t j>
- Scalar transmissibility_(Dune::index_constant<i> domainI,
- Dune::index_constant<j> domainJ,
- const Scalar insideDistance,
- const Scalar outsideDistance,
- const Scalar avgQuantityI,
- const Scalar avgQuantityJ,
- const DiffusionCoefficientAveragingType diffCoeffAvgType) const
- {
- const Scalar totalDistance = insideDistance + outsideDistance;
- if(diffCoeffAvgType == DiffusionCoefficientAveragingType::harmonic)
- {
- return harmonicMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance)
- / totalDistance;
- }
- else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::arithmetic)
- {
- return arithmeticMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance)
- / totalDistance;
- }
- else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::ffOnly)
- return domainI == freeFlowMassIndex
- ? avgQuantityI / totalDistance
- : avgQuantityJ / totalDistance;
-
- else // diffCoeffAvgType == DiffusionCoefficientAveragingType::pmOnly)
- return domainI == porousMediumIndex
- ? avgQuantityI / totalDistance
- : avgQuantityJ / totalDistance;
- }
-
- /*!
- * \brief Returns the distance between an scvf and the corresponding scv center.
- */
- template<class Scv, class Scvf>
- Scalar getDistance_(const Scv& scv, const Scvf& scvf) const
- {
- return (scv.dofPosition() - scvf.ipGlobal()).two_norm();
- }
-
- /*!
- * \brief Returns the conductive energy flux across the interface.
- */
- template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
- Scalar conductiveEnergyFlux_(Dune::index_constant<i> domainI,
- Dune::index_constant<j> domainJ,
- const FVElementGeometry<i>& fvGeometryI,
- const FVElementGeometry<j>& fvGeometryJ,
- const SubControlVolumeFace<i>& scvfI,
- const SubControlVolume<i>& scvI,
- const SubControlVolume<j>& scvJ,
- const VolumeVariables<i>& volVarsI,
- const VolumeVariables<j>& volVarsJ,
- const DiffusionCoefficientAveragingType diffCoeffAvgType) const
- {
- const Scalar insideDistance = getDistance_(scvI, scvfI);
- const Scalar outsideDistance = getDistance_(scvJ, scvfI);
-
- const Scalar deltaT = volVarsJ.temperature() - volVarsI.temperature();
- const Scalar tij = transmissibility_(
- domainI, domainJ,
- insideDistance, outsideDistance,
- volVarsI.effectiveThermalConductivity(), volVarsJ.effectiveThermalConductivity(),
- diffCoeffAvgType
- );
-
- return -tij * deltaT;
- }
-
- /*!
- * \brief Returns the pressure at the interface
- */
- template<class CouplingContext>
- static Scalar pressureAtInterface_(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
- const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars,
- const CouplingContext& context)
- {
- GlobalPosition<freeFlowMomentumIndex> velocity(0.0);
- const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
- velocity[scv.dofAxis()] = elemVolVars[scv].velocity();
- const auto& pmScvf = context.fvGeometry.scvf(context.porousMediumScvfIdx);
- return computeCouplingPhasePressureAtInterface_(context.fvGeometry, pmScvf, context.volVars, context.gravity, velocity, AdvectionType());
- }
-
- /*!
- * \brief Returns the pressure at the interface using Forchheimers's law for reconstruction
- */
- static Scalar computeCouplingPhasePressureAtInterface_(const FVElementGeometry<porousMediumIndex>& fvGeometry,
- const SubControlVolumeFace<porousMediumIndex>& scvf,
- const VolumeVariables<porousMediumIndex>& volVars,
- const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& couplingPhaseVelocity,
- ForchheimersLaw)
- {
- DUNE_THROW(Dune::NotImplemented, "Forchheimer's law pressure reconstruction");
- }
-
- /*!
- * \brief Returns the pressure at the interface using Darcy's law for reconstruction
- */
- static Scalar computeCouplingPhasePressureAtInterface_(const FVElementGeometry<porousMediumIndex>& fvGeometry,
- const SubControlVolumeFace<porousMediumIndex>& scvf,
- const VolumeVariables<porousMediumIndex>& volVars,
- const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& gravity,
- const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& couplingPhaseVelocity,
- DarcysLaw)
- {
- const auto pmPhaseIdx = couplingPhaseIdx(porousMediumIndex);
- const Scalar couplingPhaseCellCenterPressure = volVars.pressure(pmPhaseIdx);
- const Scalar couplingPhaseMobility = volVars.mobility(pmPhaseIdx);
- const Scalar couplingPhaseDensity = volVars.density(pmPhaseIdx);
- const auto K = volVars.permeability();
-
- // A tpfa approximation yields (works if mobility != 0)
- // v*n = -kr/mu*K * (gradP - rho*g)*n = mobility*(ti*(p_center - p_interface) + rho*n^TKg)
- // -> p_interface = (1/mobility * (-v*n) + rho*n^TKg)/ti + p_center
- // where v is the free-flow velocity (couplingPhaseVelocity)
- const auto alpha = vtmv(scvf.unitOuterNormal(), K, gravity);
-
- const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
- const auto ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, K, 1.0);
-
- return (-1/couplingPhaseMobility * (scvf.unitOuterNormal() * couplingPhaseVelocity) + couplingPhaseDensity * alpha)/ti
- + couplingPhaseCellCenterPressure;
- }
-};
-
-/*!
- * \ingroup FreeFlowPorousMediumCoupling
- * \brief Coupling data specialization for non-compositional models.
- */
-template<class MDTraits, class CouplingManager, bool enableEnergyBalance>
-class FreeFlowPorousMediumCouplingConditionsImplementation<MDTraits, CouplingManager, enableEnergyBalance, false>
-: public FreeFlowPorousMediumCouplingConditionsImplementationBase<MDTraits, CouplingManager>
-{
- using ParentType = FreeFlowPorousMediumCouplingConditionsImplementationBase<MDTraits, CouplingManager>;
- using Scalar = typename MDTraits::Scalar;
-
- // the sub domain type tags
- template<std::size_t id>
- using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
-
- template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
- template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
- template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
- template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
- template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
- template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
- template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
- template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
-
- static_assert(GetPropType<SubDomainTypeTag<ParentType::porousMediumIndex>, Properties::ModelTraits>::numFluidComponents() == GetPropType<SubDomainTypeTag<ParentType::porousMediumIndex>, Properties::ModelTraits>::numFluidPhases(),
- "Darcy Model must not be compositional");
-
- using DiffusionCoefficientAveragingType = typename FreeFlowPorousMediumCouplingOptions::DiffusionCoefficientAveragingType;
-
-public:
- using ParentType::ParentType;
- using ParentType::couplingPhaseIdx;
-
- /*!
- * \brief Returns the mass flux across the coupling boundary as seen from the Darcy domain.
- */
- template<std::size_t i, std::size_t j, class CouplingContext>
- static Scalar massCouplingCondition(Dune::index_constant<i> domainI,
- Dune::index_constant<j> domainJ,
- const FVElementGeometry<i>& fvGeometry,
- const SubControlVolumeFace<i>& scvf,
- const ElementVolumeVariables<i>& insideVolVars,
- const CouplingContext& context)
- {
- const Scalar normalVelocity = context.velocity * scvf.unitOuterNormal();
- const Scalar darcyDensity = insideVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(ParentType::porousMediumIndex));
- const Scalar stokesDensity = context.volVars.density();
- const bool insideIsUpstream = normalVelocity > 0.0;
- return ParentType::advectiveFlux(darcyDensity, stokesDensity, normalVelocity, insideIsUpstream);
- }
-
- /*!
- * \brief Returns the energy flux across the coupling boundary as seen from the Darcy domain.
- */
- template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
- Scalar energyCouplingCondition(const Element<ParentType::porousMediumIndex>& element,
- const FVElementGeometry<ParentType::porousMediumIndex>& fvGeometry,
- const ElementVolumeVariables<ParentType::porousMediumIndex>& darcyElemVolVars,
- const SubControlVolumeFace<ParentType::porousMediumIndex>& scvf,
- const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
- {
- DUNE_THROW(Dune::NotImplemented, "Energy coupling condition");
- }
-
- /*!
- * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain.
- */
- template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
- Scalar energyCouplingCondition(const Element<ParentType::freeFlowMassIndex>& element,
- const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry,
- const ElementVolumeVariables<ParentType::freeFlowMassIndex>& stokesElemVolVars,
- const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
- const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
- {
- DUNE_THROW(Dune::NotImplemented, "Energy coupling condition");
- }
-
-private:
-
-
-
- /*!
- * \brief Evaluate the energy flux across the interface.
- */
- template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
- Scalar energyFlux_(Dune::index_constant<i> domainI,
- Dune::index_constant<j> domainJ,
- const FVElementGeometry<i>& insideFvGeometry,
- const FVElementGeometry<j>& outsideFvGeometry,
- const SubControlVolumeFace<i>& scvf,
- const VolumeVariables<i>& insideVolVars,
- const VolumeVariables<j>& outsideVolVars,
- const Scalar velocity,
- const bool insideIsUpstream,
- const DiffusionCoefficientAveragingType diffCoeffAvgType) const
- {
- Scalar flux(0.0);
-
- const auto& insideScv = (*scvs(insideFvGeometry).begin());
- const auto& outsideScv = (*scvs(outsideFvGeometry).begin());
-
- // convective fluxes
- const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI));
- const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ));
-
- flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream);
-
- flux += this->conductiveEnergyFlux_(domainI, domainJ, insideFvGeometry, outsideFvGeometry, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType);
-
- return flux;
- }
-
-};
+template<class MDTraits, class CouplingManager>
+using FreeFlowPorousMediumCouplingConditions = typename FreeFlowPorousMediumDetail::FreeFlowPorousMediumCouplingConditionsSelector<MDTraits,CouplingManager>::type;
} // end namespace Dumux
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions_staggered_cctpfa.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions_staggered_cctpfa.hh
new file mode 100644
index 0000000000000000000000000000000000000000..91627504b6cbe3ab4489269ae3ea8c8c3af539fa
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions_staggered_cctpfa.hh
@@ -0,0 +1,455 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \copydoc Dumux::FFPMCouplingConditionsStaggeredCCTpfa
+ */
+
+#ifndef DUMUX_MD_FREEFLOW_POROUSMEDIUM_COUPLINGCONDITIONS_STAGGERED_TPFA_HH
+#define DUMUX_MD_FREEFLOW_POROUSMEDIUM_COUPLINGCONDITIONS_STAGGERED_TPFA_HH
+
+#include <numeric>
+
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/math.hh>
+
+#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
+
+#include <dumux/flux/darcyslaw_fwd.hh>
+#include <dumux/flux/fickslaw_fwd.hh>
+#include <dumux/flux/forchheimerslaw_fwd.hh>
+
+#include <dumux/multidomain/boundary/freeflowporousmedium/traits.hh>
+#include <dumux/multidomain/boundary/freeflowporousmedium/indexhelper.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief This structs holds a set of options which allow to modify the Stokes-Darcy
+ * coupling mechanism during runtime.
+ */
+struct FreeFlowPorousMediumCouplingOptions
+{
+ /*!
+ * \brief Defines which kind of averanging of diffusion coefficients
+ * (moleculat diffusion or thermal conductance) at the interface
+ * between free flow and porous medium shall be used.
+ */
+ enum class DiffusionCoefficientAveragingType
+ {
+ harmonic, arithmetic, ffOnly, pmOnly
+ };
+
+ /*!
+ * \brief Convenience function to convert user input given as std::string to the corresponding enum class used for choosing the type
+ * of averaging of the diffusion/conduction parameter at the interface between the two domains.
+ */
+ static DiffusionCoefficientAveragingType stringToEnum(DiffusionCoefficientAveragingType, const std::string& diffusionCoefficientAveragingType)
+ {
+ if (diffusionCoefficientAveragingType == "Harmonic")
+ return DiffusionCoefficientAveragingType::harmonic;
+ else if (diffusionCoefficientAveragingType == "Arithmetic")
+ return DiffusionCoefficientAveragingType::arithmetic;
+ else if (diffusionCoefficientAveragingType == "FreeFlowOnly")
+ return DiffusionCoefficientAveragingType::ffOnly;
+ else if (diffusionCoefficientAveragingType == "PorousMediumOnly")
+ return DiffusionCoefficientAveragingType::pmOnly;
+ else
+ DUNE_THROW(Dune::IOError, "Unknown DiffusionCoefficientAveragingType");
+ }
+
+};
+
+template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional>
+class FFPMCouplingConditionsStaggeredCCTpfaImpl;
+
+/*!
+* \ingroup FreeFlowPorousMediumCoupling
+* \brief Data for the coupling of a Darcy model (cell-centered finite volume)
+* with a (Navier-)Stokes model (staggerd grid).
+*/
+template<class MDTraits, class CouplingManager>
+using FFPMCouplingConditionsStaggeredCCTpfa
+ = FFPMCouplingConditionsStaggeredCCTpfaImpl<
+ MDTraits, CouplingManager,
+ GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(),
+ (GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::numFluidComponents() > 1)
+ >;
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief A base class which provides some common methods used for Stokes-Darcy coupling.
+ */
+template<class MDTraits, class CouplingManager>
+class FFPMCouplingConditionsStaggeredCCTpfaImplBase
+{
+ using Scalar = typename MDTraits::Scalar;
+
+ template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
+ template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
+ template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
+ template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
+ template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
+ template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
+ template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
+ template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+ template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>;
+ template<std::size_t id> using ModelTraits = GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>;
+ template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate;
+ template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector;
+
+public:
+ static constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex;
+ static constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex;
+ static constexpr auto porousMediumIndex = CouplingManager::porousMediumIndex;
+private:
+
+ using AdvectionType = GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::AdvectionType>;
+ using DarcysLaw = Dumux::DarcysLaw<SubDomainTypeTag<porousMediumIndex>>;
+ using ForchheimersLaw = Dumux::ForchheimersLaw<SubDomainTypeTag<porousMediumIndex>>;
+
+ static constexpr bool adapterUsed = ModelTraits<porousMediumIndex>::numFluidPhases() > 1;
+ using IndexHelper = FreeFlowPorousMediumCoupling::IndexHelper<freeFlowMassIndex, porousMediumIndex, FluidSystem<freeFlowMassIndex>, adapterUsed>;
+
+ static constexpr int enableEnergyBalance = GetPropType<SubDomainTypeTag<freeFlowMassIndex>, Properties::ModelTraits>::enableEnergyBalance();
+ static_assert(GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::ModelTraits>::enableEnergyBalance() == enableEnergyBalance,
+ "All submodels must both be either isothermal or non-isothermal");
+
+ static_assert(FreeFlowPorousMediumCoupling::IsSameFluidSystem<FluidSystem<freeFlowMassIndex>,
+ FluidSystem<porousMediumIndex>>::value,
+ "All submodels must use the same fluid system");
+
+ using DiffusionCoefficientAveragingType = typename FreeFlowPorousMediumCouplingOptions::DiffusionCoefficientAveragingType;
+
+public:
+ /*!
+ * \brief Returns the corresponding phase index needed for coupling.
+ */
+ template<std::size_t i>
+ static constexpr auto couplingPhaseIdx(Dune::index_constant<i> id, int coupledPhaseIdx = 0)
+ { return IndexHelper::couplingPhaseIdx(id, coupledPhaseIdx); }
+
+ /*!
+ * \brief Returns the corresponding component index needed for coupling.
+ */
+ template<std::size_t i>
+ static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompIdx)
+ { return IndexHelper::couplingCompIdx(id, coupledCompIdx); }
+
+ /*!
+ * \brief Returns the intrinsic permeability of the coupled Darcy element.
+ */
+ template<class Context>
+ static auto darcyPermeability(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
+ const Context& context)
+ { return context.volVars.permeability(); }
+
+ /*!
+ * \brief Returns the momentum flux across the coupling boundary.
+ *
+ * For the normal momentum coupling, the porous medium side of the coupling condition
+ * is evaluated, i.e. -[p n]^pm.
+ *
+ */
+ template<class Context>
+ static NumEqVector<freeFlowMomentumIndex> momentumCouplingCondition(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
+ const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars,
+ const Context& context)
+ {
+ static constexpr auto numPhasesDarcy = GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::ModelTraits>::numFluidPhases();
+ NumEqVector<freeFlowMomentumIndex> momentumFlux(0.0);
+
+ if (!scvf.isFrontal())
+ return momentumFlux;
+
+ const auto pmPhaseIdx = couplingPhaseIdx(porousMediumIndex);
+
+ if(numPhasesDarcy > 1)
+ momentumFlux[scvf.normalAxis()] = context.volVars.pressure(pmPhaseIdx);
+ else // use pressure reconstruction for single phase models
+ momentumFlux[scvf.normalAxis()] = pressureAtInterface_(fvGeometry, scvf, elemVolVars, context);
+
+ // TODO: generalize for permeability tensors
+
+ // normalize pressure
+ momentumFlux[scvf.normalAxis()] -= elemVolVars.gridVolVars().problem().referencePressure(fvGeometry.element(), fvGeometry, scvf);
+
+ momentumFlux[scvf.normalAxis()] *= scvf.directionSign();
+
+ return momentumFlux;
+ }
+
+ /*!
+ * \brief Evaluate an advective flux across the interface and consider upwinding.
+ */
+ static Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream)
+ {
+ const Scalar upwindWeight = 1.0; //TODO use Flux.UpwindWeight or something like Coupling.UpwindWeight?
+
+ if(insideIsUpstream)
+ return (upwindWeight * insideQuantity + (1.0 - upwindWeight) * outsideQuantity) * volumeFlow;
+ else
+ return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow;
+ }
+
+protected:
+
+ /*!
+ * \brief Returns the transmissibility used for either molecular diffusion or thermal conductivity.
+ */
+ template<std::size_t i, std::size_t j>
+ Scalar transmissibility_(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const Scalar insideDistance,
+ const Scalar outsideDistance,
+ const Scalar avgQuantityI,
+ const Scalar avgQuantityJ,
+ const DiffusionCoefficientAveragingType diffCoeffAvgType) const
+ {
+ const Scalar totalDistance = insideDistance + outsideDistance;
+ if(diffCoeffAvgType == DiffusionCoefficientAveragingType::harmonic)
+ {
+ return harmonicMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance)
+ / totalDistance;
+ }
+ else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::arithmetic)
+ {
+ return arithmeticMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance)
+ / totalDistance;
+ }
+ else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::ffOnly)
+ return domainI == freeFlowMassIndex
+ ? avgQuantityI / totalDistance
+ : avgQuantityJ / totalDistance;
+
+ else // diffCoeffAvgType == DiffusionCoefficientAveragingType::pmOnly)
+ return domainI == porousMediumIndex
+ ? avgQuantityI / totalDistance
+ : avgQuantityJ / totalDistance;
+ }
+
+ /*!
+ * \brief Returns the distance between an scvf and the corresponding scv center.
+ */
+ template<class Scv, class Scvf>
+ Scalar getDistance_(const Scv& scv, const Scvf& scvf) const
+ {
+ return (scv.dofPosition() - scvf.ipGlobal()).two_norm();
+ }
+
+ /*!
+ * \brief Returns the conductive energy flux across the interface.
+ */
+ template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ Scalar conductiveEnergyFlux_(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const FVElementGeometry<i>& fvGeometryI,
+ const FVElementGeometry<j>& fvGeometryJ,
+ const SubControlVolumeFace<i>& scvfI,
+ const SubControlVolume<i>& scvI,
+ const SubControlVolume<j>& scvJ,
+ const VolumeVariables<i>& volVarsI,
+ const VolumeVariables<j>& volVarsJ,
+ const DiffusionCoefficientAveragingType diffCoeffAvgType) const
+ {
+ const Scalar insideDistance = getDistance_(scvI, scvfI);
+ const Scalar outsideDistance = getDistance_(scvJ, scvfI);
+
+ const Scalar deltaT = volVarsJ.temperature() - volVarsI.temperature();
+ const Scalar tij = transmissibility_(
+ domainI, domainJ,
+ insideDistance, outsideDistance,
+ volVarsI.effectiveThermalConductivity(), volVarsJ.effectiveThermalConductivity(),
+ diffCoeffAvgType
+ );
+
+ return -tij * deltaT;
+ }
+
+ /*!
+ * \brief Returns the pressure at the interface
+ */
+ template<class CouplingContext>
+ static Scalar pressureAtInterface_(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
+ const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars,
+ const CouplingContext& context)
+ {
+ GlobalPosition<freeFlowMomentumIndex> velocity(0.0);
+ const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
+ velocity[scv.dofAxis()] = elemVolVars[scv].velocity();
+ const auto& pmScvf = context.fvGeometry.scvf(context.porousMediumScvfIdx);
+ return computeCouplingPhasePressureAtInterface_(context.fvGeometry, pmScvf, context.volVars, context.gravity, velocity, AdvectionType());
+ }
+
+ /*!
+ * \brief Returns the pressure at the interface using Forchheimers's law for reconstruction
+ */
+ static Scalar computeCouplingPhasePressureAtInterface_(const FVElementGeometry<porousMediumIndex>& fvGeometry,
+ const SubControlVolumeFace<porousMediumIndex>& scvf,
+ const VolumeVariables<porousMediumIndex>& volVars,
+ const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& couplingPhaseVelocity,
+ ForchheimersLaw)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Forchheimer's law pressure reconstruction");
+ }
+
+ /*!
+ * \brief Returns the pressure at the interface using Darcy's law for reconstruction
+ */
+ static Scalar computeCouplingPhasePressureAtInterface_(const FVElementGeometry<porousMediumIndex>& fvGeometry,
+ const SubControlVolumeFace<porousMediumIndex>& scvf,
+ const VolumeVariables<porousMediumIndex>& volVars,
+ const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& gravity,
+ const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& couplingPhaseVelocity,
+ DarcysLaw)
+ {
+ const auto pmPhaseIdx = couplingPhaseIdx(porousMediumIndex);
+ const Scalar couplingPhaseCellCenterPressure = volVars.pressure(pmPhaseIdx);
+ const Scalar couplingPhaseMobility = volVars.mobility(pmPhaseIdx);
+ const Scalar couplingPhaseDensity = volVars.density(pmPhaseIdx);
+ const auto K = volVars.permeability();
+
+ // A tpfa approximation yields (works if mobility != 0)
+ // v*n = -kr/mu*K * (gradP - rho*g)*n = mobility*(ti*(p_center - p_interface) + rho*n^TKg)
+ // -> p_interface = (1/mobility * (-v*n) + rho*n^TKg)/ti + p_center
+ // where v is the free-flow velocity (couplingPhaseVelocity)
+ const auto alpha = vtmv(scvf.unitOuterNormal(), K, gravity);
+
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, K, 1.0);
+
+ return (-1/couplingPhaseMobility * (scvf.unitOuterNormal() * couplingPhaseVelocity) + couplingPhaseDensity * alpha)/ti
+ + couplingPhaseCellCenterPressure;
+ }
+};
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Coupling data specialization for non-compositional models.
+ */
+template<class MDTraits, class CouplingManager, bool enableEnergyBalance>
+class FFPMCouplingConditionsStaggeredCCTpfaImpl<MDTraits, CouplingManager, enableEnergyBalance, false>
+: public FFPMCouplingConditionsStaggeredCCTpfaImplBase<MDTraits, CouplingManager>
+{
+ using ParentType = FFPMCouplingConditionsStaggeredCCTpfaImplBase<MDTraits, CouplingManager>;
+ using Scalar = typename MDTraits::Scalar;
+
+ // the sub domain type tags
+ template<std::size_t id>
+ using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
+ template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
+ template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
+ template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
+ template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
+ template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
+ template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
+
+ static_assert(GetPropType<SubDomainTypeTag<ParentType::porousMediumIndex>, Properties::ModelTraits>::numFluidComponents() == GetPropType<SubDomainTypeTag<ParentType::porousMediumIndex>, Properties::ModelTraits>::numFluidPhases(),
+ "Darcy Model must not be compositional");
+
+ using DiffusionCoefficientAveragingType = typename FreeFlowPorousMediumCouplingOptions::DiffusionCoefficientAveragingType;
+
+public:
+ using ParentType::ParentType;
+ using ParentType::couplingPhaseIdx;
+
+ /*!
+ * \brief Returns the mass flux across the coupling boundary as seen from the Darcy domain.
+ */
+ template<std::size_t i, std::size_t j, class CouplingContext>
+ static Scalar massCouplingCondition(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const FVElementGeometry<i>& fvGeometry,
+ const SubControlVolumeFace<i>& scvf,
+ const ElementVolumeVariables<i>& insideVolVars,
+ const CouplingContext& context)
+ {
+ const Scalar normalVelocity = context.velocity * scvf.unitOuterNormal();
+ const Scalar darcyDensity = insideVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(ParentType::porousMediumIndex));
+ const Scalar stokesDensity = context.volVars.density();
+ const bool insideIsUpstream = normalVelocity > 0.0;
+ return ParentType::advectiveFlux(darcyDensity, stokesDensity, normalVelocity, insideIsUpstream);
+ }
+
+ /*!
+ * \brief Returns the energy flux across the coupling boundary as seen from the Darcy domain.
+ */
+ template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ Scalar energyCouplingCondition(const Element<ParentType::porousMediumIndex>& element,
+ const FVElementGeometry<ParentType::porousMediumIndex>& fvGeometry,
+ const ElementVolumeVariables<ParentType::porousMediumIndex>& darcyElemVolVars,
+ const SubControlVolumeFace<ParentType::porousMediumIndex>& scvf,
+ const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "Energy coupling condition");
+ }
+
+ /*!
+ * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain.
+ */
+ template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ Scalar energyCouplingCondition(const Element<ParentType::freeFlowMassIndex>& element,
+ const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry,
+ const ElementVolumeVariables<ParentType::freeFlowMassIndex>& stokesElemVolVars,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "Energy coupling condition");
+ }
+
+private:
+
+
+
+ /*!
+ * \brief Evaluate the energy flux across the interface.
+ */
+ template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ Scalar energyFlux_(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const FVElementGeometry<i>& insideFvGeometry,
+ const FVElementGeometry<j>& outsideFvGeometry,
+ const SubControlVolumeFace<i>& scvf,
+ const VolumeVariables<i>& insideVolVars,
+ const VolumeVariables<j>& outsideVolVars,
+ const Scalar velocity,
+ const bool insideIsUpstream,
+ const DiffusionCoefficientAveragingType diffCoeffAvgType) const
+ {
+ Scalar flux(0.0);
+
+ const auto& insideScv = (*scvs(insideFvGeometry).begin());
+ const auto& outsideScv = (*scvs(outsideFvGeometry).begin());
+
+ // convective fluxes
+ const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI));
+ const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ));
+
+ flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream);
+
+ flux += this->conductiveEnergyFlux_(domainI, domainJ, insideFvGeometry, outsideFvGeometry, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType);
+
+ return flux;
+ }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh
index 0ad552da909e7d4d3722e3f90e7899bee9625915..b4d6b62d68bd1dff2c41cd95b80329c54ed28658 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh
@@ -7,459 +7,39 @@
/*!
* \file
* \ingroup FreeFlowPorousMediumCoupling
- * \copydoc Dumux::FreeFlowPorousMediumCouplingManager
+ * \brief Coupling managers for free and porous medium flow, specialized for different discretization schemes.
*/
#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_POROUSMEDIUM_COUPLINGMANAGER_HH
#define DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_POROUSMEDIUM_COUPLINGMANAGER_HH
-#include <utility>
-#include <memory>
+#include <dumux/discretization/method.hh>
-#include <dune/common/indices.hh>
-#include <dumux/common/properties.hh>
-#include <dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh>
-#include <dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh>
-#include <dumux/multidomain/freeflow/couplingmanager.hh>
-#include <dumux/multidomain/multibinarycouplingmanager.hh>
-
-#include "couplingconditions.hh"
+#include "couplingmanager_base.hh"
+#include "couplingmanager_staggered_cctpfa.hh"
namespace Dumux {
+#ifndef DOXYGEN
namespace FreeFlowPorousMediumDetail {
-// global subdomain indices
-static constexpr auto freeFlowMomentumIndex = Dune::index_constant<0>();
-static constexpr auto freeFlowMassIndex = Dune::index_constant<1>();
-static constexpr auto porousMediumIndex = Dune::index_constant<2>();
-
-// coupling indices
-static constexpr auto freeFlowMassToFreeFlowMomentumIndex = Dune::index_constant<0>();
-static constexpr auto freeFlowMomentumToPorousMediumIndex = Dune::index_constant<1>();
-static constexpr auto freeFlowMassToPorousMediumIndex = Dune::index_constant<2>();
-static constexpr auto noCouplingIdx = Dune::index_constant<99>();
-
-constexpr auto makeCouplingManagerMap()
-{
- auto map = std::array<std::array<std::size_t, 3>, 3>{};
-
- // free flow (momentum-mass)
- map[freeFlowMomentumIndex][freeFlowMassIndex] = freeFlowMassToFreeFlowMomentumIndex;
- map[freeFlowMassIndex][freeFlowMomentumIndex] = freeFlowMassToFreeFlowMomentumIndex;
-
- // free flow momentum - porous medium
- map[freeFlowMomentumIndex][porousMediumIndex] = freeFlowMomentumToPorousMediumIndex;
- map[porousMediumIndex][freeFlowMomentumIndex] = freeFlowMomentumToPorousMediumIndex;
-
- // free flow mass - porous medium
- map[freeFlowMassIndex][porousMediumIndex] = freeFlowMassToPorousMediumIndex;
- map[porousMediumIndex][freeFlowMassIndex] = freeFlowMassToPorousMediumIndex;
-
- return map;
-}
-
-template<std::size_t i>
-constexpr auto coupledDomains(Dune::index_constant<i> domainI)
-{
- if constexpr (i == freeFlowMomentumIndex)
- return std::make_tuple(freeFlowMassIndex, porousMediumIndex);
- else if constexpr (i == freeFlowMassIndex)
- return std::make_tuple(freeFlowMomentumIndex, porousMediumIndex);
- else // i == porousMediumIndex
- return std::make_tuple(freeFlowMomentumIndex, freeFlowMassIndex);
-}
-
-template<std::size_t i, std::size_t j>
-constexpr auto globalToLocalDomainIndices(Dune::index_constant<i>, Dune::index_constant<j>)
-{
- static_assert(i <= 2 && j <= 2);
- static_assert(i != j);
-
- if constexpr (i < j)
- return std::pair<Dune::index_constant<0>, Dune::index_constant<1>>{};
- else
- return std::pair<Dune::index_constant<1>, Dune::index_constant<0>>{};
-}
-
-struct CouplingMaps
-{
- static constexpr auto managerMap()
- {
- return FreeFlowPorousMediumDetail::makeCouplingManagerMap();
- }
-
- template<std::size_t i, std::size_t j>
- static constexpr auto globalToLocal(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
- {
- return FreeFlowPorousMediumDetail::globalToLocalDomainIndices(domainI, domainJ);
- }
-
- template<std::size_t i>
- static constexpr auto coupledDomains(Dune::index_constant<i> domainI)
- {
- return FreeFlowPorousMediumDetail::coupledDomains(domainI);
- }
-};
+// declaration (specialize for different discretization types)
+template<class MDTraits,
+ class DiscFFMomentum = typename MDTraits::template SubDomain<FreeFlowPorousMediumDetail::freeFlowMomentumIndex>::GridGeometry::DiscretizationMethod,
+ class DiscFFMass = typename MDTraits::template SubDomain<FreeFlowPorousMediumDetail::freeFlowMassIndex>::GridGeometry::DiscretizationMethod,
+ class DiscPM = typename MDTraits::template SubDomain<FreeFlowPorousMediumDetail::porousMediumIndex>::GridGeometry::DiscretizationMethod
+ >
+struct FreeFlowPorousMediumCouplingManagerSelector;
template<class MDTraits>
-struct CouplingManagers
-{
- template<std::size_t id>
- using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
-
- using FreeFlowTraits = MultiDomainTraits<
- SubDomainTypeTag<freeFlowMomentumIndex>, SubDomainTypeTag<freeFlowMassIndex>
- >;
-
- using FreeFlowMomentumPorousMediumTraits = MultiDomainTraits<
- SubDomainTypeTag<freeFlowMomentumIndex>, SubDomainTypeTag<porousMediumIndex>
- >;
-
- using FreeFlowMassPorousMediumTraits = MultiDomainTraits<
- SubDomainTypeTag<freeFlowMassIndex>, SubDomainTypeTag<porousMediumIndex>
- >;
+struct FreeFlowPorousMediumCouplingManagerSelector<MDTraits, DiscretizationMethods::FCStaggered, DiscretizationMethods::CCTpfa, DiscretizationMethods::CCTpfa>
+{ using type = FreeFlowPorousMediumCouplingManagerStaggeredCCTpfa<MDTraits>; };
- using FreeFlowCouplingManager
- = Dumux::FreeFlowCouplingManager<FreeFlowTraits>;
- using FreeFlowMomentumPorousMediumCouplingManager
- = Dumux::FreeFlowMomentumPorousMediumCouplingManager<FreeFlowMomentumPorousMediumTraits>;
- using FreeFlowMassPorousMediumCouplingManager
- = Dumux::FreeFlowMassPorousMediumCouplingManager<FreeFlowMassPorousMediumTraits>;
-};
+} // end namespace FreeFlowPorousMediumDetail
+#endif // DOXYGEN
-} // end namespace Detail
-
-
-/*!
- * \ingroup FreeFlowPorousMediumCoupling
- * \brief Coupling manager for coupling freeflow and porous medium flow models
- */
template<class MDTraits>
-class FreeFlowPorousMediumCouplingManager
-: public MultiBinaryCouplingManager<
- MDTraits,
- FreeFlowPorousMediumDetail::CouplingMaps,
- typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowCouplingManager,
- typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMomentumPorousMediumCouplingManager,
- typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMassPorousMediumCouplingManager
->
-{
- using ParentType = MultiBinaryCouplingManager<
- MDTraits,
- FreeFlowPorousMediumDetail::CouplingMaps,
- typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowCouplingManager,
- typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMomentumPorousMediumCouplingManager,
- typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMassPorousMediumCouplingManager
- >;
-
- using Scalar = typename MDTraits::Scalar;
-
- // the sub domain type tags
- template<std::size_t id>
- using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
-
- template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
- template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
- template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
- template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
- template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
- template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
- template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector;
-
- template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
- template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
- using SolutionVector = typename MDTraits::SolutionVector;
-
- using CouplingConditions = FreeFlowPorousMediumCouplingConditions<MDTraits, FreeFlowPorousMediumCouplingManager<MDTraits>>;
-
-public:
-
- template<std::size_t i, std::size_t j>
- using SubCouplingManager = typename ParentType::template SubCouplingManager<i, j>;
-
- static constexpr auto freeFlowMomentumIndex = FreeFlowPorousMediumDetail::freeFlowMomentumIndex;
- static constexpr auto freeFlowMassIndex = FreeFlowPorousMediumDetail::freeFlowMassIndex;
- static constexpr auto porousMediumIndex = FreeFlowPorousMediumDetail::porousMediumIndex;
-
-public:
- using ParentType::ParentType;
-
- template<class GridVarsTuple>
- void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
- std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
- std::shared_ptr<Problem<porousMediumIndex>> porousMediumProblem,
- GridVarsTuple&& gridVarsTuple,
- const SolutionVector& curSol)
- {
- this->updateSolution(curSol); // generic coupling manager stores tuple of shared_ptr
-
- // initialize the binary sub coupling managers
- typename SubCouplingManager<freeFlowMomentumIndex, freeFlowMassIndex>::SolutionVectorStorage ffSolVecTuple;
- std::get<0>(ffSolVecTuple) = std::get<freeFlowMomentumIndex>(this->curSol());
- std::get<1>(ffSolVecTuple) = std::get<freeFlowMassIndex>(this->curSol());
- this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).init(
- freeFlowMomentumProblem, freeFlowMassProblem,
- std::make_tuple(std::get<freeFlowMomentumIndex>(gridVarsTuple), std::get<freeFlowMassIndex>(gridVarsTuple)),
- ffSolVecTuple
- );
-
- typename SubCouplingManager<freeFlowMassIndex, porousMediumIndex>::SolutionVectorStorage ffMassPmSolVecTuple;
- std::get<0>(ffMassPmSolVecTuple) = std::get<freeFlowMassIndex>(this->curSol());
- std::get<1>(ffMassPmSolVecTuple) = std::get<porousMediumIndex>(this->curSol());
- this->subCouplingManager(freeFlowMassIndex, porousMediumIndex).init(
- freeFlowMassProblem, porousMediumProblem, ffMassPmSolVecTuple
- );
-
- typename SubCouplingManager<freeFlowMomentumIndex, porousMediumIndex>::SolutionVectorStorage ffMomentumPmSolVecTuple;
- std::get<0>(ffMomentumPmSolVecTuple) = std::get<freeFlowMomentumIndex>(this->curSol());
- std::get<1>(ffMomentumPmSolVecTuple) = std::get<porousMediumIndex>(this->curSol());
- this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).init(
- freeFlowMomentumProblem, porousMediumProblem, ffMomentumPmSolVecTuple
- );
- }
-
- /*!
- * \brief Returns the mass flux across the coupling boundary.
- */
- template<std::size_t i, std::size_t j>
- auto massCouplingCondition(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ,
- const FVElementGeometry<i>& fvGeometry,
- const typename FVElementGeometry<i>::SubControlVolumeFace& scvf,
- const ElementVolumeVariables<i>& elemVolVars) const
- {
- static_assert(domainI != freeFlowMomentumIndex && domainJ != freeFlowMomentumIndex);
-
- const auto& couplingContext = this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
- return cm.couplingContext(ii, fvGeometry, scvf);
- });
-
- const auto& freeFlowElement = [&]
- {
- if constexpr (domainI == freeFlowMassIndex)
- return fvGeometry.element();
- else
- return couplingContext.fvGeometry.element();
- }();
-
- const auto& freeFlowScvf = [&]
- {
- if constexpr (domainI == freeFlowMassIndex)
- return scvf;
- else
- return couplingContext.fvGeometry.scvf(couplingContext.freeFlowMassScvfIdx);
-
- }();
-
- // todo revise velocity (see ff mom pm mgr)
-
- couplingContext.velocity = this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).faceVelocity(freeFlowElement, freeFlowScvf);
- return CouplingConditions::massCouplingCondition(domainI, domainJ, fvGeometry, scvf, elemVolVars, couplingContext);
- }
-
-
- //////////////////////// Conditions for FreeFlowMomentum - PorousMedium coupling //////////
- ///////////////////////////////////////////////////////////////////////////////////////////
-
- NumEqVector<freeFlowMomentumIndex> momentumCouplingCondition(Dune::index_constant<freeFlowMomentumIndex> domainI,
- Dune::index_constant<porousMediumIndex> domainJ,
- const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const typename FVElementGeometry<freeFlowMomentumIndex>::SubControlVolumeFace& scvf,
- const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars) const
- {
- if (scvf.isLateral())
- return NumEqVector<freeFlowMomentumIndex>(0.0);
-
- const auto& context = this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).couplingContext(
- domainI, fvGeometry, scvf
- );
-
- return CouplingConditions::momentumCouplingCondition(fvGeometry, scvf, elemVolVars, context);
- }
-
- /*!
- * \brief Returns the intrinsic permeability of the coupled Darcy element.
- */
- auto darcyPermeability(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
- {
- if (scvf.isFrontal())
- {
- const auto& context = this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).couplingContext(
- Dune::index_constant<freeFlowMomentumIndex>(), fvGeometry, scvf
- );
-
- return CouplingConditions::darcyPermeability(fvGeometry, scvf, context);
- }
- else
- {
- const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf);
- const auto& orthogonalScv = fvGeometry.scv(orthogonalScvf.insideScvIdx());
- const auto& frontalScvfOnBoundary = fvGeometry.frontalScvfOnBoundary(orthogonalScv);
- const auto& context = this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).couplingContext(
- Dune::index_constant<freeFlowMomentumIndex>(), fvGeometry, frontalScvfOnBoundary
- );
-
- return CouplingConditions::darcyPermeability(fvGeometry, frontalScvfOnBoundary, context);
- }
- }
-
- //////////////////////// Conditions for FreeFlowMomentum - FreeFlowMass coupling //////////
- ///////////////////////////////////////////////////////////////////////////////////////////
-
- /*!
- * \brief Returns the pressure at a given sub control volume face.
- */
- Scalar pressure(const Element<freeFlowMomentumIndex>& element,
- const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
- {
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).pressure(
- element, fvGeometry, scvf
- );
- }
-
- /*!
- * \brief Returns the pressure at the center of a sub control volume corresponding to a given sub control volume face.
- * This is used for setting a Dirichlet pressure for the mass model when a fixed pressure for the momentum balance is set at another
- * boundary. Since the the pressure at the given scvf is solution-dependent and thus unknown a priori, we just use the value
- * of the interior cell here.
- */
- Scalar cellPressure(const Element<freeFlowMomentumIndex>& element,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
- {
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).cellPressure(
- element, scvf
- );
- }
-
- /*!
- * \brief Returns the density at a given sub control volume face.
- */
- Scalar density(const Element<freeFlowMomentumIndex>& element,
- const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
- const bool considerPreviousTimeStep = false) const
- {
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).density(
- element, fvGeometry, scvf, considerPreviousTimeStep
- );
- }
-
- auto insideAndOutsideDensity(const Element<freeFlowMomentumIndex>& element,
- const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
- const bool considerPreviousTimeStep = false) const
- {
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).insideAndOutsideDensity(
- element, fvGeometry, scvf, considerPreviousTimeStep
- );
- }
-
- /*!
- * \brief Returns the density at a given sub control volume.
- */
- Scalar density(const Element<freeFlowMomentumIndex>& element,
- const SubControlVolume<freeFlowMomentumIndex>& scv,
- const bool considerPreviousTimeStep = false) const
- {
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).density(
- element, scv, considerPreviousTimeStep
- );
- }
-
- /*!
- * \brief Returns the pressure at a given sub control volume face.
- */
- Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element,
- const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
- {
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).effectiveViscosity(
- element, fvGeometry, scvf
- );
- }
-
- /*!
- * \brief Returns the velocity at a given sub control volume face.
- */
- auto faceVelocity(const Element<freeFlowMassIndex>& element,
- const SubControlVolumeFace<freeFlowMassIndex>& scvf) const
- {
- return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).faceVelocity(
- element, scvf
- );
- }
-
- template<std::size_t i>
- const Problem<i>& problem(Dune::index_constant<i> domainI) const
- {
- return this->subApply(domainI, [&](const auto& cm, auto&& ii) -> const auto& {
- return cm.problem(ii);
- });
- }
-
- template<std::size_t i, std::size_t j>
- bool isCoupled(Dune::index_constant<i> domainI,
- Dune::index_constant<j> domainJ,
- const SubControlVolumeFace<i>& scvf) const
- {
- return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
- return cm.isCoupled(ii, scvf);
- });
- }
-
- /*!
- * \brief If the boundary entity is on a coupling boundary
- * \param domainI the domain index of domain i for which to compute the flux
- * \param domainJ the domain index of domain j for which to compute the flux
- * \param scv the sub control volume
- */
- template<std::size_t i, std::size_t j>
- bool isCoupled(Dune::index_constant<i> domainI,
- Dune::index_constant<j> domainJ,
- const SubControlVolume<i>& scv) const
- {
- return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
- return cm.isCoupled(ii, scv);
- });
- }
-
- /*!
- * \brief Returns whether a given scvf is coupled to the other domain
- */
- bool isCoupledLateralScvf(Dune::index_constant<freeFlowMomentumIndex> domainI,
- Dune::index_constant<porousMediumIndex> domainJ,
- const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
- {
- return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
- return cm.isCoupledLateralScvf(ii, scvf);
- });
- }
-
-
- using ParentType::couplingStencil;
- /*!
- * \brief returns an iterable container of all indices of degrees of freedom of domain j
- * that couple with / influence the residual of the given sub-control volume of domain i
- *
- * \param domainI the domain index of domain i
- * \param elementI the coupled element of domain Ã
- * \param scvI the sub-control volume of domain i
- * \param domainJ the domain index of domain j
- */
- template<std::size_t j>
- const auto& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI,
- const Element<freeFlowMomentumIndex>& elementI,
- const SubControlVolume<freeFlowMomentumIndex>& scvI,
- Dune::index_constant<j> domainJ) const
- {
- static_assert(freeFlowMomentumIndex != j);
- return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
- return cm.couplingStencil(ii, elementI, scvI, jj);
- });
- }
-};
+using FreeFlowPorousMediumCouplingManager = typename FreeFlowPorousMediumDetail::FreeFlowPorousMediumCouplingManagerSelector<MDTraits>::type;
} // end namespace Dumux
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_base.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_base.hh
new file mode 100644
index 0000000000000000000000000000000000000000..664d051fcc93102d795b7406fa5d5b9b1e9e98de
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_base.hh
@@ -0,0 +1,278 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Base class for coupling freeflow and porous medium flow models.
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_POROUSMEDIUM_COUPLINGMANAGER_BASE_HH
+#define DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_POROUSMEDIUM_COUPLINGMANAGER_BASE_HH
+
+#include <utility>
+#include <memory>
+
+#include <dune/common/indices.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh>
+#include <dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh>
+#include <dumux/multidomain/freeflow/couplingmanager.hh>
+#include <dumux/multidomain/multibinarycouplingmanager.hh>
+
+namespace Dumux {
+
+#ifndef DOXYGEN
+namespace FreeFlowPorousMediumDetail {
+
+// global subdomain indices
+static constexpr auto freeFlowMomentumIndex = Dune::index_constant<0>();
+static constexpr auto freeFlowMassIndex = Dune::index_constant<1>();
+static constexpr auto porousMediumIndex = Dune::index_constant<2>();
+
+// coupling indices
+static constexpr auto freeFlowMassToFreeFlowMomentumIndex = Dune::index_constant<0>();
+static constexpr auto freeFlowMomentumToPorousMediumIndex = Dune::index_constant<1>();
+static constexpr auto freeFlowMassToPorousMediumIndex = Dune::index_constant<2>();
+static constexpr auto noCouplingIdx = Dune::index_constant<99>();
+
+constexpr auto makeCouplingManagerMap()
+{
+ auto map = std::array<std::array<std::size_t, 3>, 3>{};
+
+ // free flow (momentum-mass)
+ map[freeFlowMomentumIndex][freeFlowMassIndex] = freeFlowMassToFreeFlowMomentumIndex;
+ map[freeFlowMassIndex][freeFlowMomentumIndex] = freeFlowMassToFreeFlowMomentumIndex;
+
+ // free flow momentum - porous medium
+ map[freeFlowMomentumIndex][porousMediumIndex] = freeFlowMomentumToPorousMediumIndex;
+ map[porousMediumIndex][freeFlowMomentumIndex] = freeFlowMomentumToPorousMediumIndex;
+
+ // free flow mass - porous medium
+ map[freeFlowMassIndex][porousMediumIndex] = freeFlowMassToPorousMediumIndex;
+ map[porousMediumIndex][freeFlowMassIndex] = freeFlowMassToPorousMediumIndex;
+
+ return map;
+}
+
+template<std::size_t i>
+constexpr auto coupledDomains(Dune::index_constant<i> domainI)
+{
+ if constexpr (i == freeFlowMomentumIndex)
+ return std::make_tuple(freeFlowMassIndex, porousMediumIndex);
+ else if constexpr (i == freeFlowMassIndex)
+ return std::make_tuple(freeFlowMomentumIndex, porousMediumIndex);
+ else // i == porousMediumIndex
+ return std::make_tuple(freeFlowMomentumIndex, freeFlowMassIndex);
+}
+
+template<std::size_t i, std::size_t j>
+constexpr auto globalToLocalDomainIndices(Dune::index_constant<i>, Dune::index_constant<j>)
+{
+ static_assert(i <= 2 && j <= 2);
+ static_assert(i != j);
+
+ if constexpr (i < j)
+ return std::pair<Dune::index_constant<0>, Dune::index_constant<1>>{};
+ else
+ return std::pair<Dune::index_constant<1>, Dune::index_constant<0>>{};
+}
+
+struct CouplingMaps
+{
+ static constexpr auto managerMap()
+ {
+ return FreeFlowPorousMediumDetail::makeCouplingManagerMap();
+ }
+
+ template<std::size_t i, std::size_t j>
+ static constexpr auto globalToLocal(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
+ {
+ return FreeFlowPorousMediumDetail::globalToLocalDomainIndices(domainI, domainJ);
+ }
+
+ template<std::size_t i>
+ static constexpr auto coupledDomains(Dune::index_constant<i> domainI)
+ {
+ return FreeFlowPorousMediumDetail::coupledDomains(domainI);
+ }
+};
+
+template<class MDTraits>
+struct CouplingManagers
+{
+ template<std::size_t id>
+ using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+
+ using FreeFlowTraits = MultiDomainTraits<
+ SubDomainTypeTag<freeFlowMomentumIndex>, SubDomainTypeTag<freeFlowMassIndex>
+ >;
+
+ using FreeFlowMomentumPorousMediumTraits = MultiDomainTraits<
+ SubDomainTypeTag<freeFlowMomentumIndex>, SubDomainTypeTag<porousMediumIndex>
+ >;
+
+ using FreeFlowMassPorousMediumTraits = MultiDomainTraits<
+ SubDomainTypeTag<freeFlowMassIndex>, SubDomainTypeTag<porousMediumIndex>
+ >;
+
+ using FreeFlowCouplingManager
+ = Dumux::FreeFlowCouplingManager<FreeFlowTraits>;
+ using FreeFlowMomentumPorousMediumCouplingManager
+ = Dumux::FreeFlowMomentumPorousMediumCouplingManager<FreeFlowMomentumPorousMediumTraits>;
+ using FreeFlowMassPorousMediumCouplingManager
+ = Dumux::FreeFlowMassPorousMediumCouplingManager<FreeFlowMassPorousMediumTraits>;
+};
+
+} // end namespace FreeFlowPorousMediumDetail
+#endif // DOXYGEN
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Base coupling manager for coupling freeflow and porous medium flow models
+ */
+template<class MDTraits>
+class FreeFlowPorousMediumCouplingManagerBase
+: public MultiBinaryCouplingManager<
+ MDTraits,
+ FreeFlowPorousMediumDetail::CouplingMaps,
+ typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowCouplingManager,
+ typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMomentumPorousMediumCouplingManager,
+ typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMassPorousMediumCouplingManager
+>
+{
+ using ParentType = MultiBinaryCouplingManager<
+ MDTraits,
+ FreeFlowPorousMediumDetail::CouplingMaps,
+ typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowCouplingManager,
+ typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMomentumPorousMediumCouplingManager,
+ typename FreeFlowPorousMediumDetail::CouplingManagers<MDTraits>::FreeFlowMassPorousMediumCouplingManager
+ >;
+
+ using Scalar = typename MDTraits::Scalar;
+
+ // the sub domain type tags
+ template<std::size_t id>
+ using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+
+ template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
+ template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
+ template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
+ template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
+ template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector;
+
+ template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
+ template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ using SolutionVector = typename MDTraits::SolutionVector;
+
+public:
+
+ template<std::size_t i, std::size_t j>
+ using SubCouplingManager = typename ParentType::template SubCouplingManager<i, j>;
+
+ static constexpr auto freeFlowMomentumIndex = FreeFlowPorousMediumDetail::freeFlowMomentumIndex;
+ static constexpr auto freeFlowMassIndex = FreeFlowPorousMediumDetail::freeFlowMassIndex;
+ static constexpr auto porousMediumIndex = FreeFlowPorousMediumDetail::porousMediumIndex;
+
+public:
+ using ParentType::ParentType;
+
+ template<class GridVarsTuple>
+ void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
+ std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
+ std::shared_ptr<Problem<porousMediumIndex>> porousMediumProblem,
+ GridVarsTuple&& gridVarsTuple,
+ const SolutionVector& curSol)
+ {
+ this->updateSolution(curSol); // generic coupling manager stores tuple of shared_ptr
+
+ // initialize the binary sub coupling managers
+ typename SubCouplingManager<freeFlowMomentumIndex, freeFlowMassIndex>::SolutionVectorStorage ffSolVecTuple;
+ std::get<0>(ffSolVecTuple) = std::get<freeFlowMomentumIndex>(this->curSol());
+ std::get<1>(ffSolVecTuple) = std::get<freeFlowMassIndex>(this->curSol());
+ this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).init(
+ freeFlowMomentumProblem, freeFlowMassProblem,
+ std::make_tuple(std::get<freeFlowMomentumIndex>(gridVarsTuple), std::get<freeFlowMassIndex>(gridVarsTuple)),
+ ffSolVecTuple
+ );
+
+ typename SubCouplingManager<freeFlowMassIndex, porousMediumIndex>::SolutionVectorStorage ffMassPmSolVecTuple;
+ std::get<0>(ffMassPmSolVecTuple) = std::get<freeFlowMassIndex>(this->curSol());
+ std::get<1>(ffMassPmSolVecTuple) = std::get<porousMediumIndex>(this->curSol());
+ this->subCouplingManager(freeFlowMassIndex, porousMediumIndex).init(
+ freeFlowMassProblem, porousMediumProblem, ffMassPmSolVecTuple
+ );
+
+ typename SubCouplingManager<freeFlowMomentumIndex, porousMediumIndex>::SolutionVectorStorage ffMomentumPmSolVecTuple;
+ std::get<0>(ffMomentumPmSolVecTuple) = std::get<freeFlowMomentumIndex>(this->curSol());
+ std::get<1>(ffMomentumPmSolVecTuple) = std::get<porousMediumIndex>(this->curSol());
+ this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).init(
+ freeFlowMomentumProblem, porousMediumProblem, ffMomentumPmSolVecTuple
+ );
+ }
+
+ template<std::size_t i>
+ const Problem<i>& problem(Dune::index_constant<i> domainI) const
+ {
+ return this->subApply(domainI, [&](const auto& cm, auto&& ii) -> const auto& {
+ return cm.problem(ii);
+ });
+ }
+
+ template<std::size_t i, std::size_t j>
+ bool isCoupled(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
+ return cm.isCoupled(ii, scvf);
+ });
+ }
+
+ /*!
+ * \brief If the boundary entity is on a coupling boundary
+ * \param domainI the domain index of domain i for which to compute the flux
+ * \param domainJ the domain index of domain j for which to compute the flux
+ * \param scv the sub control volume
+ */
+ template<std::size_t i, std::size_t j>
+ bool isCoupled(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const SubControlVolume<i>& scv) const
+ {
+ return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
+ return cm.isCoupled(ii, scv);
+ });
+ }
+
+ using ParentType::couplingStencil;
+ /*!
+ * \brief returns an iterable container of all indices of degrees of freedom of domain j
+ * that couple with / influence the residual of the given sub-control volume of domain i
+ *
+ * \param domainI the domain index of domain i
+ * \param elementI the coupled element of domain Ã
+ * \param scvI the sub-control volume of domain i
+ * \param domainJ the domain index of domain j
+ */
+ template<std::size_t j>
+ const auto& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI,
+ const Element<freeFlowMomentumIndex>& elementI,
+ const SubControlVolume<freeFlowMomentumIndex>& scvI,
+ Dune::index_constant<j> domainJ) const
+ {
+ static_assert(freeFlowMomentumIndex != j);
+ return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.couplingStencil(ii, elementI, scvI, jj);
+ });
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_staggered_cctpfa.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_staggered_cctpfa.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5ca8cc1f15ab1b8eb8184b71ce89052f7ebdca1e
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_staggered_cctpfa.hh
@@ -0,0 +1,245 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \copydoc Dumux::FreeFlowPorousMediumCouplingManagerStaggeredCCTpfa
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_POROUSMEDIUM_COUPLINGMANAGER_STAGGERED_CCTPFA_HH
+#define DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_POROUSMEDIUM_COUPLINGMANAGER_STAGGERED_CCTPFA_HH
+
+#include "couplingmanager_base.hh"
+#include "couplingconditions_staggered_cctpfa.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Coupling manager for coupling freeflow and porous medium flow models
+ * specialization for staggered-cctpfa coupling.
+ */
+template<class MDTraits>
+class FreeFlowPorousMediumCouplingManagerStaggeredCCTpfa
+: public FreeFlowPorousMediumCouplingManagerBase<MDTraits>
+{
+ using ParentType = FreeFlowPorousMediumCouplingManagerBase<MDTraits>;
+
+ using Scalar = typename MDTraits::Scalar;
+
+ // the sub domain type tags
+ template<std::size_t id>
+ using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+
+ template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
+ template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
+ template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
+ template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
+ template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector;
+
+ template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
+ template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ using SolutionVector = typename MDTraits::SolutionVector;
+
+ using CouplingConditions = FFPMCouplingConditionsStaggeredCCTpfa<MDTraits, FreeFlowPorousMediumCouplingManagerStaggeredCCTpfa<MDTraits>>;
+
+public:
+ static constexpr auto freeFlowMomentumIndex = ParentType::freeFlowMomentumIndex;
+ static constexpr auto freeFlowMassIndex = ParentType::freeFlowMassIndex;
+ static constexpr auto porousMediumIndex = ParentType::porousMediumIndex;
+
+public:
+ /*!
+ * \brief Returns the mass flux across the coupling boundary.
+ */
+ template<std::size_t i, std::size_t j>
+ auto massCouplingCondition(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ,
+ const FVElementGeometry<i>& fvGeometry,
+ const typename FVElementGeometry<i>::SubControlVolumeFace& scvf,
+ const ElementVolumeVariables<i>& elemVolVars) const
+ {
+ static_assert(domainI != freeFlowMomentumIndex && domainJ != freeFlowMomentumIndex);
+
+ const auto& couplingContext = this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.couplingContext(ii, fvGeometry, scvf);
+ });
+
+ const auto& freeFlowElement = [&]
+ {
+ if constexpr (domainI == freeFlowMassIndex)
+ return fvGeometry.element();
+ else
+ return couplingContext.fvGeometry.element();
+ }();
+
+ const auto& freeFlowScvf = [&]
+ {
+ if constexpr (domainI == freeFlowMassIndex)
+ return scvf;
+ else
+ return couplingContext.fvGeometry.scvf(couplingContext.freeFlowMassScvfIdx);
+
+ }();
+
+ // todo revise velocity (see ff mom pm mgr)
+
+ couplingContext.velocity = this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).faceVelocity(freeFlowElement, freeFlowScvf);
+ return CouplingConditions::massCouplingCondition(domainI, domainJ, fvGeometry, scvf, elemVolVars, couplingContext);
+ }
+
+
+ //////////////////////// Conditions for FreeFlowMomentum - PorousMedium coupling //////////
+ ///////////////////////////////////////////////////////////////////////////////////////////
+
+ NumEqVector<freeFlowMomentumIndex> momentumCouplingCondition(Dune::index_constant<freeFlowMomentumIndex> domainI,
+ Dune::index_constant<porousMediumIndex> domainJ,
+ const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const typename FVElementGeometry<freeFlowMomentumIndex>::SubControlVolumeFace& scvf,
+ const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars) const
+ {
+ if (scvf.isLateral())
+ return NumEqVector<freeFlowMomentumIndex>(0.0);
+
+ const auto& context = this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).couplingContext(
+ domainI, fvGeometry, scvf
+ );
+
+ return CouplingConditions::momentumCouplingCondition(fvGeometry, scvf, elemVolVars, context);
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability of the coupled Darcy element.
+ */
+ auto darcyPermeability(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
+ {
+ if (scvf.isFrontal())
+ {
+ const auto& context = this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).couplingContext(
+ Dune::index_constant<freeFlowMomentumIndex>(), fvGeometry, scvf
+ );
+
+ return CouplingConditions::darcyPermeability(fvGeometry, scvf, context);
+ }
+ else
+ {
+ const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf);
+ const auto& orthogonalScv = fvGeometry.scv(orthogonalScvf.insideScvIdx());
+ const auto& frontalScvfOnBoundary = fvGeometry.frontalScvfOnBoundary(orthogonalScv);
+ const auto& context = this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).couplingContext(
+ Dune::index_constant<freeFlowMomentumIndex>(), fvGeometry, frontalScvfOnBoundary
+ );
+
+ return CouplingConditions::darcyPermeability(fvGeometry, frontalScvfOnBoundary, context);
+ }
+ }
+
+ //////////////////////// Conditions for FreeFlowMomentum - FreeFlowMass coupling //////////
+ ///////////////////////////////////////////////////////////////////////////////////////////
+
+ /*!
+ * \brief Returns the pressure at a given sub control volume face.
+ */
+ Scalar pressure(const Element<freeFlowMomentumIndex>& element,
+ const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
+ {
+ return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).pressure(
+ element, fvGeometry, scvf
+ );
+ }
+
+ /*!
+ * \brief Returns the pressure at the center of a sub control volume corresponding to a given sub control volume face.
+ * This is used for setting a Dirichlet pressure for the mass model when a fixed pressure for the momentum balance is set at another
+ * boundary. Since the the pressure at the given scvf is solution-dependent and thus unknown a priori, we just use the value
+ * of the interior cell here.
+ */
+ Scalar cellPressure(const Element<freeFlowMomentumIndex>& element,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
+ {
+ return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).cellPressure(
+ element, scvf
+ );
+ }
+
+ /*!
+ * \brief Returns the density at a given sub control volume face.
+ */
+ Scalar density(const Element<freeFlowMomentumIndex>& element,
+ const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
+ const bool considerPreviousTimeStep = false) const
+ {
+ return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).density(
+ element, fvGeometry, scvf, considerPreviousTimeStep
+ );
+ }
+
+ auto insideAndOutsideDensity(const Element<freeFlowMomentumIndex>& element,
+ const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
+ const bool considerPreviousTimeStep = false) const
+ {
+ return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).insideAndOutsideDensity(
+ element, fvGeometry, scvf, considerPreviousTimeStep
+ );
+ }
+
+ /*!
+ * \brief Returns the density at a given sub control volume.
+ */
+ Scalar density(const Element<freeFlowMomentumIndex>& element,
+ const SubControlVolume<freeFlowMomentumIndex>& scv,
+ const bool considerPreviousTimeStep = false) const
+ {
+ return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).density(
+ element, scv, considerPreviousTimeStep
+ );
+ }
+
+ /*!
+ * \brief Returns the pressure at a given sub control volume face.
+ */
+ Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element,
+ const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
+ {
+ return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).effectiveViscosity(
+ element, fvGeometry, scvf
+ );
+ }
+
+ /*!
+ * \brief Returns the velocity at a given sub control volume face.
+ */
+ auto faceVelocity(const Element<freeFlowMassIndex>& element,
+ const SubControlVolumeFace<freeFlowMassIndex>& scvf) const
+ {
+ return this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).faceVelocity(
+ element, scvf
+ );
+ }
+
+ /*!
+ * \brief Returns whether a given scvf is coupled to the other domain
+ */
+ bool isCoupledLateralScvf(Dune::index_constant<freeFlowMomentumIndex> domainI,
+ Dune::index_constant<porousMediumIndex> domainJ,
+ const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
+ {
+ return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
+ return cm.isCoupledLateralScvf(ii, scvf);
+ });
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh
index 9f8c8493a5b976b6ff687e2dc39e530409834a41..c59221cc89224cbc0c5b73e5b64f6a9c7fca5ccb 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh
@@ -7,285 +7,37 @@
/*!
* \file
* \ingroup FreeFlowPorousMediumCoupling
- * \copydoc Dumux::FreeFlowMassPorousMediumCouplingManager
+ * \brief Coupling managers specialized for different discretization schemes for mass coupling
*/
#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMASSPM_COUPLINGMANAGER_HH
#define DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMASSPM_COUPLINGMANAGER_HH
-#include <utility>
-#include <memory>
+#include <dumux/discretization/method.hh>
-#include <dune/common/float_cmp.hh>
-#include <dune/common/exceptions.hh>
-#include <dumux/common/properties.hh>
-#include <dumux/discretization/staggered/elementsolution.hh>
-#include <dumux/multidomain/couplingmanager.hh>
-#include <dumux/multidomain/boundary/darcydarcy/couplingmapper.hh>
+#include "couplingmanager_staggered_cctpfa.hh"
namespace Dumux {
-/*!
- * \ingroup FreeFlowPorousMediumCoupling
- * \brief Coupling manager for Stokes and Darcy domains with equal dimension.
- */
-template<class MDTraits>
-class FreeFlowMassPorousMediumCouplingManager
-: public CouplingManager<MDTraits>
-{
- using Scalar = typename MDTraits::Scalar;
- using ParentType = CouplingManager<MDTraits>;
-
-public:
- static constexpr auto freeFlowMassIndex = typename MDTraits::template SubDomain<0>::Index();
- static constexpr auto porousMediumIndex = typename MDTraits::template SubDomain<1>::Index();
-
- using SolutionVectorStorage = typename ParentType::SolutionVectorStorage;
-private:
-
- // obtain the type tags of the sub problems
- using FreeFlowMassTypeTag = typename MDTraits::template SubDomain<freeFlowMassIndex>::TypeTag;
- using PorousMediumTypeTag = typename MDTraits::template SubDomain<porousMediumIndex>::TypeTag;
-
- using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >;
- using CouplingStencil = CouplingStencils::mapped_type;
-
- // the sub domain type tags
- template<std::size_t id>
- using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
-
- template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView;
- template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
- template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
- template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>;
- template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
- template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
- template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
- template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
- template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
- template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
- template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
- template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
-
- using VelocityVector = typename Element<freeFlowMassIndex>::Geometry::GlobalCoordinate;
-
- struct FreeFlowMassCouplingContext
- {
- Element<porousMediumIndex> element;
- VolumeVariables<porousMediumIndex> volVars;
- FVElementGeometry<porousMediumIndex> fvGeometry;
- std::size_t dofIdx;
- std::size_t freeFlowMassScvfIdx;
- std::size_t porousMediumScvfIdx;
- mutable VelocityVector velocity; // velocity needs to be set externally, not available in this class
- };
-
- struct PorousMediumCouplingContext
- {
- Element<freeFlowMassIndex> element;
- VolumeVariables<freeFlowMassIndex> volVars;
- FVElementGeometry<freeFlowMassIndex> fvGeometry;
- std::size_t dofIdx;
- std::size_t porousMediumScvfIdx;
- std::size_t freeFlowMassScvfIdx;
- mutable VelocityVector velocity; // velocity needs to be set externally, not available in this class
- };
-
- using CouplingMapper = DarcyDarcyBoundaryCouplingMapper<MDTraits>; // TODO rename/generalize class
-
-public:
-
- /*!
- * \brief Methods to be accessed by main
- */
- // \{
-
- //! Initialize the coupling manager
- void init(std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
- std::shared_ptr<Problem<porousMediumIndex>> darcyProblem,
- SolutionVectorStorage& curSol)
- {
- this->setSubProblems(std::make_tuple(freeFlowMassProblem, darcyProblem));
- this->attachSolution(curSol);
-
- couplingMapper_.update(*this);
- }
-
- // \}
-
- /*!
- * \brief Methods to be accessed by the assembly
- */
- // \{
-
- /*!
- * \brief prepares all data and variables that are necessary to evaluate the residual (called from the local assembler)
- */
- template<std::size_t i, class Assembler>
- void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler) const
- {
- bindCouplingContext_(domainI, element);
- }
+#ifndef DOXYGEN
+namespace FreeFlowMassPorousMediumDetail {
- /*!
- * \brief Update the coupling context
- */
- template<std::size_t i, std::size_t j, class LocalAssemblerI>
- void updateCouplingContext(Dune::index_constant<i> domainI,
- const LocalAssemblerI& localAssemblerI,
- Dune::index_constant<j> domainJ,
- std::size_t dofIdxGlobalJ,
- const PrimaryVariables<j>& priVarsJ,
- int pvIdxJ)
- {
- this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
+// declaration (specialize for different discretization types)
+template<class MDTraits,
+ class DiscFFMass = typename MDTraits::template SubDomain<0>::GridGeometry::DiscretizationMethod,
+ class DiscPM = typename MDTraits::template SubDomain<1>::GridGeometry::DiscretizationMethod
+ >
+struct FreeFlowMassPorousMediumCouplingManagerSelector;
- // we need to update all solution-depenent components of the coupling context
- // the dof of domain J has been deflected
- // if domainJ == freeFlowMassIndex: update volvars in the PorousMediumCouplingContext
- // if domainJ == porousMediumIndex: update volvars in the FreeFlowMassCouplingContext
- // as the update is symmetric we only need to write this once
- auto& context = std::get<1-j>(couplingContext_);
- for (auto& c : context)
- {
- if (c.dofIdx == dofIdxGlobalJ)
- {
- const auto elemSol = elementSolution(c.element, this->curSol(domainJ), this->problem(domainJ).gridGeometry());
- const auto& scv = *scvs(c.fvGeometry).begin();
- c.volVars.update(elemSol, this->problem(domainJ), c.element, scv);
- }
- }
- }
-
- // \}
-
- /*!
- * \brief Access the coupling context needed for the Stokes domain
- */
- template<std::size_t i>
- const auto& couplingContext(Dune::index_constant<i> domainI,
- const FVElementGeometry<i>& fvGeometry,
- const SubControlVolumeFace<i> scvf) const
- {
- auto& contexts = std::get<i>(couplingContext_);
-
- if (contexts.empty() || couplingContextBoundForElement_[i] != scvf.insideScvIdx())
- bindCouplingContext_(domainI, fvGeometry);
-
- for (const auto& context : contexts)
- {
- const auto expectedScvfIdx = domainI == freeFlowMassIndex ? context.freeFlowMassScvfIdx : context.porousMediumScvfIdx;
- if (scvf.index() == expectedScvfIdx)
- return context;
- }
-
- DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
- }
-
- /*!
- * \brief The coupling stencils
- */
- // \{
-
- /*!
- * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs
- */
- template<std::size_t i, std::size_t j>
- const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
- const Element<i>& element,
- Dune::index_constant<j> domainJ) const
- {
- const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
- return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
- }
-
- // \}
-
- /*!
- * \brief Returns whether a given scvf is coupled to the other domain
- */
- template<std::size_t i>
- bool isCoupled(Dune::index_constant<i> domainI, const SubControlVolumeFace<i>& scvf) const
- { return couplingMapper_.isCoupled(domainI, scvf); }
-
-private:
- /*!
- * \brief Returns whether a given scvf is coupled to the other domain
- */
- template<std::size_t i>
- bool isCoupledElement_(Dune::index_constant<i> domainI, std::size_t eIdx) const
- { return couplingMapper_.isCoupledElement(domainI, eIdx); }
-
- //! Return the volume variables of domain i for a given element and scv
- template<std::size_t i>
- VolumeVariables<i> volVars_(Dune::index_constant<i> domainI,
- const Element<i>& element,
- const SubControlVolume<i>& scv) const
- {
- VolumeVariables<i> volVars;
- const auto elemSol = elementSolution(
- element, this->curSol(domainI), this->problem(domainI).gridGeometry()
- );
- volVars.update(elemSol, this->problem(domainI), element, scv);
- return volVars;
- }
-
- /*!
- * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
- */
- template<std::size_t i>
- void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const
- {
- const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element);
- bindCouplingContext_(domainI, fvGeometry);
- }
-
- /*!
- * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
- */
- template<std::size_t i>
- void bindCouplingContext_(Dune::index_constant<i> domainI, const FVElementGeometry<i>& fvGeometry) const
- {
- auto& context = std::get<domainI>(couplingContext_);
- context.clear();
-
- const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element());
-
- // do nothing if the element is not coupled to the other domain
- if (!isCoupledElement_(domainI, eIdx))
- return;
-
- couplingContextBoundForElement_[domainI] = eIdx;
-
- for (const auto& scvf : scvfs(fvGeometry))
- {
- if (isCoupled(domainI, scvf))
- {
- const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
- constexpr auto domainJ = Dune::index_constant<1-domainI>();
- const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
- const auto& otherElement = otherGridGeometry.element(otherElementIdx);
- auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
-
- // there is only one scv for TPFA
- context.push_back({
- otherElement,
- volVars_(domainJ, otherElement, *std::begin(scvs(otherFvGeometry))),
- std::move(otherFvGeometry),
- otherElementIdx,
- scvf.index(),
- couplingMapper_.flipScvfIndex(domainI, scvf),
- VelocityVector{}
- });
- }
- }
- }
+template<class MDTraits>
+struct FreeFlowMassPorousMediumCouplingManagerSelector<MDTraits, DiscretizationMethods::CCTpfa, DiscretizationMethods::CCTpfa>
+{ using type = FFMassPMCouplingManagerStaggeredCCTpfa<MDTraits>; };
- mutable std::tuple<std::vector<FreeFlowMassCouplingContext>, std::vector<PorousMediumCouplingContext>> couplingContext_;
- mutable std::array<std::size_t, 2> couplingContextBoundForElement_;
+} // end namespace FreeFlowMassPorousMediumDetail
+#endif // DOXYGEN
- CouplingMapper couplingMapper_;
-};
+template<class MDTraits>
+using FreeFlowMassPorousMediumCouplingManager = typename FreeFlowMassPorousMediumDetail::FreeFlowMassPorousMediumCouplingManagerSelector<MDTraits>::type;
} // end namespace Dumux
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager_staggered_cctpfa.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager_staggered_cctpfa.hh
new file mode 100644
index 0000000000000000000000000000000000000000..16fb3d0751f34a5db40540d09809883cd89b6f6f
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager_staggered_cctpfa.hh
@@ -0,0 +1,293 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \copydoc Dumux::FFMassPMCouplingManagerStaggeredCCTpfa
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMASSPM_COUPLINGMANAGER_STAGGERED_TPFA_HH
+#define DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMASSPM_COUPLINGMANAGER_STAGGERED_TPFA_HH
+
+#include <utility>
+#include <memory>
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/exceptions.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/staggered/elementsolution.hh>
+#include <dumux/multidomain/couplingmanager.hh>
+#include <dumux/multidomain/boundary/darcydarcy/couplingmapper.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Coupling manager for Stokes and Darcy domains with equal dimension.
+ * Specialization for staggered-cctpfa coupling.
+ */
+template<class MDTraits>
+class FFMassPMCouplingManagerStaggeredCCTpfa
+: public CouplingManager<MDTraits>
+{
+ using Scalar = typename MDTraits::Scalar;
+ using ParentType = CouplingManager<MDTraits>;
+
+public:
+ static constexpr auto freeFlowMassIndex = typename MDTraits::template SubDomain<0>::Index();
+ static constexpr auto porousMediumIndex = typename MDTraits::template SubDomain<1>::Index();
+
+ using SolutionVectorStorage = typename ParentType::SolutionVectorStorage;
+private:
+
+ // obtain the type tags of the sub problems
+ using FreeFlowMassTypeTag = typename MDTraits::template SubDomain<freeFlowMassIndex>::TypeTag;
+ using PorousMediumTypeTag = typename MDTraits::template SubDomain<porousMediumIndex>::TypeTag;
+
+ using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >;
+ using CouplingStencil = CouplingStencils::mapped_type;
+
+ // the sub domain type tags
+ template<std::size_t id>
+ using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+
+ template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView;
+ template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+ template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
+ template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>;
+ template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
+ template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
+ template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
+ template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
+ template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
+
+ using VelocityVector = typename Element<freeFlowMassIndex>::Geometry::GlobalCoordinate;
+
+ struct FreeFlowMassCouplingContext
+ {
+ Element<porousMediumIndex> element;
+ VolumeVariables<porousMediumIndex> volVars;
+ FVElementGeometry<porousMediumIndex> fvGeometry;
+ std::size_t dofIdx;
+ std::size_t freeFlowMassScvfIdx;
+ std::size_t porousMediumScvfIdx;
+ mutable VelocityVector velocity; // velocity needs to be set externally, not available in this class
+ };
+
+ struct PorousMediumCouplingContext
+ {
+ Element<freeFlowMassIndex> element;
+ VolumeVariables<freeFlowMassIndex> volVars;
+ FVElementGeometry<freeFlowMassIndex> fvGeometry;
+ std::size_t dofIdx;
+ std::size_t porousMediumScvfIdx;
+ std::size_t freeFlowMassScvfIdx;
+ mutable VelocityVector velocity; // velocity needs to be set externally, not available in this class
+ };
+
+ using CouplingMapper = DarcyDarcyBoundaryCouplingMapper<MDTraits>; // TODO rename/generalize class
+
+public:
+
+ /*!
+ * \brief Methods to be accessed by main
+ */
+ // \{
+
+ //! Initialize the coupling manager
+ void init(std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
+ std::shared_ptr<Problem<porousMediumIndex>> darcyProblem,
+ SolutionVectorStorage& curSol)
+ {
+ this->setSubProblems(std::make_tuple(freeFlowMassProblem, darcyProblem));
+ this->attachSolution(curSol);
+
+ couplingMapper_.update(*this);
+ }
+
+ // \}
+
+ /*!
+ * \brief Methods to be accessed by the assembly
+ */
+ // \{
+
+ /*!
+ * \brief prepares all data and variables that are necessary to evaluate the residual (called from the local assembler)
+ */
+ template<std::size_t i, class Assembler>
+ void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler) const
+ {
+ bindCouplingContext_(domainI, element);
+ }
+
+ /*!
+ * \brief Update the coupling context
+ */
+ template<std::size_t i, std::size_t j, class LocalAssemblerI>
+ void updateCouplingContext(Dune::index_constant<i> domainI,
+ const LocalAssemblerI& localAssemblerI,
+ Dune::index_constant<j> domainJ,
+ std::size_t dofIdxGlobalJ,
+ const PrimaryVariables<j>& priVarsJ,
+ int pvIdxJ)
+ {
+ this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
+
+ // we need to update all solution-depenent components of the coupling context
+ // the dof of domain J has been deflected
+ // if domainJ == freeFlowMassIndex: update volvars in the PorousMediumCouplingContext
+ // if domainJ == porousMediumIndex: update volvars in the FreeFlowMassCouplingContext
+ // as the update is symmetric we only need to write this once
+ auto& context = std::get<1-j>(couplingContext_);
+ for (auto& c : context)
+ {
+ if (c.dofIdx == dofIdxGlobalJ)
+ {
+ const auto elemSol = elementSolution(c.element, this->curSol(domainJ), this->problem(domainJ).gridGeometry());
+ const auto& scv = *scvs(c.fvGeometry).begin();
+ c.volVars.update(elemSol, this->problem(domainJ), c.element, scv);
+ }
+ }
+ }
+
+ // \}
+
+ /*!
+ * \brief Access the coupling context needed for the Stokes domain
+ */
+ template<std::size_t i>
+ const auto& couplingContext(Dune::index_constant<i> domainI,
+ const FVElementGeometry<i>& fvGeometry,
+ const SubControlVolumeFace<i> scvf) const
+ {
+ auto& contexts = std::get<i>(couplingContext_);
+
+ if (contexts.empty() || couplingContextBoundForElement_[i] != scvf.insideScvIdx())
+ bindCouplingContext_(domainI, fvGeometry);
+
+ for (const auto& context : contexts)
+ {
+ const auto expectedScvfIdx = domainI == freeFlowMassIndex ? context.freeFlowMassScvfIdx : context.porousMediumScvfIdx;
+ if (scvf.index() == expectedScvfIdx)
+ return context;
+ }
+
+ DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
+ }
+
+ /*!
+ * \brief The coupling stencils
+ */
+ // \{
+
+ /*!
+ * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs
+ */
+ template<std::size_t i, std::size_t j>
+ const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
+ const Element<i>& element,
+ Dune::index_constant<j> domainJ) const
+ {
+ const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
+ return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
+ }
+
+ // \}
+
+ /*!
+ * \brief Returns whether a given scvf is coupled to the other domain
+ */
+ template<std::size_t i>
+ bool isCoupled(Dune::index_constant<i> domainI, const SubControlVolumeFace<i>& scvf) const
+ { return couplingMapper_.isCoupled(domainI, scvf); }
+
+private:
+ /*!
+ * \brief Returns whether a given scvf is coupled to the other domain
+ */
+ template<std::size_t i>
+ bool isCoupledElement_(Dune::index_constant<i> domainI, std::size_t eIdx) const
+ { return couplingMapper_.isCoupledElement(domainI, eIdx); }
+
+ //! Return the volume variables of domain i for a given element and scv
+ template<std::size_t i>
+ VolumeVariables<i> volVars_(Dune::index_constant<i> domainI,
+ const Element<i>& element,
+ const SubControlVolume<i>& scv) const
+ {
+ VolumeVariables<i> volVars;
+ const auto elemSol = elementSolution(
+ element, this->curSol(domainI), this->problem(domainI).gridGeometry()
+ );
+ volVars.update(elemSol, this->problem(domainI), element, scv);
+ return volVars;
+ }
+
+ /*!
+ * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
+ */
+ template<std::size_t i>
+ void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const
+ {
+ const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element);
+ bindCouplingContext_(domainI, fvGeometry);
+ }
+
+ /*!
+ * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
+ */
+ template<std::size_t i>
+ void bindCouplingContext_(Dune::index_constant<i> domainI, const FVElementGeometry<i>& fvGeometry) const
+ {
+ auto& context = std::get<domainI>(couplingContext_);
+ context.clear();
+
+ const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element());
+
+ // do nothing if the element is not coupled to the other domain
+ if (!isCoupledElement_(domainI, eIdx))
+ return;
+
+ couplingContextBoundForElement_[domainI] = eIdx;
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ if (isCoupled(domainI, scvf))
+ {
+ const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
+ constexpr auto domainJ = Dune::index_constant<1-domainI>();
+ const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
+ const auto& otherElement = otherGridGeometry.element(otherElementIdx);
+ auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
+
+ // there is only one scv for TPFA
+ context.push_back({
+ otherElement,
+ volVars_(domainJ, otherElement, *std::begin(scvs(otherFvGeometry))),
+ std::move(otherFvGeometry),
+ otherElementIdx,
+ scvf.index(),
+ couplingMapper_.flipScvfIndex(domainI, scvf),
+ VelocityVector{}
+ });
+ }
+ }
+ }
+
+ mutable std::tuple<std::vector<FreeFlowMassCouplingContext>, std::vector<PorousMediumCouplingContext>> couplingContext_;
+ mutable std::array<std::size_t, 2> couplingContextBoundForElement_;
+
+ CouplingMapper couplingMapper_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh
index d07569c68a0f4558847397a43abc1d269210f4cd..52156ab9c40688e449392c68bc708627669ddfd9 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh
@@ -7,403 +7,37 @@
/*!
* \file
* \ingroup FreeFlowPorousMediumCoupling
- * \copydoc Dumux::FreeFlowMomentumPorousMediumCouplingManager
+ * \brief Coupling managers specialized for different discretization schemes for momentum coupling
*/
#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMMOMENTUMPM_COUPLINGMANAGER_HH
#define DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMMOMENTUMPM_COUPLINGMANAGER_HH
-#include <utility>
-#include <memory>
+#include <dumux/discretization/method.hh>
-#include <dune/common/float_cmp.hh>
-#include <dune/common/exceptions.hh>
-#include <dumux/common/properties.hh>
-#include <dumux/discretization/staggered/elementsolution.hh>
-#include <dumux/multidomain/couplingmanager.hh>
-
-#include "couplingmapper.hh"
+#include "couplingmanager_staggered_cctpfa.hh"
namespace Dumux {
-/*!
- * \ingroup FreeFlowPorousMediumCoupling
- * \brief Coupling manager for Stokes and Darcy domains with equal dimension.
- */
-template<class MDTraits>
-class FreeFlowMomentumPorousMediumCouplingManager
-: public CouplingManager<MDTraits>
-{
- using Scalar = typename MDTraits::Scalar;
- using ParentType = CouplingManager<MDTraits>;
-
-public:
- static constexpr auto freeFlowMomentumIndex = typename MDTraits::template SubDomain<0>::Index();
- static constexpr auto porousMediumIndex = typename MDTraits::template SubDomain<1>::Index();
-
- using SolutionVectorStorage = typename ParentType::SolutionVectorStorage;
-private:
- // obtain the type tags of the sub problems
- using FreeFlowMomentumTypeTag = typename MDTraits::template SubDomain<freeFlowMomentumIndex>::TypeTag;
- using PorousMediumTypeTag = typename MDTraits::template SubDomain<porousMediumIndex>::TypeTag;
-
- using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >;
- using CouplingStencil = CouplingStencils::mapped_type;
-
- // the sub domain type tags
- template<std::size_t id>
- using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
-
- template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView;
- template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
- template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
- template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>;
- template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
- template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
- template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
- template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
- template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
- template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
- template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
- template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
-
- using VelocityVector = typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate;
-
- struct FreeFlowMomentumCouplingContext
- {
- Element<porousMediumIndex> element;
- VolumeVariables<porousMediumIndex> volVars;
- FVElementGeometry<porousMediumIndex> fvGeometry;
- std::size_t freeFlowMomentumScvfIdx;
- std::size_t porousMediumScvfIdx;
- std::size_t porousMediumDofIdx;
- VelocityVector gravity;
- };
-
- struct PorousMediumCouplingContext
- {
- Element<freeFlowMomentumIndex> element;
- FVElementGeometry<freeFlowMomentumIndex> fvGeometry;
- std::size_t porousMediumScvfIdx;
- std::size_t freeFlowMomentumScvfIdx;
- std::size_t freeFlowMomentumDofIdx;
- VelocityVector faceVelocity;
- };
-
- using CouplingMapper = FreeFlowMomentumPorousMediumCouplingMapper<MDTraits, FreeFlowMomentumPorousMediumCouplingManager<MDTraits>>;
-
-public:
-
- /*!
- * \brief Methods to be accessed by main
- */
- // \{
-
- //! Initialize the coupling manager
- void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
- std::shared_ptr<Problem<porousMediumIndex>> porousMediumProblem,
- SolutionVectorStorage& curSol)
- {
- this->setSubProblems(std::make_tuple(freeFlowMomentumProblem, porousMediumProblem));
- this->attachSolution(curSol);
-
- couplingMapper_.update(*this);
- }
-
- // \}
-
-
- /*!
- * \brief Methods to be accessed by the assembly
- */
- // \{
-
- /*!
- * \brief prepares all data and variables that are necessary to evaluate the residual (called from the local assembler)
- */
- template<std::size_t i, class Assembler>
- void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler) const
- {
- bindCouplingContext_(domainI, element);
- }
-
- /*!
- * \brief Update the coupling context
- */
- template<std::size_t i, std::size_t j, class LocalAssemblerI>
- void updateCouplingContext(Dune::index_constant<i> domainI,
- const LocalAssemblerI& localAssemblerI,
- Dune::index_constant<j> domainJ,
- std::size_t dofIdxGlobalJ,
- const PrimaryVariables<j>& priVarsJ,
- int pvIdxJ)
- {
- this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
-
- // we need to update all solution-depenent components of the coupling context
- // the dof of domain J has been deflected
-
- // update the faceVelocity in the PorousMediumCouplingContext
- if constexpr (domainJ == freeFlowMomentumIndex)
- {
- // we only need to update if we are assembling the porous medium domain
- // since the freeflow domain will not use the velocity from the context
- if constexpr (domainI == porousMediumIndex)
- {
- auto& context = std::get<porousMediumIndex>(couplingContext_);
- for (auto& c : context)
- {
- if (c.freeFlowMomentumDofIdx == dofIdxGlobalJ)
- {
- const auto& scvf = localAssemblerI.fvGeometry().scvf(c.porousMediumScvfIdx);
- c.faceVelocity = faceVelocity(localAssemblerI.element(), scvf);
- }
- }
- }
- }
-
- // update volVars in the FreeFlowMomentumCouplingContext
- else if (domainJ == porousMediumIndex)
- {
- auto& context = std::get<freeFlowMomentumIndex>(couplingContext_);
- for (auto& c : context)
- {
- if (c.porousMediumDofIdx == dofIdxGlobalJ)
- {
- const auto& ggJ = c.fvGeometry.gridGeometry();
- const auto& scv = *scvs(c.fvGeometry).begin();
- const auto elemSol = elementSolution(c.element, this->curSol(domainJ), ggJ);
- c.volVars.update(elemSol, this->problem(domainJ), c.element, scv);
- }
- }
- }
- }
-
- // \}
-
- /*!
- * \brief Access the coupling context needed for the Stokes domain
- */
- template<std::size_t i>
- const auto& couplingContext(Dune::index_constant<i> domainI,
- const FVElementGeometry<i>& fvGeometry,
- const SubControlVolumeFace<i> scvf) const
- {
- auto& contexts = std::get<i>(couplingContext_);
+#ifndef DOXYGEN
+namespace FreeFlowMomentumPorousMediumDetail {
- if (contexts.empty() || couplingContextBoundForElement_[i] != scvf.insideScvIdx())
- bindCouplingContext_(domainI, fvGeometry);
+// declaration (specialize for different discretization types)
+template<class MDTraits,
+ class DiscFFMomentum = typename MDTraits::template SubDomain<0>::GridGeometry::DiscretizationMethod,
+ class DiscPM = typename MDTraits::template SubDomain<1>::GridGeometry::DiscretizationMethod
+ >
+struct FreeFlowMomentumPorousMediumCouplingManagerSelector;
- for (const auto& context : contexts)
- {
- const auto expectedScvfIdx = domainI == freeFlowMomentumIndex ? context.freeFlowMomentumScvfIdx : context.porousMediumScvfIdx;
- if (scvf.index() == expectedScvfIdx)
- return context;
- }
-
- DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
- }
-
- /*!
- * \brief The coupling stencils
- */
- // \{
-
- /*!
- * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs
- */
- const CouplingStencil& couplingStencil(Dune::index_constant<porousMediumIndex> domainI,
- const Element<porousMediumIndex>& element,
- Dune::index_constant<freeFlowMomentumIndex> domainJ) const
- {
- const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
- return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
- }
-
- /*!
- * \brief returns an iterable container of all indices of degrees of freedom of domain j
- * that couple with / influence the residual of the given sub-control volume of domain i
- *
- * \param domainI the domain index of domain i
- * \param elementI the coupled element of domain Ã
- * \param scvI the sub-control volume of domain i
- * \param domainJ the domain index of domain j
- */
- const CouplingStencil& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI,
- const Element<freeFlowMomentumIndex>& elementI,
- const SubControlVolume<freeFlowMomentumIndex>& scvI,
- Dune::index_constant<porousMediumIndex> domainJ) const
- {
- return couplingMapper_.couplingStencil(domainI, elementI, scvI, domainJ);
- }
-
- using ParentType::evalCouplingResidual;
-
- /*!
- * \brief evaluate the coupling residual
- * special interface for fcstaggered methods
- */
- template<class LocalAssemblerI, std::size_t j>
- decltype(auto) evalCouplingResidual(Dune::index_constant<freeFlowMomentumIndex> domainI,
- const LocalAssemblerI& localAssemblerI,
- const SubControlVolume<freeFlowMomentumIndex>& scvI,
- Dune::index_constant<j> domainJ,
- std::size_t dofIdxGlobalJ) const
- {
- return localAssemblerI.evalLocalResidual();
- }
-
- // \}
-
- /*!
- * \brief Returns whether a given scvf is coupled to the other domain
- */
- bool isCoupledLateralScvf(Dune::index_constant<freeFlowMomentumIndex> domainI, const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
- { return couplingMapper_.isCoupledLateralScvf(domainI, scvf); }
-
- /*!
- * \brief Returns whether a given scvf is coupled to the other domain
- */
- template<std::size_t i>
- bool isCoupled(Dune::index_constant<i> domainI, const SubControlVolumeFace<i>& scvf) const
- {
- if constexpr (i == freeFlowMomentumIndex)
- return couplingMapper_.isCoupled(domainI, scvf) || couplingMapper_.isCoupledLateralScvf(domainI, scvf);
- else
- return couplingMapper_.isCoupled(domainI, scvf);
- }
-
- /*!
- * \brief If the boundary entity is on a coupling boundary
- * \param domainI the domain index for which to compute the flux
- * \param scv the sub control volume
- */
- bool isCoupled(Dune::index_constant<freeFlowMomentumIndex> domainI,
- const SubControlVolume<freeFlowMomentumIndex>& scv) const
- { return couplingMapper_.isCoupled(domainI, scv); }
-
- /*!
- * \brief Returns the velocity at a given sub control volume face.
- */
- auto faceVelocity(const Element<porousMediumIndex>& element,
- const SubControlVolumeFace<porousMediumIndex>& scvf) const
- {
- // create a unit normal vector oriented in positive coordinate direction
- auto velocity = scvf.unitOuterNormal();
- using std::abs;
- std::for_each(velocity.begin(), velocity.end(), [](auto& v){ v = abs(v); });
-
- // create the actual velocity vector
- velocity *= this->curSol(freeFlowMomentumIndex)[
- couplingMapper_.outsideDofIndex(Dune::index_constant<porousMediumIndex>(), scvf)
- ];
-
- return velocity;
- }
-
-private:
- //! Return the volume variables of domain i for a given element and scv
- template<std::size_t i>
- VolumeVariables<i> volVars_(Dune::index_constant<i> domainI,
- const Element<i>& element,
- const SubControlVolume<i>& scv) const
- {
- VolumeVariables<i> volVars;
- const auto elemSol = elementSolution(
- element, this->curSol(domainI), this->problem(domainI).gridGeometry()
- );
- volVars.update(elemSol, this->problem(domainI), element, scv);
- return volVars;
- }
-
- /*!
- * \brief Returns whether a given scvf is coupled to the other domain
- */
- template<std::size_t i>
- bool isCoupledElement_(Dune::index_constant<i> domainI, std::size_t eIdx) const
- { return couplingMapper_.isCoupledElement(domainI, eIdx); }
-
- /*!
- * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
- */
- template<std::size_t i>
- void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const
- {
- const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element);
- bindCouplingContext_(domainI, fvGeometry);
- }
-
- /*!
- * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
- */
- template<std::size_t i>
- void bindCouplingContext_(Dune::index_constant<i> domainI, const FVElementGeometry<i>& fvGeometry) const
- {
- auto& context = std::get<domainI>(couplingContext_);
- context.clear();
-
- const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element());
-
- // do nothing if the element is not coupled to the other domain
- if (!isCoupledElement_(domainI, eIdx))
- return;
-
- couplingContextBoundForElement_[domainI] = eIdx;
-
- for (const auto& scvf : scvfs(fvGeometry))
- {
- if (couplingMapper_.isCoupled(domainI, scvf))
- {
- if constexpr (domainI == freeFlowMomentumIndex)
- {
- const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
- constexpr auto domainJ = Dune::index_constant<1-i>();
- const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
- const auto& otherElement = otherGridGeometry.element(otherElementIdx);
- auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
-
- // there is only one scv for TPFA
- context.push_back({
- otherElement,
- volVars_(domainJ, otherElement, *std::begin(scvs(otherFvGeometry))),
- std::move(otherFvGeometry),
- scvf.index(),
- couplingMapper_.flipScvfIndex(domainI, scvf),
- otherElementIdx,
- this->problem(domainJ).spatialParams().gravity(scvf.center())
- });
- }
-
- else if constexpr (domainI == porousMediumIndex)
- {
- const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
- constexpr auto domainJ = Dune::index_constant<1-i>();
- const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
- const auto& otherElement = otherGridGeometry.element(otherElementIdx);
- auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
-
- const auto otherScvfIdx = couplingMapper_.flipScvfIndex(domainI, scvf);
- const auto& otherScvf = otherFvGeometry.scvf(otherScvfIdx);
- const auto& otherScv = otherFvGeometry.scv(otherScvf.insideScvIdx());
-
- context.push_back({
- otherElement,
- std::move(otherFvGeometry),
- scvf.index(),
- otherScvfIdx,
- otherScv.dofIndex(),
- faceVelocity(fvGeometry.element(), scvf)
- });
- }
- }
- }
- }
+template<class MDTraits>
+struct FreeFlowMomentumPorousMediumCouplingManagerSelector<MDTraits, DiscretizationMethods::FCStaggered, DiscretizationMethods::CCTpfa>
+{ using type = FFMomentumPMCouplingManagerStaggeredCCTpfa<MDTraits>; };
- mutable std::tuple<std::vector<FreeFlowMomentumCouplingContext>, std::vector<PorousMediumCouplingContext>> couplingContext_;
- mutable std::array<std::size_t, 2> couplingContextBoundForElement_;
+} // end namespace FreeFlowMomentumPorousMediumDetail
+#endif // DOXYGEN
- CouplingMapper couplingMapper_;
-};
+template<class MDTraits>
+using FreeFlowMomentumPorousMediumCouplingManager = typename FreeFlowMomentumPorousMediumDetail::FreeFlowMomentumPorousMediumCouplingManagerSelector<MDTraits>::type;
} // end namespace Dumux
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager_staggered_cctpfa.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager_staggered_cctpfa.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6e0c3c0cbdc136e0d77377261dc2123d395e9f5a
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager_staggered_cctpfa.hh
@@ -0,0 +1,411 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \copydoc Dumux::FFMomentumPMCouplingManagerStaggeredCCTpfa
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMMOMENTUMPM_COUPLINGMANAGER_STAGGERED_TPFA_HH
+#define DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMMOMENTUMPM_COUPLINGMANAGER_STAGGERED_TPFA_HH
+
+#include <utility>
+#include <memory>
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/exceptions.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/staggered/elementsolution.hh>
+#include <dumux/multidomain/couplingmanager.hh>
+
+#include "couplingmapper_staggered_cctpfa.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Coupling manager for Stokes and Darcy domains with equal dimension
+ * specialization for staggered-cctpfa coupling.
+ */
+template<class MDTraits>
+class FFMomentumPMCouplingManagerStaggeredCCTpfa
+: public CouplingManager<MDTraits>
+{
+ using Scalar = typename MDTraits::Scalar;
+ using ParentType = CouplingManager<MDTraits>;
+
+public:
+ static constexpr auto freeFlowMomentumIndex = typename MDTraits::template SubDomain<0>::Index();
+ static constexpr auto porousMediumIndex = typename MDTraits::template SubDomain<1>::Index();
+
+ using SolutionVectorStorage = typename ParentType::SolutionVectorStorage;
+private:
+ // obtain the type tags of the sub problems
+ using FreeFlowMomentumTypeTag = typename MDTraits::template SubDomain<freeFlowMomentumIndex>::TypeTag;
+ using PorousMediumTypeTag = typename MDTraits::template SubDomain<porousMediumIndex>::TypeTag;
+
+ using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >;
+ using CouplingStencil = CouplingStencils::mapped_type;
+
+ // the sub domain type tags
+ template<std::size_t id>
+ using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+
+ template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView;
+ template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+ template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
+ template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>;
+ template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
+ template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
+ template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
+ template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
+ template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
+
+ using VelocityVector = typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate;
+
+ struct FreeFlowMomentumCouplingContext
+ {
+ Element<porousMediumIndex> element;
+ VolumeVariables<porousMediumIndex> volVars;
+ FVElementGeometry<porousMediumIndex> fvGeometry;
+ std::size_t freeFlowMomentumScvfIdx;
+ std::size_t porousMediumScvfIdx;
+ std::size_t porousMediumDofIdx;
+ VelocityVector gravity;
+ };
+
+ struct PorousMediumCouplingContext
+ {
+ Element<freeFlowMomentumIndex> element;
+ FVElementGeometry<freeFlowMomentumIndex> fvGeometry;
+ std::size_t porousMediumScvfIdx;
+ std::size_t freeFlowMomentumScvfIdx;
+ std::size_t freeFlowMomentumDofIdx;
+ VelocityVector faceVelocity;
+ };
+
+ using CouplingMapper = FFMomentumPMCouplingMapperStaggeredCCTpfa<MDTraits, FFMomentumPMCouplingManagerStaggeredCCTpfa<MDTraits>>;
+
+public:
+
+ /*!
+ * \brief Methods to be accessed by main
+ */
+ // \{
+
+ //! Initialize the coupling manager
+ void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
+ std::shared_ptr<Problem<porousMediumIndex>> porousMediumProblem,
+ SolutionVectorStorage& curSol)
+ {
+ this->setSubProblems(std::make_tuple(freeFlowMomentumProblem, porousMediumProblem));
+ this->attachSolution(curSol);
+
+ couplingMapper_.update(*this);
+ }
+
+ // \}
+
+
+ /*!
+ * \brief Methods to be accessed by the assembly
+ */
+ // \{
+
+ /*!
+ * \brief prepares all data and variables that are necessary to evaluate the residual (called from the local assembler)
+ */
+ template<std::size_t i, class Assembler>
+ void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler) const
+ {
+ bindCouplingContext_(domainI, element);
+ }
+
+ /*!
+ * \brief Update the coupling context
+ */
+ template<std::size_t i, std::size_t j, class LocalAssemblerI>
+ void updateCouplingContext(Dune::index_constant<i> domainI,
+ const LocalAssemblerI& localAssemblerI,
+ Dune::index_constant<j> domainJ,
+ std::size_t dofIdxGlobalJ,
+ const PrimaryVariables<j>& priVarsJ,
+ int pvIdxJ)
+ {
+ this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
+
+ // we need to update all solution-depenent components of the coupling context
+ // the dof of domain J has been deflected
+
+ // update the faceVelocity in the PorousMediumCouplingContext
+ if constexpr (domainJ == freeFlowMomentumIndex)
+ {
+ // we only need to update if we are assembling the porous medium domain
+ // since the freeflow domain will not use the velocity from the context
+ if constexpr (domainI == porousMediumIndex)
+ {
+ auto& context = std::get<porousMediumIndex>(couplingContext_);
+ for (auto& c : context)
+ {
+ if (c.freeFlowMomentumDofIdx == dofIdxGlobalJ)
+ {
+ const auto& scvf = localAssemblerI.fvGeometry().scvf(c.porousMediumScvfIdx);
+ c.faceVelocity = faceVelocity(localAssemblerI.element(), scvf);
+ }
+ }
+ }
+ }
+
+ // update volVars in the FreeFlowMomentumCouplingContext
+ else if (domainJ == porousMediumIndex)
+ {
+ auto& context = std::get<freeFlowMomentumIndex>(couplingContext_);
+ for (auto& c : context)
+ {
+ if (c.porousMediumDofIdx == dofIdxGlobalJ)
+ {
+ const auto& ggJ = c.fvGeometry.gridGeometry();
+ const auto& scv = *scvs(c.fvGeometry).begin();
+ const auto elemSol = elementSolution(c.element, this->curSol(domainJ), ggJ);
+ c.volVars.update(elemSol, this->problem(domainJ), c.element, scv);
+ }
+ }
+ }
+ }
+
+ // \}
+
+ /*!
+ * \brief Access the coupling context needed for the Stokes domain
+ */
+ template<std::size_t i>
+ const auto& couplingContext(Dune::index_constant<i> domainI,
+ const FVElementGeometry<i>& fvGeometry,
+ const SubControlVolumeFace<i> scvf) const
+ {
+ auto& contexts = std::get<i>(couplingContext_);
+
+ if (contexts.empty() || couplingContextBoundForElement_[i] != scvf.insideScvIdx())
+ bindCouplingContext_(domainI, fvGeometry);
+
+ for (const auto& context : contexts)
+ {
+ const auto expectedScvfIdx = domainI == freeFlowMomentumIndex ? context.freeFlowMomentumScvfIdx : context.porousMediumScvfIdx;
+ if (scvf.index() == expectedScvfIdx)
+ return context;
+ }
+
+ DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
+ }
+
+ /*!
+ * \brief The coupling stencils
+ */
+ // \{
+
+ /*!
+ * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs
+ */
+ const CouplingStencil& couplingStencil(Dune::index_constant<porousMediumIndex> domainI,
+ const Element<porousMediumIndex>& element,
+ Dune::index_constant<freeFlowMomentumIndex> domainJ) const
+ {
+ const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
+ return couplingMapper_.couplingStencil(domainI, eIdx, domainJ);
+ }
+
+ /*!
+ * \brief returns an iterable container of all indices of degrees of freedom of domain j
+ * that couple with / influence the residual of the given sub-control volume of domain i
+ *
+ * \param domainI the domain index of domain i
+ * \param elementI the coupled element of domain Ã
+ * \param scvI the sub-control volume of domain i
+ * \param domainJ the domain index of domain j
+ */
+ const CouplingStencil& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI,
+ const Element<freeFlowMomentumIndex>& elementI,
+ const SubControlVolume<freeFlowMomentumIndex>& scvI,
+ Dune::index_constant<porousMediumIndex> domainJ) const
+ {
+ return couplingMapper_.couplingStencil(domainI, elementI, scvI, domainJ);
+ }
+
+ using ParentType::evalCouplingResidual;
+
+ /*!
+ * \brief evaluate the coupling residual
+ * special interface for fcstaggered methods
+ */
+ template<class LocalAssemblerI, std::size_t j>
+ decltype(auto) evalCouplingResidual(Dune::index_constant<freeFlowMomentumIndex> domainI,
+ const LocalAssemblerI& localAssemblerI,
+ const SubControlVolume<freeFlowMomentumIndex>& scvI,
+ Dune::index_constant<j> domainJ,
+ std::size_t dofIdxGlobalJ) const
+ {
+ return localAssemblerI.evalLocalResidual();
+ }
+
+ // \}
+
+ /*!
+ * \brief Returns whether a given scvf is coupled to the other domain
+ */
+ bool isCoupledLateralScvf(Dune::index_constant<freeFlowMomentumIndex> domainI, const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
+ { return couplingMapper_.isCoupledLateralScvf(domainI, scvf); }
+
+ /*!
+ * \brief Returns whether a given scvf is coupled to the other domain
+ */
+ template<std::size_t i>
+ bool isCoupled(Dune::index_constant<i> domainI, const SubControlVolumeFace<i>& scvf) const
+ {
+ if constexpr (i == freeFlowMomentumIndex)
+ return couplingMapper_.isCoupled(domainI, scvf) || couplingMapper_.isCoupledLateralScvf(domainI, scvf);
+ else
+ return couplingMapper_.isCoupled(domainI, scvf);
+ }
+
+ /*!
+ * \brief If the boundary entity is on a coupling boundary
+ * \param domainI the domain index for which to compute the flux
+ * \param scv the sub control volume
+ */
+ bool isCoupled(Dune::index_constant<freeFlowMomentumIndex> domainI,
+ const SubControlVolume<freeFlowMomentumIndex>& scv) const
+ { return couplingMapper_.isCoupled(domainI, scv); }
+
+ /*!
+ * \brief Returns the velocity at a given sub control volume face.
+ */
+ auto faceVelocity(const Element<porousMediumIndex>& element,
+ const SubControlVolumeFace<porousMediumIndex>& scvf) const
+ {
+ // create a unit normal vector oriented in positive coordinate direction
+ auto velocity = scvf.unitOuterNormal();
+ using std::abs;
+ std::for_each(velocity.begin(), velocity.end(), [](auto& v){ v = abs(v); });
+
+ // create the actual velocity vector
+ velocity *= this->curSol(freeFlowMomentumIndex)[
+ couplingMapper_.outsideDofIndex(Dune::index_constant<porousMediumIndex>(), scvf)
+ ];
+
+ return velocity;
+ }
+
+private:
+ //! Return the volume variables of domain i for a given element and scv
+ template<std::size_t i>
+ VolumeVariables<i> volVars_(Dune::index_constant<i> domainI,
+ const Element<i>& element,
+ const SubControlVolume<i>& scv) const
+ {
+ VolumeVariables<i> volVars;
+ const auto elemSol = elementSolution(
+ element, this->curSol(domainI), this->problem(domainI).gridGeometry()
+ );
+ volVars.update(elemSol, this->problem(domainI), element, scv);
+ return volVars;
+ }
+
+ /*!
+ * \brief Returns whether a given scvf is coupled to the other domain
+ */
+ template<std::size_t i>
+ bool isCoupledElement_(Dune::index_constant<i> domainI, std::size_t eIdx) const
+ { return couplingMapper_.isCoupledElement(domainI, eIdx); }
+
+ /*!
+ * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
+ */
+ template<std::size_t i>
+ void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const
+ {
+ const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element);
+ bindCouplingContext_(domainI, fvGeometry);
+ }
+
+ /*!
+ * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information)
+ */
+ template<std::size_t i>
+ void bindCouplingContext_(Dune::index_constant<i> domainI, const FVElementGeometry<i>& fvGeometry) const
+ {
+ auto& context = std::get<domainI>(couplingContext_);
+ context.clear();
+
+ const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element());
+
+ // do nothing if the element is not coupled to the other domain
+ if (!isCoupledElement_(domainI, eIdx))
+ return;
+
+ couplingContextBoundForElement_[domainI] = eIdx;
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ if (couplingMapper_.isCoupled(domainI, scvf))
+ {
+ if constexpr (domainI == freeFlowMomentumIndex)
+ {
+ const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
+ constexpr auto domainJ = Dune::index_constant<1-i>();
+ const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
+ const auto& otherElement = otherGridGeometry.element(otherElementIdx);
+ auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
+
+ // there is only one scv for TPFA
+ context.push_back({
+ otherElement,
+ volVars_(domainJ, otherElement, *std::begin(scvs(otherFvGeometry))),
+ std::move(otherFvGeometry),
+ scvf.index(),
+ couplingMapper_.flipScvfIndex(domainI, scvf),
+ otherElementIdx,
+ this->problem(domainJ).spatialParams().gravity(scvf.center())
+ });
+ }
+
+ else if constexpr (domainI == porousMediumIndex)
+ {
+ const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf);
+ constexpr auto domainJ = Dune::index_constant<1-i>();
+ const auto& otherGridGeometry = this->problem(domainJ).gridGeometry();
+ const auto& otherElement = otherGridGeometry.element(otherElementIdx);
+ auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement);
+
+ const auto otherScvfIdx = couplingMapper_.flipScvfIndex(domainI, scvf);
+ const auto& otherScvf = otherFvGeometry.scvf(otherScvfIdx);
+ const auto& otherScv = otherFvGeometry.scv(otherScvf.insideScvIdx());
+
+ context.push_back({
+ otherElement,
+ std::move(otherFvGeometry),
+ scvf.index(),
+ otherScvfIdx,
+ otherScv.dofIndex(),
+ faceVelocity(fvGeometry.element(), scvf)
+ });
+ }
+ }
+ }
+ }
+
+ mutable std::tuple<std::vector<FreeFlowMomentumCouplingContext>, std::vector<PorousMediumCouplingContext>> couplingContext_;
+ mutable std::array<std::size_t, 2> couplingContextBoundForElement_;
+
+ CouplingMapper couplingMapper_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper_staggered_cctpfa.hh
similarity index 98%
rename from dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh
rename to dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper_staggered_cctpfa.hh
index 6e02da1a72cc6ca5c1369166fb391284c085230a..56ecf1a5cec28b86d0822a2306af6d1dbc2e002c 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper_staggered_cctpfa.hh
@@ -7,11 +7,11 @@
/*!
* \file
* \ingroup FreeFlowPorousMediumCoupling
- * \copydoc Dumux::FreeFlowMomentumPorousMediumCouplingMapper
+ * \copydoc Dumux::FFMomentumPMCouplingMapperStaggeredCCTpfa
*/
-#ifndef DUMUX_MULTIDOMAIN_FREEFLOWMOMENTUM_POROUSMEDIUM_COUPLINGMAPPER_HH
-#define DUMUX_MULTIDOMAIN_FREEFLOWMOMENTUM_POROUSMEDIUM_COUPLINGMAPPER_HH
+#ifndef DUMUX_MULTIDOMAIN_FREEFLOWMOMENTUM_POROUSMEDIUM_COUPLINGMAPPER_STAGGERED_TPFA_HH
+#define DUMUX_MULTIDOMAIN_FREEFLOWMOMENTUM_POROUSMEDIUM_COUPLINGMAPPER_STAGGERED_TPFA_HH
#include <iostream>
#include <unordered_map>
@@ -33,7 +33,7 @@ namespace Dumux {
* \todo how to extend to arbitrary number of domains?
*/
template<class MDTraits, class CouplingManager>
-class FreeFlowMomentumPorousMediumCouplingMapper
+class FFMomentumPMCouplingMapperStaggeredCCTpfa
{
using Scalar = typename MDTraits::Scalar;
@@ -80,7 +80,6 @@ public:
*/
void update(const CouplingManager& couplingManager)
{
- // TODO: Box and multiple domains
static_assert(numSD == 2, "More than two subdomains not implemented!");
static_assert(isFcStaggered<0>() && isCCTpfa<1>(), "Only coupling between fcstaggered and cctpfa implemented!");