diff --git a/exercises/exercise-coupling-ff-pm/1pspatialparams.hh b/exercises/exercise-coupling-ff-pm/1pspatialparams.hh
index 64887e33307e8a45908dabc131e919a99f95dfa1..e95c24917d05382c77b53c769357a61ee41a9994 100644
--- a/exercises/exercise-coupling-ff-pm/1pspatialparams.hh
+++ b/exercises/exercise-coupling-ff-pm/1pspatialparams.hh
@@ -38,14 +38,14 @@ namespace Dumux
  */
 template<class TypeTag>
 class OnePSpatialParams
-: public FVSpatialParamsOneP<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
-                             typename GET_PROP_TYPE(TypeTag, Scalar),
+: public FVSpatialParamsOneP<GetPropType<TypeTag, Properties::FVGridGeometry>,
+                             GetPropType<TypeTag, Properties::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 FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
     using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<TypeTag>>;
 
     using Element = typename GridView::template Codim<0>::Entity;
diff --git a/exercises/exercise-coupling-ff-pm/2pspatialparams.hh b/exercises/exercise-coupling-ff-pm/2pspatialparams.hh
index 1830f3cbf2381280ca1213b8a07a17890846203c..948b7a7b7ce4e00bbba83f88b399b8889eee7177 100644
--- a/exercises/exercise-coupling-ff-pm/2pspatialparams.hh
+++ b/exercises/exercise-coupling-ff-pm/2pspatialparams.hh
@@ -40,15 +40,15 @@ namespace Dumux
  */
 template<class TypeTag>
 class TwoPSpatialParams
-: public FVSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
-                         typename GET_PROP_TYPE(TypeTag, Scalar),
+: public FVSpatialParams<GetPropType<TypeTag, Properties::FVGridGeometry>,
+                         GetPropType<TypeTag, Properties::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 Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
     using Element = typename GridView::template Codim<0>::Entity;
-    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using ParentType = FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<TypeTag>>;
diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
index 58b75459bbae1fe406dca15062d96d27b8adff8f..d37eabf37f241a0d6c493c51bd17348df8ac623c 100644
--- a/exercises/exercise-coupling-ff-pm/README.md
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -83,6 +83,7 @@ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
 {
     values.setCouplingNeumann(Indices::conti0EqIdx);
     values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+    values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
 }
 ```
 
@@ -123,7 +124,11 @@ $`\frac{\partial v_x}{\partial y} = \frac{\alpha}{\sqrt K} (v_x - q_{pm})\quad`$
 
 with  $`\quad q_{pm}=0`$.
 
-To include this, just set a Beavers-Joseph-Saffman (BJS) boundary condition for the respective momentum balance:
+To include this, just replace the no-slip condition at the interface
+``` cpp
+values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
+```
+with a Beavers-Joseph-Saffman (BJS) boundary condition for the respective momentum balance:
 ``` cpp
 values.setBJS(Indices::momentumXBalanceIdx);
 ```
@@ -223,7 +228,8 @@ liquid saturation as primary variables (p_g-S_l -> `p1s0`) or vice versa.
 * Set the property
 
 ```
-SET_PROP(DarcyTypeTag, Formulation)
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyTypeTag>
 { static constexpr auto value = TwoPFormulation::p1s0; };
 ```
   in the Properties section in the problem file.
@@ -263,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.
@@ -289,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.
 
@@ -321,7 +323,7 @@ For using the compositional zero equation turbulence model, the following header
 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)
+* Change the entry in the `ZeroEqTypeTag` 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.
@@ -330,9 +332,9 @@ Here, the turbulent free flow is wall bounded which means that the main reason f
 of turbulent flow is the presence of walls.
 The porous medium at the bottom and also the channel wall are such walls.
 Inside the problem you have to tell the turbulence model, where these walls are located
-by providing an `isOnWall()` function:
+by providing an `isOnWallAtPos()` function:
 ```c++
-bool isOnWall(const GlobalPosition& globalPos) const
+bool isOnWallAtPos(const GlobalPosition& globalPos) const
 {
     return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
 }
@@ -349,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`).
@@ -369,7 +371,7 @@ Instead of computing the whole cross-section of a channel, you can use symmetric
 values.setAllSymmetry();
 ```
 
