diff --git a/exercises/CMakeLists.txt b/exercises/CMakeLists.txt
index 55370342d6b61fa363689ea63e83d70e14ce11f8..82031e6c316e1cec7f6f58a2021df810acb1f14f 100644
--- a/exercises/CMakeLists.txt
+++ b/exercises/CMakeLists.txt
@@ -2,6 +2,7 @@
add_custom_target(test_exercises)
add_subdirectory(exercise-basic)
+add_subdirectory(exercise-coupling-ff-pm)
add_subdirectory(exercise-runtimeparams)
add_subdirectory(exercise-grids)
add_subdirectory(exercise-properties)
diff --git a/exercises/exercise-coupling-ff-pm/1pspatialparams.hh b/exercises/exercise-coupling-ff-pm/1pspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..64887e33307e8a45908dabc131e919a99f95dfa1
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/1pspatialparams.hh
@@ -0,0 +1,98 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup OnePTests
+ * \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 OnePModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief The spatial parameters class for the test problem using the
+ * 1p cc model
+ */
+template<class TypeTag>
+class OnePSpatialParams
+: public FVSpatialParamsOneP<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar),
+ OnePSpatialParams<TypeTag>>
+{
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<TypeTag>>;
+
+ 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 FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
+ 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 Define the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return porosity_; }
+
+ /*! \brief Define the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+
+private:
+ Scalar permeability_;
+ Scalar porosity_;
+ Scalar alphaBJ_;
+};
+
+} // end namespace
+
+#endif
diff --git a/exercises/exercise-coupling-ff-pm/2pspatialparams.hh b/exercises/exercise-coupling-ff-pm/2pspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1830f3cbf2381280ca1213b8a07a17890846203c
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/2pspatialparams.hh
@@ -0,0 +1,139 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPTests
+ * \brief The spatial parameters class for the test problem using the 2p cc model
+ */
+#ifndef DUMUX_TWOPHASE_SPATIAL_PARAMS_HH
+#define DUMUX_TWOPHASE_SPATIAL_PARAMS_HH
+
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief The spatial parameters class for the test problem using the 2p cc model
+ */
+template<class TypeTag>
+class TwoPSpatialParams
+: public FVSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar),
+ TwoPSpatialParams<TypeTag>>
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using ParentType = FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<TypeTag>>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
+
+public:
+ using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
+ using MaterialLawParams = typename MaterialLaw::Params;
+ using PermeabilityType = Scalar;
+
+ TwoPSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
+ alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+
+ // residual saturations
+ params_.setSwr(getParam<Scalar>("Darcy.SpatialParams.Swr"));
+ params_.setSnr(getParam<Scalar>("Darcy.SpatialParams.Snr"));
+ // parameters for the vanGenuchten law
+ params_.setVgAlpha(getParam<Scalar>("Darcy.SpatialParams.VgAlpha"));
+ params_.setVgn(getParam<Scalar>("Darcy.SpatialParams.VgN"));
+ params_.setPcLowSw(params_.swr()*5.0);
+ params_.setPcHighSw(1.0-params_.snr()*5.0);
+ }
+
+ /*!
+ * \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 Define the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return porosity_; }
+
+ /*! \brief Define the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+ /*!
+ * \brief Returns the parameter object for the Brooks-Corey material law.
+ * In this test, we use element-wise distributed material parameters.
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \return the material parameters object
+ */
+ template<class ElementSolutionVector>
+ const MaterialLawParams& materialLawParams(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return params_; }
+
+ /*!
+ * \brief Function for defining which phase is to be considered as the wetting phase.
+ *
+ * \return the wetting phase index
+ * \param globalPos The global position
+ */
+ template<class FluidSystem>
+ int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+ { return FluidSystem::phase0Idx; }
+
+private:
+ Scalar permeability_;
+ Scalar porosity_;
+ Scalar alphaBJ_;
+ MaterialLawParams params_;
+ static constexpr Scalar eps_ = 1.0e-7;
+};
+
+} // end namespace
+
+#endif
diff --git a/exercises/exercise-coupling-ff-pm/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bae4c263a6bf7ae6a06d90af66a2e45d224f9b2e
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_subdirectory(interface)
+add_subdirectory(models)
+add_subdirectory(turbulence)
diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..6b0b300811c751b29fd28c65d916f334dac1a66b
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -0,0 +1,388 @@
+# Exercise Coupling free flow/porous medium flow
+
+The aim of this exercise is to get familiar with setting up coupled free flow/porous medium flow problems.
+
+## Problem set-up
+
+The model domain consists of two non-overlapping subdomains.
+Free flow is modeled in the upper subdomain, while the lower subdomain models a flow in a porous medium. Both single-phase flow and two-phase flow will be considered in the porous domain.
+
+
+### 0. Getting familiar with the code
+
+* Navigate to the directory `exercises/exercise-coupling-ff-pm`
+
+There are three sub folders: `interface`, `models` and `turbulence`.
+
+The problem-related files for this exercise are:
+* Three __main files__ for the three sub-tasks :`ex_interface_coupling_ff-pm.cc`, `ex_models_coupling_ff-pm.cc`, `ex_turbulence_coupling_ff-pm.cc`,
+* Three __free flow problem files__: `ex_interface_ffproblem.hh`, `ex_models_ffproblem.hh`, `ex_turbulence__ffproblem.hh`
+* Three __porous medium flow problem files__: `ex_interface_pmproblem.hh`, `ex_models_pmproblem.hh`, `ex_turbulence_pmproblem.hh`
+* The __input files__: `ex_interface_coupling_ff-pm.input`, `ex_models_coupling_ff-pm.input`, `ex_turbulence_coupling_ff-pm.input`,
+* The __spatial parameters files__: `1pspatialparams.hh`, `2pspatialparams.hh`
+
+
+In the main file, `TypeTags` for both submodels are defined.
+The same applies for types such as `GridManager`, `FVGridGeometry`, `Problem`, etc...
+Since we use a monolithic coupling scheme, there is only one `Assembler` and one `NewtonSolver`.
+
+The problem files very much look like "regular", uncoupled ones with the exception that they hold a pointer to the `CouplingManager` which allows to evaluate the coupling conditions and to exchange information between the coupled models.
+The coupling conditions are realized technically in terms of boundary condition. For instance, in lines 171 and 172
+in `ex_interface_ffproblem.hh`, `couplingNeumann` boundary conditions are set, which means that the free flow models evaluates the
+mass and momentum fluxes coming from the porous domain and uses these values as boundary conditions at the interface.
+
+Note the certain checks are performed when combining different models, e.g., the fluid system has to be the same for both domains
+and the sub-control-volume faces at the interface have to match.
+
+
+We will use a staggered grid for the free flow and a cell-centered finite volume method for the porous medium.
+Keep in mind that the staggered grid implementation distinguishes between face variables (velocity components) and
+cell center variables (all other variables), therefore in some cases either the `stokesCellCenterIdx`
+or the `stokesFaceIdx` is used respectively, while for the porous medium all variables can be accessed with `darcyIdx`.
+
+__Task__:
+Take a closer look at the Stokes/Darcy coupling files before moving to the next part of the exercise:
+
+
+### 1. Changing the interface
+
+In this part of the exercise, a simple coupled system consisting of a one-phase (1p) free flow and a one-phase flow in a porous medium is set up. Both subproblems have no-flow boundaries at the sides.
+Currently, a velocity profile is set on the upper free flow boundary, which leads to a vertical flow into the porous medium:
+
+
+
+* We will first change the flow direction such that the free flow is parallel to the porous medium.
+* Afterwards, the Beavers-Joseph-Saffman condition will be used as an interface condition for the tangential momentum transfer.
+* Last, we change the flat interface between the two domains to a wave-shaped one.
+
+__Task A: Change the flow direction__
+
+Open the file `ex_interface_ffproblem.hh` and navigate to line 148, where the types of boundary condition are set.
+Instead of applying a fixed velocity profile at the top of the domain, we want to use fixed pressure boundary conditions
+at the left and right side of the free flow domain, while the top represents an impermeable wall.
+
+Set a Dirichlet boundary condition for the pressure at the left and right side of the domain:
+``` cpp
+if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
+ values.setDirichlet(Indices::pressureIdx);
+
+```
+
+Set a Dirichlet boundary condition for the velocities at the top:
+``` cpp
+if(onUpperBoundary_(globalPos))
+{
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+}
+```
+
+Keep the coupling boundary condition:
+``` cpp
+if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+{
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+}
+```
+
+Having changed the types of boundary conditions, we must now assign the correct values for them.
+
+Set a no-slip, no-flow condition for the velocity at the top:
+``` cpp
+values[Indices::velocityXIdx] = 0.0;
+values[Indices::velocityYIdx] = 0.0;
+```
+Apply a fixed pressure difference between the inlet and outlet, e.g.:
+``` cpp
+if(onLeftBoundary_(globalPos))
+ values[Indices::pressureIdx] = deltaP_;
+if(onRightBoundary_(globalPos))
+ values[Indices::pressureIdx] = 0.0;
+```
+
+For changing the flow direction, the boundary conditions for the porous medium have to be changed as well.
+
+Use Neumann no-flow boundaries everywhere, keep the coupling conditions.
+``` cpp
+values.setAllNeumann();
+
+if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+```
+
+This should make the flow go from left to right.
+
+__Task B: Include slip-condition__
+
+However, we are still missing one important feature:
+at the moment, the velocity component tangential to the interface gets a no-slip condition.
+In the next step we want to implement the Beavers-Joseph-Saffman slip condition at the interface:
+
+$`\frac{\partial v_x}{\partial y} = \frac{\alpha}{\sqrt K} (v_x - q_{pm})\quad`$ at $`\quad y=0`$
+
+with $`\quad q_{pm}=0`$.
+
+To include this, just set a Beavers-Joseph-Saffman (BJS) boundary condition for the respective momentum balance:
+``` cpp
+values.setBJS(Indices::momentumXBalanceIdx);
+```
+
+at the position where the coupling boundary conditions are set in `ex_interface_ffproblem.hh`.
+
+To check if the simulation behaves as expected, we can compare the velocity profile $`v_x(y)`$ with the analytical solution provided by [Beavers and Joseph (1967)](https://doi.org/10.1017/S0022112067001375).
+For doing so, we uncomment line 212 in `ex_interface_coupling_ff-pm.cc`.
+```cpp
+stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+```
+
+After re-compiling and re-running the executable, we should be able to see also
+the analytical solution of $`v_x`$ on the free flow domain. Play around with the grid resolution to see how that affects the velocity profile.
+
+__Task C: Cange shape of interface__
+
+Now we want to include a non-flat interface between the two domains. We use `dune-subgrid` to construct
+two grids for the two domains from one common host grid. Comment out lines 93-103 in `ex_interface_coupling_ff-pm.cc` and comment lines 112-147 in the same file. This will instantiate a host grid and define two helper lambda functions that are used to choose elements from to host grid for the respective sub grid. In the given case,
+the domain is split in two haves, separated by a sinusoidal interface.
+
+```cpp
+auto elementSelectorStokes = [&](const auto& element)
+{
+ double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ return element.geometry().center()[1] > interface;
+};
+
+auto elementSelectorDarcy = [&](const auto& element)
+{
+ double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ return element.geometry().center()[1] < interface;
+};
+```
+
+Make sure, that you have uncommented the lines including the gridcreators in both problem files
+```cpp
+#include <dumux/io/grid/subgridgridcreator.hh>
+```
+
+and do the changes in the respective lines for the `Grid` property.
+
+The problem should now compile and run. However, an error occurs due to the coupling conditions.
+So far, we assumed a flat interface, therefore the normal momentum coupling condition
+
+ $`[\sigma \cdot \mathbf{n}]^{FF} = p^{PM}`$
+
+ was always set for a fixed $`\mathbf{n} = (0,1)^T`$. We need to account for the curvature of the interface and thus replace
+ ```cpp
+values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ ```
+ with
+ ```cpp
+values.setCouplingNeumann(scvf.directionIndex());
+ ```
+The same if true for the BJS condition, however, here we need to consider the tangential direction:
+```cpp
+values.setCouplingNeumann(1 - scvf.directionIndex());
+```
+
+The final result should look something like this:
+
+
+
+
+### 2. Changing the porous medium model
+
+In this part of the exercise, the coupled system will be extended such that the transport of components
+in both parts is included and the presence of a second phase in the porous medium is considered.
+This enables the simulation of the drying of a porous medium (however no energy balance is included yet).
+
+We start with an example in which the transport of water vapor in the gas phase is considered.
+The porous medium is filled with gas, initially the mole fraction of water vapor is $`x^w_g = 0.1`$.
+Above the porous medium, a dry gas flows and by diffusion, the porous medium dries out.
+
+
+__Task A: Change the model__:
+
+In the first task, the porous-medium model will be changed from a 1p2c system to a 2p2c system.
+Although a 2p2c system is plugged in, we still want to simulate the same situation as before, i.e., air with water vapor in both domains.
+The following changes have to be made in the porous-medium model (`ex_models_pmproblem.hh`):
+* Include the 2pnc model: include the respective headers and inherit from the new model `TwoPNC`
+* Exchange the spatial parameters for the 1-phase system by those for a 2-phase system (hint: two occurrences).
+* Since two phases are involved now, we do not need to use the OnePAdapter anymore. Change to property of the FluidSystem such that `H2OAir` is used directly.
+ Afterwards, set the `transportCompIdx` to `Indices::switchIdx`.
+
+One big difference between the 1p and 2p model is, that the primary variables for which
+the problem is solved, are not fixed.
+It is possible to use different formulations, e.g. with the gas pressure and the
+liquid saturation as primary variables (p_g-S_l -> `p1s0`) or vice versa.
+* Set the property
+
+```
+SET_PROP(DarcyTypeTag, Formulation)
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+```
+ in the Properties section in the problem file.
+
+In contrast to the formulation, which stays the same during one simulation, the meaning of
+the primary variables may vary during one simulation.
+In the case we are investigating, we want to use the gas pressure and the liquid saturation as primary variables.
+However, if only the gas phase is present, the liquid saturation is always zero.
+In this case, the chosen formulation will set the given value as the mole fraction of water vapor in the gas phase.
+* To tell to program which phases are present in which parts of the domain at the beginning of the simulation,
+ you have to call `values.setState(MY_PHASE_PRESENCE);` in `initialAtPos()`. Have a look at the `indices.hh`
+ in the `2p2c` model to figure out which is the correct value of `MY_PHASE_PRESENCE` for the presence of
+ a gas-phase only (hint: the numbering of phase indices begins with 0, the numbering of the phase presence states begins
+ with 1. Take a look at your formulation to find out which phase index to use for the gas phase.)
+
+Compile and run the new simulation.
+(Don't mind the compiler warning, we will deal with it in the next task.)
+
+__Task B: Add output__:
+
+In the next step, we want to add some output to the simulation.
+First we want to know the water mass in the (porous-medium) system.
+Therefore, we evaluate the storage term of the water component.
+* Have a look at the function `evaluateWaterMassStorageTerm()` in the porous medium subproblem.
+ Then implement a calculation of the total water mass:
+ $`\sum_{\alpha \in \textrm{g,l}} \left( \phi S_\alpha X^\text{w}_\alpha \varrho_\alpha V_\textrm{scv} \right)`$.
+ Afterwards, adapt the method `init()` such that the variable `initialWaterContent_` is initialized correctly using the `evaluateWaterMassStorageTerm()` method and assign that value also to the variable `lastWaterMass_`.
+
+We also want to investigate the temporal evolution of the water mass.
+The following steps need to be done to do so.
+Check if all instructions are implemented accordingly:
+* Calculate the initial water mass at the beginning of the simulation and add the water mass loss to the output.
+ Based on the water mass loss you can derive the evaporation rate. Because despite at the interface, all
+ boundaries have Neumann no-flow conditions and no sources are present, the water can only leave the system
+ via the porous-medium free-flow interface. The evaporation in [mm/d] can be calculated by:
+ $`e = \frac{\textrm{d}\, M^\text{w}}{\textrm{d}\, t} / A_\textrm{interface}`$.
+
+Finally we want to know the distribution of the water mass fluxes across the interface.
+* Have a look at the function `evaluateInterfaceFluxes()` in the porous medium problem.
+ Use the facilities therein to return the values of ...massCouplingCondition... from the `couplingManager`
+ for each coupling scvf. Then the fluxes are visualized with gnuplot, when setting `Problem.PlotFluxes = true`.
+ If the simulation is too fast, you can have a look at the flux*.png files after the simulation.
+* You can use the property `Problem.PlotStorage = true` to see the temporal evolution of the evaporation rate
+ and the cumulative water mass loss.
+
+Compile and run the simulation and take a look at the results.
+
+__Task C: Let it dry out__:
+
+In this exercise we want to completely dry an initial water-filled porous medium.
+- Change the initial condition such that it is started with a two-phase system.
+ Set the initial saturation to $`S_\text{l} = 0.1`$.
+
+If one wants to simulate the complete drying of a porous medium, the standard capillary pressure--saturation
+relationships have to be regularized. If no regularization is used, then the capillary pressure would
+approach infinity when the liquid saturation goes to zero. This means that water can flow to the
+interface from any region of the porous medium and of course this poses numerical problems.
+Therefore, the capillary pressure--saturation relationship has to be regularized for low saturations.
+The regularization has a strong influence on how long liquid water is present at the interface, see
+[Mosthaf (2014)](http://dx.doi.org/10.18419/opus-519).
+* Take a look at how the regularization is set in the `2pspatialparams.hh` and see how adapating
+ the parameter for `pcLowSw` and `pcHighSw` affects the results.
+
+The porous-medium model can now reach a liquid saturation of zero. As already explained above,
+for this case the liquid saturation cannot serve as primary variable anymore. However,
+manually adapting the primary variable states and values is inconvenient.
+[Class et al. (2002)](http://dx.doi.org/10.1016/S0309-1708(02)00014-3)
+describe an algorithm to switch the primary variables, if phases should appear or disappear during a simulation.
+- Replace the current implementation of the Newton solver with the version which can handle
+ primary variable switches: `dumux/multidomain/privarswitchnewtonsolver.hh`.
+ You also have to uncomment the line containing the `PriVarSwitchTuple` and to overwrite
+ the last argument with the `PrimaryVariableSwitch` property from the Darcy model.
+
+Now you are able to simulate a complete drying of the porous medium.
+
+
+### 4. Use a turbulence model in the free flow domain
+
+Several RANS turbulence models are implemented in DuMuX.
+This part of the exercise consists of the following steps:
+* replacing the Navier-Stokes model by the zero equation turbulence model,
+* switching to a symmetry boundary condition,
+* applying a grid refinement towards the interface,
+* subsequently refining the grid (convergence study).
+
+We will work with a `1p2cni/2p2cni` coupled problem, where `ni` stands for non-isothermal, and take the inertial forces into account.
+All the prepared files can be found in the subfolder `exercise-coupling-ff-pm/turbulence`.
+
+__Task A__:
+
+The file `ex_turbulence_ffproblem.hh` is your free flow problem file within this exercise.
+
+For using the compositional zero equation turbulence model, the following header files need to be included:
+```
+#include <dumux/freeflow/compositional/zeroeqncmodel.hh>
+#include <dumux/freeflow/rans/zeroeq/problem.hh>
+```
+The includes for the NavierStokesNC model and the NavierStokesProblem are no longer needed and can be removed.
+
+Make sure your free flow problem inherits from the correct parent type:
+* Change the last entry in the `NEW_TYPE_TAG` definition accordingly (non-isothermal zero equation model)
+* Adapt the inheritance of the problem class (hint: two occurrences)
+
+Take a look into the two headers included above to see how the correct TypeTag and the Problem class the inherit from are called.
+
+Additionally, you have to update the static and dynamic wall properties with
+```c++
+stokesProblem->updateStaticWallProperties();
+stokesProblem->updateDynamicWallProperties(stokesSol);
+```
+after applying the initial solution to the Stokes problem in the main file.
+
+In order to update these properties, the Stokes problem file needs to be extended to provide an isOnWall() method:
+```c++
+bool isOnWall(const GlobalPosition& globalPos) const
+{
+ return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+}
+```
+
+Since the dynamic wall properties can change during the simulation and depend on the solution, they have to be updated in each time step.
+Include the following lines into your main file (after `// update dynamic wall properties`):
+```c++
+stokesProblem->updateDynamicWallProperties(stokesSol);
+```
+
+Compile and run your new coupled problem and take a look at the results in Paraview.
+In addition to the standard variables and parameters, you can now choose turbulence model specific quantities (e.g. the turbulent viscosity `nu_t`) for the free flow domain.
+The result for the turbulent viscosity should look like this:
+
+
+
+__Task B__:
+
+Instead of computing the whole cross-section of a channel, you can use symmetric boundary conditions at the top boundary of your free flow domain by replacing all previous boundary conditions with
+```c++
+values.setAllSymmetry();
+```
+
+In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `isOnWall(globalPos)` method.
+
+__Task C__:
+
+Choose `Surface With Edges` instead of `Surface` in the Paraview toolbar to see the discretization grid.
+We will refine the grid now in order to better resolve the processes at the coupling interface.
+Since not much is happening at the upper and lower boundaries of the whole domain, we want to keep the resolution low in these areas to save some computation time.
+
+A grid refinement towards the interface is called _grading_.
+Try different gradings by changing the values in the respective sections in the input file:
+```c++
+Grading0 = 1.0
+Grading1 = 1.0
+```
+
+__Task D__:
+
+For the grid convergence study, run various simulations with the following grading parameters:
+```c++
+* [Stokes.Grid] Grading1 = 1.2, [Darcy.Grid] Grading1 = -1.2
+* [Stokes.Grid] Grading1 = 1.3, [Darcy.Grid] Grading1 = -1.3
+* [Stokes.Grid] Grading1 = 1.4, [Darcy.Grid] Grading1 = -1.4
+* [Stokes.Grid] Grading1 = 1.5, [Darcy.Grid] Grading1 = -1.5
+* [Stokes.Grid] Grading1 = 1.6, [Darcy.Grid] Grading1 = -1.6
+```
+
+By changing the parameter `Problem.Name` for each grading factor, you avoid loosing the `.vtu` and `.pvd` files of the previous simulation runs.
+Check the first lines of the output to see how the grading factors change the height of your grid cells.
+Compare the velocity fields with Paraview.
diff --git a/exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..249f34aa90d7853fc35384d9953d4c0284939ca8
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_input_file_links()
+
+# executables for ex_interface_coupling_ff-pm
+dune_add_test(NAME ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ CMD_ARGS ex_interface_coupling_ff-pm.input)
+
+# add tutorial to the common target
+add_dependencies(test_exercises ex_interface_coupling_ff-pm)
diff --git a/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.cc b/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..85226b998887b5ea030ed889ae598573ce250da3
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.cc
@@ -0,0 +1,281 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the coupled Stokes/Darcy problem (1p)
+ */
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/istl/io.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/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_interface_pmproblem.hh"
+#include "ex_interface_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+SET_PROP(StokesOnePTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyOnePTypeTag)>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+SET_PROP(DarcyOnePTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TTAG(StokesOnePTypeTag), TTAG(StokesOnePTypeTag), TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ 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 StokesTypeTag = TTAG(StokesOnePTypeTag);
+ using DarcyTypeTag = TTAG(DarcyOnePTypeTag);
+
+
+
+ // ******************** comment-out this section for the last exercise **************** //
+
+ // create two individual grids (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+
+ // ************************************************************************************ //
+
+
+ // ******************** uncomment this section for the last exercise ****************** //
+
+ // // use dune-subgrid to create the individual grids
+ // static constexpr int dim = 2;
+ // using HostGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >;
+ // using HostGridManager = Dumux::GridManager<HostGrid>;
+ // HostGridManager hostGridManager;
+ // hostGridManager.init();
+ // auto& hostGrid = hostGridManager.grid();
+ //
+ // struct Params
+ // {
+ // double amplitude = getParam<double>("Grid.Amplitude");
+ // double baseline = getParam<double>("Grid.Baseline");
+ // double offset = getParam<double>("Grid.Offset");
+ // double scaling = getParam<double>("Grid.Scaling");
+ // };
+ //
+ // Params params;
+ //
+ // auto elementSelectorStokes = [&](const auto& element)
+ // {
+ // double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ // return element.geometry().center()[1] > interface;
+ // };
+ //
+ // auto elementSelectorDarcy = [&](const auto& element)
+ // {
+ // double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ // return element.geometry().center()[1] < interface;
+ // };
+ //
+ // // subgrid Pointer
+ // auto stokesGridPtr = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorStokes, "Stokes");
+ // auto darcyGridPtr = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorDarcy, "Darcy");
+ //
+ // // we compute on the leaf grid view
+ // const auto& darcyGridView = darcyGridPtr->leafGridView();
+ // const auto& stokesGridView = stokesGridPtr->leafGridView();
+
+ // ************************************************************************************ //
+
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ const auto& cellCenterSol = sol[stokesCellCenterIdx];
+ const auto& faceSol = sol[stokesFaceIdx];
+
+ // apply initial solution for instationary problems
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ std::get<0>(stokesSol) = cellCenterSol;
+ std::get<1>(stokesSol) = faceSol;
+ stokesProblem->applyInitialSolution(stokesSol);
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // intialize the vtk output module
+ const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+ StaggeredVtkOutputModule<StokesTypeTag> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+
+ //****** uncomment the add analytical solution of v_x *****//
+ // stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ 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(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ 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);
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // write vtk output
+ stokesVtkWriter.write(1.0);
+ darcyVtkWriter.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
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input b/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
new file mode 100644
index 0000000000000000000000000000000000000000..e3813950725ba4837e96b79d76dc9ccc74dc0d4f
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
@@ -0,0 +1,45 @@
+ # for dune-subgrid
+ [Grid]
+Positions0 = 0 1
+Positions1 = 0 0.2 0.3 0.65
+Cells0 = 100
+Cells1 = 10 50 18
+Baseline = 0.25 # [m]
+Amplitude = 0.04 # [m]
+Offset = 0.5 # [m]
+Scaling = 0.2 #[m]
+
+[Stokes.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 1.0 2.0
+Cells0 = 20
+Cells1 = 100
+Grading1 = 1
+
+[Darcy.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 0.0 1.0
+Cells0 = 20
+Cells1 = 20
+Grading1 = 1
+
+[Stokes.Problem]
+Name = stokes
+PressureDifference = 1e-9
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+Permeability = 1e-6 # m^2
+Porosity = 0.4
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+Name = ex_ff-pm-interface
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1
diff --git a/exercises/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh b/exercises/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..48d29d667983db35a315d3417554f0c1880784db
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
@@ -0,0 +1,323 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief The free flow sub problem
+ */
+#ifndef DUMUX_STOKES_SUBPROBLEM_HH
+#define DUMUX_STOKES_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+//****** uncomment for the last exercise *****//
+// #include <dumux/io/grid/subgridgridcreator.hh>
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/simpleh2o.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(StokesOnePTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
+
+// the fluid system
+SET_PROP(StokesOnePTypeTag, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+SET_PROP(StokesOnePTypeTag, Grid)
+{
+ static constexpr auto dim = 2;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+//****** comment out for the last exercise *****//
+ using type = TensorGrid;
+
+//****** uncomment for the last exercise *****//
+ // using HostGrid = TensorGrid;
+ // using type = Dune::SubGrid<dim, HostGrid>;
+};
+
+// Set the problem property
+SET_TYPE_PROP(StokesOnePTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
+
+SET_BOOL_PROP(StokesOnePTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(StokesOnePTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(StokesOnePTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(StokesOnePTypeTag, EnableInertiaTerms, false);
+}
+
+/*!
+ * \brief The free flow sub problem
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+
+public:
+ StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return 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
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.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.dofPosition();
+
+ if(onUpperBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+
+ // left/right wall
+ if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+
+ // coupling interface
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(globalPos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry 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 ElementFaceVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ }
+
+ return values;
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]);
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ /*!
+ * \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967)
+ */
+ void calculateAnalyticalVelocityX() const
+ {
+ analyticalVelocityX_.resize(this->fvGridGeometry().gridView().size(0));
+
+ using std::sqrt;
+ const Scalar dPdX = -deltaP_ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]);
+ static const Scalar mu = FluidSystem::viscosity(temperature(), 1e5);
+ 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 = (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().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(this->fvGridGeometry().gridView()))
+ {
+ const auto eIdx = this->fvGridGeometry().gridView().indexSet().index(element);
+ const Scalar y = element.geometry().center()[1] - this->fvGridGeometry().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;
+ }
+ }
+
+ /*!
+ * \brief Get the analytical velocity in x direction
+ */
+ const std::vector<Scalar>& getAnalyticalVelocityX() const
+ {
+ if(analyticalVelocityX_.empty())
+ calculateAnalyticalVelocityX();
+ return analyticalVelocityX_;
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+ Scalar deltaP_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ mutable std::vector<Scalar> analyticalVelocityX_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES_SUBPROBLEM_HH
diff --git a/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh b/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7dd08bc3d41b34bbe7609a386d1bce637a24a5fc
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
@@ -0,0 +1,251 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+* \file
+*
+* \brief The porous medium flow sub problem
+*/
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+//****** uncomment for the last exercise *****//
+// #include <dumux/io/grid/subgridgridcreator.hh>
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "../1pspatialparams.hh"
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(DarcyOnePTypeTag, INHERITS_FROM(CCTpfaModel, OneP));
+
+// Set the problem property
+SET_TYPE_PROP(DarcyOnePTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// the fluid system
+SET_PROP(DarcyOnePTypeTag, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+SET_PROP(DarcyOnePTypeTag, Grid)
+{
+ static constexpr auto dim = 2;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+//****** comment out for the last exercise *****//
+ using type = TensorGrid;
+
+//****** uncomment for the last exercise *****//
+ // using HostGrid = TensorGrid;
+ // using type = Dune::SubGrid<dim, HostGrid>;
+};
+
+SET_TYPE_PROP(DarcyOnePTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+}
+
+/*!
+ * \brief The porous medium flow sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {}
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Return 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 (onLowerBoundary_(scvf.center()))
+ values.setAllDirichlet();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ PrimaryVariables values(0.0);
+ values = initial(element);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ * For this method, the \a values variable stores primary variables.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param element The element
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initial(const Element &element) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1776059b237e7dd90b48d577bdf01b3f463e7c34
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_input_file_links()
+
+# executables for ex_interface_coupling_ff-pm
+dune_add_test(NAME ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+# add tutorial to the common target
+add_dependencies(test_exercises ex_models_coupling_ff-pm)
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b54ae19ed5f17fa1d5bcc47a0d3e01f1c5f502ed
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
@@ -0,0 +1,296 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the coupled Stokes/Darcy problem (1p)
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+#include <fstream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/geometry/diameter.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_models_pmproblem.hh"
+#include "ex_models_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+SET_PROP(StokesTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTypeTag)>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+SET_PROP(DarcyTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TTAG(StokesTypeTag), TTAG(StokesTypeTag), TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ 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 StokesTypeTag = TTAG(StokesTypeTag);
+ using DarcyTypeTag = TTAG(DarcyTypeTag);
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+ stokesProblem->setTimeLoop(timeLoop);
+ darcyProblem->setTimeLoop(timeLoop);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ const auto& cellCenterSol = sol[stokesCellCenterIdx];
+ const auto& faceSol = sol[stokesFaceIdx];
+
+ // apply initial solution for instationary problems
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ std::get<0>(stokesSol) = cellCenterSol;
+ std::get<1>(stokesSol) = faceSol;
+ stokesProblem->applyInitialSolution(stokesSol);
+ auto solStokesOld = stokesSol;
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+ auto solDarcyOld = sol[darcyIdx];
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol, solStokesOld);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+
+ // intialize the vtk output module
+ const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+ StaggeredVtkOutputModule<StokesTypeTag> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // intialize the subproblems
+ darcyProblem->init(sol[darcyIdx], *darcyGridVariables);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the primary variable switches used by the sub models and the non-linear solver
+// using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, NoPrimaryVariableSwitch>;
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ const auto episodeLength = getParam<Scalar>("TimeLoop.EpisodeLength");
+ if (episodeLength > 0.0)
+ timeLoop->setPeriodicCheckPoint(episodeLength);
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(solOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(sol, *timeLoop);
+
+ // make the new solution the old solution
+ solOld = sol;
+ stokesGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // call the postTimeStep routine for output
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables);
+
+ // write vtk output
+ if (timeLoop->isCheckPoint() || timeLoop->finished() || episodeLength < 0.0)
+ {
+ stokesVtkWriter.write(timeLoop->time());
+ darcyVtkWriter.write(timeLoop->time());
+ }
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(stokesGridView.comm());
+ timeLoop->finalize(darcyGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // 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
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
new file mode 100644
index 0000000000000000000000000000000000000000..d3f68518fa9a38504ee876af35f0c9e79f31cf17
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
@@ -0,0 +1,49 @@
+[TimeLoop]
+DtInitial = 100 # s
+EpisodeLength = -360 # s # 0.25 days
+TEnd = 256000 # s # 2 days
+
+[Stokes.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 16 16
+
+[Darcy.Grid]
+UpperRight = 1 1
+Cells = 16 16
+
+[Stokes.Problem]
+Name = stokes
+EnableGravity = false
+MoleFraction = 0.0 # -
+Pressure = 1e5 # Pa
+Velocity = 1e-3 # m/s
+
+[Darcy.Problem]
+Name = darcy
+EnableGravity = true
+Saturation = 0.1 # -
+MoleFraction = 0.1 # -
+Pressure = 1e5 # Pa
+
+[Darcy.SpatialParams]
+Permeability = 2.65e-10 # m^2
+Porosity = 0.4 # -
+AlphaBeaversJoseph = 1.0 # -
+# EXNUMBER >= 1
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.5e-4
+VgN = 8.0
+
+[Problem]
+Name = models_coupling
+PlotFluxes = false
+PlotStorage = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = 1
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh b/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1fe98c085805a45045474cece286db0f4561b399
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
@@ -0,0 +1,341 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesTests
+ * \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model.
+ */
+#ifndef DUMUX_STOKES1P2C_SUBPROBLEM_HH
+#define DUMUX_STOKES1P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1padapter.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(StokesTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+
+// Set the grid type
+SET_TYPE_PROP(StokesTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// The fluid system
+SET_PROP(StokesTypeTag, FluidSystem)
+{
+ using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
+};
+
+// Do not replace one equation with a total mass balance
+SET_INT_PROP(StokesTypeTag, ReplaceCompEqIdx, 3);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(StokesTypeTag, UseMoles, true);
+
+// Set the problem property
+SET_TYPE_PROP(StokesTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
+
+SET_BOOL_PROP(StokesTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(StokesTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(StokesTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(StokesTypeTag, EnableInertiaTerms, false);
+}
+
+/*!
+ * \ingroup NavierStokesTests
+ * \brief Test problem for the one-phase compositional (Navier-) Stokes problem.
+ *
+ * Horizontal flow from left to right with a parabolic velocity profile.
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+
+public:
+ StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ velocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+ moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return 293.15; }
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.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();
+
+ if(onLeftBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::conti0EqIdx + 1);
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+ else if(onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx + 1);
+ }
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element
+ * \param scvf The subcontrolvolume face
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(pos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ }
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, pressure_);
+ const Scalar density = FluidSystem::density(fluidState, 0);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = pressure_ + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[Indices::conti0EqIdx + 1] = moleFraction_;
+ values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1])
+ * (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
+ / (height_() * height_());
+
+ return values;
+ }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ //! \brief updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState, const Scalar pressure) const
+ {
+ fluidState.setTemperature(temperature());
+ fluidState.setPressure(0, pressure);
+ fluidState.setSaturation(0, 1.0);
+ fluidState.setMoleFraction(0, 1, moleFraction_);
+ fluidState.setMoleFraction(0, 0, 1.0 - moleFraction_);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, 0);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+ fluidState.setDensity(0, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+ fluidState.setMolarDensity(0, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+ fluidState.setEnthalpy(0, enthalpy);
+ }
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar velocity_;
+ Scalar pressure_;
+ Scalar moleFraction_;
+
+ TimeLoopPtr timeLoop_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh b/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0e48d23560788095b46c641d55c64606f9ecf6d6
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
@@ -0,0 +1,404 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A simple Darcy test problem (cell-centered finite volume method).
+ */
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/material/fluidsystems/1padapter.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/porousmediumflow/1pnc/model.hh>
+#include "../1pspatialparams.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
+
+// Set the problem property
+SET_TYPE_PROP(DarcyTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// The fluid system
+SET_PROP(DarcyTypeTag, FluidSystem)
+{
+ using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
+};
+
+// Use moles
+SET_BOOL_PROP(DarcyTypeTag, UseMoles, true);
+
+// Do not replace one equation with a total mass balance
+SET_INT_PROP(DarcyTypeTag, ReplaceCompEqIdx, 3);
+
+//! Use a model with constant tortuosity for the effective diffusivity
+SET_TYPE_PROP(DarcyTypeTag, EffectiveDiffusivityModel,
+ DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set the grid type
+SET_TYPE_PROP(DarcyTypeTag, Grid, Dune::YaspGrid<2>);
+
+// Set the spatial paramaters type
+SET_TYPE_PROP(DarcyTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+}
+
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ enum {
+ // grid and world dimension
+ dim = GridView::dimension,
+ dimworld = GridView::dimensionworld,
+
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ pressureIdx = Indices::pressureIdx,
+ phaseIdx = 0,
+ transportCompIdx = 1
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+ moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+
+ // initialize output file
+ plotFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotFluxes", false);
+ plotStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotStorage", false);
+ storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+ storageFile_.open(storageFileName_);
+ storageFile_ << "#Time[s]" << ";"
+ << "WaterMass[kg]" << ";"
+ << "WaterMassLoss[kg]" << ";"
+ << "EvaporationRate[mm/d]"
+ << std::endl;
+ }
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Initialize the problem.
+ */
+ template<class SolutionVector, class GridVariables>
+ void init(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+ { }
+
+ template<class SolutionVector, class GridVariables>
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ evaluateWaterMassStorageTerm(curSol, gridVariables);
+ evaluateInterfaceFluxes(curSol, gridVariables);
+
+ gnuplotStorage_.resetPlot();
+ gnuplotStorage_.setDatafileSeparator(';');
+ gnuplotStorage_.setXlabel("time [d]");
+ gnuplotStorage_.setXRange(0.0, getParam<Scalar>("TimeLoop.TEnd"));
+ gnuplotStorage_.setYlabel("evaporation rate [mm/d]");
+ gnuplotStorage_.setOption("set yrange [0.0:]");
+ gnuplotStorage_.setOption("set y2label 'cumulative mass loss'");
+ gnuplotStorage_.setOption("set y2range [0.0:0.5]");
+ gnuplotStorage_.setOption("set y2range [0.0:0.5]");
+ gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:4 with lines title 'evaporation rate'");
+ gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:3 axes x1y2 with lines title 'cumulative mass loss'");
+ if (plotStorage_)
+ gnuplotStorage_.plot("temp");
+ }
+
+ template<class SolutionVector, class GridVariables>
+ Scalar evaluateWaterMassStorageTerm(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ // compute the mass in the entire domain
+ Scalar waterMass = 0.0;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ // insert calculation of the water mass here
+ waterMass += 0.0;
+ }
+ }
+
+ Scalar cumMassLoss = initialWaterContent_ - waterMass;
+ Scalar evaporationRate = (lastWaterMass_ - waterMass) * 86400
+ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0])
+ / timeLoop_->timeStepSize();
+ lastWaterMass_ = waterMass;
+
+ storageFile_ << timeLoop_->time() << ";"
+ << waterMass << ";"
+ << cumMassLoss << ";"
+ << evaporationRate
+ << std::endl;
+
+ return waterMass;
+ }
+
+ template<class SolutionVector, class GridVariables>
+ void evaluateInterfaceFluxes(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ std::vector<Scalar> x;
+ std::vector<Scalar> y;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ continue;
+
+ // NOTE: binding the coupling context is necessary
+ couplingManager_->bindCouplingContext(CouplingManager::darcyIdx, element);
+ NumEqVector flux(0.0); // use "massCouplingCondition" from the couplingManager here
+
+ x.push_back(scvf.center()[0]);
+ y.push_back(flux[transportCompIdx]);
+ }
+ }
+
+ gnuplotInterfaceFluxes_.resetPlot();
+ gnuplotInterfaceFluxes_.setXlabel("x-position [m]");
+ gnuplotInterfaceFluxes_.setXRange(this->fvGridGeometry().bBoxMin()[0], this->fvGridGeometry().bBoxMax()[0]);
+ gnuplotInterfaceFluxes_.setYlabel("flux [kg/(m^2 s)]");
+ gnuplotInterfaceFluxes_.setYRange(-5e-4, 0.0);
+ gnuplotInterfaceFluxes_.setOption("set label 'time: " + std::to_string(timeLoop_->time()/86400.) + "d' at graph 0.8,0.8 ");
+ std::string fluxFileName = "flux_" + std::to_string(timeLoop_->timeStepIndex()) +
+ "_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+ gnuplotInterfaceFluxes_.addDataSetToPlot(x, y, fluxFileName, "with lines title 'water mass flux'");
+ if (plotFluxes_)
+ gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex()));
+ }
+
+ /*!
+ * \brief Returns true if a restart file should be written to
+ * disk.
+ */
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ bool shouldWriteOutput() const // define output
+ { return true; }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 293.15; }
+ // \}
+
+ /*!
+ * \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().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ * For this method, the \a values variable stores primary variables.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param element The element
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values[transportCompIdx] = moleFraction_;
+ values[pressureIdx] = pressure_;
+ return values;
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+ Scalar moleFraction_;
+ Scalar pressure_;
+
+ Scalar initialWaterContent_ = 0.0;
+ Scalar lastWaterMass_ = 0.0;
+
+ TimeLoopPtr timeLoop_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ std::string storageFileName_;
+ std::ofstream storageFile_;
+ bool plotFluxes_;
+ bool plotStorage_;
+ Dumux::GnuplotInterface<Scalar> gnuplotInterfaceFluxes_;
+ Dumux::GnuplotInterface<Scalar> gnuplotStorage_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3fb93b14e46940ca5110f1b2ab3f69143b089593
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_input_file_links()
+
+# executables for ex_interface_coupling_ff-pm
+dune_add_test(NAME ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ CMD_ARGS ex_turbulence_coupling_ff-pm.input)
+
+# add tutorial to the common target
+add_dependencies(test_exercises ex_turbulence_coupling_ff-pm)
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..80a5b9be520e8573a020b083c213fa95c2dfee9e
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
@@ -0,0 +1,300 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the isothermal coupled Stokes/Darcy problem (1p2c/2p2c)
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+#include <fstream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/geometry/diameter.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_turbulence_pmproblem.hh"
+#include "ex_turbulence_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+SET_PROP(ZeroEqTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTwoPTwoCTypeTag)>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+SET_PROP(DarcyTwoPTwoCTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TTAG(ZeroEqTypeTag), TTAG(ZeroEqTypeTag), TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ 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 StokesTypeTag = TTAG(ZeroEqTypeTag);
+ using DarcyTypeTag = TTAG(DarcyTwoPTwoCTypeTag);
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // set timeloop for the subproblems, needed for boundary value variations
+ stokesProblem->setTimeLoop(timeLoop);
+ darcyProblem->setTimeLoop(timeLoop);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ // apply initial solution for instationary problems
+ // auxiliary free flow solution vector
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ stokesSol[stokesCellCenterIdx].resize(sol[stokesCellCenterIdx].size());
+ stokesSol[stokesFaceIdx].resize(sol[stokesFaceIdx].size());
+ stokesProblem->applyInitialSolution(stokesSol);
+ auto solStokesOld = stokesSol;
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+
+ // TODO: update static wall properties
+ // TODO: update dynamic wall properties
+
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+ auto solDarcyOld = sol[darcyIdx];
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol, solStokesOld);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+
+ // intialize the vtk output module
+ const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+ StaggeredVtkOutputModule<StokesTypeTag, GET_PROP_VALUE(StokesTypeTag, PhaseIdx)> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the primary variable switches used by the sub models
+ using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(solOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(sol, *timeLoop);
+
+ // make the new solution the old solution
+ solOld = sol;
+
+ // update the auxiliary free flow solution vector
+ stokesSol[stokesCellCenterIdx] = sol[stokesCellCenterIdx];
+ stokesSol[stokesFaceIdx] = sol[stokesFaceIdx];
+
+ // TODO: update dynamic wall properties
+
+ // post time step treatment of Darcy problem
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables, timeLoop->timeStepSize());
+
+ // advance grid variables to the next time step
+ stokesGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ stokesVtkWriter.write(timeLoop->time());
+ darcyVtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(stokesGridView.comm());
+ timeLoop->finalize(darcyGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // 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
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input
new file mode 100644
index 0000000000000000000000000000000000000000..514cf156d55211777c2ec0980f37d1a95124867e
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input
@@ -0,0 +1,62 @@
+[TimeLoop]
+DtInitial = 1e-1 # [s]
+MaxTimeStepSize = 43200 # [s] (12 hours)
+TEnd = 864000 # [s] (6 days)
+
+[Stokes.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.25 0.5
+Grading0 = 1.0
+Grading1 = 1.0
+Cells0 = 15
+Cells1 = 20
+Verbosity = true
+
+[Darcy.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.0 0.25
+Cells0 = 15
+Cells1 = 10
+Grading0 = 1.0
+Grading1 = 1.0
+Verbosity = true
+
+[Stokes.Problem]
+Name = stokes
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+refMoleFrac = 0 # [-]
+RefTemperature = 298.15 # [K]
+
+[Darcy.Problem]
+Name = darcy
+Pressure = 1e5
+Saturation = 0.5 # initial Sw
+Temperature = 298.15 # [K]
+InitPhasePresence = 3 # bothPhases
+
+[Darcy.SpatialParams]
+Porosity = 0.41
+Permeability = 2.65e-10
+AlphaBeaversJoseph = 1.0
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.371e-4
+VgN = 6.9
+
+[Problem]
+Name = ex_coupling_turbulence_ff-pm
+EnableGravity = true
+InterfaceDiffusionCoefficientAvg = Harmonic
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Assembly]
+NumericDifferenceMethod = 0
+NumericDifference.BaseEpsilon = 1e-8
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..861e86ce9ff366b24376846708bb88fec1a7e4c3
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
@@ -0,0 +1,378 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief The free-flow sub problem
+ */
+#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1padapter.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+#include <dumux/freeflow/navierstokes/problem.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class FreeFlowSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
+
+// Set the grid type
+SET_TYPE_PROP(ZeroEqTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// The fluid system
+SET_PROP(ZeroEqTypeTag, FluidSystem)
+{
+ using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
+ using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
+
+SET_INT_PROP(ZeroEqTypeTag, ReplaceCompEqIdx, 3);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(ZeroEqTypeTag, UseMoles, true);
+
+// Set the problem property
+SET_TYPE_PROP(ZeroEqTypeTag, Problem, Dumux::FreeFlowSubProblem<TypeTag> );
+
+SET_BOOL_PROP(ZeroEqTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(ZeroEqTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(ZeroEqTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(ZeroEqTypeTag, EnableInertiaTerms, true);
+}
+
+/*!
+ * \brief The free-flow sub problem
+ */
+template <class TypeTag>
+class FreeFlowSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+ static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+
+public:
+ FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ refVelocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefVelocity");
+ refPressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefPressure");
+ refMoleFrac_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.refMoleFrac");
+ refTemperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefTemperature");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return refTemperature_; }
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.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();
+
+ if (onLeftBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setDirichlet(Indices::conti0EqIdx + 1);
+ values.setDirichlet(Indices::energyBalanceIdx);
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyBalanceIdx);
+ }
+
+ if (onUpperBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyBalanceIdx);
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx + 1);
+ values.setOutflow(Indices::energyBalanceIdx);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ values.setCouplingNeumann(Indices::energyBalanceIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element
+ * \param scvf The subcontrolvolume face
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(pos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ values[Indices::energyBalanceIdx] = couplingManager().couplingData().energyCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ }
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, refTemperature(), refPressure(), refMoleFrac());
+
+ const Scalar density = FluidSystem::density(fluidState, 0);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = refPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[Indices::conti0EqIdx + 1] = refMoleFrac();
+ values[Indices::velocityXIdx] = refVelocity();
+ values[Indices::temperatureIdx] = refTemperature();
+
+ if(onUpperBoundary_(globalPos))
+ values[Indices::velocityXIdx] = 0.0;
+
+ return values;
+ }
+
+ //! \brief Returns the reference velocity.
+ const Scalar refVelocity() const
+ { return refVelocity_ ;}
+
+ //! \brief Returns the reference pressure.
+ const Scalar refPressure() const
+ { return refPressure_; }
+
+ //! \brief Returns the reference mass fraction.
+ const Scalar refMoleFrac() const
+ { return refMoleFrac_; }
+
+ //! \brief Returns the reference temperature.
+ const Scalar refTemperature() const
+ { return refTemperature_; }
+
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ //! \brief updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState, const Scalar temperature,
+ const Scalar pressure, const Scalar moleFraction) const
+ {
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(0, pressure);
+ fluidState.setSaturation(0, 1.0);
+ fluidState.setMoleFraction(0, 1, moleFraction);
+ fluidState.setMoleFraction(0, 0, 1.0 - moleFraction);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, 0);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+ fluidState.setDensity(0, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+ fluidState.setMolarDensity(0, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+ fluidState.setEnthalpy(0, enthalpy);
+ }
+
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar refVelocity_;
+ Scalar refPressure_;
+ Scalar refMoleFrac_;
+ Scalar refTemperature_;
+
+ TimeLoopPtr timeLoop_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2d3a84d25c4884441720cc2ad10b0ccb8ab7a4ab
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
@@ -0,0 +1,324 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+ /*!
+ * \file
+ *
+ * \brief The porous medium sub problem
+ */
+#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH
+#define DUMUX_DARCY2P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+
+#include <dumux/porousmediumflow/2p2c/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+#include "../2pspatialparams.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(DarcyTwoPTwoCTypeTag, INHERITS_FROM(CCTpfaModel, TwoPTwoCNI));
+
+// Set the problem property
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// the fluid system
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+SET_PROP(DarcyTwoPTwoCTypeTag, Formulation)
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+
+// The gas component balance (air) is replaced by the total mass balance
+SET_INT_PROP(DarcyTwoPTwoCTypeTag, ReplaceCompEqIdx, 3);
+
+// Set the grid type
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+SET_BOOL_PROP(DarcyTwoPTwoCTypeTag, UseMoles, true);
+
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+}
+
+/*!
+ * \brief The porous medium sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ enum {
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ contiWEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
+ contiNEqIdx = Indices::conti0EqIdx + FluidSystem::AirIdx,
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+ initialSw_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
+ temperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Temperature");
+ initialPhasePresence_ = getParamFromGroup<int>(this->paramGroup(), "Problem.InitPhasePresence");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ }
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ template<class SolutionVector, class GridVariables>
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables,
+ const Scalar timeStepSize)
+
+ {
+ // compute the mass in the entire domain
+ Scalar massWater = 0.0;
+
+ // bulk elements
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
+ {
+ massWater += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx)*volVars.density(phaseIdx)
+ * scv.volume() * volVars.saturation(phaseIdx) * volVars.porosity() * volVars.extrusionFactor();
+ }
+ }
+ }
+
+ std::cout << "mass of water is: " << massWater << std::endl;
+ }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return temperature_; }
+ // \}
+
+ /*!
+ * \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().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(scvf.center());
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann
+ * control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ {
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+
+ for(int i = 0; i< massFlux.size(); ++i)
+ values[i] = massFlux[i];
+
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+ }
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ *
+ * For this method, the \a values variable stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values.setState(initialPhasePresence_);
+
+ values[pressureIdx] = pressure_ + 1. * this->gravity()[1] * (globalPos[1] - this->fvGridGeometry().bBoxMax()[1]);
+ values[switchIdx] = initialSw_;
+ values[Indices::temperatureIdx] = temperature_;
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar pressure_;
+ Scalar initialSw_;
+ Scalar temperature_;
+ int initialPhasePresence_;
+
+ TimeLoopPtr timeLoop_;
+
+ Scalar eps_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY2P2C_SUBPROBLEM_HH
diff --git a/exercises/extradoc/ex_ff-pm-turb_viscosity.png b/exercises/extradoc/ex_ff-pm-turb_viscosity.png
new file mode 100644
index 0000000000000000000000000000000000000000..1e7ddab16db7d8a718490ffbf5c2a57ddb3aeacf
Binary files /dev/null and b/exercises/extradoc/ex_ff-pm-turb_viscosity.png differ
diff --git a/exercises/extradoc/ex_ff-pm-vertical-flow.png b/exercises/extradoc/ex_ff-pm-vertical-flow.png
new file mode 100644
index 0000000000000000000000000000000000000000..2c771210ad0f7d1c2d37a6f077642a32925483b4
Binary files /dev/null and b/exercises/extradoc/ex_ff-pm-vertical-flow.png differ
diff --git a/exercises/extradoc/ex_ff-pm-wave-interface.png b/exercises/extradoc/ex_ff-pm-wave-interface.png
new file mode 100644
index 0000000000000000000000000000000000000000..15c3a01caebfdd30885afd462f9eff4313436bd5
Binary files /dev/null and b/exercises/extradoc/ex_ff-pm-wave-interface.png differ
diff --git a/exercises/solution/CMakeLists.txt b/exercises/solution/CMakeLists.txt
index 50ac3c623f92d870c3ae0bafb54df8d4b3863476..b55b7ace7f2cd022db71b2769f2548e5bb55136f 100644
--- a/exercises/solution/CMakeLists.txt
+++ b/exercises/solution/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory(exercise-basic)
+add_subdirectory(exercise-coupling-ff-pm)
add_subdirectory(exercise-fluidsystem)
add_subdirectory(exercise-grids)
add_subdirectory(exercise-mainfile)
diff --git a/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh b/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..64887e33307e8a45908dabc131e919a99f95dfa1
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh
@@ -0,0 +1,98 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup OnePTests
+ * \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 OnePModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief The spatial parameters class for the test problem using the
+ * 1p cc model
+ */
+template<class TypeTag>
+class OnePSpatialParams
+: public FVSpatialParamsOneP<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar),
+ OnePSpatialParams<TypeTag>>
+{
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<TypeTag>>;
+
+ 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 FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
+ 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 Define the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return porosity_; }
+
+ /*! \brief Define the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+
+private:
+ Scalar permeability_;
+ Scalar porosity_;
+ Scalar alphaBJ_;
+};
+
+} // end namespace
+
+#endif
diff --git a/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh b/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1830f3cbf2381280ca1213b8a07a17890846203c
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh
@@ -0,0 +1,139 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPTests
+ * \brief The spatial parameters class for the test problem using the 2p cc model
+ */
+#ifndef DUMUX_TWOPHASE_SPATIAL_PARAMS_HH
+#define DUMUX_TWOPHASE_SPATIAL_PARAMS_HH
+
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief The spatial parameters class for the test problem using the 2p cc model
+ */
+template<class TypeTag>
+class TwoPSpatialParams
+: public FVSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar),
+ TwoPSpatialParams<TypeTag>>
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using ParentType = FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<TypeTag>>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
+
+public:
+ using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
+ using MaterialLawParams = typename MaterialLaw::Params;
+ using PermeabilityType = Scalar;
+
+ TwoPSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
+ alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+
+ // residual saturations
+ params_.setSwr(getParam<Scalar>("Darcy.SpatialParams.Swr"));
+ params_.setSnr(getParam<Scalar>("Darcy.SpatialParams.Snr"));
+ // parameters for the vanGenuchten law
+ params_.setVgAlpha(getParam<Scalar>("Darcy.SpatialParams.VgAlpha"));
+ params_.setVgn(getParam<Scalar>("Darcy.SpatialParams.VgN"));
+ params_.setPcLowSw(params_.swr()*5.0);
+ params_.setPcHighSw(1.0-params_.snr()*5.0);
+ }
+
+ /*!
+ * \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 Define the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return porosity_; }
+
+ /*! \brief Define the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+ /*!
+ * \brief Returns the parameter object for the Brooks-Corey material law.
+ * In this test, we use element-wise distributed material parameters.
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \return the material parameters object
+ */
+ template<class ElementSolutionVector>
+ const MaterialLawParams& materialLawParams(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return params_; }
+
+ /*!
+ * \brief Function for defining which phase is to be considered as the wetting phase.
+ *
+ * \return the wetting phase index
+ * \param globalPos The global position
+ */
+ template<class FluidSystem>
+ int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+ { return FluidSystem::phase0Idx; }
+
+private:
+ Scalar permeability_;
+ Scalar porosity_;
+ Scalar alphaBJ_;
+ MaterialLawParams params_;
+ static constexpr Scalar eps_ = 1.0e-7;
+};
+
+} // end namespace
+
+#endif
diff --git a/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bae4c263a6bf7ae6a06d90af66a2e45d224f9b2e
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_subdirectory(interface)
+add_subdirectory(models)
+add_subdirectory(turbulence)
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..24d0d12ed5870a60bd54491d14d04c9993359fc4
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt
@@ -0,0 +1,25 @@
+add_input_file_links()
+
+dune_add_test(NAME orig_ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=0
+ CMD_ARGS ex_interface_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_a_ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=1
+ CMD_ARGS ex_interface_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_b_ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=2
+ CMD_ARGS ex_interface_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_c_ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=3
+ CMD_ARGS ex_interface_coupling_ff-pm.input)
+
+
+# add exercise to the common target
+add_dependencies(test_exercises sol_a_ex_interface_coupling_ff-pm sol_b_ex_interface_coupling_ff-pm sol_b_ex_interface_coupling_ff-pm)
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.cc b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..116ed7e2d8e4469f44bde197b8b71974422264b1
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.cc
@@ -0,0 +1,274 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the coupled Stokes/Darcy problem (1p)
+ */
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/istl/io.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/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_interface_pmproblem.hh"
+#include "ex_interface_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+SET_PROP(StokesOnePTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyOnePTypeTag)>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+SET_PROP(DarcyOnePTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TTAG(StokesOnePTypeTag), TTAG(StokesOnePTypeTag), TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ 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 StokesTypeTag = TTAG(StokesOnePTypeTag);
+ using DarcyTypeTag = TTAG(DarcyOnePTypeTag);
+
+#if EXNUMBER < 3
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+#else
+ // use dune-subgrid to create the individual grids
+ static constexpr int dim = 2;
+ using HostGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >;
+ using HostGridManager = Dumux::GridManager<HostGrid>;
+ HostGridManager hostGridManager;
+ hostGridManager.init();
+ auto& hostGrid = hostGridManager.grid();
+
+ struct Params
+ {
+ double amplitude = getParam<double>("Grid.Amplitude");
+ double baseline = getParam<double>("Grid.Baseline");
+ double offset = getParam<double>("Grid.Offset");
+ double scaling = getParam<double>("Grid.Scaling");
+ };
+
+ Params params;
+
+ auto elementSelectorStokes = [&](const auto& element)
+ {
+ double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ return element.geometry().center()[1] > interface;
+ };
+
+ auto elementSelectorDarcy = [&](const auto& element)
+ {
+ double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ return element.geometry().center()[1] < interface;
+ };
+
+ // subgrid Pointer
+ auto stokesGridPtr = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorStokes, "Stokes");
+ auto darcyGridPtr = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorDarcy, "Darcy");
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridPtr->leafGridView();
+ const auto& stokesGridView = stokesGridPtr->leafGridView();
+
+#endif
+
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ const auto& cellCenterSol = sol[stokesCellCenterIdx];
+ const auto& faceSol = sol[stokesFaceIdx];
+
+ // apply initial solution for instationary problems
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ std::get<0>(stokesSol) = cellCenterSol;
+ std::get<1>(stokesSol) = faceSol;
+ stokesProblem->applyInitialSolution(stokesSol);
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // intialize the vtk output module
+ const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+ StaggeredVtkOutputModule<StokesTypeTag> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+
+#if EXNUMBER >= 2
+ stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+#endif
+
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ 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(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ 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);
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // write vtk output
+ stokesVtkWriter.write(1.0);
+ darcyVtkWriter.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
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
new file mode 100644
index 0000000000000000000000000000000000000000..e3813950725ba4837e96b79d76dc9ccc74dc0d4f
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
@@ -0,0 +1,45 @@
+ # for dune-subgrid
+ [Grid]
+Positions0 = 0 1
+Positions1 = 0 0.2 0.3 0.65
+Cells0 = 100
+Cells1 = 10 50 18
+Baseline = 0.25 # [m]
+Amplitude = 0.04 # [m]
+Offset = 0.5 # [m]
+Scaling = 0.2 #[m]
+
+[Stokes.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 1.0 2.0
+Cells0 = 20
+Cells1 = 100
+Grading1 = 1
+
+[Darcy.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 0.0 1.0
+Cells0 = 20
+Cells1 = 20
+Grading1 = 1
+
+[Stokes.Problem]
+Name = stokes
+PressureDifference = 1e-9
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+Permeability = 1e-6 # m^2
+Porosity = 0.4
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+Name = ex_ff-pm-interface
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..369a7bd2b9238272a36759a1856d4f7369f55a37
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
@@ -0,0 +1,360 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief The free flow sub problem
+ */
+#ifndef DUMUX_STOKES_SUBPROBLEM_HH
+#define DUMUX_STOKES_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#if EXNUMBER >= 3
+#include <dumux/io/grid/subgridgridcreator.hh>
+#endif
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/simpleh2o.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(StokesOnePTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
+
+// the fluid system
+SET_PROP(StokesOnePTypeTag, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+SET_PROP(StokesOnePTypeTag, Grid)
+{
+ static constexpr auto dim = 2;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+#if EXNUMBER < 3 // use "normal" grid
+ using type = TensorGrid;
+#else // use dune-subgrid
+ using HostGrid = TensorGrid;
+ using type = Dune::SubGrid<dim, HostGrid>;
+#endif
+};
+
+// Set the problem property
+SET_TYPE_PROP(StokesOnePTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
+
+SET_BOOL_PROP(StokesOnePTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(StokesOnePTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(StokesOnePTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(StokesOnePTypeTag, EnableInertiaTerms, false);
+}
+
+/*!
+ * \brief The free flow sub problem
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+
+public:
+ StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return 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
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.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.dofPosition();
+
+#if EXNUMBER == 0 // flow from top to bottom
+ if(onUpperBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+
+ if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+#else // flow flom left to right
+ if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
+ values.setDirichlet(Indices::pressureIdx);
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+#endif
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+#if EXNUMBER < 3
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+#else
+ //consider orientation of face automatically
+ values.setCouplingNeumann(scvf.directionIndex());
+#endif
+
+#if EXNUMBER == 2
+ // set the Beaver-Joseph-Saffman slip condition for the tangential momentum balance equation
+ values.setBJS(Indices::momentumXBalanceIdx);
+#endif
+
+#if EXNUMBER == 3
+ // set the Beaver-Joseph-Saffman slip condition for the tangential momentum balance equation,
+ // consider orientation of face automatically
+ values.setBJS(1 - scvf.directionIndex());
+#endif
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(globalPos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry 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 ElementFaceVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+#if EXNUMBER < 3
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+#else
+ values[scvf.directionIndex()] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+#endif
+
+ }
+ return values;
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+#if EXNUMBER == 0
+ values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]);
+#else
+ // set fixed pressures on the left and right boundary
+ if(onLeftBoundary_(globalPos))
+ values[Indices::pressureIdx] = deltaP_;
+ if(onRightBoundary_(globalPos))
+ values[Indices::pressureIdx] = 0.0;
+#endif
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ /*!
+ * \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967)
+ */
+ void calculateAnalyticalVelocityX() const
+ {
+ analyticalVelocityX_.resize(this->fvGridGeometry().gridView().size(0));
+
+ using std::sqrt;
+ const Scalar dPdX = -deltaP_ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]);
+ static const Scalar mu = FluidSystem::viscosity(temperature(), 1e5);
+ 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 = (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().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(this->fvGridGeometry().gridView()))
+ {
+ const auto eIdx = this->fvGridGeometry().gridView().indexSet().index(element);
+ const Scalar y = element.geometry().center()[1] - this->fvGridGeometry().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;
+ }
+ }
+
+ /*!
+ * \brief Get the analytical velocity in x direction
+ */
+ const std::vector<Scalar>& getAnalyticalVelocityX() const
+ {
+ if(analyticalVelocityX_.empty())
+ calculateAnalyticalVelocityX();
+ return analyticalVelocityX_;
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+ Scalar deltaP_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ mutable std::vector<Scalar> analyticalVelocityX_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES_SUBPROBLEM_HH
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7a5c4baaa4cdcc2997ed77781f1ce5308f8718ef
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
@@ -0,0 +1,257 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief The porous medium flow sub problem
+ */
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#if EXNUMBER >= 3
+#include <dumux/io/grid/subgridgridcreator.hh>
+#endif
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "../1pspatialparams.hh"
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(DarcyOnePTypeTag, INHERITS_FROM(CCTpfaModel, OneP));
+
+// Set the problem property
+SET_TYPE_PROP(DarcyOnePTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// the fluid system
+SET_PROP(DarcyOnePTypeTag, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+SET_PROP(DarcyOnePTypeTag, Grid)
+{
+ static constexpr auto dim = 2;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+#if EXNUMBER < 3 // use "normal" grid
+ using type = TensorGrid;
+#else // use dune-subgrid
+ using HostGrid = TensorGrid;
+ using type = Dune::SubGrid<dim, HostGrid>;
+#endif
+};
+
+SET_TYPE_PROP(DarcyOnePTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+}
+
+/*!
+ * \brief The porous medium flow sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {}
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Return 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 EXNUMBER == 0 // flow from top to bottom
+ if (onLowerBoundary_(scvf.center()))
+ values.setAllDirichlet();
+#endif
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ PrimaryVariables values(0.0);
+ values = initial(element);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ * For this method, the \a values variable stores primary variables.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param element The element
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initial(const Element &element) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f7e7af5fdf7667b674cdf987c20a7f23985487e5
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
@@ -0,0 +1,24 @@
+add_input_file_links()
+
+dune_add_test(NAME orig_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=0
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_a_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=1
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_b_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=2
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_c_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=3
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+# add exercise to the common target
+add_dependencies(test_exercises sol_a_ex_models_coupling_ff-pm sol_b_ex_models_coupling_ff-pm sol_b_ex_models_coupling_ff-pm)
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cdb7abdd800c9994ade2dab51adda2076ed7a51d
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
@@ -0,0 +1,311 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the coupled Stokes/Darcy problem (1p)
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+#include <fstream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/geometry/diameter.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#if EXNUMBER >= 3
+#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+#else
+#include <dumux/multidomain/newtonsolver.hh>
+#endif
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_models_pmproblem.hh"
+#include "ex_models_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+SET_PROP(StokesTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTypeTag)>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+SET_PROP(DarcyTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TTAG(StokesTypeTag), TTAG(StokesTypeTag), TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ 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 StokesTypeTag = TTAG(StokesTypeTag);
+ using DarcyTypeTag = TTAG(DarcyTypeTag);
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+ stokesProblem->setTimeLoop(timeLoop);
+ darcyProblem->setTimeLoop(timeLoop);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ const auto& cellCenterSol = sol[stokesCellCenterIdx];
+ const auto& faceSol = sol[stokesFaceIdx];
+
+ // apply initial solution for instationary problems
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ std::get<0>(stokesSol) = cellCenterSol;
+ std::get<1>(stokesSol) = faceSol;
+ stokesProblem->applyInitialSolution(stokesSol);
+ auto solStokesOld = stokesSol;
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+ auto solDarcyOld = sol[darcyIdx];
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol, solStokesOld);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+
+ // intialize the vtk output module
+#if EXNUMBER >= 1
+ const std::array<std::string, 3> part = {"a", "b", "c"};
+ const auto stokesName = "sol_" + part[EXNUMBER-1] + "_" + getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = "sol_" + part[EXNUMBER-1] + "_" + getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+#else
+ const auto stokesName = "orig_" + getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = "orig_" + getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+#endif
+
+ StaggeredVtkOutputModule<StokesTypeTag> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // intialize the subproblems
+ darcyProblem->init(sol[darcyIdx], *darcyGridVariables);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the primary variable switches used by the sub models and the non-linear solver
+#if EXNUMBER >= 3
+ using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
+ using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+#else
+// using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, NoPrimaryVariableSwitch>;
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+#endif
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ const auto episodeLength = getParam<Scalar>("TimeLoop.EpisodeLength");
+ if (episodeLength > 0.0)
+ timeLoop->setPeriodicCheckPoint(episodeLength);
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(solOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(sol, *timeLoop);
+
+ // make the new solution the old solution
+ solOld = sol;
+ stokesGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // call the postTimeStep routine for output
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables);
+
+ // write vtk output
+ if (timeLoop->isCheckPoint() || timeLoop->finished() || episodeLength < 0.0)
+ {
+ stokesVtkWriter.write(timeLoop->time());
+ darcyVtkWriter.write(timeLoop->time());
+ }
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(stokesGridView.comm());
+ timeLoop->finalize(darcyGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // 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
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
new file mode 100644
index 0000000000000000000000000000000000000000..8d08fe836a2682e36007bbafa51c6438a5d15c6e
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
@@ -0,0 +1,49 @@
+[TimeLoop]
+DtInitial = 100 # s
+EpisodeLength = -360 # s # 0.25 days
+TEnd = 256000 # s # 2 days
+
+[Stokes.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 16 16
+
+[Darcy.Grid]
+UpperRight = 1 1
+Cells = 16 16
+
+[Stokes.Problem]
+Name = stokes
+EnableGravity = false
+MoleFraction = 0.0 # -
+Pressure = 1e5 # Pa
+Velocity = 1e-3 # m/s
+
+[Darcy.Problem]
+Name = darcy
+EnableGravity = true
+Saturation = 0.1 # -
+MoleFraction = 0.1 # -
+Pressure = 1e5 # Pa
+
+[Darcy.SpatialParams]
+Permeability = 2.65e-10 # m^2
+Porosity = 0.4 # -
+AlphaBeaversJoseph = 1.0 # -
+# EXNUMBER >= 1
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.5e-4
+VgN = 8.0
+
+[Problem]
+Name = ex_models_coupling
+PlotFluxes = false
+PlotStorage = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = 1
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3d59b25a071ed0e71a0c3baac902b4847a2ce840
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
@@ -0,0 +1,344 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesTests
+ * \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model.
+ */
+#ifndef DUMUX_STOKES1P2C_SUBPROBLEM_HH
+#define DUMUX_STOKES1P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1padapter.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(StokesTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+
+// Set the grid type
+SET_TYPE_PROP(StokesTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// The fluid system
+SET_PROP(StokesTypeTag, FluidSystem)
+{
+ using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
+};
+
+// Do not replace one equation with a total mass balance
+SET_INT_PROP(StokesTypeTag, ReplaceCompEqIdx, 3);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(StokesTypeTag, UseMoles, true);
+
+// Set the problem property
+SET_TYPE_PROP(StokesTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
+
+SET_BOOL_PROP(StokesTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(StokesTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(StokesTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(StokesTypeTag, EnableInertiaTerms, false);
+}
+
+/*!
+ * \ingroup NavierStokesTests
+ * \brief Test problem for the one-phase compositional (Navier-) Stokes problem.
+ *
+ * Horizontal flow from left to right with a parabolic velocity profile.
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+
+ static constexpr auto dim = GET_PROP_TYPE(TypeTag, ModelTraits)::dim();
+ static constexpr auto transportCompIdx = 1;
+
+public:
+ StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ velocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+ moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return 293.15; }
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.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();
+
+ if(onLeftBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::conti0EqIdx + 1);
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+ else if(onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx + 1);
+ }
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element
+ * \param scvf The subcontrolvolume face
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(pos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ }
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, pressure_);
+ const Scalar density = FluidSystem::density(fluidState, 0);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = pressure_ + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[Indices::conti0EqIdx + 1] = moleFraction_;
+ values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1])
+ * (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
+ / (height_() * height_());
+
+ return values;
+ }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ //! \brief updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState, const Scalar pressure) const
+ {
+ fluidState.setTemperature(temperature());
+ fluidState.setPressure(0, pressure);
+ fluidState.setSaturation(0, 1.0);
+ fluidState.setMoleFraction(0, 1, moleFraction_);
+ fluidState.setMoleFraction(0, 0, 1.0 - moleFraction_);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, 0);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+ fluidState.setDensity(0, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+ fluidState.setMolarDensity(0, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+ fluidState.setEnthalpy(0, enthalpy);
+ }
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar velocity_;
+ Scalar pressure_;
+ Scalar moleFraction_;
+
+ TimeLoopPtr timeLoop_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..18ad042dfdb65cbfc9a61389f8a7084a21e34eed
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
@@ -0,0 +1,471 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A simple Darcy test problem (cell-centered finite volume method).
+ */
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1padapter.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+
+#if EXNUMBER >= 1
+#include <dumux/porousmediumflow/2pnc/model.hh>
+#include "../2pspatialparams.hh"
+#else
+#include <dumux/porousmediumflow/1pnc/model.hh>
+#include "../1pspatialparams.hh"
+#endif
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+#if EXNUMBER >= 1
+NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, TwoPNC));
+#else
+NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
+#endif
+
+// Set the problem property
+SET_TYPE_PROP(DarcyTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// The fluid system
+SET_PROP(DarcyTypeTag, FluidSystem)
+{
+ using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+#if EXNUMBER == 0
+ using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
+#else
+ using type = H2OAir;
+#endif
+};
+
+// Use moles
+SET_BOOL_PROP(DarcyTypeTag, UseMoles, true);
+
+// Do not replace one equation with a total mass balance
+SET_INT_PROP(DarcyTypeTag, ReplaceCompEqIdx, 3);
+
+//! Use a model with constant tortuosity for the effective diffusivity
+SET_TYPE_PROP(DarcyTypeTag, EffectiveDiffusivityModel,
+ DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set the grid type
+SET_TYPE_PROP(DarcyTypeTag, Grid, Dune::YaspGrid<2>);
+
+#if EXNUMBER >= 1
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+SET_PROP(DarcyTypeTag, Formulation)
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+#endif
+
+// Set the spatial paramaters type
+#if EXNUMBER >= 1
+SET_TYPE_PROP(DarcyTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+#else
+SET_TYPE_PROP(DarcyTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+#endif
+}
+
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ enum {
+ // grid and world dimension
+ dim = GridView::dimension,
+ dimworld = GridView::dimensionworld,
+
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ pressureIdx = Indices::pressureIdx,
+#if EXNUMBER >= 3
+ saturationIdx = Indices::switchIdx,
+ transportCompIdx = Indices::switchIdx
+#elif EXNUMBER >= 1
+ transportCompIdx = Indices::switchIdx
+#else
+ phaseIdx = 0,
+ transportCompIdx = 1
+#endif
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+#if EXNUMBER >= 3
+ saturation_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
+#else
+ moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
+#endif
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+
+ // initialize output file
+ plotFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotFluxes", false);
+ plotStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotStorage", false);
+ storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+ storageFile_.open(storageFileName_);
+ storageFile_ << "#Time[s]" << ";"
+ << "WaterMass[kg]" << ";"
+ << "WaterMassLoss[kg]" << ";"
+ << "EvaporationRate[mm/d]"
+ << std::endl;
+ }
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Initialize the problem.
+ */
+ template<class SolutionVector, class GridVariables>
+ void init(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+ {
+#if EXNUMBER >= 2
+ initialWaterContent_ = evaluateWaterMassStorageTerm(curSol, gridVariables);
+ lastWaterMass_ = initialWaterContent_;
+#endif
+ }
+
+ template<class SolutionVector, class GridVariables>
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ evaluateWaterMassStorageTerm(curSol, gridVariables);
+ evaluateInterfaceFluxes(curSol, gridVariables);
+
+ gnuplotStorage_.resetPlot();
+ gnuplotStorage_.setDatafileSeparator(';');
+ gnuplotStorage_.setXlabel("time [d]");
+ gnuplotStorage_.setXRange(0.0, getParam<Scalar>("TimeLoop.TEnd"));
+ gnuplotStorage_.setYlabel("evaporation rate [mm/d]");
+ gnuplotStorage_.setOption("set yrange [0.0:]");
+ gnuplotStorage_.setOption("set y2label 'cumulative mass loss'");
+ gnuplotStorage_.setOption("set y2range [0.0:0.5]");
+ gnuplotStorage_.setOption("set y2range [0.0:0.5]");
+ gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:4 with lines title 'evaporation rate'");
+ gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:3 axes x1y2 with lines title 'cumulative mass loss'");
+ if (plotStorage_)
+ gnuplotStorage_.plot("temp");
+ }
+
+ template<class SolutionVector, class GridVariables>
+ Scalar evaluateWaterMassStorageTerm(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ // compute the mass in the entire domain
+ Scalar waterMass = 0.0;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
+ {
+ // insert calculation of the water mass here
+#if EXNUMBER >= 2
+ waterMass += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx) * volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx) * volVars.porosity()
+ * scv.volume() * volVars.extrusionFactor();
+#else
+ waterMass += 0.0;
+#endif
+ }
+ }
+ }
+#if EXNUMBER >= 2
+ std::cout << "Mass of water is: " << waterMass << std::endl;
+#endif
+
+ Scalar cumMassLoss = initialWaterContent_ - waterMass;
+ Scalar evaporationRate = (lastWaterMass_ - waterMass) * 86400
+ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0])
+ / timeLoop_->timeStepSize();
+ lastWaterMass_ = waterMass;
+
+ storageFile_ << timeLoop_->time() << ";"
+ << waterMass << ";"
+ << cumMassLoss << ";"
+ << evaporationRate
+ << std::endl;
+
+ return waterMass;
+ }
+
+ template<class SolutionVector, class GridVariables>
+ void evaluateInterfaceFluxes(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ std::vector<Scalar> x;
+ std::vector<Scalar> y;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ continue;
+
+ // NOTE: binding the coupling context is necessary
+ couplingManager_->bindCouplingContext(CouplingManager::darcyIdx, element);
+#if EXNUMBER >= 2
+ NumEqVector flux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+#else
+ NumEqVector flux(0.0); // add "massCouplingCondition" from the couplingManager here
+#endif
+
+ x.push_back(scvf.center()[0]);
+ y.push_back(flux[transportCompIdx]);
+ }
+ }
+
+ gnuplotInterfaceFluxes_.resetPlot();
+ gnuplotInterfaceFluxes_.setXlabel("x-position [m]");
+ gnuplotInterfaceFluxes_.setXRange(this->fvGridGeometry().bBoxMin()[0], this->fvGridGeometry().bBoxMax()[0]);
+ gnuplotInterfaceFluxes_.setYlabel("flux [kg/(m^2 s)]");
+ gnuplotInterfaceFluxes_.setYRange(-5e-4, 0.0);
+ gnuplotInterfaceFluxes_.setOption("set label 'time: " + std::to_string(timeLoop_->time()/86400.) + "d' at graph 0.8,0.8 ");
+ std::string fluxFileName = "flux_" + std::to_string(timeLoop_->timeStepIndex()) +
+ "_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+ gnuplotInterfaceFluxes_.addDataSetToPlot(x, y, fluxFileName, "with lines title 'water mass flux'");
+ if (plotFluxes_)
+ gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex()));
+ }
+
+ /*!
+ * \brief Returns true if a restart file should be written to
+ * disk.
+ */
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ bool shouldWriteOutput() const // define output
+ { return true; }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 293.15; }
+ // \}
+
+ /*!
+ * \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().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ * For this method, the \a values variable stores primary variables.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param element The element
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+#if EXNUMBER >= 3
+ values.setState(3/*bothPhases*/);
+ values[saturationIdx] = saturation_;
+#elif EXNUMBER >= 1
+ values.setState(2/*secondPhaseOnly*/);
+ values[transportCompIdx] = moleFraction_;
+#else
+ values[transportCompIdx] = moleFraction_;
+#endif
+ values[pressureIdx] = pressure_;
+ return values;
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+#if EXNUMBER >= 3
+ Scalar saturation_;
+#else
+ Scalar moleFraction_;
+#endif
+ Scalar pressure_;
+
+ Scalar initialWaterContent_ = 0.0;
+ Scalar lastWaterMass_ = 0.0;
+
+ TimeLoopPtr timeLoop_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ std::string storageFileName_;
+ std::ofstream storageFile_;
+ bool plotFluxes_;
+ bool plotStorage_;
+ Dumux::GnuplotInterface<Scalar> gnuplotInterfaceFluxes_;
+ Dumux::GnuplotInterface<Scalar> gnuplotStorage_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e19d81c519756532c8f6fd2b7f28926723f56952
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
@@ -0,0 +1,19 @@
+add_input_file_links()
+
+dune_add_test(NAME orig_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=0
+ CMD_ARGS ex_turbulence_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_a_ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=1
+ CMD_ARGS ex_turbulence_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_b_ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=2
+ CMD_ARGS ex_turbulence_coupling_ff-pm.input)
+
+# add exercise to the common target
+add_dependencies(test_exercises sol_a_ex_turbulence_coupling_ff-pm sol_b_ex_turbulence_coupling_ff-pm sol_b_ex_turbulence_coupling_ff-pm)
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9b02336e65a9ac918d99a2cd34c6353d66e258d7
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
@@ -0,0 +1,306 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the isothermal coupled Stokes/Darcy problem (1p2c/2p2c)
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+#include <fstream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/geometry/diameter.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_turbulence_pmproblem.hh"
+#include "ex_turbulence_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+SET_PROP(ZeroEqTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTwoPTwoCTypeTag)>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+SET_PROP(DarcyTwoPTwoCTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TTAG(ZeroEqTypeTag), TTAG(ZeroEqTypeTag), TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ 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 StokesTypeTag = TTAG(ZeroEqTypeTag);
+ using DarcyTypeTag = TTAG(DarcyTwoPTwoCTypeTag);
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // set timeloop for the subproblems, needed for boundary value variations
+ stokesProblem->setTimeLoop(timeLoop);
+ darcyProblem->setTimeLoop(timeLoop);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ // apply initial solution for instationary problems
+ // auxiliary free flow solution vector
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ stokesSol[stokesCellCenterIdx].resize(sol[stokesCellCenterIdx].size());
+ stokesSol[stokesFaceIdx].resize(sol[stokesFaceIdx].size());
+ stokesProblem->applyInitialSolution(stokesSol);
+ auto solStokesOld = stokesSol;
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+ // TODO: update static wall properties
+ // TODO: update dynamic wall properties
+#if EXNUMBER >= 1
+ stokesProblem->updateStaticWallProperties();
+ stokesProblem->updateDynamicWallProperties(stokesSol);
+#endif
+
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+ auto solDarcyOld = sol[darcyIdx];
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol, solStokesOld);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+
+ // intialize the vtk output module
+ const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+ StaggeredVtkOutputModule<StokesTypeTag, GET_PROP_VALUE(StokesTypeTag, PhaseIdx)> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the primary variable switches used by the sub models
+ using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(solOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(sol, *timeLoop);
+
+ // make the new solution the old solution
+ solOld = sol;
+
+ // update the auxiliary free flow solution vector
+ stokesSol[stokesCellCenterIdx] = sol[stokesCellCenterIdx];
+ stokesSol[stokesFaceIdx] = sol[stokesFaceIdx];
+
+#if EXNUMBER >= 1
+ // TODO: update dynamic wall properties
+ stokesProblem->updateDynamicWallProperties(stokesSol);
+#endif
+
+ // post time step treatment of Darcy problem
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables, timeLoop->timeStepSize());
+
+ // advance grid variables to the next time step
+ stokesGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ stokesVtkWriter.write(timeLoop->time());
+ darcyVtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(stokesGridView.comm());
+ timeLoop->finalize(darcyGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // 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
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input
new file mode 100644
index 0000000000000000000000000000000000000000..514cf156d55211777c2ec0980f37d1a95124867e
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input
@@ -0,0 +1,62 @@
+[TimeLoop]
+DtInitial = 1e-1 # [s]
+MaxTimeStepSize = 43200 # [s] (12 hours)
+TEnd = 864000 # [s] (6 days)
+
+[Stokes.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.25 0.5
+Grading0 = 1.0
+Grading1 = 1.0
+Cells0 = 15
+Cells1 = 20
+Verbosity = true
+
+[Darcy.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.0 0.25
+Cells0 = 15
+Cells1 = 10
+Grading0 = 1.0
+Grading1 = 1.0
+Verbosity = true
+
+[Stokes.Problem]
+Name = stokes
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+refMoleFrac = 0 # [-]
+RefTemperature = 298.15 # [K]
+
+[Darcy.Problem]
+Name = darcy
+Pressure = 1e5
+Saturation = 0.5 # initial Sw
+Temperature = 298.15 # [K]
+InitPhasePresence = 3 # bothPhases
+
+[Darcy.SpatialParams]
+Porosity = 0.41
+Permeability = 2.65e-10
+AlphaBeaversJoseph = 1.0
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.371e-4
+VgN = 6.9
+
+[Problem]
+Name = ex_coupling_turbulence_ff-pm
+EnableGravity = true
+InterfaceDiffusionCoefficientAvg = Harmonic
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Assembly]
+NumericDifferenceMethod = 0
+NumericDifference.BaseEpsilon = 1e-8
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b02930b9ba255a3cd6982de4dc102e21c2635780
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
@@ -0,0 +1,409 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+ /*!
+ * \file
+ * \brief The free-flow sub problem
+ */
+#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1padapter.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+
+#if EXNUMBER >= 1
+#include <dumux/freeflow/compositional/zeroeqncmodel.hh>
+#include <dumux/freeflow/rans/zeroeq/problem.hh>
+#else
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+#include <dumux/freeflow/navierstokes/problem.hh>
+#endif
+
+namespace Dumux
+{
+template <class TypeTag>
+class FreeFlowSubProblem;
+
+namespace Properties
+{
+#if EXNUMBER >= 1
+NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNCNI));
+#else
+NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
+#endif
+
+// Set the grid type
+SET_TYPE_PROP(ZeroEqTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// The fluid system
+SET_PROP(ZeroEqTypeTag, FluidSystem)
+{
+ using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
+ using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
+
+SET_INT_PROP(ZeroEqTypeTag, ReplaceCompEqIdx, 3);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(ZeroEqTypeTag, UseMoles, true);
+
+// Set the problem property
+SET_TYPE_PROP(ZeroEqTypeTag, Problem, Dumux::FreeFlowSubProblem<TypeTag> );
+
+SET_BOOL_PROP(ZeroEqTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(ZeroEqTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(ZeroEqTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(ZeroEqTypeTag, EnableInertiaTerms, true);
+}
+
+/*!
+ * \brief The free-flow sub problem
+ */
+template <class TypeTag>
+#if EXNUMBER >= 1
+class FreeFlowSubProblem : public ZeroEqProblem<TypeTag>
+{
+ using ParentType = ZeroEqProblem<TypeTag>;
+#else
+class FreeFlowSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+#endif
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+ static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+
+public:
+ FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ refVelocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefVelocity");
+ refPressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefPressure");
+ refMoleFrac_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.refMoleFrac");
+ refTemperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefTemperature");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return refTemperature_; }
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.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();
+
+ if (onLeftBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setDirichlet(Indices::conti0EqIdx + 1);
+ values.setDirichlet(Indices::energyBalanceIdx);
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyBalanceIdx);
+ }
+
+ if (onUpperBoundary_(globalPos))
+ {
+#if EXNUMBER >=2
+ values.setAllSymmetry();
+#else
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyBalanceIdx);
+#endif
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx + 1);
+ values.setOutflow(Indices::energyBalanceIdx);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ values.setCouplingNeumann(Indices::energyBalanceIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element
+ * \param scvf The subcontrolvolume face
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(pos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ values[Indices::energyBalanceIdx] = couplingManager().couplingData().energyCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ }
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+#if EXNUMBER >= 2
+ bool isOnWall(const GlobalPosition& globalPos) const
+ {
+ return (onLowerBoundary_(globalPos));
+ }
+#elif EXNUMBER >= 1
+ bool isOnWall(const GlobalPosition& globalPos) const
+ {
+ return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+ }
+#endif
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, refTemperature(), refPressure(), refMoleFrac());
+
+ const Scalar density = FluidSystem::density(fluidState, 0);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = refPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[Indices::conti0EqIdx + 1] = refMoleFrac();
+ values[Indices::velocityXIdx] = refVelocity();
+ values[Indices::temperatureIdx] = refTemperature();
+
+ if(onUpperBoundary_(globalPos))
+ values[Indices::velocityXIdx] = 0.0;
+
+ return values;
+ }
+
+ //! \brief Returns the reference velocity.
+ const Scalar refVelocity() const
+ { return refVelocity_ ;}
+
+ //! \brief Returns the reference pressure.
+ const Scalar refPressure() const
+ { return refPressure_; }
+
+ //! \brief Returns the reference mass fraction.
+ const Scalar refMoleFrac() const
+ { return refMoleFrac_; }
+
+ //! \brief Returns the reference temperature.
+ const Scalar refTemperature() const
+ { return refTemperature_; }
+
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ //! \brief updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState, const Scalar temperature,
+ const Scalar pressure, const Scalar moleFraction) const
+ {
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(0, pressure);
+ fluidState.setSaturation(0, 1.0);
+ fluidState.setMoleFraction(0, 1, moleFraction);
+ fluidState.setMoleFraction(0, 0, 1.0 - moleFraction);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, 0);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+ fluidState.setDensity(0, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+ fluidState.setMolarDensity(0, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+ fluidState.setEnthalpy(0, enthalpy);
+ }
+
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar refVelocity_;
+ Scalar refPressure_;
+ Scalar refMoleFrac_;
+ Scalar refTemperature_;
+
+ TimeLoopPtr timeLoop_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bd589d6303f95b7c9ba96b94102174af24c8a129
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
@@ -0,0 +1,324 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief The porous medium sub problem
+ */
+#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH
+#define DUMUX_DARCY2P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+
+#include <dumux/porousmediumflow/2p2c/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+#include "../2pspatialparams.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(DarcyTwoPTwoCTypeTag, INHERITS_FROM(CCTpfaModel, TwoPTwoCNI));
+
+// Set the problem property
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// the fluid system
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+SET_PROP(DarcyTwoPTwoCTypeTag, Formulation)
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+
+// The gas component balance (air) is replaced by the total mass balance
+SET_INT_PROP(DarcyTwoPTwoCTypeTag, ReplaceCompEqIdx, 3);
+
+// Set the grid type
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+SET_BOOL_PROP(DarcyTwoPTwoCTypeTag, UseMoles, true);
+
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+}
+
+/*!
+ * \brief The porous medium sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ enum {
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ contiWEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
+ contiNEqIdx = Indices::conti0EqIdx + FluidSystem::AirIdx,
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+ initialSw_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
+ temperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Temperature");
+ initialPhasePresence_ = getParamFromGroup<int>(this->paramGroup(), "Problem.InitPhasePresence");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ }
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ template<class SolutionVector, class GridVariables>
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables,
+ const Scalar timeStepSize)
+
+ {
+ // compute the mass in the entire domain
+ Scalar massWater = 0.0;
+
+ // bulk elements
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
+ {
+ massWater += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx)*volVars.density(phaseIdx)
+ * scv.volume() * volVars.saturation(phaseIdx) * volVars.porosity() * volVars.extrusionFactor();
+ }
+ }
+ }
+
+ std::cout << "mass of water is: " << massWater << std::endl;
+ }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return temperature_; }
+ // \}
+
+ /*!
+ * \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().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(scvf.center());
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann
+ * control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ {
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+
+ for(int i = 0; i< massFlux.size(); ++i)
+ values[i] = massFlux[i];
+
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+ }
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ *
+ * For this method, the \a values variable stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values.setState(initialPhasePresence_);
+
+ values[pressureIdx] = pressure_ + 1. * this->gravity()[1] * (globalPos[1] - this->fvGridGeometry().bBoxMax()[1]);
+ values[switchIdx] = initialSw_;
+ values[Indices::temperatureIdx] = temperature_;
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar pressure_;
+ Scalar initialSw_;
+ Scalar temperature_;
+ int initialPhasePresence_;
+
+ TimeLoopPtr timeLoop_;
+
+ Scalar eps_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY2P2C_SUBPROBLEM_HH