diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
index f3d56c7e2e9e292510d60db415d9233447f88259..d37eabf37f241a0d6c493c51bd17348df8ac623c 100644
--- a/exercises/exercise-coupling-ff-pm/README.md
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -323,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.
@@ -332,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));
 }
@@ -371,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/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/models/ex_models_coupling_ff-pm.input b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
index 1b0a62ad6015111620d85d017969309b93322fb5..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
@@ -39,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/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;