-In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `isOnWall(globalPos)`
+In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `isOnWallAtPos(globalPos)`
 and `initialAtPos(globalPos)` method.
 
 __Task C__:
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
index 85226b998887b5ea030ed889ae598573ce250da3..62673d5a81469edb1b67376c506b54a2dd840507 100644
--- 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
@@ -31,10 +31,11 @@
 #include <dumux/common/properties.hh>
 #include <dumux/common/parameters.hh>
 #include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/partial.hh>
 #include <dumux/linear/seqsolverbackend.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/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>
@@ -51,15 +52,17 @@
 namespace Dumux {
 namespace Properties {
 
-SET_PROP(StokesOnePTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesOnePTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyOnePTypeTag)>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOnePTypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
-SET_PROP(DarcyOnePTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOnePTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TTAG(StokesOnePTypeTag), TTAG(StokesOnePTypeTag), TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOnePTypeTag, Properties::TTag::StokesOnePTypeTag, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
@@ -81,8 +84,8 @@ int main(int argc, char** argv) try
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = TTAG(StokesOnePTypeTag);
-    using DarcyTypeTag = TTAG(DarcyOnePTypeTag);
+    using StokesTypeTag = Properties::TTag::StokesOnePTypeTag;
+    using DarcyTypeTag = Properties::TTag::DarcyOnePTypeTag;
 
 
 
@@ -90,11 +93,11 @@ int main(int argc, char** argv) try
 
     // 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)>;
+    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
 
@@ -149,10 +152,10 @@ int main(int argc, char** argv) try
 
 
     // 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();
 
@@ -168,9 +171,9 @@ 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);
 
     // the solution vector
@@ -179,24 +182,18 @@ 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);
-    sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
-    sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+    darcyProblem->applyInitialSolution(sol[darcyIdx]);
 
     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);
-    using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
     auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
     darcyGridVariables->init(sol[darcyIdx]);
 
@@ -204,16 +201,18 @@ int main(int argc, char** argv) try
     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);
 
     //****** 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);
+    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);
 
     // the assembler for a stationary problem
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
index e3813950725ba4837e96b79d76dc9ccc74dc0d4f..4117b7f17c501817b738dfdb6836cc2c3223e4a8 100644
--- 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
@@ -28,6 +28,7 @@ Grading1 = 1
 [Stokes.Problem]
 Name = stokes
 PressureDifference = 1e-9
+EnableInertiaTerms = false
 
 [Darcy.Problem]
 Name = darcy
