diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
index ec4e22aa1ec9038d3aeaa6a7bd6fa0309f14be20..f3d56c7e2e9e292510d60db415d9233447f88259 100644
--- a/exercises/exercise-coupling-ff-pm/README.md
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -269,7 +269,7 @@ Finally we want to know the distribution of the water mass fluxes across the int
   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
+* You can use the parameter `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.
@@ -295,10 +295,6 @@ for this case the liquid saturation cannot serve as primary variable anymore. Ho
 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.
 
@@ -355,7 +351,7 @@ However, there is also a solution-dependent component of these interactions, e.g
 damping of the eddy viscosity toward the wall, the velocity gradient at the wall and inside the
 cells is needed.
 These dynamic interactions are to be updated by calling  
-```cpp 
+```cpp
 stokesProblem->updateDynamicWallProperties(stokesSol)
 ```
 in the time loop (after `// update dynamic wall properties`).
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
index b54ae19ed5f17fa1d5bcc47a0d3e01f1c5f502ed..8dc210c527a0f39bc03a19b311e228b4c272a70d 100644
--- 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
@@ -34,11 +34,11 @@
 #include <dumux/common/properties.hh>
 #include <dumux/common/parameters.hh>
 #include <dumux/common/dumuxmessage.hh>
-#include <dumux/common/geometry/diameter.hh>
+#include <dumux/common/partial.hh>
 #include <dumux/linear/seqsolverbackend.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
-#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/method.hh>
 #include <dumux/io/vtkoutputmodule.hh>
 #include <dumux/io/staggeredvtkoutputmodule.hh>
 #include <dumux/io/grid/gridmanager.hh>
@@ -55,15 +55,17 @@
 namespace Dumux {
 namespace Properties {
 
-SET_PROP(StokesTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTypeTag)>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
-SET_PROP(DarcyTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TTAG(StokesTypeTag), TTAG(StokesTypeTag), TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesTypeTag, Properties::TTag::StokesTypeTag, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
@@ -85,16 +87,16 @@ int main(int argc, char** argv) try
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = TTAG(StokesTypeTag);
-    using DarcyTypeTag = TTAG(DarcyTypeTag);
+    using StokesTypeTag = Properties::TTag::StokesTypeTag;
+    using DarcyTypeTag = Properties::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)>;
+    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
     DarcyGridManager darcyGridManager;
     darcyGridManager.init("Darcy"); // pass parameter group
 
-    using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+    using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
     StokesGridManager stokesGridManager;
     stokesGridManager.init("Stokes"); // pass parameter group
 
@@ -103,10 +105,10 @@ int main(int argc, char** argv) try
     const auto& stokesGridView = stokesGridManager.grid().leafGridView();
 
     // create the finite volume grid geometry
-    using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+    using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
     auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
     stokesFvGridGeometry->update();
-    using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>;
     auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
     darcyFvGridGeometry->update();
 
@@ -122,27 +124,22 @@ int main(int argc, char** argv) try
     constexpr auto darcyIdx = CouplingManager::darcyIdx;
 
     // the problem (initial and boundary conditions)
-    using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+    using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
     auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
-    using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
     auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
 
     // initialize the fluidsystem (tabulation)
-    GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+    GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
 
     // get some time loop parameters
-    using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+    using Scalar = GetPropType<StokesTypeTag, Properties::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);
+    auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
     stokesProblem->setTimeLoop(timeLoop);
     darcyProblem->setTimeLoop(timeLoop);
@@ -153,43 +150,35 @@ int main(int argc, char** argv) try
     sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
     sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
 
-    const auto& cellCenterSol = sol[stokesCellCenterIdx];
-    const auto& faceSol = sol[stokesFaceIdx];
+    auto stokesSol = partial(sol, stokesCellCenterIdx, 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);
+    using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
     auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
-    stokesGridVariables->init(stokesSol, solStokesOld);
-    using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+    stokesGridVariables->init(stokesSol);
+    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
     auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
-    darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+    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);
+    StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+    GetPropType<StokesTypeTag, Properties::VtkOutputFields>::initOutputModule(stokesVtkWriter);
     stokesVtkWriter.write(0.0);
 
-    VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
-    GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
     darcyVtkWriter.write(0.0);
 
     // intialize the subproblems
@@ -211,8 +200,6 @@ int main(int argc, char** argv) try
     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);
 
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
index d3f68518fa9a38504ee876af35f0c9e79f31cf17..07bb4dd0e86001358836ac361a1399f344add5d5 100644
--- 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
@@ -1,5 +1,5 @@
 [TimeLoop]
