diff --git a/dumux/discretization/facecentered/staggered/fvgridgeometry.hh b/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
index c53cd6e1bad9de339b46b21f81cd461a707fd457..18af7159e96a1d63d70189f4e7ae20ffaff2c0fb 100644
--- a/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/facecentered/staggered/fvgridgeometry.hh
@@ -822,6 +822,9 @@ private:
}
}
+ // set number of subcontrolvolume faces
+ numScvf_ = scvfIdx;
+
connectivityMap_.update(*this);
}
diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh
index 092603f0324a6b1c326f92ede9dd6bb86291e7b8..d6832349bd3de307d22f03878c7a455517d8c5d8 100644
--- a/dumux/freeflow/navierstokes/problem.hh
+++ b/dumux/freeflow/navierstokes/problem.hh
@@ -471,7 +471,7 @@ public:
const VelocityVector beaversJosephVelocity(const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
- const Scalar tangentialVelocityGradient) const
+ Scalar tangentialVelocityGradient) const
{
assert(scvf.isLateral());
assert(scvf.boundary());
@@ -487,6 +487,10 @@ public:
const Scalar betaBJ = asImp_().betaBJ(fvGeometry, scvf, orientation);
const Scalar distanceNormalToBoundary = (scvf.ipGlobal() - scv.dofPosition()).two_norm();
+ static const bool onlyNormalGradient = getParamFromGroup<bool>(this->paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
+ if (onlyNormalGradient)
+ tangentialVelocityGradient = 0.0;
+
const Scalar scalarSlipVelocity = (tangentialVelocityGradient*distanceNormalToBoundary
+ asImp_().porousMediumVelocity(fvGeometry, scvf) * orientation * betaBJ * distanceNormalToBoundary
+ elemVolVars[scv].velocity()) / (betaBJ*distanceNormalToBoundary + 1.0);
diff --git a/dumux/multidomain/boundary/CMakeLists.txt b/dumux/multidomain/boundary/CMakeLists.txt
index addd460feb4d2b502e6a92a0630ad445b11020d9..5ffc5873b26785265ec0c9ad64e4cc7aa93cce1b 100644
--- a/dumux/multidomain/boundary/CMakeLists.txt
+++ b/dumux/multidomain/boundary/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory(darcydarcy)
+add_subdirectory(freeflowporousmedium)
add_subdirectory(stokesdarcy)
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/CMakeLists.txt b/dumux/multidomain/boundary/freeflowporousmedium/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..063a18b3b6887225bcbd24d05c411231071873f7
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/CMakeLists.txt
@@ -0,0 +1,6 @@
+add_subdirectory(ffmasspm)
+add_subdirectory(ffmomentumpm)
+
+file(GLOB DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOWPOROUSMEDIUM_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOWPOROUSMEDIUM_HEADERS}
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/multidomain/boundary/freeflowporousmedium)
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d4d9dceec66947ac3b9967ff2c4f1e60f671e476
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingconditions.hh
@@ -0,0 +1,580 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \copydoc Dumux::FreeFlowPorousMediumCouplingData
+ */
+
+#ifndef DUMUX_MD_FREEFLOW_POROUSMEDIUM_COUPLINGCONDITIONS_HH
+#define DUMUX_MD_FREEFLOW_POROUSMEDIUM_COUPLINGCONDITIONS_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>
+
+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, arithmethic, ffOnly, pmOnly
+ };
+
+ /*!
+ * \brief Convenience function to convert user input given as std::string to the corresponding enum class used for chosing 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 == "Arithmethic")
+ return DiffusionCoefficientAveragingType::arithmethic;
+ else if (diffusionCoefficientAveragingType == "FreeFlowOnly")
+ return DiffusionCoefficientAveragingType::ffOnly;
+ else if (diffusionCoefficientAveragingType == "PorousMediumOnly")
+ return DiffusionCoefficientAveragingType::pmOnly;
+ else
+ DUNE_THROW(Dune::IOError, "Unknown DiffusionCoefficientAveragingType");
+ }
+
+};
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief This structs helps to check if the two sub models use the same fluidsystem.
+ * Specialization for the case of using an adapter only for the free-flow model.
+ * \tparam FFFS The free-flow fluidsystem
+ * \tparam PMFS The porous-medium flow fluidsystem
+ */
+template<class FFFS, class PMFS>
+struct IsSameFluidSystem
+{
+ static_assert(FFFS::numPhases == 1, "Only single-phase fluidsystems may be used for free flow.");
+ static constexpr bool value = std::is_same<typename FFFS::MultiPhaseFluidSystem, PMFS>::value;
+};
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief This structs helps to check if the two sub models use the same fluidsystem.
+ * \tparam FS The fluidsystem
+ */
+template<class FS>
+struct IsSameFluidSystem<FS, FS>
+{
+ static_assert(FS::numPhases == 1, "Only single-phase fluidsystems may be used for free flow.");
+ static constexpr bool value = std::is_same<FS, FS>::value; // always true
+};
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief This structs indicates that Fick's law is not used for diffusion.
+ * \tparam DiffLaw The diffusion law.
+ */
+template<class DiffLaw>
+struct IsFicksLaw : public std::false_type {};
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief This structs indicates that Fick's law is used for diffusion.
+ * \tparam DiffLaw The diffusion law.
+ */
+template<class T>
+struct IsFicksLaw<Dumux::FicksLaw<T>> : public std::true_type {};
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Helper struct to choose the correct index for phases and components. This is need if the porous-medium-flow model
+ features more fluid phases than the free-flow model.
+ * \tparam stokesIdx The domain index of the free-flow model.
+ * \tparam porousMediumIndex The domain index of the porous-medium-flow model.
+ * \tparam FFFS The free-flow fluidsystem.
+ * \tparam hasAdapter Specifies whether an adapter class for the fluidsystem is used.
+ */
+template<std::size_t stokesIdx, std::size_t porousMediumIndex, class FFFS, bool hasAdapter>
+struct IndexHelper;
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Helper struct to choose the correct index for phases and components. This is need if the porous-medium-flow model
+ features more fluid phases than the free-flow model. Specialization for the case that no adapter is used.
+ * \tparam stokesIdx The domain index of the free-flow model.
+ * \tparam porousMediumIndex The domain index of the porous-medium-flow model.
+ * \tparam FFFS The free-flow fluidsystem.
+ */
+template<std::size_t stokesIdx, std::size_t porousMediumIndex, class FFFS>
+struct IndexHelper<stokesIdx, porousMediumIndex, FFFS, false>
+{
+ /*!
+ * \brief No adapter is used, just return the input index.
+ */
+ template<std::size_t i>
+ static constexpr auto couplingPhaseIdx(Dune::index_constant<i>, int coupledPhaseIdx = 0)
+ { return coupledPhaseIdx; }
+
+ /*!
+ * \brief No adapter is used, just return the input index.
+ */
+ template<std::size_t i>
+ static constexpr auto couplingCompIdx(Dune::index_constant<i>, int coupledCompdIdx)
+ { return coupledCompdIdx; }
+};
+
+/*!
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \brief Helper struct to choose the correct index for phases and components. This is need if the porous-medium-flow model
+ features more fluid phases than the free-flow model. Specialization for the case that a adapter is used.
+ * \tparam stokesIdx The domain index of the free-flow model.
+ * \tparam porousMediumIndex The domain index of the porous-medium-flow model.
+ * \tparam FFFS The free-flow fluidsystem.
+ */
+template<std::size_t stokesIdx, std::size_t porousMediumIndex, class FFFS>
+struct IndexHelper<stokesIdx, porousMediumIndex, FFFS, true>
+{
+ /*!
+ * \brief The free-flow model always uses phase index 0.
+ */
+ static constexpr int couplingPhaseIdx(Dune::index_constant<stokesIdx>, int coupledPhaseIdx = 0)
+ { return 0; }
+
+ /*!
+ * \brief The phase index of the porous-medium-flow model is given by the adapter fluidsytem (i.e., user input).
+ */
+ static constexpr auto couplingPhaseIdx(Dune::index_constant<porousMediumIndex>, int coupledPhaseIdx = 0)
+ { return FFFS::multiphaseFluidsystemPhaseIdx; }
+
+ /*!
+ * \brief The free-flow model does not need any change of the component index.
+ */
+ static constexpr auto couplingCompIdx(Dune::index_constant<stokesIdx>, int coupledCompdIdx)
+ { return coupledCompdIdx; }
+
+ /*!
+ * \brief The component index of the porous-medium-flow model is mapped by the adapter fluidsytem.
+ */
+ static constexpr auto couplingCompIdx(Dune::index_constant<porousMediumIndex>, int coupledCompdIdx)
+ { return FFFS::compIdx(coupledCompdIdx); }
+};
+
+template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional>
+class FreeFlowPorousMediumCouplingConditionsImplementation;
+
+/*!
+* \ingroup BoundaryCoupling
+* \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 = Dumux::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(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 coupledCompdIdx)
+ { return IndexHelper::couplingCompIdx(id, coupledCompdIdx); }
+
+ /*!
+ * \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::arithmethic)
+ {
+ 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 acorss 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(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;
+ }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cb3b2b57ade83aa1c4952885cf3f9b9d1ad79bec
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh
@@ -0,0 +1,477 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryCoupling
+ * \copydoc Dumux::FreeFlowPorousMediumCouplingManager
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_POROUSMEDIUM_COUPLINGMANAGER_HH
+#define DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_POROUSMEDIUM_COUPLINGMANAGER_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/staggeredfreeflow/couplingmanager.hh>
+#include <dumux/multidomain/multibinarycouplingmanager.hh>
+
+#include "couplingconditions.hh"
+
+namespace Dumux {
+
+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::StaggeredFreeFlowCouplingManager<FreeFlowTraits>;
+ using FreeFlowMomentumPorousMediumCouplingManager
+ = Dumux::FreeFlowMomentumPorousMediumCouplingManager<FreeFlowMomentumPorousMediumTraits>;
+ using FreeFlowMassPorousMediumCouplingManager
+ = Dumux::FreeFlowMassPorousMediumCouplingManager<FreeFlowMassPorousMediumTraits>;
+};
+
+} // end namespace Detail
+
+
+/*!
+ * \ingroup BoundaryCoupling
+ * \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 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 iteratable 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/ffmasspm/CMakeLists.txt b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2c7f3e04baf43525e71f8e78ad20d4eb23e64c53
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOWPOROUSMEDIUM_FFMASSPM_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOWPOROUSMEDIUM_FFMASSPM_HEADERS}
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm)
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e8156b97747aff358f6f33891fcb029db0f2380d
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmasspm/couplingmanager.hh
@@ -0,0 +1,304 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \copydoc Dumux::FreeFlowMassPorousMediumCouplingManager
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMASSPM_COUPLINGMANAGER_HH
+#define DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMASSPM_COUPLINGMANAGER_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 StokesDarcyCoupling
+ * \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 availabe 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 availabe 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/CMakeLists.txt b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b79ca4a8b025f3fdad01e0495885236acbc424a0
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOWPOROUSMEDIUM_FFMOMENTUMPM_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOWPOROUSMEDIUM_FFMOMENTUMPM_HEADERS}
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm)
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f957e9f65f15d9c92dac5e867681b8f63f7d14f3
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmanager.hh
@@ -0,0 +1,422 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup FreeFlowPorousMediumCoupling
+ * \copydoc Dumux::FreeFlowMomentumPorousMediumCouplingManager
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMMOMENTUMPM_COUPLINGMANAGER_HH
+#define DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMMOMENTUMPM_COUPLINGMANAGER_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.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup StokesDarcyCoupling
+ * \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_);
+
+ 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 iteratable 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.hh
new file mode 100644
index 0000000000000000000000000000000000000000..59f15ea1385dbc547c749a1c8bb77f76b9d6c162
--- /dev/null
+++ b/dumux/multidomain/boundary/freeflowporousmedium/ffmomentumpm/couplingmapper.hh
@@ -0,0 +1,347 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup DarcyDarcyCoupling
+ * \copydoc Dumux::FreeFlowMomentumPorousMediumCouplingMapper
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_FREEFLOWMOMENTUM_POROUSMEDIUM_COUPLINGMAPPER_HH
+#define DUMUX_MULTIDOMAIN_FREEFLOWMOMENTUM_POROUSMEDIUM_COUPLINGMAPPER_HH
+
+#include <iostream>
+#include <unordered_map>
+#include <tuple>
+#include <vector>
+
+#include <dune/common/timer.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/indices.hh>
+
+#include <dumux/geometry/intersectingentities.hh>
+#include <dumux/discretization/method.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup DarcyDarcyCoupling
+ * \brief the default mapper for conforming equal dimension boundary coupling between two domains (box or cc)
+ * \todo how to extend to arbitrary number of domains?
+ */
+template<class MDTraits, class CouplingManager>
+class FreeFlowMomentumPorousMediumCouplingMapper
+{
+ using Scalar = typename MDTraits::Scalar;
+
+ template<std::size_t i> using GridGeometry = typename MDTraits::template SubDomain<i>::GridGeometry;
+ template<std::size_t i> using SubControlVolume = typename GridGeometry<i>::SubControlVolume;
+ template<std::size_t i> using SubControlVolumeFace = typename GridGeometry<i>::SubControlVolumeFace;
+ template<std::size_t i> using GridView = typename GridGeometry<i>::GridView;
+ template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity;
+
+ template<std::size_t i>
+ static constexpr auto domainIdx()
+ { return typename MDTraits::template SubDomain<i>::Index{}; }
+
+ template<std::size_t i>
+ static constexpr bool isFcStaggered()
+ { return GridGeometry<i>::discMethod == DiscretizationMethod::fcstaggered; }
+
+ template<std::size_t i>
+ static constexpr bool isCCTpfa()
+ { return GridGeometry<i>::discMethod == DiscretizationMethod::cctpfa; }
+
+ struct ScvfInfoPM
+ {
+ std::size_t eIdxOutside;
+ std::size_t flipScvfIdx;
+ };
+
+ struct ScvfInfoFF
+ {
+ std::size_t eIdxOutside;
+ std::size_t flipScvfIdx;
+ std::size_t dofIdxOutside;
+ };
+
+ using FlipScvfMapTypePM = std::unordered_map<std::size_t, ScvfInfoPM>;
+ using FlipScvfMapTypeFF = std::unordered_map<std::size_t, ScvfInfoFF>;
+ using MapType = std::unordered_map<std::size_t, std::vector<std::size_t>>;
+
+ static constexpr std::size_t numSD = MDTraits::numSubDomains;
+
+public:
+ /*!
+ * \brief Main update routine
+ */
+ 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!");
+
+ Dune::Timer watch;
+ std::cout << "Initializing the coupling map..." << std::endl;
+
+ for (std::size_t domIdx = 0; domIdx < numSD; ++domIdx)
+ stencils_[domIdx].clear();
+
+ std::get<CouplingManager::freeFlowMomentumIndex>(scvfInfo_).clear();
+ std::get<CouplingManager::porousMediumIndex>(scvfInfo_).clear();
+
+ const auto& freeFlowMomentumProblem = couplingManager.problem(CouplingManager::freeFlowMomentumIndex);
+ const auto& porousMediumProblem = couplingManager.problem(CouplingManager::porousMediumIndex);
+ const auto& freeFlowMomentumGG = freeFlowMomentumProblem.gridGeometry();
+ const auto& porousMediumGG = porousMediumProblem.gridGeometry();
+
+ isCoupledFFDof_.resize(freeFlowMomentumGG.numScvf(), false);
+ isCoupledFFElement_.resize(freeFlowMomentumGG.gridView().size(0), false);
+ isCoupledScvf_[CouplingManager::freeFlowMomentumIndex].resize(freeFlowMomentumGG.numScvf(), false);
+ isCoupledScvf_[CouplingManager::porousMediumIndex].resize(porousMediumGG.numScvf(), false);
+
+ auto pmFvGeometry = localView(porousMediumGG);
+ auto ffFvGeometry = localView(freeFlowMomentumGG);
+
+ for (const auto& pmElement : elements(porousMediumGG.gridView()))
+ {
+ pmFvGeometry.bindElement(pmElement);
+
+ for (const auto& pmScvf : scvfs(pmFvGeometry))
+ {
+ // skip all non-boundaries
+ if (!pmScvf.boundary())
+ continue;
+
+ // get elements intersecting with the scvf center
+ // for robustness add epsilon in unit outer normal direction
+ const auto eps = (pmScvf.ipGlobal() - pmScvf.geometry().corner(0)).two_norm()*1e-8;
+ auto globalPos = pmScvf.ipGlobal(); globalPos.axpy(eps, pmScvf.unitOuterNormal());
+ const auto indices = intersectingEntities(globalPos, freeFlowMomentumGG.boundingBoxTree());
+
+ // skip if no intersection was found
+ if (indices.empty())
+ continue;
+
+ // sanity check
+ if (indices.size() > 1)
+ DUNE_THROW(Dune::InvalidStateException, "Are you sure your sub-domain grids are conformingly discretized on the common interface?");
+
+ // add the pair to the multimap
+ const auto pmElemIdx = porousMediumGG.elementMapper().index(pmElement);
+ const auto ffElemIdx = indices[0];
+ const auto& ffElement = freeFlowMomentumGG.element(ffElemIdx);
+ ffFvGeometry.bindElement(ffElement);
+
+ for (const auto& ffScvf : scvfs(ffFvGeometry)) // TODO: maybe better iterate over all scvs
+ {
+ // TODO this only takes the scv directly at the interface, maybe extend
+ if (!ffScvf.boundary() || !ffScvf.isFrontal())
+ continue;
+
+ const auto dist = (pmScvf.ipGlobal() - ffScvf.ipGlobal()).two_norm();
+ if (dist > eps)
+ continue;
+
+ const auto& ffScv = ffFvGeometry.scv(ffScvf.insideScvIdx());
+ stencils_[CouplingManager::porousMediumIndex][pmElemIdx].push_back(ffScv.dofIndex());
+ stencils_[CouplingManager::freeFlowMomentumIndex][ffScv.dofIndex()].push_back(pmElemIdx);
+
+ // mark the scvf and find and mark the flip scvf
+ isCoupledScvf_[CouplingManager::porousMediumIndex][pmScvf.index()] = true;
+ isCoupledScvf_[CouplingManager::freeFlowMomentumIndex][ffScvf.index()] = true;
+
+ // add all lateral free-flow scvfs touching the coupling interface ("standing" or "lying")
+ for (const auto& otherFfScvf : scvfs(ffFvGeometry, ffScv))
+ {
+ if (otherFfScvf.isLateral())
+ {
+ // the orthogonal scvf has to lie on the coupling interface
+ const auto& lateralOrthogonalScvf = ffFvGeometry.lateralOrthogonalScvf(otherFfScvf);
+ isCoupledLateralScvf_[lateralOrthogonalScvf.index()] = true;
+
+ // the "standing" scvf is marked as coupled if it does not lie on a boundary or if that boundary itself is a coupling interface
+ if (!otherFfScvf.boundary())
+ isCoupledLateralScvf_[otherFfScvf.index()] = true;
+ else
+ {
+ // for robustness add epsilon in unit outer normal direction
+ const auto otherScvfEps = (otherFfScvf.ipGlobal() - otherFfScvf.geometry().corner(0)).two_norm()*1e-8;
+ auto otherScvfGlobalPos = otherFfScvf.center(); otherScvfGlobalPos.axpy(otherScvfEps, otherFfScvf.unitOuterNormal());
+
+ if (!intersectingEntities(otherScvfGlobalPos, porousMediumGG.boundingBoxTree()).empty())
+ isCoupledLateralScvf_[otherFfScvf.index()] = true;
+ }
+ }
+ }
+ isCoupledFFDof_[ffScv.dofIndex()] = true;
+ isCoupledFFElement_[ffElemIdx] = true;
+
+ std::get<CouplingManager::porousMediumIndex>(scvfInfo_)[pmScvf.index()] = ScvfInfoPM{ffElemIdx, ffScvf.index()};
+ std::get<CouplingManager::freeFlowMomentumIndex>(scvfInfo_)[ffScvf.index()] = ScvfInfoFF{pmElemIdx, pmScvf.index(), ffScv.dofIndex()};
+ }
+ }
+ }
+
+ std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
+ }
+
+ /*!
+ * \brief returns an iteratable container of all indices of degrees of freedom of domain j
+ * that couple with / influence the element residual of the given element of domain i
+ *
+ * \param domainI the domain index of domain i
+ * \param eIdxI the index of the coupled element of domain Ã
+ * \param domainJ the domain index of domain j
+ *
+ * \note The element residual definition depends on the discretization scheme of domain i
+ * box: a container of the residuals of all sub control volumes
+ * cc : the residual of the (sub) control volume
+ * fem: the residual of the element
+ * \note This function has to be implemented by all coupling managers for all combinations of i and j
+ */
+ const std::vector<std::size_t>& couplingStencil(Dune::index_constant<CouplingManager::porousMediumIndex> domainI,
+ const std::size_t eIdxI,
+ Dune::index_constant<CouplingManager::freeFlowMomentumIndex> domainJ) const
+ {
+ if (isCoupledElement(domainI, eIdxI))
+ return stencils_[CouplingManager::porousMediumIndex].at(eIdxI);
+ else
+ return emptyStencil_;
+ }
+
+ /*!
+ * \brief returns an iteratable container of all indices of degrees of freedom of domain j
+ * that couple with / influence the element residual of the given element of domain i
+ *
+ * \param domainI the domain index of domain i
+ * \param eIdxI the index of the coupled element of domain Ã
+ * \param domainJ the domain index of domain j
+ *
+ * \note The element residual definition depends on the discretization scheme of domain i
+ * box: a container of the residuals of all sub control volumes
+ * cc : the residual of the (sub) control volume
+ * fem: the residual of the element
+ * \note This function has to be implemented by all coupling managers for all combinations of i and j
+ */
+ const std::vector<std::size_t>& couplingStencil(Dune::index_constant<CouplingManager::freeFlowMomentumIndex> domainI,
+ const Element<CouplingManager::freeFlowMomentumIndex>& elementI,
+ const SubControlVolume<CouplingManager::freeFlowMomentumIndex>& scvI,
+ Dune::index_constant<CouplingManager::porousMediumIndex> domainJ) const
+ {
+ if (isCoupled(domainI, scvI))
+ return stencils_[CouplingManager::freeFlowMomentumIndex].at(scvI.dofIndex());
+ else
+ return emptyStencil_;
+ }
+
+ /*!
+ * \brief Return if an element residual with index eIdx of domain i is coupled to domain j
+ */
+ template<std::size_t i>
+ bool isCoupledElement(Dune::index_constant<i>, std::size_t eIdx) const
+ {
+ if constexpr (i == CouplingManager::porousMediumIndex)
+ return static_cast<bool>(stencils_[i].count(eIdx));
+ else
+ return isCoupledFFElement_[eIdx];
+ }
+
+ /*!
+ * \brief If the boundary entity is on a coupling boundary
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ template<std::size_t i>
+ bool isCoupled(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ return isCoupledScvf_[i].at(scvf.index());
+ }
+
+ /*!
+ * \brief If the boundary entity is on a coupling boundary
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ bool isCoupledLateralScvf(Dune::index_constant<CouplingManager::freeFlowMomentumIndex> domainI,
+ const SubControlVolumeFace<CouplingManager::freeFlowMomentumIndex>& scvf) const
+ { return isCoupledLateralScvf_.count(scvf.index()); }
+
+ /*!
+ * \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<CouplingManager::freeFlowMomentumIndex> domainI,
+ const SubControlVolume<CouplingManager::freeFlowMomentumIndex>& scv) const
+ { return isCoupledFFDof_[scv.dofIndex()]; }
+
+ /*!
+ * \brief Return the scvf index of the flipped scvf in the other domain
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ template<std::size_t i>
+ std::size_t flipScvfIndex(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ return std::get<i>(scvfInfo_).at(scvf.index()).flipScvfIdx;
+ }
+
+ /*!
+ * \brief Return the outside element index (the element index of the other domain)
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ template<std::size_t i>
+ std::size_t outsideElementIndex(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ return std::get<i>(scvfInfo_).at(scvf.index()).eIdxOutside;
+ }
+
+ /*!
+ * \brief Return the outside element index (the element index of the other domain)
+ * \param domainI the domain index for which to compute the flux
+ * \param scvf the sub control volume face
+ */
+ template<std::size_t i>
+ std::size_t outsideDofIndex(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvf) const
+ {
+ if constexpr (i == CouplingManager::porousMediumIndex)
+ return outsideElementIndex(domainI, scvf);
+ else
+ return std::get<i>(scvfInfo_).at(scvf.index()).dofIdxOutside;
+ }
+
+private:
+ std::array<MapType, numSD> stencils_;
+ std::vector<std::size_t> emptyStencil_;
+ std::array<std::vector<bool>, numSD> isCoupledScvf_;
+ std::unordered_map<std::size_t, bool> isCoupledLateralScvf_;
+ std::vector<bool> isCoupledFFDof_;
+ std::vector<bool> isCoupledFFElement_;
+ std::tuple<FlipScvfMapTypeFF, FlipScvfMapTypePM> scvfInfo_;
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/multibinarycouplingmanager.hh b/dumux/multidomain/multibinarycouplingmanager.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cd335ffb5bc48361b169ae87d43f27e9512f2df7
--- /dev/null
+++ b/dumux/multidomain/multibinarycouplingmanager.hh
@@ -0,0 +1,397 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup MultiDomain
+ * \copydoc Dumux::MultiBinaryCouplingManager
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_MULTIBINARY_COUPLINGMANAGER_HH
+#define DUMUX_MULTIDOMAIN_MULTIBINARY_COUPLINGMANAGER_HH
+
+#include <utility>
+#include <memory>
+#include <dune/common/hybridutilities.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/multidomain/traits.hh>
+
+namespace Dumux {
+
+namespace Detail {
+
+template <std::size_t, typename Tuple>
+struct HasIndex;
+
+template <std::size_t i, typename... Indices>
+struct HasIndex<i, std::tuple<Indices...>>
+: std::disjunction<std::is_same<Dune::index_constant<i>, Indices>...>
+{};
+
+} // end namespace Detail
+
+/*!
+ * \ingroup MultiDomain
+ * \brief Coupling manager that combines an arbitrary number of binary coupling manager (coupling two domains each)
+ * \tparam MDTraits the multidomain traits
+ * \tparam CouplingMap a coupling policy class
+ * \tparam CouplingMgrs the binary sub-coupling manager types
+ *
+ * The coupling policy has to provide the interfaces
+ * - CouplingMap::coupledDomains(i): returns a tuple of Dune::index_constants with the coupled domains
+ * - CouplingMap::globalToLocal(i, j): maps the indices i, j to the local index pair of the responsible sub coupling manager
+ * - CouplingMap::managerMap(): returns a two-dimensional array mapping two indices to the coupling manager index
+ */
+template<class MDTraits, class CouplingMap, class ...CouplingMgrs>
+class MultiBinaryCouplingManager
+{
+ // 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 GridView = typename GridGeometry<id>::GridView;
+ template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ 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;
+ using CouplingStencil = std::vector<std::size_t>;
+
+ template<std::size_t id>
+ using SubCouplingManagerT = typename std::tuple_element_t<id, std::tuple<CouplingMgrs...>>;
+
+ using CMIndices = std::make_index_sequence<sizeof...(CouplingMgrs)>;
+ using CouplingManagers = typename Detail::MultiDomainTupleSharedPtr<SubCouplingManagerT, CMIndices>::type;
+
+ template<std::size_t id>
+ using SubSolutionVector = std::decay_t<decltype(std::declval<typename MDTraits::SolutionVector>()[Dune::index_constant<id>()])>;
+ using SolutionVectors = typename MDTraits::template TupleOfSharedPtr<SubSolutionVector>;
+
+ static constexpr auto couplingManagerMap_ = CouplingMap::managerMap();
+
+ //! Returns the local domain indices used within the binary sub coupling managers.
+ template<std::size_t i, std::size_t j>
+ static constexpr auto globalToLocal_(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
+ {
+ static_assert(i <= MDTraits::numSubDomains && j <= MDTraits::numSubDomains);
+ return CouplingMap::globalToLocal(domainI, domainJ);
+ }
+
+ //! If two domain are coupled
+ template<class Map, std::size_t i, std::size_t j>
+ static constexpr bool isCoupled_(Dune::index_constant<i> domainI, Dune::index_constant<j>)
+ { return Detail::HasIndex<j, std::decay_t<decltype(Map::coupledDomains(domainI))>>::value; }
+
+ //! Returns the coupling manager index for a given domain combination
+ template<std::size_t i, std::size_t j>
+ static constexpr auto subCouplingManagerIndex_(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
+ {
+ static_assert(
+ isCoupled_<CouplingMap>(domainI, domainJ),
+ "Sub-coupling manager only exists for coupled domains."
+ );
+ return couplingManagerMap_[i][j];
+ }
+
+public:
+ template<std::size_t i, std::size_t j>
+ using SubCouplingManager = SubCouplingManagerT<couplingManagerMap_[i][j]>;
+
+ MultiBinaryCouplingManager()
+ {
+ using namespace Dune::Hybrid;
+ forEach(couplingManagers_, [&](auto&& couplingManager)
+ {
+ couplingManager = std::make_shared<typename std::decay_t<decltype(couplingManager)>::element_type>();
+ });
+
+ forEach(solutionVectors_, [&](auto&& solutionVector)
+ {
+ solutionVector = std::make_shared<typename std::decay_t<decltype(solutionVector)>::element_type>();
+ });
+ }
+
+ //! return the binary sub-coupling manager
+ template<std::size_t i, std::size_t j>
+ auto& subCouplingManager(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ)
+ {
+ constexpr auto idx = subCouplingManagerIndex_(domainI, domainJ);
+ return *std::get<idx>(couplingManagers_);
+ }
+
+ //! return the binary sub-coupling manager
+ template<std::size_t i, std::size_t j>
+ const auto& subCouplingManager(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ) const
+ {
+ constexpr auto idx = subCouplingManagerIndex_(domainI, domainJ);
+ return *std::get<idx>(couplingManagers_);
+ }
+
+ //! apply a function to the domainI-domainJ sub coupling manager using its local indices
+ template<std::size_t i, std::size_t j, class Apply>
+ decltype(auto) subApply(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ Apply&& apply)
+ {
+ constexpr auto localIndices = globalToLocal_(domainI, domainJ);
+ return apply(subCouplingManager(domainI, domainJ), localIndices.first, localIndices.second);
+ }
+
+ //! apply a function to the domainI-domainJ sub coupling manager using its local indices
+ template<std::size_t i, std::size_t j, class Apply>
+ decltype(auto) subApply(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const Apply& apply) const
+ {
+ constexpr auto localIndices = globalToLocal_(domainI, domainJ);
+ return apply(subCouplingManager(domainI, domainJ), localIndices.first, localIndices.second);
+ }
+
+ //! apply a function to a sub coupling manager containing this domain
+ template<std::size_t i, class Apply>
+ decltype(auto) subApply(Dune::index_constant<i> domainI, Apply&& apply)
+ {
+ constexpr auto dm = CouplingMap::coupledDomains(domainI);
+ static_assert(std::tuple_size_v<std::decay_t<decltype(dm)>> != 0, "subApply on uncoupled domain!");
+ constexpr auto domainJ = std::get<0>(dm);
+ constexpr auto localIndices = globalToLocal_(domainI, domainJ);
+ return apply(subCouplingManager(domainI, domainJ), localIndices.first);
+ }
+
+ //! apply a function to a sub coupling manager containing this domain
+ template<std::size_t i, class Apply>
+ decltype(auto) subApply(Dune::index_constant<i> domainI, const Apply& apply) const
+ {
+ constexpr auto dm = CouplingMap::coupledDomains(domainI);
+ static_assert(std::tuple_size_v<std::decay_t<decltype(dm)>> != 0, "subApply on uncoupled domain!");
+ constexpr auto domainJ = std::get<0>(dm);
+ constexpr auto localIndices = globalToLocal_(domainI, domainJ);
+ return apply(subCouplingManager(domainI, domainJ), localIndices.first);
+ }
+
+ //! Update the solution vector before assembly
+ void updateSolution(const typename MDTraits::SolutionVector& curSol)
+ {
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Dune::Hybrid::size(solutionVectors_)), [&](const auto id)
+ {
+ *std::get<id>(solutionVectors_) = curSol[id];
+ });
+ }
+
+ /*!
+ * \brief extend the jacobian pattern of the diagonal block of domain i
+ * by those entries that are not already in the uncoupled pattern
+ * \note per default we do not add such additional dependencies
+ * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
+ * term depends on a non-local average of a quantity of the same domain
+ * \warning if you overload this also implement evalAdditionalDomainDerivatives
+ */
+ template<std::size_t id, class JacobianPattern>
+ void extendJacobianPattern(Dune::index_constant<id> domainI, JacobianPattern& pattern) const
+ {}
+
+ /*!
+ * \brief Return the coupling element stencil for a given bulk domain element
+ */
+ template<std::size_t i, class Entity, std::size_t j>
+ const auto& couplingStencil(Dune::index_constant<i> domainI,
+ const Entity& entity,
+ Dune::index_constant<j> domainJ) const
+ {
+ // if the domains are coupled according to the map, forward to sub-coupling manager
+ if constexpr (isCoupled_<CouplingMap>(domainI, domainJ))
+ return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.couplingStencil(ii, entity, jj);
+ });
+ else
+ return emptyStencil_;
+ }
+
+ // ! evaluate coupling residual for the derivative low dim DOF with respect to bulk DOF
+ // ! we only need to evaluate the part of the residual that will be influence by the bulk DOF
+ template<std::size_t i, class LocalAssemblerI, std::size_t j>
+ decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
+ const SubControlVolumeFace<i>& scvfI,
+ const LocalAssemblerI& localAssemblerI,
+ Dune::index_constant<j> domainJ,
+ std::size_t dofIdxGlobalJ) const
+ {
+ if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
+ return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
+ return cm.evalCouplingResidual(ii, scvfI, localAssemblerI, jj, dofIdxGlobalJ);
+ });
+ else
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
+ );
+ return localAssemblerI.evalLocalResidual();
+ }
+ }
+
+ //! evaluate coupling residual for the derivative low dim DOF with respect to bulk DOF
+ //! we only need to evaluate the part of the residual that will be influence by the bulk DOF
+ template<std::size_t i, class LocalAssemblerI, std::size_t j>
+ decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
+ const LocalAssemblerI& localAssemblerI,
+ Dune::index_constant<j> domainJ,
+ std::size_t dofIdxGlobalJ) const
+ {
+ if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
+ return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
+ return cm.evalCouplingResidual(ii, localAssemblerI, jj, dofIdxGlobalJ);
+ });
+ else
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
+ );
+ return localAssemblerI.evalLocalResidual();
+ }
+ }
+
+ //! evaluate coupling residual for the derivative low dim DOF with respect to bulk DOF
+ //! we only need to evaluate the part of the residual that will be influence by the bulk DOF
+ template<std::size_t i, class LocalAssemblerI, std::size_t j>
+ decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
+ const LocalAssemblerI& localAssemblerI,
+ const SubControlVolume<i>& scvI,
+ Dune::index_constant<j> domainJ,
+ std::size_t dofIdxGlobalJ) const
+ {
+ if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
+ return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
+ return cm.evalCouplingResidual(ii, localAssemblerI, scvI, jj, dofIdxGlobalJ);
+ });
+ else
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
+ );
+ return localAssemblerI.evalLocalResidual();
+ }
+ }
+
+ /*!
+ * \brief Update the coupling context for the bulk face residual w.r.t to the lowDim dofs
+ */
+ template<std::size_t i, class LocalAssemblerI, std::size_t j, class PrimaryVariables>
+ void updateCouplingContext(Dune::index_constant<i> domainI,
+ const LocalAssemblerI& localAssemblerI,
+ Dune::index_constant<j> domainJ,
+ const std::size_t dofIdxGlobalJ,
+ const PrimaryVariables& priVars,
+ int pvIdxJ)
+ {
+ // only one other manager needs to take action if i != j
+ if constexpr (i != j)
+ {
+ if constexpr (isCoupled_<CouplingMap>(domainI, domainJ))
+ subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
+ cm.updateCouplingContext(ii, localAssemblerI, jj, dofIdxGlobalJ, priVars, pvIdxJ);
+ });
+ }
+ else
+ {
+ // for i == j, we need to call all relevant managers where domainI is involved and make sure that the
+ // context is updated with respect to its own domain (I-I coupling context)
+ using namespace Dune::Hybrid;
+ forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
+ static_assert(domainI != domainJ);
+ subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
+ cm.updateCouplingContext(ii, localAssemblerI, ii, dofIdxGlobalJ, priVars, pvIdxJ);
+ });
+ });
+ }
+ }
+
+ //! Bind the coupling context for a low dim element TODO remove Assembler
+ template<std::size_t i, class Assembler = int>
+ void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler = 0)
+ {
+ // for the coupling blocks
+ using namespace Dune::Hybrid;
+ forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
+ subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj) -> void {
+ cm.bindCouplingContext(ii, element, assembler);
+ });
+ });
+ }
+
+ /*!
+ * \brief return the numeric epsilon used for deflecting primary variables of coupled domain i.
+ * \note specialization for free-flow schemes
+ */
+ template<std::size_t i>
+ decltype(auto) numericEpsilon(Dune::index_constant<i> domainI,
+ const std::string& paramGroup) const
+ {
+ return subApply(domainI, [&](const auto& cm, auto&& ii) -> decltype(auto) {
+ return cm.numericEpsilon(ii, paramGroup);
+ });
+ }
+
+ /*!
+ * \brief evaluate additional derivatives of the element residual of a domain with respect
+ * to dofs in the same domain that are not in the regular stencil (see CouplingManager::extendJacobianPattern)
+ * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
+ * term depends on a non-local average of a quantity of the same domain
+ */
+ template<std::size_t i, class LocalAssemblerI, class JacobianMatrixDiagBlock, class GridVariables>
+ void evalAdditionalDomainDerivatives(Dune::index_constant<i> domainI,
+ const LocalAssemblerI& localAssemblerI,
+ const typename LocalAssemblerI::LocalResidual::ElementResidualVector& origResiduals,
+ JacobianMatrixDiagBlock& A,
+ GridVariables& gridVariables)
+ {}
+
+ template<std::size_t i, class LocalAssemblerI, class UpdatableElementVolVars, class UpdatableFluxVarCache>
+ void updateCoupledVariables(Dune::index_constant<i> domainI,
+ const LocalAssemblerI& localAssemblerI,
+ UpdatableElementVolVars& elemVolVars,
+ UpdatableFluxVarCache& elemFluxVarsCache)
+ {
+ // for the coupling blocks
+ using namespace Dune::Hybrid;
+ forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
+ subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
+ cm.updateCoupledVariables(ii, localAssemblerI, elemVolVars, elemFluxVarsCache);
+ });
+ });
+ }
+
+protected:
+ SolutionVectors& curSol()
+ { return solutionVectors_; }
+
+ const SolutionVectors& curSol() const
+ { return solutionVectors_; }
+
+private:
+ CouplingManagers couplingManagers_;
+ SolutionVectors solutionVectors_;
+
+ CouplingStencil emptyStencil_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/staggeredfreeflow/couplingmanager.hh b/dumux/multidomain/staggeredfreeflow/couplingmanager.hh
index d04276a6e468bf2f991ee5d6214b8d3cbe1ff5ce..08fc15548a45d3e381cb067f36f78ef4a344ed07 100644
--- a/dumux/multidomain/staggeredfreeflow/couplingmanager.hh
+++ b/dumux/multidomain/staggeredfreeflow/couplingmanager.hh
@@ -48,9 +48,14 @@ template<class Traits>
class StaggeredFreeFlowCouplingManager
: public CouplingManager<Traits>
{
+ using ParentType = CouplingManager<Traits>;
public:
static constexpr auto freeFlowMomentumIndex = typename Traits::template SubDomain<0>::Index();
static constexpr auto freeFlowMassIndex = typename Traits::template SubDomain<1>::Index();
+
+ // this can be used if the coupling manager is used inside a meta-coupling manager (e.g. multi-binary)
+ // to manager the solution vector storage outside this class
+ using SolutionVectorStorage = typename ParentType::SolutionVectorStorage;
private:
template<std::size_t id> using SubDomainTypeTag = typename Traits::template SubDomain<id>::TypeTag;
template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
@@ -69,8 +74,6 @@ private:
using Scalar = typename Traits::Scalar;
using SolutionVector = typename Traits::SolutionVector;
- using ParentType = CouplingManager<Traits>;
-
using CouplingStencilType = std::vector<std::size_t>;
using GridVariablesTuple = typename Traits::template TupleOfSharedPtr<GridVariables>;
diff --git a/test/multidomain/boundary/CMakeLists.txt b/test/multidomain/boundary/CMakeLists.txt
index addd460feb4d2b502e6a92a0630ad445b11020d9..5ffc5873b26785265ec0c9ad64e4cc7aa93cce1b 100644
--- a/test/multidomain/boundary/CMakeLists.txt
+++ b/test/multidomain/boundary/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory(darcydarcy)
+add_subdirectory(freeflowporousmedium)
add_subdirectory(stokesdarcy)
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt b/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..617ade1705be78321c072af769414efc7b66d337
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt
@@ -0,0 +1,32 @@
+add_subdirectory(convergence)
+
+add_input_file_links()
+
+add_executable(test_md_boundary_ff1p_pm1p EXCLUDE_FROM_ALL main.cc)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_horizontal
+ LABELS multidomain multidomain_boundary stokesdarcy 1p navierstokes
+ TARGET test_md_boundary_ff1p_pm1p
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_horizontal_stokes-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p_horizontal_freeflow-00000.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_horizontal_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p_horizontal_darcy-00000.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p params_horizontalflow.input
+ -Vtk.OutputName test_md_boundary_ff1p_pm1p_horizontal")
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_vertical
+ LABELS multidomain multidomain_boundary stokesdarcy 1p navierstokes
+ TARGET test_md_boundary_ff1p_pm1p
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_vertical_stokes-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p_vertical_freeflow-00000.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_vertical_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p_vertical_darcy-00000.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p params_verticalflow.input
+ -Vtk.OutputName test_md_boundary_ff1p_pm1p_vertical"
+ --zeroThreshold {"velocity_liq \(m/s\)_0":6e-17})
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/CMakeLists.txt b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8336a0b30df4dc11d70415c1f3e81e6be49475f5
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/CMakeLists.txt
@@ -0,0 +1,50 @@
+add_input_file_links()
+dune_symlink_to_source_files(FILES "convergencetest.py")
+
+add_executable(test_md_boundary_ff1p_pm1p_convtest EXCLUDE_FROM_ALL main.cc)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_convtest_shiue1
+ TARGET test_md_boundary_ff1p_pm1p_convtest
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_ff1p_pm1p_convtest params.input
+ -Problem.TestCase ShiueExampleOne
+ -Darcy.Problem.BoundaryConditions Dirichlet
+ -FreeFlow.Problem.BoundaryConditions Dirichlet)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_convtest_shiue2
+ TARGET test_md_boundary_ff1p_pm1p_convtest
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_ff1p_pm1p_convtest params.input
+ -Problem.TestCase ShiueExampleTwo
+ -Darcy.Problem.BoundaryConditions Mixed
+ -FreeFlow.Problem.BoundaryConditions Mixed)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_convtest_rybak
+ TARGET test_md_boundary_ff1p_pm1p_convtest
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_ff1p_pm1p_convtest params.input
+ -Problem.TestCase Rybak
+ -Darcy.Problem.BoundaryConditions Neumann
+ -FreeFlow.Problem.BoundaryConditions Stress)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_convtest_navierstokes
+ TARGET test_md_boundary_ff1p_pm1p_convtest
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_ff1p_pm1p_convtest params.input
+ -Problem.TestCase Schneider
+ -FreeFlow.Problem.EnableInertiaTerms true
+ -FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph true
+ -Darcy.Problem.BoundaryConditions Mixed
+ -FreeFlow.Problem.BoundaryConditions Mixed)
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/analyticalsolutions.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/analyticalsolutions.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6df9efadd3fba9d9581d38606c584d73af590fe9
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/analyticalsolutions.hh
@@ -0,0 +1,340 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief Analytical solutions and rhs for the different test cases
+ */
+
+#ifndef DUMUX_CONVERGENCE_TEST_ANALYTICAL_SOLUTIONS_HH
+#define DUMUX_CONVERGENCE_TEST_ANALYTICAL_SOLUTIONS_HH
+
+#include <cmath>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+namespace Dumux::Solution::DarcyStokes {
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// see Rybak et al., 2015: "Multirate time integration for coupled
+// saturated/unsaturated porous medium and free flow systems"
+/////////////////////////////////////////////////////////////////////////////////////////////////
+namespace Rybak {
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[0] = -0.5*M_PI*y*y*cos(M_PI*x);
+ sol[1] = -1.0*y*sin(M_PI*x);
+ sol[2] = 0.5*y*y*sin(M_PI*x);
+ return sol;
+}
+
+// 0: mass
+Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::sin;
+ return { (0.5*M_PI*y*M_PI*y - 1)*sin(M_PI*x) };
+}
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::sin; using std::cos;
+ sol[0] = -cos(M_PI*x)*sin(M_PI*y);
+ sol[1] = sin(M_PI*x)*cos(M_PI*y);
+ sol[2] = 0.5*y*sin(M_PI*x);
+ return sol;
+}
+
+// 0: mom-x | 1: mom-y | 2: mass
+Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ Dune::FieldVector<double, 3> source(0.0);
+ using std::sin; using std::cos;
+ source[0] = 0.5*M_PI*y*cos(M_PI*x) - 2*M_PI*M_PI*cos(M_PI*x)*sin(M_PI*y);
+ source[1] = 0.5*sin(M_PI*x) + 2*M_PI*M_PI*sin(M_PI*x)*cos(M_PI*y);
+ return source;
+}
+
+// sigma_ff = -sym(grad(v)) + pI
+Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+
+ Dune::FieldMatrix<double, 2, 2> stress(0.0);
+ stress[0][0] = (0.5*y - 2*M_PI*sin(M_PI*y))*sin(M_PI*x);
+ stress[1][0] = 0.0;
+ stress[0][1] = stress[1][0]; // symmetric
+ stress[1][1] = (0.5*y - 2*M_PI*sin(M_PI*y))*sin(M_PI*x);
+ return stress;
+}
+
+} // end namespace Rybak
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+// Example 1
+/////////////////////////////////////////////////////////////////////////////////////////////////
+namespace ShiueOne {
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[0] = M_PI*(-exp(1)*y + exp(y))*sin(M_PI*x);
+ sol[1] = (exp(1) - exp(y))*cos(M_PI*x);
+ sol[2] = (exp(y) - y*exp(1)) * cos(M_PI*x);
+
+ return sol;
+}
+
+// 0: mass
+Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ return { cos(M_PI*x) * (M_PI*M_PI*(exp(y) - y*exp(1)) - exp(y)) };
+}
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[0] = -1/M_PI * exp(y) * sin(M_PI*x);
+ sol[1] = (exp(y) - exp(1)) * cos(M_PI*x);
+ sol[2] = 2*exp(y) * cos(M_PI*x);
+ return sol;
+}
+
+// 0: mom-x | 1: mom-y | 2: mass
+Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ Dune::FieldVector<double, 3> source(0.0);
+ source[0] = exp(y)*sin(M_PI*x) * (1/M_PI -3*M_PI);
+ source[1] = cos(M_PI*x) * (M_PI*M_PI*(exp(y)- exp(1)) + exp(y));
+ return source;
+}
+
+// sigma_ff = -sym(grad(v)) + pI
+Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+
+ Dune::FieldMatrix<double, 2, 2> stress(0.0);
+ stress[0][0] = 4.0*exp(y)*cos(M_PI*x);
+ stress[1][0] = (M_PI*M_PI*(exp(y) - exp(1)) + 1.0*exp(y))*sin(M_PI*x)/M_PI;
+ stress[0][1] = stress[1][0]; // symmetric
+ stress[1][1] = 0.0;
+ return stress;
+}
+
+} // end namespace ShiueExOne
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+// Example 2
+/////////////////////////////////////////////////////////////////////////////////////////////////
+namespace ShiueTwo {
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ sol[0] = x*(y - 1.0) + (x - 1.0)*(y - 1.0) - 2.0;
+ sol[1] = x*(x - 1.0) - 1.0*(y - 1.0)*(y - 1.0) - 2.0;
+ sol[2] = x*(1.0-x)*(y-1.0) + (y-1.0)*(y-1.0)*(y-1.0)/3.0 + 2.0*x + 2.0*y + 4.0;
+
+ return sol;
+}
+
+// 0: mass
+Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
+{ return { 0.0 }; }
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ sol[0] = (y-1.0)*(y-1.0) + x*(y-1.0) + 3.0*x - 1.0;
+ sol[1] = x*(x-1.0) - 0.5*(y-1.0)*(y-1.0) - 3.0*y + 1.0;
+ sol[2] = 2.0*x + y - 1.0;
+ return sol;
+}
+
+// 0: mom-x | 1: mom-y | 2: mass
+Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
+{ return { 0.0, 0.0, 0.0 }; }
+
+// sigma_ff = -sym(grad(v)) + pI
+Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ Dune::FieldMatrix<double, 2, 2> stress(0.0);
+ stress[0][0] = 2.0*x - y - 5.0;
+ stress[1][0] = -3*x - 2*y + 3.0;
+ stress[0][1] = stress[1][0]; // symmetric
+ stress[1][1] = 2.0*x + 3.0*y + 3.0;
+ return stress;
+}
+
+} // end namespace ShiueTwo
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for
+// free flow/porous-medium flow problems (with c = 0)"
+/////////////////////////////////////////////////////////////////////////////////////////////////
+namespace Schneider {
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ static constexpr double omega = M_PI;
+ static constexpr double c = 0.0;
+ using std::exp; using std::sin; using std::cos;
+ const double sinOmegaX = sin(omega*x);
+ const double cosOmegaX = cos(omega*x);
+ static const double expTwo = exp(2);
+ const double expYPlusOne = exp(y+1);
+
+ sol[0] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX
+ -omega*(expYPlusOne + 2 - expTwo)*cosOmegaX;
+ sol[1] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX
+ -(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX;
+ sol[2] = (expYPlusOne + 2 - expTwo)*sinOmegaX;
+
+ return sol;
+}
+
+// 0: mass
+Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ static constexpr double omega = M_PI;
+ static constexpr double c = 0.0;
+ const double cosOmegaX = cos(omega*x);
+ static const double expTwo = exp(2);
+ const double expYPlusOne = exp(y+1);
+
+ return { sin(omega*x)*(-(c*cosOmegaX + 1)*exp(y - 1)
+ + 1.5*c*expYPlusOne*cosOmegaX
+ + omega*omega*(expYPlusOne - expTwo + 2)) };
+}
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::sin;
+ static constexpr double omega = M_PI;
+ const double sinOmegaX = sin(omega*x);
+
+ sol[0] = y;
+ sol[1] = -y*sinOmegaX;
+ sol[2] = -y*y*sinOmegaX*sinOmegaX;
+ return sol;
+}
+
+// 0: mom-x | 1: mom-y | 2: mass
+Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ static constexpr double omega = M_PI;
+ const double sinOmegaX = sin(omega*x);
+ const double cosOmegaX = cos(omega*x);
+
+ Dune::FieldVector<double, 3> source(0.0);
+ source[0] = -2*omega*y*y*sinOmegaX*cosOmegaX
+ -2*y*sinOmegaX + omega*cosOmegaX;
+ source[1] = -omega*y*y*cosOmegaX - omega*omega*y*sinOmegaX;
+ source[2] = -sinOmegaX;
+ return source;
+}
+
+// sigma_ff = vvT + sym(grad(v)) + pI
+Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
+{
+ using std::sin; using std::cos;
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ static constexpr double omega = M_PI;
+ const double sinOmegaX = sin(omega*x);
+ const double cosOmegaX = cos(omega*x);
+
+ Dune::FieldMatrix<double, 2, 2> stress(0.0);
+ stress[0][0] = y*y * cosOmegaX*cosOmegaX;
+ stress[0][1] = -y*y*sinOmegaX + omega*y*cosOmegaX - 1.0;
+ stress[1][0] = stress[0][1]; // symmetric
+ stress[1][1] = 2*sinOmegaX;
+ return stress;
+}
+
+} // end namespace Schneider
+
+} // end namespace Dumux::Solution::DarcyStokes
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/convergencetest.py b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/convergencetest.py
new file mode 100755
index 0000000000000000000000000000000000000000..b52826c3c48cc713f2f526dffe2bfb551e720cb2
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/convergencetest.py
@@ -0,0 +1,124 @@
+#!/usr/bin/env python3
+
+from math import *
+import subprocess
+import sys
+
+if len(sys.argv) < 2:
+ sys.stderr.write('Please provide a single argument <testname> to the script\n')
+ sys.exit(1)
+
+executableName = str(sys.argv[1])
+testargs = [str(i) for i in sys.argv][2:]
+testname = testargs[testargs.index('-Problem.TestCase')+1]
+
+# remove the old log files
+subprocess.call(['rm', testname + '_freeFlow.log'])
+print("Removed old log file ({})!".format(testname + '_freeFlow.log'))
+subprocess.call(['rm', testname + '_darcy.log'])
+print("Removed old log file ({})!".format(testname + '_darcy.log'))
+
+# do the runs with different refinement
+for i in [0, 1, 2]:
+ subprocess.call(['./' + executableName] + testargs + ['-Grid.Refinement', str(i)])
+
+def checkRatesFreeFlow():
+ # check the rates and append them to the log file
+ with open(testname + '_freeFlow.log', "r+") as logfile:
+
+ errorP, errorVx, errorVy = [], [], []
+ for line in logfile:
+ line = line.strip("\n").strip("\[ConvergenceTest\]").split()
+ errorP.append(float(line[2]))
+ errorVx.append(float(line[5]))
+ errorVy.append(float(line[8]))
+
+ resultsP, resultsVx, resultsVy = [], [], []
+ logfile.truncate(0)
+ logfile.write("n\terrorP\t\trateP\t\terrorVx\t\trateVx\t\terrorVy\t\trateVy\n")
+ logfile.write("-"*50 + "\n")
+ for i in range(len(errorP)-1):
+ if isnan(errorP[i]) or isinf(errorP[i]):
+ continue
+ if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12) and (errorVx[i] < 1e-12 or errorVx[i+1] < 1e-12) and (errorVy[i] < 1e-12 or errorVy[i+1] < 1e-12)):
+ rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
+ rateVx = (log(errorVx[i])-log(errorVx[i+1]))/log(2)
+ rateVy = (log(errorVy[i])-log(errorVy[i+1]))/log(2)
+ message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP, errorVx[i], rateVx, errorVy[i], rateVy)
+ logfile.write(message)
+ resultsP.append(rateP)
+ resultsVx.append(rateVx)
+ resultsVy.append(rateVy)
+ else:
+ logfile.write("error: exact solution!?")
+ i = len(errorP)-1
+ message = "{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\n".format(i, errorP[i], "", errorVx[i], "", errorVy[i], "")
+ logfile.write(message)
+
+ print("\nComputed the following convergence rates for {}:\n".format(testname))
+
+ subprocess.call(['cat', testname + '_freeFlow.log'])
+
+ return {"p" : resultsP, "v_x" : resultsVx, "v_y" : resultsVy}
+
+def checkRatesDarcy():
+ # check the rates and append them to the log file
+ with open(testname + '_darcy.log', "r+") as logfile:
+
+ errorP = []
+ for line in logfile:
+ line = line.strip("\n")
+ line = line.strip("\[ConvergenceTest\]")
+ line = line.split()
+ errorP.append(float(line[2]))
+
+ resultsP = []
+ logfile.truncate(0)
+ logfile.write("n\terrorP\t\trateP\n")
+ logfile.write("-"*50 + "\n")
+ for i in range(len(errorP)-1):
+ if isnan(errorP[i]) or isinf(errorP[i]):
+ continue
+ if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12)):
+ rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
+ message = "{}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP)
+ logfile.write(message)
+ resultsP.append(rateP)
+ else:
+ logfile.write("error: exact solution!?")
+ i = len(errorP)-1
+ message = "{}\t{:0.4e}\n".format(i, errorP[i], "")
+ logfile.write(message)
+
+ print("\nComputed the following convergence rates for {}:\n".format(testname))
+
+ subprocess.call(['cat', testname + '_darcy.log'])
+
+ return {"p" : resultsP}
+
+def checkRatesFreeFlowAndDarcy():
+ resultsFreeFlow = checkRatesFreeFlow()
+ resultsDarcy = checkRatesDarcy()
+
+ def mean(numbers):
+ return float(sum(numbers)) / len(numbers)
+
+ # check the rates, we expect rates around 2
+ if mean(resultsFreeFlow["p"]) < 2.05 and mean(resultsFreeFlow["p"]) < 1.7:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+ if mean(resultsFreeFlow["v_x"]) < 2.05 and mean(resultsFreeFlow["v_x"]) < 1.95:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates for x-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+ if mean(resultsFreeFlow["v_y"]) < 2.05 and mean(resultsFreeFlow["v_y"]) < 1.95:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates for y-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+ if mean(resultsDarcy["p"]) < 2.05 and mean(resultsDarcy["p"]) < 1.95:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+
+checkRatesFreeFlowAndDarcy()
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/main.cc b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8a9a0044d80d02a2b5b3d8d08ee2301cf9edc6ae
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/main.cc
@@ -0,0 +1,320 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief A convergence test for the coupled FreeFlow/Darcy problem (1p).
+ */
+
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/format.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+#include <dumux/freeflow/navierstokes/velocityoutput.hh>
+
+#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
+#include <test/freeflow/navierstokes/errors.hh>
+
+#include "testcase.hh"
+#include "properties.hh"
+
+/*!
+* \brief Creates analytical solution vector for output
+* Returns a tuple of the analytical solution for the pressure and velocity at cell centers
+* \param problem the problem for which to evaluate the analytical solution
+*/
+template<class Scalar, class Problem>
+auto createDarcyAnalyticalSolution(const Problem& problem)
+{
+ const auto& gridGeometry = problem.gridGeometry();
+ using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
+
+ static constexpr auto dim = GridView::dimension;
+ static constexpr auto dimWorld = GridView::dimensionworld;
+
+ using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+
+ std::vector<Scalar> analyticalPressure;
+ std::vector<VelocityVector> analyticalVelocity;
+
+ analyticalPressure.resize(gridGeometry.numDofs());
+ analyticalVelocity.resize(gridGeometry.numDofs());
+
+ for (const auto& element : elements(gridGeometry.gridView()))
+ {
+ auto fvGeometry = localView(gridGeometry);
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto ccDofIdx = scv.dofIndex();
+ const auto ccDofPosition = scv.dofPosition();
+ const auto analyticalSolutionAtCc = problem.fullAnalyticalSolution(ccDofPosition);
+ analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[dim];
+
+ for (int dirIdx = 0; dirIdx < dim; ++dirIdx)
+ analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[dirIdx];
+ }
+ }
+
+ return std::make_tuple(analyticalPressure, analyticalVelocity);
+}
+
+template<class MassP, class MomP, class SolutionVector>
+void printFreeFlowL2Error(std::shared_ptr<MomP> momentumProblem,
+ std::shared_ptr<MassP> massProblem,
+ const SolutionVector& sol)
+{
+ using namespace Dumux;
+
+ NavierStokesTest::Errors errors(momentumProblem, massProblem, sol);
+ const int numCellCenterDofs = massProblem->gridGeometry().numDofs();
+ const int numFaceDofs = momentumProblem->gridGeometry().numDofs();
+ const auto absL2 = errors.l2Absolute();
+ const auto relL2 = errors.l2Relative();
+
+ std::cout << Fmt::format("** L2 error (abs/rel) for {} cc dofs and {} face dofs (total: {}): ",
+ numCellCenterDofs, numFaceDofs, numCellCenterDofs + numFaceDofs)
+ << Fmt::format("L2(p) = {:.8e} / {:.8e}", absL2[0], relL2[0])
+ << Fmt::format(", L2(vx) = {:.8e} / {:.8e}", absL2[1], relL2[1])
+ << Fmt::format(", L2(vy) = {:.8e} / {:.8e}", absL2[2], relL2[2])
+ << std::endl;
+
+ // write the norm into a log file
+ std::ofstream logFile(massProblem->name() + ".log", std::ios::app);
+ logFile << "[ConvergenceTest] L2(p) = " << absL2[0]
+ << " L2(vx) = " << absL2[1]
+ << " L2(vy) = " << absL2[2]
+ << std::endl;
+}
+
+template<class Problem, class SolutionVector>
+void printDarcyL2Error(std::shared_ptr<Problem> problem, const SolutionVector& x)
+{
+ using namespace Dumux;
+ using Scalar = double;
+
+ Scalar l2error = 0.0;
+ for (const auto& element : elements(problem->gridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(problem->gridGeometry());
+ fvGeometry.bindElement(element);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto dofIdx = scv.dofIndex();
+ const Scalar delta = x[dofIdx] - problem->fullAnalyticalSolution(scv.center())[2/*pressureIdx*/];
+ l2error += scv.volume()*(delta*delta);
+ }
+ }
+ using std::sqrt;
+ l2error = sqrt(l2error);
+
+ const auto numDofs = problem->gridGeometry().numDofs();
+ std::cout << Fmt::format("** L2 error (abs) for {} cc dofs", numDofs)
+ << Fmt::format("L2 error = {:.8e}", l2error)
+ << std::endl;
+
+ // write the norm into a log file
+ std::ofstream logFile(problem->name() + ".log", std::ios::app);
+ logFile << "[ConvergenceTest] L2(p) = " << l2error << std::endl;
+}
+
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // Define the sub problem type tags
+ using FreeFlowMomentumTypeTag = Properties::TTag::FreeFlowOnePMomentum;
+ using FreeFlowMassTypeTag = Properties::TTag::FreeFlowOnePMass;
+ using DarcyTypeTag = Properties::TTag::DarcyOneP;
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using FreeFlowGridManager = Dumux::GridManager<GetPropType<FreeFlowMomentumTypeTag, Properties::Grid>>;
+ FreeFlowGridManager freeFlowGridManager;
+ freeFlowGridManager.init("FreeFlow"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& freeFlowGridView = freeFlowGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FreeFlowMomentumGridGeometry = GetPropType<FreeFlowMomentumTypeTag, Properties::GridGeometry>;
+ auto freeFlowMomentumGridGeometry = std::make_shared<FreeFlowMomentumGridGeometry>(freeFlowGridView);
+ using FreeFlowMassGridGeometry = GetPropType<FreeFlowMassTypeTag, Properties::GridGeometry>;
+ auto freeFlowMassGridGeometry = std::make_shared<FreeFlowMassGridGeometry>(freeFlowGridView);
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
+ auto darcyGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+
+ using Traits = MultiDomainTraits<FreeFlowMomentumTypeTag, FreeFlowMassTypeTag, DarcyTypeTag>;
+ using CouplingManager = FreeFlowPorousMediumCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>();
+
+ // the indices
+ constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex;
+ constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex;
+ constexpr auto porousMediumIndex = CouplingManager::porousMediumIndex;
+
+ // the problem (initial and boundary conditions)
+ const auto testCaseName = getParam<std::string>("Problem.TestCase");
+ const auto testCase = [&]()
+ {
+ if (testCaseName == "ShiueExampleOne")
+ return DarcyStokesTestCase::ShiueExampleOne;
+ else if (testCaseName == "ShiueExampleTwo")
+ return DarcyStokesTestCase::ShiueExampleTwo;
+ else if (testCaseName == "Rybak")
+ return DarcyStokesTestCase::Rybak;
+ else if (testCaseName == "Schneider")
+ return DarcyStokesTestCase::Schneider;
+ else
+ DUNE_THROW(Dune::InvalidStateException, testCaseName + " is not a valid test case");
+ }();
+
+ // the problem (initial and boundary conditions)
+ using FreeFlowMomentumProblem = GetPropType<FreeFlowMomentumTypeTag, Properties::Problem>;
+ auto freeFlowMomentumProblem = std::make_shared<FreeFlowMomentumProblem>(freeFlowMomentumGridGeometry, couplingManager, testCase, testCaseName);
+ using FreeFlowMassProblem = GetPropType<FreeFlowMassTypeTag, Properties::Problem>;
+ auto freeFlowMassProblem = std::make_shared<FreeFlowMassProblem>(freeFlowMassGridGeometry, couplingManager, testCase, testCaseName);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto spatialParams = std::make_shared<typename DarcyProblem::SpatialParams>(darcyGridGeometry, testCase);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyGridGeometry, couplingManager, spatialParams, testCase, testCaseName);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[freeFlowMomentumIndex].resize(freeFlowMomentumGridGeometry->numDofs());
+ sol[freeFlowMassIndex].resize(freeFlowMassGridGeometry->numDofs());
+ sol[porousMediumIndex].resize(darcyGridGeometry->numDofs());
+
+ // the grid variables
+ using FreeFlowMomentumGridVariables = GetPropType<FreeFlowMomentumTypeTag, Properties::GridVariables>;
+ auto freeFlowMomentumGridVariables = std::make_shared<FreeFlowMomentumGridVariables>(freeFlowMomentumProblem, freeFlowMomentumGridGeometry);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ using FreeFlowMassGridVariables = GetPropType<FreeFlowMassTypeTag, Properties::GridVariables>;
+ auto freeFlowMassGridVariables = std::make_shared<FreeFlowMassGridVariables>(freeFlowMassProblem, freeFlowMassGridGeometry);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyGridGeometry);
+
+ couplingManager->init(freeFlowMomentumProblem, freeFlowMassProblem, darcyProblem,
+ std::make_tuple(freeFlowMomentumGridVariables, freeFlowMassGridVariables, darcyGridVariables),
+ sol);
+
+ freeFlowMomentumGridVariables->init(sol[freeFlowMomentumIndex]);
+ freeFlowMassGridVariables->init(sol[freeFlowMassIndex]);
+ darcyGridVariables->init(sol[porousMediumIndex]);
+
+ // initialize the vtk output module
+ VtkOutputModule freeflowVtkWriter(*freeFlowMassGridVariables, sol[freeFlowMassIndex], freeFlowMassProblem->name());
+ GetPropType<FreeFlowMassTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter); // Add model specific output fields
+ freeflowVtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<FreeFlowMassGridVariables>>());
+ NavierStokesTest::AnalyticalSolutionVectors freeFlowAnalyticSol(freeFlowMomentumProblem, freeFlowMassProblem);
+ freeflowVtkWriter.addField(freeFlowAnalyticSol.analyticalPressureSolution(), "pressureExact");
+ freeflowVtkWriter.addField(freeFlowAnalyticSol.analyticalVelocitySolution(), "velocityExact");
+ freeflowVtkWriter.write(0.0);
+
+ VtkOutputModule darcyVtkWriter(*darcyGridVariables, sol[porousMediumIndex], darcyProblem->name());
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+ darcyVtkWriter.addVelocityOutput(std::make_shared<GetPropType<DarcyTypeTag, Properties::VelocityOutput>>(*darcyGridVariables));
+ const auto darcyAnalyticalSolution = createDarcyAnalyticalSolution<double>(*darcyProblem);
+ darcyVtkWriter.addField(std::get<0>(darcyAnalyticalSolution), "pressureExact");
+ darcyVtkWriter.addField(std::get<1>(darcyAnalyticalSolution), "velocityExact");
+ darcyVtkWriter.write(0.0);
+
+ // the assembler for a stationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(freeFlowMomentumProblem, freeFlowMassProblem, darcyProblem),
+ std::make_tuple(freeFlowMomentumGridGeometry,
+ freeFlowMassGridGeometry,
+ darcyGridGeometry),
+ std::make_tuple(freeFlowMomentumGridVariables,
+ freeFlowMassGridVariables,
+ darcyGridVariables),
+ couplingManager);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ const auto dofsVS = freeFlowMomentumGridGeometry->numDofs();
+ const auto dofsPS = freeFlowMassGridGeometry->numDofs();
+ const auto dofsPD = darcyGridGeometry->numDofs();
+ std::cout << "Stokes velocity dofs: " << dofsVS << std::endl;
+ std::cout << "Stokes pressure dofs: " << dofsPS << std::endl;
+ std::cout << "Darcy pressure dofs: " << dofsPD << std::endl;
+ std::cout << "Total number of dofs: " << dofsVS + dofsPS + dofsPD << std::endl;
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // write vtk output
+ freeflowVtkWriter.write(1.0);
+ darcyVtkWriter.write(1.0);
+
+ // print L2 errors
+ printFreeFlowL2Error(freeFlowMomentumProblem, freeFlowMassProblem, partial(sol, freeFlowMomentumIndex, freeFlowMassIndex));
+ printDarcyL2Error(darcyProblem, sol[porousMediumIndex]);
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+}
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/manufactured_solution.py b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/manufactured_solution.py
new file mode 100644
index 0000000000000000000000000000000000000000..7bc744363ba71de33bd5c423e56dbab10e873221
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/manufactured_solution.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+
+"""
+This script determines the source terms needed for analytical solutions for coupled Stokes/Darcy problems.
+Given an analytical solution, it evaluates the momentum and mass balance equations and outputs source terms
+such that the balance equations are fulfilled. It furthermore calculates the Darcy velocity from a given pressure solution.
+"""
+
+import sympy as sp
+import numpy as np
+from sympy import *
+
+# case = "Rybak"
+#case = "Shiue_1"
+case = "Schneider"
+
+# divergence of a matrix
+def div(M):
+ return np.array([sp.diff(M[0,0],x) + sp.diff(M[0,1],y) , sp.diff(M[1,0],x) + sp.diff(M[1,1],y) ])
+
+# gradient of a vector
+def grad(v):
+ return np.array([[sp.diff(v[0],x), sp.diff(v[0],y)], [sp.diff(v[1],x), sp.diff(v[1],y)]])
+
+# coordinates
+x, y = sp.symbols("x y")
+
+# analytical solutions for v and p in free-flow domain (see reference given in problems)
+def analyticalSolutionStokes(case):
+ if case == "Rybak":
+ vFF = np.array([-sp.cos(sp.pi*x)*sp.sin(sp.pi*y), sp.sin(sp.pi*x)*sp.cos(sp.pi*y)])
+ pFF = 0.5*y*sp.sin(sp.pi*x)
+ elif case == "Shiue_1":
+ vFF = np.array([-1.0/sp.pi * sp.exp(y) * sp.sin(sp.pi*x), (sp.exp(y) - sp.exp(1)) * sp.cos(sp.pi*x)])
+ pFF = 2.0*sp.exp(y) * sp.cos(sp.pi*x)
+ elif case == "Shiue_2":
+ vFF = np.array([(y-1.0)*(y-1.0) + x*(y-1.0) + 3.0*x - 1.0, x*(x-1.0) - 0.5*(y-1.0)*(y-1.0) - 3.0*y + 1.0])
+ pFF = 2.0*x + y - 1.0
+ else: # Schneider
+ vFF = np.array([y, -y*sp.sin(sp.pi*x)])
+ pFF = -y*y*sp.sin(sp.pi*x)*sp.sin(sp.pi*x)
+ return [vFF, pFF]
+
+vFF, pFF = analyticalSolutionStokes(case)
+
+# individual terms of the Navier-Stokes eq.
+vvT = np.outer(vFF, vFF)
+gradV = grad(vFF)
+gradVT = grad(vFF).T
+pI = np.array([[pFF, 0], [0, pFF]])
+
+# complete momentum flux and its divergence
+if case == "Schneider":
+ momentumFlux = vvT - (gradV + gradVT) +pI
+else:
+ momentumFlux = - (gradV + gradVT) +pI # only Stokes
+divMomentumFlux = div(momentumFlux)
+
+print("\n\nStress tensor flux for case", case)
+print("-- [0][0]: ", sp.simplify(momentumFlux[0][0]))
+print("-- [1][1]: ", sp.simplify(momentumFlux[1][1]))
+print("-- [0][1] == [1][0]: ", sp.simplify(momentumFlux[0][1]))
+
+print("\n\nSource terms for case", case)
+
+# solution for source term
+print(" \nStokes:")
+print("-- Source term mass:", sp.diff(vFF[0],x) + sp.diff(vFF[1], y))
+print("-- Source term momentum x:", divMomentumFlux[0])
+print("-- Source term momentum y:", divMomentumFlux[1])
+
+# analytical solutions for p in Darcy domain (see reference given in problems)
+def analyticalSolutionDarcy(case):
+ if case == "Rybak":
+ pD = 0.5*y*y*sp.sin(sp.pi*x)
+ elif case == "Shiue_1":
+ pD = (sp.exp(y) - y*sp.exp(1)) * sp.cos(sp.pi*x)
+ elif case == "Shiue_2":
+ pD = x*(1.0-x)*(y-1.0) + pow(y-1.0, 3)/3.0 + 2.0*x + 2.0*y + 4.0
+ else: # Schneider
+ pD = (sp.exp(y+1) + 2.0 - sp.exp(2))*sp.sin(sp.pi*x)
+ return pD
+
+pD = analyticalSolutionDarcy(case)
+
+gradPdK = np.array([sp.diff(pD,x), sp.diff(pD,y)])
+vD = -gradPdK
+divVd = sp.diff(vD[0],x) + sp.diff(vD[1],y)
+
+print("\nDarcy:")
+print("-- Source term mass:", simplify(divVd))
+print("-- v x:", simplify(vD[0]))
+print("-- v_y", simplify(vD[1]))
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/params.input b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/params.input
new file mode 100644
index 0000000000000000000000000000000000000000..6d656742ed3c074b442d523783964a7dc81bf8b9
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/params.input
@@ -0,0 +1,34 @@
+[Darcy.Grid]
+UpperRight = 1 1
+Cells = 40 40
+
+[FreeFlow.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 40 40
+
+[FreeFlow.Problem]
+Name = freeFlow
+EnableInertiaTerms = false
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1
+
+[Component]
+LiquidDensity = 1.0
+LiquidKinematicViscosity = 1.0
+
+[Assembly]
+NumericDifference.BaseEpsilon = 1e-4
+
+[FreeFlow.Flux]
+UpwindWeight = 0.5
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_darcy.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_darcy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a575666dbe796d2a1c744c75bdfbd6d1362813d0
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_darcy.hh
@@ -0,0 +1,271 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The Darcy sub-problem of coupled Stokes-Darcy convergence test
+ */
+
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/common/numeqvector.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "testcase.hh"
+#include "analyticalsolutions.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The Darcy sub-problem of coupled Stokes-Darcy convergence test
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+ using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using GridView = typename GridGeometry::GridView;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+ enum class BC {
+ dirichlet, neumann, mixed
+ };
+
+public:
+ //! export the Indices
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ DarcySubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager,
+ std::shared_ptr<typename ParentType::SpatialParams> spatialParams,
+ const DarcyStokesTestCase testCase,
+ const std::string& name)
+ : ParentType(gridGeometry, spatialParams, "Darcy")
+ , couplingManager_(couplingManager)
+ , testCase_(testCase)
+ {
+ problemName_ = name + "_"
+ + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+
+ auto bc = getParamFromGroup<std::string>(this->paramGroup(), "Problem.BoundaryConditions", "Dirichlet");
+ if (bc == "Dirichlet")
+ boundaryConditions_ = BC::dirichlet;
+ else if (bc == "Neumann")
+ boundaryConditions_ = BC::neumann;
+ else if (bc == "Mixed")
+ boundaryConditions_ = BC::mixed;
+ else
+ DUNE_THROW(Dune::Exception, "Wrong BC type choose: Dirichlet, Neumann or Mixed");
+
+ std::cout << "Porous medium domain: Using " << bc << " boundary conditions" << std::endl;
+ }
+
+ const std::string& name() const
+ { return problemName_; }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param element The element
+ * \param scvf The boundary sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ BoundaryTypes values;
+
+ if (couplingManager().isCoupled(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, scvf))
+ values.setAllCouplingNeumann();
+ else
+ {
+ if (boundaryConditions_ == BC::dirichlet)
+ values.setAllDirichlet();
+ else if (boundaryConditions_ == BC::neumann)
+ values.setAllNeumann();
+ else
+ {
+ if (onLeftBoundary_(scvf.center()))
+ values.setAllNeumann();
+ else
+ values.setAllDirichlet();
+ }
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary sub-control-volume-face
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ const auto p = fullAnalyticalSolution(scvf.center())[2];
+ return PrimaryVariables(p);
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeometry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFluxVarsCache Flux variables caches for all faces in stencil
+ * \param scvf The boundary sub control volume face
+ */
+ template<class ElementVolumeVariables, class ElementFluxVarsCache>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVarsCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupled(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, scvf))
+ values[Indices::conti0EqIdx] = couplingManager().massCouplingCondition(
+ CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+
+ else
+ {
+ const auto sol = fullAnalyticalSolution(scvf.center());
+ const auto n = scvf.unitOuterNormal();
+ auto v = n; v[0] = sol[0]; v[1] = sol[1];
+ values[Indices::conti0EqIdx] = v*n;
+ }
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluates the source term for all phases within a given
+ * sub control volume.
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
+ {
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return ShiueOne::darcyRHS(globalPos);
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return ShiueTwo::darcyRHS(globalPos);
+ case DarcyStokesTestCase::Rybak:
+ return Rybak::darcyRHS(globalPos);
+ case DarcyStokesTestCase::Schneider:
+ return Schneider::darcyRHS(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ // \}
+
+ /*!
+ * \brief Evaluates the initial value for a control volume.
+ * \param element The element
+ */
+ PrimaryVariables initial(const Element &element) const
+ { return PrimaryVariables(0.0); }
+
+ /*!
+ * \brief Returns the analytical solution of the problem at a given position.
+ * \param globalPos The global position
+ * Returns vector with entries: (velocity-x | velocity-y | pressure)
+ */
+ Dune::FieldVector<Scalar, 3> fullAnalyticalSolution(const GlobalPosition& globalPos) const
+ {
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return ShiueOne::darcy(globalPos);
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return ShiueTwo::darcy(globalPos);
+ case DarcyStokesTestCase::Rybak:
+ return Rybak::darcy(globalPos);
+ case DarcyStokesTestCase::Schneider:
+ return Schneider::darcy(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ // \}
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition& globalPos) const
+ { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
+
+ static constexpr Scalar eps_ = 1e-7;
+ std::shared_ptr<CouplingManager> couplingManager_;
+ std::string problemName_;
+ DarcyStokesTestCase testCase_;
+ BC boundaryConditions_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9ef9fd3db9a051d18c4c65faf43449b7cc663c7a
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
@@ -0,0 +1,391 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The free-flow sub-problem of coupled FreeFlow/Darcy convergence test
+ */
+
+#ifndef DUMUX_FREEFLOW_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW_SUBPROBLEM_HH
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
+#include <dune/common/fmatrix.hh>
+
+#include "testcase.hh"
+#include "analyticalsolutions.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The Stokes sub-problem of coupled Stokes-Darcy convergence test
+ */
+template <class TypeTag>
+class FreeFlowSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using GridView = typename GridGeometry::GridView;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using PrimaryVariables = typename ParentType::PrimaryVariables;
+ using NumEqVector = typename ParentType::NumEqVector;
+ using BoundaryTypes = typename ParentType::BoundaryTypes;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+ enum class BC {
+ dirichlet, stress, mixed
+ };
+
+public:
+ //! export the Indices
+ using Indices = typename ModelTraits::Indices;
+
+ FreeFlowSubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager,
+ const DarcyStokesTestCase testCase,
+ const std::string& name)
+ : ParentType(gridGeometry, couplingManager, "FreeFlow")
+ , couplingManager_(couplingManager)
+ , testCase_(testCase)
+ {
+ problemName_ = name + "_"
+ + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+
+ auto bc = getParamFromGroup<std::string>(this->paramGroup(), "Problem.BoundaryConditions", "Dirichlet");
+ if (bc == "Dirichlet")
+ boundaryConditions_ = BC::dirichlet;
+ else if (bc == "Stress")
+ boundaryConditions_ = BC::stress;
+ else if (bc == "Mixed")
+ boundaryConditions_ = BC::mixed;
+ else
+ DUNE_THROW(Dune::Exception, "Wrong BC type choose: Dirichlet, Stress or Mixed");
+
+ std::cout << "Free flow domain: Using " << bc << " boundary conditions" << std::endl;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ const std::string& name() const
+ { return problemName_; }
+
+ /*!
+ * \brief Returns the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param element The finite element
+ * \param scvf The sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ if (couplingManager().isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
+ {
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setCouplingNeumann(Indices::momentumXBalanceIdx);
+ }
+ else
+ {
+ if (boundaryConditions_ == BC::dirichlet)
+ values.setAllDirichlet();
+ else if (boundaryConditions_ == BC::stress)
+ values.setAllNeumann();
+ else
+ {
+ if (onLeftBoundary_(scvf.center()))
+ values.setAllNeumann();
+ else
+ values.setAllDirichlet();
+ }
+ }
+ }
+ else
+ {
+ if (couplingManager().isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
+ values.setAllCouplingNeumann();
+ else
+ values.setAllNeumann();
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ { return analyticalSolution(globalPos); }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeometry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ template<class ElementVolumeVariables, class ElementFluxVariablesCache>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+ const auto& globalPos = scvf.ipGlobal();
+ using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
+
+ // momentum boundary conditions
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ // We need to take care here: If the integration point of an scvf lies both on the coupling interface and the
+ // domain boundary (i.e., the leftmost and rightmost points of the interface), do not evaluate the coupling or slip condition
+ // because they would require data from the boundary not available in this case (velocities for evaluating gradients,
+ // those would only be available for Dirichlet BCs). Instead, directly use the given Neumann boundary stress.
+ // TODO: Maybe couplingManager().isCoupled(...) could return false for these scvfs.
+ //
+ // | <-- left Neumann boundary
+ // |
+ // integration point --> o##### <-- scvf lying on coupling interface with integration point touching left Neumann boundary
+ //
+ if (couplingManager().isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf)
+ && scvf.ipGlobal()[0] > this->gridGeometry().bBoxMin()[0] + eps_
+ && scvf.ipGlobal()[0] < this->gridGeometry().bBoxMax()[0] - eps_)
+ {
+ values += couplingManager().momentumCouplingCondition(
+ CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+
+ values += FluxHelper::slipVelocityMomentumFlux(
+ *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache
+ );
+ }
+ else
+ {
+ const auto stress = stressTensorAtPos(globalPos);
+ const auto n = scvf.unitOuterNormal();
+ stress.mv(n, values);
+ }
+ }
+
+ // mass boundary conditions
+ else
+ {
+ if (couplingManager().isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
+ {
+ values = couplingManager().massCouplingCondition(
+ CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+ }
+ else
+ {
+ using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
+ values = FluxHelper::scalarOutflowFlux(
+ *this, element, fvGeometry, scvf, elemVolVars
+ );
+ }
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns true if the scvf lies on a porous slip boundary
+ */
+ bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ {
+ return scvf.boundary()
+ && couplingManager().isCoupled(
+ CouplingManager::freeFlowMomentumIndex,
+ CouplingManager::porousMediumIndex,
+ scvf
+ );
+ }
+
+ // \}
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Returns the sources within the domain.
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ const auto select = [&](const auto& rhs) -> NumEqVector {
+ if constexpr (ParentType::isMomentumProblem())
+ return { rhs[0], rhs[1] };
+ else
+ return { rhs[2] };
+ };
+
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return select(ShiueOne::stokesRHS(globalPos));
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return select(ShiueTwo::stokesRHS(globalPos));
+ case DarcyStokesTestCase::Rybak:
+ return select(Rybak::stokesRHS(globalPos));
+ case DarcyStokesTestCase::Schneider:
+ return select(Schneider::stokesRHS(globalPos));
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ /*!
+ * \brief Evaluates the initial value for a control volume.
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+ { return PrimaryVariables(0.0); }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter
+ for the Beavers-Joseph-Saffman boundary condition
+ */
+ auto permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ { return couplingManager().darcyPermeability(fvGeometry, scvf); }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the
+ Beavers-Joseph-Saffman boundary condition.
+ */
+ Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ { return couplingManager().problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal()); }
+
+ /*!
+ * \brief Returns the analytical solution of the problem at a given position.
+ * \param globalPos The global position
+ * Returns vector with entries: (velocity-x | velocity-y | pressure)
+ */
+ Dune::FieldVector<Scalar, 3> fullAnalyticalSolution(const GlobalPosition& globalPos) const
+ {
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return ShiueOne::stokes(globalPos);
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return ShiueTwo::stokes(globalPos);
+ case DarcyStokesTestCase::Rybak:
+ return Rybak::stokes(globalPos);
+ case DarcyStokesTestCase::Schneider:
+ return Schneider::stokes(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ /*!
+ * \brief Returns the analytical solution of the problem at a given position.
+ * \param globalPos The global position
+ * Returns analytical solution depending on the type of sub-problem
+ */
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
+ {
+ using namespace Solution::DarcyStokes;
+ const auto sol = fullAnalyticalSolution(globalPos);
+ if constexpr (ParentType::isMomentumProblem())
+ return { sol[0], sol[1] };
+ else
+ return { sol[2] };
+ }
+
+ // \}
+
+private:
+ /*!
+ * \brief Returns the stress tensor within the domain.
+ * \param globalPos The global position
+ */
+ Dune::FieldMatrix<double, 2, 2> stressTensorAtPos(const GlobalPosition &globalPos) const
+ {
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return ShiueOne::stokesStress(globalPos);
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return ShiueTwo::stokesStress(globalPos);
+ case DarcyStokesTestCase::Rybak:
+ return Rybak::stokesStress(globalPos);
+ case DarcyStokesTestCase::Schneider:
+ return Schneider::stokesStress(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ bool onLeftBoundary_(const GlobalPosition& globalPos) const
+ { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
+
+ static constexpr Scalar eps_ = 1e-7;
+ std::string problemName_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+ DarcyStokesTestCase testCase_;
+ BC boundaryConditions_;
+};
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/properties.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..70c62ded2fe14d3e218baeaee6951bc20f4f72a7
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/properties.hh
@@ -0,0 +1,129 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The properties for a simple Darcy test (cell-centered finite volume method)
+ */
+#ifndef DUMUX_DARCYSTOKES_PROPERTIES_HH
+#define DUMUX_DARCYSTOKES_PROPERTIES_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/constant.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
+#include <dumux/freeflow/navierstokes/momentum/model.hh>
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/discretization/fcstaggered.hh>
+
+#include <dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh>
+#include <dumux/multidomain/traits.hh>
+
+#include "spatialparams.hh"
+#include "problem_darcy.hh"
+#include "problem_freeflow.hh"
+
+namespace Dumux::Properties {
+
+// Create new type tags
+namespace TTag {
+struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOneP>
+{ using type = DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOneP>
+{ using type = Dune::YaspGrid<2>; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOneP>
+{
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = ConvergenceTestSpatialParams<GridGeometry, Scalar>;
+};
+
+// Create new type tags
+namespace TTag {
+struct FreeFlowOneP {};
+struct FreeFlowOnePMass { using InheritsFrom = std::tuple<FreeFlowOneP, NavierStokesMassOneP, CCTpfaModel>; };
+struct FreeFlowOnePMomentum { using InheritsFrom = std::tuple<FreeFlowOneP, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
+} // end namespace TTag
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::FreeFlowOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::FreeFlowOneP>
+{
+ using Coords = Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2>;
+ using type = Dune::YaspGrid<2, Coords>;
+};
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::FreeFlowOneP>
+{ using type = FreeFlowSubProblem<TypeTag>; };
+
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOneP>
+{
+ using Traits = MultiDomainTraits<TTag::FreeFlowOnePMomentum, TTag::FreeFlowOnePMass, TTag::DarcyOneP>;
+ using type = FreeFlowPorousMediumCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::FreeFlowOneP>
+{
+ using Traits = MultiDomainTraits<TTag::FreeFlowOnePMomentum, TTag::FreeFlowOnePMass, TTag::DarcyOneP>;
+ using type = FreeFlowPorousMediumCouplingManager<Traits>;
+};
+
+} // end namespace Dumux::Properties
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/spatialparams.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..82eb797c14cf2134c53b701a9524a32a39056cd3
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/spatialparams.hh
@@ -0,0 +1,129 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The spatial parameters class for the test problem using the 1p cc model.
+ */
+
+#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH
+#define DUMUX_1P_TEST_SPATIALPARAMS_HH
+
+#include <dumux/material/spatialparams/fv1p.hh>
+#include "testcase.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The spatial parameters class for the test problem using the
+ * 1p cc model.
+ */
+template<class GridGeometry, class Scalar>
+class ConvergenceTestSpatialParams
+: public FVSpatialParamsOneP<GridGeometry, Scalar,
+ ConvergenceTestSpatialParams<GridGeometry, Scalar>>
+{
+ using GridView = typename GridGeometry::GridView;
+ using ParentType = FVSpatialParamsOneP<
+ GridGeometry, Scalar, ConvergenceTestSpatialParams<GridGeometry, Scalar>
+ >;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+
+
+public:
+ // export permeability type
+ using PermeabilityType = DimWorldMatrix;
+
+ ConvergenceTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry,
+ const DarcyStokesTestCase testCase)
+ : ParentType(gridGeometry)
+ , testCase_(testCase)
+ {
+ alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ *
+ * \param element The element
+ * \param scv The sub control volume
+ * \param elemSol The element solution vector
+ * \return the intrinsic permeability
+ */
+ template<class SubControlVolume, class ElementSolution>
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ PermeabilityType K(0.0);
+
+ if (testCase_ == DarcyStokesTestCase::Schneider)
+ {
+ using std::cos;
+ using std::sin;
+ using std::exp;
+
+ static constexpr Scalar c = 0.0;
+ static constexpr Scalar omega = M_PI;
+
+ const Scalar x = scv.center()[0];
+ K[0][0] = 1.0;
+ K[0][1] = -c/(2*omega) * sin(omega*x);
+ K[1][0] = K[0][1];
+ K[1][1] = exp(-2)*(1 + c*cos(omega*x));
+ }
+ else
+ {
+ const static Scalar permeability = getParam<Scalar>("Darcy.SpatialParams.Permeability", 1.0);
+ K[0][0] = permeability;
+ K[1][1] = permeability;
+ }
+
+ return K;
+ }
+
+ /*! \brief Defines the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return 0.4; }
+
+ /*! \brief Defines the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+private:
+ DarcyStokesTestCase testCase_;
+ Scalar permeability_;
+ Scalar alphaBJ_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/testcase.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/testcase.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9496ffe6ae773a194afc28d6ef1d6ca0e070d4b2
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/testcase.hh
@@ -0,0 +1,37 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The different test cases
+ */
+
+#ifndef DUMUX_CONVERGENCE_TEST_TESTCASE_HH
+#define DUMUX_CONVERGENCE_TEST_TESTCASE_HH
+
+namespace Dumux {
+
+enum class DarcyStokesTestCase
+{
+ ShiueExampleOne, ShiueExampleTwo, Rybak, Schneider
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/main.cc b/test/multidomain/boundary/freeflowporousmedium/1p_1p/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..17abbaac06ee5cc9eb75793adddd40c1e1433c10
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/main.cc
@@ -0,0 +1,213 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief A test problem for the coupled Stokes/Darcy problem (1p).
+ */
+
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+#include <dumux/freeflow/navierstokes/velocityoutput.hh>
+
+#include "properties.hh"
+
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // Define the sub problem type tags
+ using FreeFlowMomentumTypeTag = Properties::TTag::FreeFlowOnePMomentum;
+ using FreeFlowMassTypeTag = Properties::TTag::FreeFlowOnePMass;
+ using DarcyTypeTag = Properties::TTag::DarcyOneP;
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using FreeFlowGridManager = Dumux::GridManager<GetPropType<FreeFlowMomentumTypeTag, Properties::Grid>>;
+ FreeFlowGridManager freeFlowGridManager;
+ freeFlowGridManager.init("FreeFlow"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& freeFlowGridView = freeFlowGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FreeFlowMomentumGridGeometry = GetPropType<FreeFlowMomentumTypeTag, Properties::GridGeometry>;
+ auto freeFlowMomentumGridGeometry = std::make_shared<FreeFlowMomentumGridGeometry>(freeFlowGridView);
+ using FreeFlowMassGridGeometry = GetPropType<FreeFlowMassTypeTag, Properties::GridGeometry>;
+ auto freeFlowMassGridGeometry = std::make_shared<FreeFlowMassGridGeometry>(freeFlowGridView);
+ using DarcyGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
+ auto darcyGridGeometry = std::make_shared<DarcyGridGeometry>(darcyGridView);
+
+ using Traits = MultiDomainTraits<FreeFlowMomentumTypeTag, FreeFlowMassTypeTag, DarcyTypeTag>;
+ using CouplingManager = FreeFlowPorousMediumCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>();
+
+ // the indices
+ constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex;
+ constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex;
+ constexpr auto porousMediumIndex = CouplingManager::porousMediumIndex;
+
+ // the problem (initial and boundary conditions)
+ using FreeFlowMomentumProblem = GetPropType<FreeFlowMomentumTypeTag, Properties::Problem>;
+ auto freeFlowMomentumProblem = std::make_shared<FreeFlowMomentumProblem>(freeFlowMomentumGridGeometry, couplingManager);
+ using FreeFlowMassProblem = GetPropType<FreeFlowMassTypeTag, Properties::Problem>;
+ auto freeFlowMassProblem = std::make_shared<FreeFlowMassProblem>(freeFlowMassGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyGridGeometry, couplingManager);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[freeFlowMomentumIndex].resize(freeFlowMomentumGridGeometry->numDofs());
+ sol[freeFlowMassIndex].resize(freeFlowMassGridGeometry->numDofs());
+ sol[porousMediumIndex].resize(darcyGridGeometry->numDofs());
+
+ // the grid variables
+ using FreeFlowMomentumGridVariables = GetPropType<FreeFlowMomentumTypeTag, Properties::GridVariables>;
+ auto freeFlowMomentumGridVariables = std::make_shared<FreeFlowMomentumGridVariables>(freeFlowMomentumProblem, freeFlowMomentumGridGeometry);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ using FreeFlowMassGridVariables = GetPropType<FreeFlowMassTypeTag, Properties::GridVariables>;
+ auto freeFlowMassGridVariables = std::make_shared<FreeFlowMassGridVariables>(freeFlowMassProblem, freeFlowMassGridGeometry);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyGridGeometry);
+
+ couplingManager->init(freeFlowMomentumProblem, freeFlowMassProblem, darcyProblem,
+ std::make_tuple(freeFlowMomentumGridVariables, freeFlowMassGridVariables, darcyGridVariables),
+ sol);
+
+ freeFlowMomentumGridVariables->init(sol[freeFlowMomentumIndex]);
+ freeFlowMassGridVariables->init(sol[freeFlowMassIndex]);
+ darcyGridVariables->init(sol[porousMediumIndex]);
+
+ // initialize the vtk output module
+ VtkOutputModule vtkWriterFF(*freeFlowMassGridVariables, sol[freeFlowMassIndex], freeFlowMassProblem->name());
+ GetPropType<FreeFlowMassTypeTag, Properties::IOFields>::initOutputModule(vtkWriterFF); // Add model specific output fields
+ vtkWriterFF.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<FreeFlowMassGridVariables>>());
+
+ const bool addAnalyticSolutionForV = getParamFromGroup<bool>(freeFlowMassProblem->paramGroup(), "Problem.AddAnalyticSolutionForV", false);
+
+ if (!freeFlowMassProblem->verticalFlow() && addAnalyticSolutionForV)
+ {
+ using Scalar = typename Traits::Scalar;
+ using std::sqrt;
+ std::vector<Scalar> analyticalVelocityX(freeFlowMassGridGeometry->gridView().size(0));
+ const Scalar dPdX = -freeFlowMomentumProblem->pressureDifference() / (freeFlowMassGridGeometry->bBoxMax()[0] - freeFlowMassGridGeometry->bBoxMin()[0]);
+ static const Scalar mu = getParam<Scalar>("Component.LiquidKinematicViscosity") * getParam<Scalar>("Component.LiquidDensity");
+ static const Scalar alpha = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ static const Scalar K = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ static const Scalar sqrtK = sqrt(K);
+ const Scalar sigma = (freeFlowMassGridGeometry->bBoxMax()[1] - freeFlowMassGridGeometry->bBoxMin()[1])/sqrtK;
+
+ const Scalar uB = -K/(2.0*mu) * ((sigma*sigma + 2.0*alpha*sigma) / (1.0 + alpha*sigma)) * dPdX;
+
+ for (const auto& element : elements(freeFlowMassGridGeometry->gridView()))
+ {
+ const auto eIdx = freeFlowMassGridGeometry->gridView().indexSet().index(element);
+ const Scalar y = element.geometry().center()[1] - freeFlowMassGridGeometry->bBoxMin()[1];
+
+ const Scalar u = uB*(1.0 + alpha/sqrtK*y) + 1.0/(2.0*mu) * (y*y + 2*alpha*y*sqrtK) * dPdX;
+ analyticalVelocityX[eIdx] = u;
+ }
+
+ vtkWriterFF.addField(analyticalVelocityX, "analyticalV_x");
+ }
+
+ VtkOutputModule pmVtkWriter(*darcyGridVariables, sol[porousMediumIndex], darcyProblem->name());
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(pmVtkWriter);
+ pmVtkWriter.addVelocityOutput(std::make_shared<GetPropType<DarcyTypeTag, Properties::VelocityOutput>>(*darcyGridVariables));
+
+ // the assembler for a stationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(freeFlowMomentumProblem, freeFlowMassProblem, darcyProblem),
+ std::make_tuple(freeFlowMomentumGridGeometry,
+ freeFlowMassGridGeometry,
+ darcyGridGeometry),
+ std::make_tuple(freeFlowMomentumGridVariables,
+ freeFlowMassGridVariables,
+ darcyGridVariables),
+ couplingManager);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ const auto dofsVS = freeFlowMomentumGridGeometry->numDofs();
+ const auto dofsPS = freeFlowMassGridGeometry->numDofs();
+ const auto dofsPD = darcyGridGeometry->numDofs();
+ std::cout << "Stokes velocity dofs: " << dofsVS << std::endl;
+ std::cout << "Stokes pressure dofs: " << dofsPS << std::endl;
+ std::cout << "Darcy pressure dofs: " << dofsPD << std::endl;
+ std::cout << "Total number of dofs: " << dofsVS + dofsPS + dofsPD << std::endl;
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // write vtk output
+ vtkWriterFF.write(1.0);
+ pmVtkWriter.write(1.0);
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_horizontalflow.input b/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_horizontalflow.input
new file mode 100644
index 0000000000000000000000000000000000000000..3da8b5649dfa203a9fc834ded0109be31b547c6c
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_horizontalflow.input
@@ -0,0 +1,33 @@
+[Darcy.Grid]
+UpperRight = 1 1
+Cells = 20 20
+
+[FreeFlow.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 20 20
+
+[FreeFlow.Problem]
+Name = freeflow
+PressureDifference = 1e-9
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+Permeability = 1e-6 # m^2
+AlphaBeaversJoseph = 1.0
+
+[Component]
+LiquidKinematicViscosity = 1e-6
+LiquidDensity = 1e3
+
+[Vtk]
+OutputName = test_md_boundary_freeflow1p_darcy1p_horizontal
+
+[Problem]
+EnableGravity = false
+EnableInertiaTerms = false
+
+[Vtk]
+AddVelocity = 1
\ No newline at end of file
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_verticalflow.input b/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_verticalflow.input
new file mode 100644
index 0000000000000000000000000000000000000000..54358b805bc91bb253dab9fdf5e0e3b6d783cd6a
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_verticalflow.input
@@ -0,0 +1,37 @@
+[Darcy.Grid]
+UpperRight = 2.0 2.0
+Cells = 20 20
+
+[FreeFlow.Grid]
+LowerLeft = 0.0 2.0
+UpperRight = 2.0 4.0
+Cells = 20 20
+
+[FreeFlow.Problem]
+Name = freeflow
+Velocity = -1e-6
+
+[Darcy.Problem]
+Name = darcy
+Pressure = 0.0
+
+[Darcy.SpatialParams]
+Permeability = 1e-10 # m^2
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+EnableGravity = false
+EnableInertiaTerms = false
+
+[Component]
+LiquidKinematicViscosity = 1e-6
+LiquidDensity = 1e3
+
+[Vtk]
+OutputName = test_md_boundary_freeflow1p_darcy1p_vertical
+
+[Vtk]
+AddVelocity = 1
+
+[Assembly.NumericDifference]
+BaseEpsilon = 1e-6
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_darcy.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_darcy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..28490814a2bca3dff391e65a0d61298f0dc7323c
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_darcy.hh
@@ -0,0 +1,183 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief A simple Darcy test problem (cell-centered finite volume method).
+ */
+
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/common/numeqvector.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+
+namespace Dumux {
+
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+ using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+public:
+ DarcySubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(gridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+ problemName_ = getParam<std::string>("Vtk.OutputName")
+ + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+
+ // determine whether to simulate a vertical or horizontal flow configuration
+ verticalFlow_ = problemName_.find("vertical") != std::string::npos;
+ }
+
+ /*!
+ * \brief The problem name.
+ */
+ const std::string& name() const
+ { return problemName_; }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param element The element
+ * \param scvf The boundary sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+
+ if (couplingManager().isCoupled(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, scvf))
+ values.setAllCouplingNeumann();
+
+ if (verticalFlow_)
+ {
+ if (onLowerBoundary_(scvf.center()))
+ values.setAllDirichlet();
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary sub-control-volume-face
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ { return initial(element); }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeometry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFluxVarsCache Flux variables caches for all faces in stencil
+ * \param scvf The boundary sub control volume face
+ */
+ template<class ElementVolumeVariables, class ElementFluxVarsCache>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVarsCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupled(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, scvf))
+ values[Indices::conti0EqIdx] = couplingManager().massCouplingCondition(
+ CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Evaluates the initial value for a control volume.
+ * \param element The element
+ */
+ PrimaryVariables initial(const Element &element) const
+ { return PrimaryVariables(0.0); }
+
+ // \}
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+private:
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
+
+ Scalar eps_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+ std::string problemName_;
+ bool verticalFlow_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d184aa12247975602496061572dea14df8dbb729
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh
@@ -0,0 +1,307 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesTests
+ * \brief Test for the staggered grid Navier-Stokes model with analytical solution (Kovasznay 1948, \cite Kovasznay1948)
+ */
+
+#ifndef DUMUX_KOVASZNAY_TEST_PROBLEM_NEW_HH
+#define DUMUX_KOVASZNAY_TEST_PROBLEM_NEW_HH
+
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesTests
+ * \brief Test problem for the staggered grid (Kovasznay 1948, \cite Kovasznay1948)
+ *
+ * A two-dimensional Navier-Stokes flow with a periodicity in one direction
+ * is considered. The set-up represents a wake behind a two-dimensional grid
+ * and is chosen in a way such that an exact solution is available.
+ */
+template <class TypeTag>
+class FreeFlowOnePTestProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using BoundaryTypes = typename ParentType::BoundaryTypes;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using NumEqVector = typename ParentType::NumEqVector;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using PrimaryVariables = typename ParentType::PrimaryVariables;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+
+ static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
+ using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+ // static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
+
+public:
+ FreeFlowOnePTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(gridGeometry, couplingManager, "FreeFlow")
+ , couplingManager_(couplingManager)
+ {
+ problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+ deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference", 0.0);
+
+ // determine whether to simulate a vertical or horizontal flow configuration
+ verticalFlow_ = problemName_.find("vertical") != std::string::npos;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ const std::string& name() const
+ { return problemName_; }
+
+ /*!
+ * \brief Returns the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 298.0; }
+
+ // \}
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param element The finite element
+ * \param scvf The sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+ const auto& globalPos = scvf.center(); //avoid ambiguities at corners
+
+ if (verticalFlow_)
+ {
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
+ {
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setCouplingNeumann(Indices::momentumXBalanceIdx);
+ }
+ else
+ values.setAllDirichlet();
+ }
+ else
+ {
+ if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
+ values.setAllCouplingNeumann();
+ else
+ values.setAllNeumann();
+ }
+ }
+ else // horizontal flow
+ {
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
+ values.setAllNeumann();
+ else if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
+ {
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setCouplingNeumann(Indices::momentumXBalanceIdx);
+ }
+ else
+ values.setAllDirichlet();
+ }
+ else
+ {
+ if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
+ values.setAllCouplingNeumann();
+ else
+ values.setAllNeumann();
+ }
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns Dirichlet boundary values at a given position.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ if (verticalFlow_)
+ {
+ static const Scalar inletVelocity = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
+ values[Indices::velocityYIdx] = inletVelocity * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]);
+ }
+ }
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeometry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ template<class ElementVolumeVariables, class ElementFluxVariablesCache>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+ const auto& globalPos = scvf.ipGlobal();
+ using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
+
+ if constexpr (ParentType::isMomentumProblem())
+ {
+
+ if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
+ {
+ values += couplingManager_->momentumCouplingCondition(
+ CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+
+ values += FluxHelper::slipVelocityMomentumFlux(
+ *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache
+ );
+ }
+ else
+ {
+ const Scalar pressure = onLeftBoundary_(globalPos) ? deltaP_ : 0.0;
+ values = FluxHelper::fixedPressureMomentumFlux(
+ *this, fvGeometry, scvf, elemVolVars,
+ elemFluxVarsCache, pressure, true /*zeroNormalVelocityGradient*/
+ );
+ }
+ }
+ else
+ {
+ if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
+ {
+ values = couplingManager_->massCouplingCondition(
+ CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+ }
+ else
+ {
+ using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
+ values = FluxHelper::scalarOutflowFlux(
+ *this, element, fvGeometry, scvf, elemVolVars
+ );
+ }
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns true if the scvf lies on a porous slip boundary
+ */
+ bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ {
+ assert(scvf.isFrontal());
+ return scvf.boundary() && couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf);
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ { return couplingManager_->darcyPermeability(fvGeometry, scvf); }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ { return couplingManager_->problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal()); }
+
+ /*!
+ * \brief Returns the pressure difference between inlet and outlet for the horizontal flow scenario
+ */
+ Scalar pressureDifference() const
+ {
+ assert(!verticalFlow_);
+ return deltaP_;
+ }
+
+ /*!
+ * \brief Returns true if a vertical flow scenario is considered
+ */
+ bool verticalFlow() const
+ { return verticalFlow_; }
+
+ // \}
+
+private:
+
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
+
+ std::string problemName_;
+ static constexpr Scalar eps_ = 1e-6;
+ Scalar deltaP_;
+ bool verticalFlow_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/properties.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4bf95acba7c79323df243bd349f4f2677dd736f0
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/properties.hh
@@ -0,0 +1,125 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The properties for a simple Darcy test (cell-centered finite volume method)
+ */
+#ifndef DUMUX_DARCYSTOKES_PROPERTIES_HH
+#define DUMUX_DARCYSTOKES_PROPERTIES_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/constant.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
+#include <dumux/freeflow/navierstokes/momentum/model.hh>
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/discretization/fcstaggered.hh>
+
+#include <dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh>
+
+#include "spatialparams.hh"
+#include "problem_darcy.hh"
+#include "problem_freeflow.hh"
+
+namespace Dumux::Properties {
+
+// Create new type tags
+namespace TTag {
+struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOneP> { using type = Dune::YaspGrid<2>; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOneP>
+{
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = OnePSpatialParams<GridGeometry, Scalar>;
+};
+
+// Create new type tags
+namespace TTag {
+struct FreeFlowOneP {};
+struct FreeFlowOnePMass { using InheritsFrom = std::tuple<FreeFlowOneP, NavierStokesMassOneP, CCTpfaModel>; };
+struct FreeFlowOnePMomentum { using InheritsFrom = std::tuple<FreeFlowOneP, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
+} // end namespace TTag
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::FreeFlowOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::FreeFlowOneP>
+{
+ using Coords = Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2>;
+ using type = Dune::YaspGrid<2, Coords>;
+};
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::FreeFlowOneP> { using type = Dumux::FreeFlowOnePTestProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOneP>
+{
+ using Traits = MultiDomainTraits<Properties::TTag::FreeFlowOnePMomentum, Properties::TTag::FreeFlowOnePMass, TTag::DarcyOneP>;
+ using type = Dumux::FreeFlowPorousMediumCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::FreeFlowOneP>
+{
+ using Traits = MultiDomainTraits<Properties::TTag::FreeFlowOnePMomentum, Properties::TTag::FreeFlowOnePMass, TTag::DarcyOneP>;
+ using type = Dumux::FreeFlowPorousMediumCouplingManager<Traits>;
+};
+
+} // end namespace Dumux::Properties
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/spatialparams.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..40ef1710e2229ddb9fc273f9ea636a80384b44e3
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/spatialparams.hh
@@ -0,0 +1,92 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The spatial parameters class for the test problem using the 1p cc model.
+ */
+
+#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH
+#define DUMUX_1P_TEST_SPATIALPARAMS_HH
+
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The spatial parameters class for the test problem using the
+ * 1p cc model.
+ */
+template<class GridGeometry, class Scalar>
+class OnePSpatialParams
+: public FVSpatialParamsOneP<GridGeometry, Scalar,
+ OnePSpatialParams<GridGeometry, Scalar>>
+{
+ using GridView = typename GridGeometry::GridView;
+ using ParentType = FVSpatialParamsOneP<
+ GridGeometry, Scalar, OnePSpatialParams<GridGeometry, Scalar>
+ >;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ // export permeability type
+ using PermeabilityType = Scalar;
+
+ OnePSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ *
+ * \param globalPos The global position
+ * \return the intrinsic permeability
+ */
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+ { return permeability_; }
+
+ /*! \brief Defines the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return 0.4; }
+
+ /*! \brief Defines the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+
+private:
+ Scalar permeability_;
+ Scalar alphaBJ_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/CMakeLists.txt b/test/multidomain/boundary/freeflowporousmedium/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f62f49df141e549022cf0cb7aac8d6f5c4d7f83f
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(1p_1p)