diff --git a/exercises/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh b/exercises/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
index 48d29d667983db35a315d3417554f0c1880784db..c5649a77e27284fe946d8821ca4c526e4f4bf17c 100644
--- a/exercises/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
@@ -42,20 +42,25 @@ class StokesSubProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(StokesOnePTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
+// Create new type tags
+namespace TTag {
+struct StokesOnePTypeTag { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+} // end namespace TTag
 
 // the fluid system
-SET_PROP(StokesOnePTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesOnePTypeTag>
 {
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
 };
 
 // Set the grid type
-SET_PROP(StokesOnePTypeTag, Grid)
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesOnePTypeTag>
 {
     static constexpr auto dim = 2;
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
 
 //****** comment out for the last exercise *****//
@@ -67,13 +72,15 @@ SET_PROP(StokesOnePTypeTag, Grid)
 };
 
 // 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);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesOnePTypeTag> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesOnePTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOnePTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOnePTypeTag> { static constexpr bool value = true; };
 }
 
 /*!
@@ -84,25 +91,25 @@ 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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
-    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    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 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 PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
-    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
 
 public:
     StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
@@ -170,6 +177,7 @@ public:
         {
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+            values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
         }
 
         return values;
@@ -208,8 +216,8 @@ public:
 
         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);
+            values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
         }
 
         return values;
@@ -246,9 +254,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/interface/ex_interface_pmproblem.hh b/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
index 8b9f74954fb20bb8163788262cbae5f820f85d82..e447100ca1c7e3d2852acf622f1b3570ab0772a7 100644
--- a/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
@@ -46,23 +46,29 @@ class DarcySubProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(DarcyOnePTypeTag, INHERITS_FROM(CCTpfaModel, OneP));
+// Create new type tags
+namespace TTag {
+struct DarcyOnePTypeTag { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+} // end namespace TTag
 
 // Set the problem property
-SET_TYPE_PROP(DarcyOnePTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOnePTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
 
 // the fluid system
-SET_PROP(DarcyOnePTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOnePTypeTag>
 {
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
 };
 
 // Set the grid type
-SET_PROP(DarcyOnePTypeTag, Grid)
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOnePTypeTag>
 {
     static constexpr auto dim = 2;
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
 
 //****** comment out for the last exercise *****//
@@ -73,7 +79,8 @@ SET_PROP(DarcyOnePTypeTag, Grid)
     // using type = Dune::SubGrid<dim, HostGrid>;
 };
 
-SET_TYPE_PROP(DarcyOnePTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOnePTypeTag> { using type = OnePSpatialParams<TypeTag>; };
 }
 
 /*!
@@ -83,23 +90,23 @@ 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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+    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>;
 
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
-    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
 
 public:
     DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
@@ -188,7 +195,7 @@ public:
 
         // ... except at the coupling interface
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
-            values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+            values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
 
         return values;
     }
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/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
index 233cdd3fe0dff13f2afd3052f56dd7bce30ded28..523a33706876711eb7bce77fd58758ad54b8c239 100644
--- 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
@@ -34,18 +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>
-#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+#include <dumux/multidomain/newtonsolver.hh>
 
 #include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
 
@@ -55,15 +55,17 @@
 namespace Dumux {
 namespace Properties {
 
-SET_PROP(ZeroEqTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::ZeroEqTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTwoPTwoCTypeTag)>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTwoPTwoCTypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
-SET_PROP(DarcyTwoPTwoCTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyTwoPTwoCTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TTAG(ZeroEqTypeTag), TTAG(ZeroEqTypeTag), TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::ZeroEqTypeTag, Properties::TTag::ZeroEqTypeTag, 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(ZeroEqTypeTag);
-    using DarcyTypeTag = TTAG(DarcyTwoPTwoCTypeTag);
+    using StokesTypeTag = Properties::TTag::ZeroEqTypeTag;
+    using DarcyTypeTag = Properties::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)>;
+    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<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
 
     // set timeloop for the subproblems, needed for boundary value variations
@@ -155,18 +152,10 @@ int main(int argc, char** argv) try
     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];
+    auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
 
+    stokesProblem->applyInitialSolution(stokesSol);
     darcyProblem->applyInitialSolution(sol[darcyIdx]);
-    auto solDarcyOld = sol[darcyIdx];
 
     auto solOld = sol;
 
@@ -176,23 +165,25 @@ int main(int argc, char** argv) try
     // TODO: update dynamic wall properties
 
     // 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, GET_PROP_VALUE(StokesTypeTag, PhaseIdx)> 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);
 
     // the assembler with time loop for instationary problem
@@ -211,11 +202,8 @@ 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
-    using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
-
     // the non-linear solver
-    using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+    using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
     NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
 
     // time loop
@@ -230,10 +218,6 @@ int main(int argc, char** argv) try
         // 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
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
index 514cf156d55211777c2ec0980f37d1a95124867e..00e02465ce811bb20125e851dbe57964dc023e52 100644
--- 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
@@ -27,6 +27,7 @@ RefVelocity = 3.5 # [m/s]
 RefPressure = 1e5 # [Pa]
 refMoleFrac = 0 # [-]
 RefTemperature = 298.15 # [K]
+EnableInertiaTerms = true
 
 [Darcy.Problem]
 Name = darcy
@@ -60,3 +61,8 @@ MaxRelativeShift = 1e-5
 [Assembly]
 NumericDifferenceMethod = 0
 NumericDifference.BaseEpsilon = 1e-8
+
+[Component]
+SolidDensity = 2700
+SolidThermalConductivity = 2.8
+SolidHeatCapacity = 790
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
index 6c89102d08587d591a7d4f1125feede34dd0b6da..88d9cdeea6cd9551067187e2a7ab6a7929fba8a5 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
@@ -39,32 +39,41 @@ class FreeFlowSubProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
+// Create new type tags
+namespace TTag {
+struct ZeroEqTypeTag { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
+} // end namespace TTag
 
 // Set the grid type
-SET_TYPE_PROP(ZeroEqTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::ZeroEqTypeTag> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
 // The fluid system
-SET_PROP(ZeroEqTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::ZeroEqTypeTag>
 {
-  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
   static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
   using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
 };
 
-SET_INT_PROP(ZeroEqTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::ZeroEqTypeTag> { static constexpr int value = 3; };
 
 // Use formulation based on mass fractions
-SET_BOOL_PROP(ZeroEqTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::ZeroEqTypeTag> { static constexpr bool value = 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);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::ZeroEqTypeTag> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::ZeroEqTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::ZeroEqTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::ZeroEqTypeTag> { static constexpr bool value = true; };
 }
 
 /*!
@@ -75,31 +84,31 @@ 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 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>>;
 
     using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
 
-    static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+    static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
 
 public:
     FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
@@ -158,7 +167,7 @@ public:
             values.setDirichlet(Indices::velocityXIdx);
             values.setDirichlet(Indices::velocityYIdx);
             values.setDirichlet(Indices::conti0EqIdx + 1);
-            values.setDirichlet(Indices::energyBalanceIdx);
+            values.setDirichlet(Indices::energyEqIdx);
         }
 
         if (onLowerBoundary_(globalPos))
@@ -167,7 +176,7 @@ public:
             values.setDirichlet(Indices::velocityYIdx);
             values.setNeumann(Indices::conti0EqIdx);
             values.setNeumann(Indices::conti0EqIdx + 1);
-            values.setNeumann(Indices::energyBalanceIdx);
+            values.setNeumann(Indices::energyEqIdx);
         }
 
         if (onUpperBoundary_(globalPos))
@@ -176,21 +185,21 @@ public:
             values.setDirichlet(Indices::velocityYIdx);
             values.setNeumann(Indices::conti0EqIdx);
             values.setNeumann(Indices::conti0EqIdx + 1);
-            values.setNeumann(Indices::energyBalanceIdx);
+            values.setNeumann(Indices::energyEqIdx);
         }
 
         if (onRightBoundary_(globalPos))
         {
             values.setDirichlet(Indices::pressureIdx);
             values.setOutflow(Indices::conti0EqIdx + 1);
-            values.setOutflow(Indices::energyBalanceIdx);
+            values.setOutflow(Indices::energyEqIdx);
         }
 
         if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
         {
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::conti0EqIdx + 1);
-            values.setCouplingNeumann(Indices::energyBalanceIdx);
+            values.setCouplingNeumann(Indices::energyEqIdx);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
             values.setBJS(Indices::momentumXBalanceIdx);
         }
@@ -229,12 +238,12 @@ public:
         PrimaryVariables values(0.0);
         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, diffCoeffAvgType_);
+            const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, 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_);
+            values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
         }
         return values;
     }
@@ -305,7 +314,7 @@ 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().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center());
     }
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
index 0a9cd5f54e1f81694257ddf4d8f14347343a5fcc..de4bb4e71ccca3ce12c36ba5cbb455600c15aa13 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
@@ -42,27 +42,37 @@ class DarcySubProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(DarcyTwoPTwoCTypeTag, INHERITS_FROM(CCTpfaModel, TwoPTwoCNI));
+// Create new type tags
+namespace TTag {
+struct DarcyTwoPTwoCTypeTag { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
+} // end namespace TTag
 
 // Set the problem property
-SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
 
 // the fluid system
-SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
 
 //! Set the default formulation to pw-Sn: This can be over written in the problem.
-SET_PROP(DarcyTwoPTwoCTypeTag, Formulation)
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyTwoPTwoCTypeTag>
 { static constexpr auto value = TwoPFormulation::p1s0; };
 
 // The gas component balance (air) is replaced by the total mass balance
-SET_INT_PROP(DarcyTwoPTwoCTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { static constexpr int value = 3; };
 
 // Set the grid type
-SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
-SET_BOOL_PROP(DarcyTwoPTwoCTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { static constexpr bool value = true; };
 
-SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = TwoPSpatialParams<TypeTag>; };
 }
 
 /*!
@@ -72,22 +82,22 @@ 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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+    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 ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
 
-    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
     // copy some indices for convenience
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
     enum {
         // primary variable indices
         conti0EqIdx = Indices::conti0EqIdx,
@@ -100,7 +110,7 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
-    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
     using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
 
     using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
@@ -220,12 +230,12 @@ public:
 
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
         {
-            const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+            const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, 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_);
+            values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
         }
 
         return values;
diff --git a/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh b/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh
index 64887e33307e8a45908dabc131e919a99f95dfa1..e95c24917d05382c77b53c769357a61ee41a9994 100644
--- a/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh
@@ -38,14 +38,14 @@ namespace Dumux
  */
 template<class TypeTag>
 class OnePSpatialParams
-: public FVSpatialParamsOneP<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
-                             typename GET_PROP_TYPE(TypeTag, Scalar),
+: public FVSpatialParamsOneP<GetPropType<TypeTag, Properties::FVGridGeometry>,
+                             GetPropType<TypeTag, Properties::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 FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
     using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<TypeTag>>;
 
     using Element = typename GridView::template Codim<0>::Entity;
diff --git a/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh b/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh
index 1830f3cbf2381280ca1213b8a07a17890846203c..948b7a7b7ce4e00bbba83f88b399b8889eee7177 100644
--- a/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh
@@ -40,15 +40,15 @@ namespace Dumux
  */
 template<class TypeTag>
 class TwoPSpatialParams
-: public FVSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
-                         typename GET_PROP_TYPE(TypeTag, Scalar),
+: public FVSpatialParams<GetPropType<TypeTag, Properties::FVGridGeometry>,
+                         GetPropType<TypeTag, Properties::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 Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
     using Element = typename GridView::template Codim<0>::Entity;
-    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using ParentType = FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<TypeTag>>;
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
index 116ed7e2d8e4469f44bde197b8b71974422264b1..84434750b1961a57b3d6a06be27df697eb06e601 100644
--- 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
@@ -31,10 +31,11 @@
 #include <dumux/common/properties.hh>
 #include <dumux/common/parameters.hh>
 #include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/partial.hh>
 #include <dumux/linear/seqsolverbackend.hh>
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/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>
@@ -171,16 +172,10 @@ 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);
-    sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
-    sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+    darcyProblem->applyInitialSolution(sol[darcyIdx]);
 
     couplingManager->init(stokesProblem, darcyProblem, sol);
 