-DtInitial = 100 # s
+DtInitial = 1 # s
 EpisodeLength = -360 # s # 0.25 days
 TEnd = 256000 # s # 2 days
 
@@ -14,10 +14,11 @@ Cells = 16 16
 
 [Stokes.Problem]
 Name = stokes
-EnableGravity = false
+EnableGravity = true
+EnableInertiaTerms = true
 MoleFraction = 0.0 # -
 Pressure = 1e5 # Pa
-Velocity = 1e-3 # m/s
+Velocity = 1 # m/s
 
 [Darcy.Problem]
 Name = darcy
@@ -43,7 +44,11 @@ PlotStorage = false
 
 [Newton]
 MaxSteps = 12
-MaxRelativeShift = 1e-5
+MaxRelativeShift = 1e-7
+TargetSteps = 5
 
 [Vtk]
 AddVelocity = 1
+
+[Assembly]
+NumericDifference.BaseEpsilon = 1e-8
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh b/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
index 1fe98c085805a45045474cece286db0f4561b399..6f00aa101770dcd2c1e0e4f0a5b3cd5c874f99e0 100644
--- a/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
@@ -40,32 +40,41 @@ class StokesSubProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(StokesTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+// Create new type tags
+namespace TTag {
+struct StokesTypeTag { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+} // end namespace TTag
 
 // Set the grid type
-SET_TYPE_PROP(StokesTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesTypeTag> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
 // The fluid system
-SET_PROP(StokesTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesTypeTag>
 {
-    using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+    using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
     using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
 };
 
 // Do not replace one equation with a total mass balance
-SET_INT_PROP(StokesTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::StokesTypeTag> { static constexpr int value = 3; };
 
 // Use formulation based on mass fractions
-SET_BOOL_PROP(StokesTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = 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);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesTypeTag> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
 }
 
 /*!
@@ -79,29 +88,29 @@ 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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
 
-    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVGridGeometry = GetPropType<TypeTag, Properties::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 ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+    using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
+    using FluidState = GetPropType<TypeTag, Properties::FluidState>;
 
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
-    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
-    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
     using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
 
-    static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+    static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
 
 public:
     StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
@@ -117,10 +126,6 @@ public:
      */
     // \{
 
-
-    bool shouldWriteRestartFile() const
-    { return false; }
-
    /*!
      * \brief Return the temperature within the domain in [K].
      */
@@ -218,9 +223,9 @@ public:
 
         if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
         {
-            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
 
-            const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+            const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
             values[Indices::conti0EqIdx] = massFlux[0];
             values[Indices::conti0EqIdx + 1] = massFlux[1];
         }
@@ -251,14 +256,22 @@ public:
      *
      * \param globalPos The global position
      */
-    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
     {
+        PrimaryVariables values(0.0);
+
+        // This is only an approximation of the actual hydorostatic pressure gradient.
+        // Air is compressible (the density depends on pressure), thus the actual
+        // vertical pressure profile is non-linear.
+        // This discrepancy can lead to spurious flows at the outlet if the inlet
+        // velocity is chosen too small.
         FluidState fluidState;
         updateFluidStateForBC_(fluidState, pressure_);
         const Scalar density = FluidSystem::density(fluidState, 0);
+        values[Indices::pressureIdx] = pressure_ - density*this->gravity()[1]*(this->fvGridGeometry().bBoxMax()[1] - globalPos[1]);
 
-        PrimaryVariables values(0.0);
-        values[Indices::pressureIdx] = pressure_ + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+
+        // gravity has negative sign
         values[Indices::conti0EqIdx + 1] = moleFraction_;
         values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1])
                                                         * (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
@@ -273,9 +286,9 @@ public:
     /*!
      * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
-    Scalar permeability(const SubControlVolumeFace& scvf) const
+    Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
     {
-        return couplingManager().couplingData().darcyPermeability(scvf);
+        return couplingManager().couplingData().darcyPermeability(element, scvf);
     }
 
     /*!
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh b/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
index c20e9e42963d7f778a13a38b8414606bf6e9d259..fb4f258c0b28b1dd7b6de8bc54df920cfb698f73 100644
--- a/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
@@ -44,52 +44,61 @@ class DarcySubProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
+// Create new type tags
+namespace TTag {
+struct DarcyTypeTag { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
+} // end namespace TTag
 
 // Set the problem property
-SET_TYPE_PROP(DarcyTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
 
 // The fluid system
-SET_PROP(DarcyTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyTypeTag>
 {
-    using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+    using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
     using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
 };
 
 // Use moles
-SET_BOOL_PROP(DarcyTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyTypeTag> { static constexpr bool value = true; };
 
 // Do not replace one equation with a total mass balance
-SET_INT_PROP(DarcyTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTypeTag> { static constexpr int value = 3; };
 
 //! Use a model with constant tortuosity for the effective diffusivity
 SET_TYPE_PROP(DarcyTypeTag, EffectiveDiffusivityModel,
-              DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+              DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>);
 
 // Set the grid type
-SET_TYPE_PROP(DarcyTypeTag, Grid, Dune::YaspGrid<2>);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyTypeTag> { using type = Dune::YaspGrid<2>; };
 
 // Set the spatial paramaters type
-SET_TYPE_PROP(DarcyTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTypeTag> { using type = 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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+    using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
 
     // copy some indices for convenience
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
     enum {
         // grid and world dimension
         dim = GridView::dimension,
@@ -105,7 +114,7 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
 
-    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
     using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
 
 public:
@@ -114,7 +123,6 @@ public:
     : 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);
@@ -224,8 +232,6 @@ public:
                 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]);
@@ -246,21 +252,11 @@ public:
             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].
      *
@@ -311,7 +307,7 @@ public:
         NumEqVector values(0.0);
 
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
-            values = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+            values = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
 
         return values;
     }
@@ -332,10 +328,10 @@ public:
      * \param scv The subcontrolvolume
      */
     template<class ElementVolumeVariables>
-    NumEqVector source(const Element &element,
+    NumEqVector source(const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
-                       const SubControlVolume &scv) const
+                       const SubControlVolume& scv) const
     { return NumEqVector(0.0); }
 
     // \}
@@ -348,11 +344,13 @@ public:
      * For this method, the \a priVars parameter stores primary
      * variables.
      */
-    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
     {
+        static const Scalar stokesPressure = getParamFromGroup<Scalar>("Stokes", "Problem.Pressure");
+
         PrimaryVariables values(0.0);
+        values[Indices::pressureIdx] = stokesPressure;
         values[transportCompIdx] = moleFraction_;
-        values[pressureIdx] = pressure_;
         return values;
     }
 
@@ -384,7 +382,6 @@ private:
 
     Scalar eps_;
     Scalar moleFraction_;
-    Scalar pressure_;
 
     Scalar initialWaterContent_ = 0.0;
     Scalar lastWaterMass_ = 0.0;
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
index cdb7abdd800c9994ade2dab51adda2076ed7a51d..8528aeba78ac6a54a16aec617d21ba99395b431d 100644
--- 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
@@ -34,22 +34,18 @@
 #include <dumux/common/properties.hh>
 #include <dumux/common/parameters.hh>
 #include <dumux/common/dumuxmessage.hh>
-#include <dumux/common/geometry/diameter.hh>
+#include <dumux/common/partial.hh>
 #include <dumux/linear/seqsolverbackend.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
-#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/method.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>
 
@@ -59,15 +55,17 @@
 namespace Dumux {
 namespace Properties {
 
-SET_PROP(StokesTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTypeTag)>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
-SET_PROP(DarcyTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TTAG(StokesTypeTag), TTAG(StokesTypeTag), TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesTypeTag, Properties::TTag::StokesTypeTag, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
@@ -89,16 +87,16 @@ int main(int argc, char** argv) try
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = TTAG(StokesTypeTag);
-    using DarcyTypeTag = TTAG(DarcyTypeTag);
+    using StokesTypeTag = Properties::TTag::StokesTypeTag;
+    using DarcyTypeTag = Properties::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)>;
+    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
     DarcyGridManager darcyGridManager;
     darcyGridManager.init("Darcy"); // pass parameter group
 
-    using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+    using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
     StokesGridManager stokesGridManager;
     stokesGridManager.init("Stokes"); // pass parameter group
 
@@ -107,10 +105,10 @@ int main(int argc, char** argv) try
     const auto& stokesGridView = stokesGridManager.grid().leafGridView();
 
     // create the finite volume grid geometry
-    using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+    using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
     auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
     stokesFvGridGeometry->update();
-    using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>;
     auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
     darcyFvGridGeometry->update();
 
@@ -126,27 +124,22 @@ int main(int argc, char** argv) try
     constexpr auto darcyIdx = CouplingManager::darcyIdx;
 
     // the problem (initial and boundary conditions)
-    using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+    using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
     auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
-    using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
     auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
 
     // initialize the fluidsystem (tabulation)
-    GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+    GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
 
     // get some time loop parameters
-    using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+    using Scalar = GetPropType<StokesTypeTag, Properties::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);
+    auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
     stokesProblem->setTimeLoop(timeLoop);
     darcyProblem->setTimeLoop(timeLoop);
@@ -157,32 +150,22 @@ int main(int argc, char** argv) try
     sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
     sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
 
-    const auto& cellCenterSol = sol[stokesCellCenterIdx];
-    const auto& faceSol = sol[stokesFaceIdx];
+    auto stokesSol = partial(sol, stokesCellCenterIdx, 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);
+    using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
     auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
-    stokesGridVariables->init(stokesSol, solStokesOld);
-    using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+    stokesGridVariables->init(stokesSol);
+    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
     auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
-    darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+    darcyGridVariables->init(sol[darcyIdx]);
 
     // intialize the vtk output module
 #if EXNUMBER >= 1
@@ -194,12 +177,14 @@ int main(int argc, char** argv) try
     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);
+    StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+    GetPropType<StokesTypeTag, Properties::VtkOutputFields>::initOutputModule(stokesVtkWriter);
+    stokesVtkWriter.write(0);
 
-    VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
-    GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
     darcyVtkWriter.write(0.0);
 
     // intialize the subproblems
@@ -221,14 +206,7 @@ int main(int argc, char** argv) try
     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
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
index 8d08fe836a2682e36007bbafa51c6438a5d15c6e..1b0a62ad6015111620d85d017969309b93322fb5 100644
--- 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
@@ -14,10 +14,11 @@ Cells = 16 16
 
 [Stokes.Problem]
 Name = stokes
-EnableGravity = false
+EnableGravity = true
+EnableInertiaTerms = true
 MoleFraction = 0.0 # -
 Pressure = 1e5 # Pa
-Velocity = 1e-3 # m/s
+Velocity = 1 # m/s
 
 [Darcy.Problem]
 Name = darcy
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
index 3d59b25a071ed0e71a0c3baac902b4847a2ce840..ae63df0d88f092b00544e4d3a6e38d0c976129a3 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
@@ -40,32 +40,41 @@ class StokesSubProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(StokesTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+// Create new type tags
+namespace TTag {
+struct StokesTypeTag { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+} // end namespace TTag
 
 // Set the grid type
-SET_TYPE_PROP(StokesTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesTypeTag> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
 // The fluid system
-SET_PROP(StokesTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesTypeTag>
 {
-    using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+    using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
     using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
 };
 
 // Do not replace one equation with a total mass balance
-SET_INT_PROP(StokesTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::StokesTypeTag> { static constexpr int value = 3; };
 
 // Use formulation based on mass fractions
-SET_BOOL_PROP(StokesTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = 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);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesTypeTag> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
 }
 
 /*!
@@ -79,31 +88,31 @@ 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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
 
-    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVGridGeometry = GetPropType<TypeTag, Properties::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 ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+    using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
+    using FluidState = GetPropType<TypeTag, Properties::FluidState>;
 
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
-    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
-    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
     using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
 
-    static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+    static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
 
-    static constexpr auto dim = GET_PROP_TYPE(TypeTag, ModelTraits)::dim();
+    static constexpr auto dim = GetPropType<TypeTag, Properties::ModelTraits>::dim();
     static constexpr auto transportCompIdx = 1;
 
 public:
@@ -120,10 +129,6 @@ public:
      */
     // \{
 
-
-    bool shouldWriteRestartFile() const
-    { return false; }
-
    /*!
      * \brief Return the temperature within the domain in [K].
      */
@@ -221,9 +226,9 @@ public:
 
         if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
         {
-            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
 
-            const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+            const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
             values[Indices::conti0EqIdx] = massFlux[0];
             values[Indices::conti0EqIdx + 1] = massFlux[1];
         }
@@ -249,26 +254,34 @@ public:
      */
     // \{
 
-   /*!
-     * \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;
-    }
+    /*!
+      * \brief Evaluate the initial value for a control volume.
+      *
+      * \param globalPos The global position
+      */
+     PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+     {
+         PrimaryVariables values(0.0);
+
+         // This is only an approximation of the actual hydorostatic pressure gradient.
+         // Air is compressible (the density depends on pressure), thus the actual
+         // vertical pressure profile is non-linear.
+         // This discrepancy can lead to spurious flows at the outlet if the inlet
+         // velocity is chosen too small.
+         FluidState fluidState;
+         updateFluidStateForBC_(fluidState, pressure_);
+         const Scalar density = FluidSystem::density(fluidState, 0);
+         values[Indices::pressureIdx] = pressure_ - density*this->gravity()[1]*(this->fvGridGeometry().bBoxMax()[1] - globalPos[1]);
+
+
+         // gravity has negative sign
+         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; }
@@ -276,9 +289,9 @@ public:
     /*!
      * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
-    Scalar permeability(const SubControlVolumeFace& scvf) const
+    Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
     {
-        return couplingManager().couplingData().darcyPermeability(scvf);
+        return couplingManager().couplingData().darcyPermeability(element, scvf);
     }
 
     /*!
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
index b2d951bc2addbf1250f87ee8eba511827ed27a26..cb79f6a3ebebaaa125c5829c880bfc1da3b6d704 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
@@ -49,19 +49,24 @@ class DarcySubProblem;
 
 namespace Properties
 {
+// Create new type tags
+namespace TTag {
 #if EXNUMBER >= 1
-NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, TwoPNC));
+struct DarcyTypeTag { using InheritsFrom = std::tuple<TwoPNC, CCTpfaModel>; };
 #else
-NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
+struct DarcyTypeTag { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
 #endif
+} // end namespace TTag
 
 // Set the problem property
-SET_TYPE_PROP(DarcyTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
 
 // The fluid system
-SET_PROP(DarcyTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyTypeTag>
 {
-    using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+    using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
 #if EXNUMBER == 0
     using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
 #else
@@ -70,29 +75,35 @@ SET_PROP(DarcyTypeTag, FluidSystem)
 };
 
 // Use moles
-SET_BOOL_PROP(DarcyTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyTypeTag> { static constexpr bool value = true; };
 
 // Do not replace one equation with a total mass balance
-SET_INT_PROP(DarcyTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTypeTag> { static constexpr int value = 3; };
 
 //! Use a model with constant tortuosity for the effective diffusivity
 SET_TYPE_PROP(DarcyTypeTag, EffectiveDiffusivityModel,
-              DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+              DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>);
 
 // Set the grid type
-SET_TYPE_PROP(DarcyTypeTag, Grid, Dune::YaspGrid<2>);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyTypeTag> { using type = 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)
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyTypeTag>
 { static constexpr auto value = TwoPFormulation::p1s0; };
 #endif
 
 // Set the spatial paramaters type
 #if EXNUMBER >= 1
-SET_TYPE_PROP(DarcyTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTypeTag> { using type = TwoPSpatialParams<TypeTag>; };
 #else
-SET_TYPE_PROP(DarcyTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTypeTag> { using type = OnePSpatialParams<TypeTag>; };
 #endif
 }
 
@@ -100,19 +111,19 @@ 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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+    using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
 
     // copy some indices for convenience
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
     enum {
         // grid and world dimension
         dim = GridView::dimension,
@@ -135,7 +146,7 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
 
-    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
     using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
 
 public:
@@ -148,7 +159,6 @@ public:
 #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);
@@ -275,10 +285,8 @@ public:
                 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);
+                NumEqVector flux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
 #else
                 NumEqVector flux(0.0); // add "massCouplingCondition" from the couplingManager here
 #endif
@@ -301,21 +309,6 @@ public:
             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].
      *
@@ -366,7 +359,7 @@ public:
         NumEqVector values(0.0);
 
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
-            values = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+            values = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
 
         return values;
     }
@@ -405,6 +398,8 @@ public:
      */
     PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
     {
+        static const Scalar stokesPressure = getParamFromGroup<Scalar>("Stokes", "Problem.Pressure");
+
         PrimaryVariables values(0.0);
 #if EXNUMBER >= 3
         values.setState(3/*bothPhases*/);
@@ -415,7 +410,7 @@ public:
 #else
         values[transportCompIdx] = moleFraction_;
 #endif
-        values[pressureIdx] = pressure_;
+        values[pressureIdx] = stokesPressure;
         return values;
     }
 
@@ -451,7 +446,6 @@ private:
 #else
     Scalar moleFraction_;
 #endif
-    Scalar pressure_;
 
     Scalar initialWaterContent_ = 0.0;
     Scalar lastWaterMass_ = 0.0;