@@ -196,8 +191,8 @@ int main(int argc, char** argv) try
     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);
 
 #if EXNUMBER >= 2
     stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
@@ -205,8 +200,10 @@ int main(int argc, char** argv) try
 
     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);
 
     // the assembler for a stationary problem
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
index e3813950725ba4837e96b79d76dc9ccc74dc0d4f..4117b7f17c501817b738dfdb6836cc2c3223e4a8 100644
--- 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
@@ -28,6 +28,7 @@ Grading1 = 1
 [Stokes.Problem]
 Name = stokes
 PressureDifference = 1e-9
+EnableInertiaTerms = false
 
 [Darcy.Problem]
 Name = darcy
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
index 369a7bd2b9238272a36759a1856d4f7369f55a37..c0dc45575729802b20b101c00970934afb41e5c4 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_ffproblem.hh
@@ -74,7 +74,6 @@ SET_BOOL_PROP(StokesOnePTypeTag, EnableFVGridGeometryCache, true);
 SET_BOOL_PROP(StokesOnePTypeTag, EnableGridFluxVariablesCache, true);
 SET_BOOL_PROP(StokesOnePTypeTag, EnableGridVolumeVariablesCache, true);
 
-SET_BOOL_PROP(StokesOnePTypeTag, EnableInertiaTerms, false);
 }
 
 /*!
@@ -185,12 +184,12 @@ public:
             values.setCouplingNeumann(scvf.directionIndex());
 #endif
 
-#if EXNUMBER == 2
+#if EXNUMBER < 2
+            values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
+#elif EXNUMBER == 2
             // set the Beaver-Joseph-Saffman slip condition for the tangential momentum balance equation
             values.setBJS(Indices::momentumXBalanceIdx);
-#endif
-
-#if EXNUMBER == 3
+#else
             // set the Beaver-Joseph-Saffman slip condition for the tangential momentum balance equation,
             // consider orientation of face automatically
             values.setBJS(1 - scvf.directionIndex());
@@ -233,11 +232,11 @@ public:
 
         if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
         {
-            values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+            values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
 #if EXNUMBER < 3
-            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+            values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
 #else
-            values[scvf.directionIndex()] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+            values[scvf.directionIndex()] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
 #endif
 
         }
@@ -283,9 +282,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/interface/ex_interface_pmproblem.hh b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
index d78c76f0ed0958c319724812d1008337de2f00ea..3549ff03521c727f481ccf32af16817f2fd01e4b 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
@@ -184,7 +184,7 @@ public:
         NumEqVector values(0.0);
 
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
-            values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+            values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
 
         return values;
     }
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..7dda6b1b36292024dc4d5f410eb9c08f9b8b0675 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
@@ -38,8 +39,8 @@ VgN = 8.0
 
 [Problem]
 Name = ex_models_coupling
-PlotFluxes = false
-PlotStorage = false
+PlotFluxes = true
+PlotStorage = true
 
 [Newton]
 MaxSteps = 12
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;
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
index e1e4739c3fabc8afdc943d6318472b984b23c6c4..ac7f7584f622064966d4eefc47487889de0180d9 100644
--- 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
@@ -34,18 +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>
-#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+#include <dumux/multidomain/newtonsolver.hh>
 
 #include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
 
@@ -55,15 +55,17 @@
 namespace Dumux {
 namespace Properties {
 
-SET_PROP(ZeroEqTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::ZeroEqTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTwoPTwoCTypeTag)>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTwoPTwoCTypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
-SET_PROP(DarcyTwoPTwoCTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyTwoPTwoCTypeTag>
 {
-    using Traits = StaggeredMultiDomainTraits<TTAG(ZeroEqTypeTag), TTAG(ZeroEqTypeTag), TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::ZeroEqTypeTag, Properties::TTag::ZeroEqTypeTag, 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(ZeroEqTypeTag);
-    using DarcyTypeTag = TTAG(DarcyTwoPTwoCTypeTag);
+    using StokesTypeTag = Properties::TTag::ZeroEqTypeTag;
+    using DarcyTypeTag = Properties::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)>;
+    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<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
 
     // set timeloop for the subproblems, needed for boundary value variations
@@ -155,18 +152,10 @@ int main(int argc, char** argv) try
     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];
+    auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
 
+    stokesProblem->applyInitialSolution(stokesSol);
     darcyProblem->applyInitialSolution(sol[darcyIdx]);
-    auto solDarcyOld = sol[darcyIdx];
 
     auto solOld = sol;
 
@@ -179,23 +168,25 @@ int main(int argc, char** argv) try
 #endif
 
     // 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, GET_PROP_VALUE(StokesTypeTag, PhaseIdx)> 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);
 
     // the assembler with time loop for instationary problem
@@ -214,11 +205,8 @@ 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
-    using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
-
     // the non-linear solver
-    using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+    using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
     NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
 
     // time loop
@@ -233,10 +221,6 @@ int main(int argc, char** argv) try
         // 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);
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
index 514cf156d55211777c2ec0980f37d1a95124867e..00e02465ce811bb20125e851dbe57964dc023e52 100644
--- 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
@@ -27,6 +27,7 @@ RefVelocity = 3.5 # [m/s]
 RefPressure = 1e5 # [Pa]
 refMoleFrac = 0 # [-]
 RefTemperature = 298.15 # [K]
+EnableInertiaTerms = true
 
 [Darcy.Problem]
 Name = darcy
@@ -60,3 +61,8 @@ MaxRelativeShift = 1e-5
 [Assembly]
 NumericDifferenceMethod = 0
 NumericDifference.BaseEpsilon = 1e-8
+
+[Component]
+SolidDensity = 2700
+SolidThermalConductivity = 2.8
+SolidHeatCapacity = 790
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
index 194597f76375334fccdd7473af2054e32f481028..a78efe81c38b40db995cca48a2413720d0787e33 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
@@ -44,36 +44,45 @@ class FreeFlowSubProblem;
 
 namespace Properties
 {
+// Create new type tags
+namespace TTag {
 #if EXNUMBER >= 1
-NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNCNI));
+struct ZeroEqTypeTag { using InheritsFrom = std::tuple<ZeroEqNCNI, StaggeredFreeFlowModel>; };
 #else
-NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
+struct ZeroEqTypeTag { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
 #endif
+} // end namespace TTag
 
 // Set the grid type
-SET_TYPE_PROP(ZeroEqTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::ZeroEqTypeTag> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
 // The fluid system
-SET_PROP(ZeroEqTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::ZeroEqTypeTag>
 {
-  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
   static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
   using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
 };
 
-SET_INT_PROP(ZeroEqTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::ZeroEqTypeTag> { static constexpr int value = 3; };
 
 // Use formulation based on mass fractions
-SET_BOOL_PROP(ZeroEqTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::ZeroEqTypeTag> { static constexpr bool value = 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);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::ZeroEqTypeTag> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::ZeroEqTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::ZeroEqTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::ZeroEqTypeTag> { static constexpr bool value = true; };
 }
 
 /*!
@@ -90,31 +99,31 @@ 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 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>>;
 
     using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
 
-    static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+    static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
 
 public:
     FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
@@ -173,7 +182,7 @@ public:
             values.setDirichlet(Indices::velocityXIdx);
             values.setDirichlet(Indices::velocityYIdx);
             values.setDirichlet(Indices::conti0EqIdx + 1);
-            values.setDirichlet(Indices::energyBalanceIdx);
+            values.setDirichlet(Indices::energyEqIdx);
         }
 
         if (onLowerBoundary_(globalPos))
@@ -182,7 +191,7 @@ public:
             values.setDirichlet(Indices::velocityYIdx);
             values.setNeumann(Indices::conti0EqIdx);
             values.setNeumann(Indices::conti0EqIdx + 1);
-            values.setNeumann(Indices::energyBalanceIdx);
+            values.setNeumann(Indices::energyEqIdx);
         }
 
         if (onUpperBoundary_(globalPos))
@@ -194,7 +203,7 @@ public:
             values.setDirichlet(Indices::velocityYIdx);
             values.setNeumann(Indices::conti0EqIdx);
             values.setNeumann(Indices::conti0EqIdx + 1);
-            values.setNeumann(Indices::energyBalanceIdx);
+            values.setNeumann(Indices::energyEqIdx);
 #endif
         }
 
@@ -202,14 +211,14 @@ public:
         {
             values.setDirichlet(Indices::pressureIdx);
             values.setOutflow(Indices::conti0EqIdx + 1);
-            values.setOutflow(Indices::energyBalanceIdx);
+            values.setOutflow(Indices::energyEqIdx);
         }
 
         if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
         {
             values.setCouplingNeumann(Indices::conti0EqIdx);
             values.setCouplingNeumann(Indices::conti0EqIdx + 1);
-            values.setCouplingNeumann(Indices::energyBalanceIdx);
+            values.setCouplingNeumann(Indices::energyEqIdx);
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
             values.setBJS(Indices::momentumXBalanceIdx);
         }
@@ -248,12 +257,12 @@ public:
         PrimaryVariables values(0.0);
         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, diffCoeffAvgType_);
+            const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, 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_);
+            values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
         }
         return values;
     }
@@ -273,12 +282,12 @@ public:
     { return *couplingManager_; }
 
 #if EXNUMBER >= 2
-    bool isOnWall(const GlobalPosition& globalPos) const
+    bool isOnWallAtPos(const GlobalPosition& globalPos) const
     {
         return (onLowerBoundary_(globalPos));
     }
 #elif EXNUMBER >= 1
-    bool isOnWall(const GlobalPosition& globalPos) const
+    bool isOnWallAtPos(const GlobalPosition& globalPos) const
     {
         return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
     }
@@ -341,7 +350,7 @@ 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().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center());
     }
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
index 1b46ce77d818164c849de90f3855db0665c96775..e8d4477f21c0b02c5796b458b11e3d622677fc08 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
@@ -42,27 +42,37 @@ class DarcySubProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(DarcyTwoPTwoCTypeTag, INHERITS_FROM(CCTpfaModel, TwoPTwoCNI));
+// Create new type tags
+namespace TTag {
+struct DarcyTwoPTwoCTypeTag { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
+} // end namespace TTag
 
 // Set the problem property
-SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
 
 // the fluid system
-SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
 
 //! Set the default formulation to pw-Sn: This can be over written in the problem.
-SET_PROP(DarcyTwoPTwoCTypeTag, Formulation)
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyTwoPTwoCTypeTag>
 { static constexpr auto value = TwoPFormulation::p1s0; };
 
 // The gas component balance (air) is replaced by the total mass balance
-SET_INT_PROP(DarcyTwoPTwoCTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { static constexpr int value = 3; };
 
 // Set the grid type
-SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
-SET_BOOL_PROP(DarcyTwoPTwoCTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { static constexpr bool value = true; };
 
-SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = TwoPSpatialParams<TypeTag>; };
 }
 
 /*!
@@ -72,22 +82,22 @@ 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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+    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 ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
 
-    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
     // copy some indices for convenience
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
     enum {
         // primary variable indices
         conti0EqIdx = Indices::conti0EqIdx,
@@ -100,7 +110,7 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
-    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
     using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
 
     using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
@@ -220,12 +230,12 @@ public:
 
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
         {
-            const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+            const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, 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_);
+            values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
         }
 
         return values;