diff --git a/dumux/freeflow/navierstokes/model.hh b/dumux/freeflow/navierstokes/model.hh
index 912a3d242aa9977823c5cbaae954a8df1ef2f611..42e0b2e4bd739eaff5d136affb26cc985884d0f2 100644
--- a/dumux/freeflow/navierstokes/model.hh
+++ b/dumux/freeflow/navierstokes/model.hh
@@ -105,13 +105,13 @@ struct NavierStokesModelTraits
     using Indices = NavierStokesIndices<dim()>;
 
     //! return the names of the primary variables in cells
-    static std::string primaryVariableNameCell(int pvIdx = 0, int state = 0)
+    static std::string primaryVariableNameCell(int pvIdx = 0)
     {
         return "p";
     }
 
     //! return the names of the primary variables on faces
-    static std::string primaryVariableNameFace(int pvIdx = 0, int state = 0)
+    static std::string primaryVariableNameFace(int pvIdx = 0)
     {
         return "v";
     }
diff --git a/dumux/io/loadsolution.hh b/dumux/io/loadsolution.hh
index 4b06ce183502ba9579ba33d1a9af1357a82291e8..26be5e1282dec2ef68783968a2b76f8e88af0e41 100644
--- a/dumux/io/loadsolution.hh
+++ b/dumux/io/loadsolution.hh
@@ -28,6 +28,8 @@
 #include <iostream>
 #include <vector>
 #include <unordered_set>
+#include <unordered_map>
+#include <type_traits>
 
 #include <dune/common/exceptions.hh>
 #include <dune/common/indices.hh>
@@ -38,8 +40,8 @@
 #include <dumux/io/vtk/vtkreader.hh>
 
 namespace Dumux {
-
 namespace Detail {
+
 //! helper struct detecting if a PrimaryVariables object has a state() function
 struct hasState
 {
@@ -51,17 +53,21 @@ struct hasState
 
 } // end namespace Detail
 
-template <class FieldType, class Container, class GridView, int codim>
+/*!
+ * \ingroup InputOutput
+ * \brief a data handle to communicate the solution on ghosts and overlaps
+ *        when reading from vtk file in parallel
+ */
+template <class Container, class EntityMapper, int codim>
 class LoadSolutionDataHandle
-: public Dune::CommDataHandleIF< LoadSolutionDataHandle<FieldType, Container, GridView, codim>,
-                                 FieldType >
+: public Dune::CommDataHandleIF< LoadSolutionDataHandle<Container, EntityMapper, codim>,
+                                 std::decay_t<decltype(std::declval<Container>()[0])> >
 {
-    using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
-
+    using FieldType = std::decay_t<decltype(std::declval<Container>()[0])>;
 public:
-    LoadSolutionDataHandle(Container &container,
-                           const GridView &gridView)
-    : mapper_(gridView, Dune::mcmgLayout(Dune::Codim<codim>()))
+    LoadSolutionDataHandle(Container& container,
+                           const EntityMapper& mapper)
+    : mapper_(mapper)
     , container_(container)
     {}
 
@@ -76,101 +82,80 @@ public:
     { return 1; }
 
     template<class MessageBufferImp, class EntityType>
-    void gather(MessageBufferImp &buff, const EntityType &e) const
+    void gather(MessageBufferImp& buff, const EntityType& e) const
     {
-        int vIdx = mapper_.index(e);
+        const auto vIdx = mapper_.index(e);
         buff.write(container_[vIdx]);
     }
 
     template<class MessageBufferImp, class EntityType>
-    void scatter(MessageBufferImp &buff, const EntityType &e, size_t n)
+    void scatter(MessageBufferImp& buff, const EntityType& e, size_t n)
     {
-        int vIdx = mapper_.index(e);
-
+        const auto vIdx = mapper_.index(e);
         FieldType tmp;
         buff.read(tmp);
         container_[vIdx] = tmp;
     }
 
 private:
-    ElementMapper mapper_;
-    Container &container_;
+    EntityMapper mapper_;
+    Container& container_;
 };
 
 /*!
  * \ingroup InputOutput
- * \brief read from a sequential file into a solution vector without state
+ * \brief read from a vtk file into a solution vector with primary variables without state
  */
-template <class SolutionVector, class PvNamesFunc>
-auto loadSolutionFromSequentialVtkFile(const std::string fileName,
-                                       const VTKReader::DataType& dataType,
-                                       PvNamesFunc&& pvNamesFunc,
-                                       SolutionVector& sol)
+template <class SolutionVector, class PvNameFunc, class FVGridGeometry>
+auto loadSolutionFromVtkFile(SolutionVector& sol,
+                             const std::string fileName,
+                             PvNameFunc&& pvNameFunc,
+                             const FVGridGeometry& fvGridGeometry,
+                             const VTKReader::DataType& dataType)
 -> typename std::enable_if_t<!decltype(isValid(Detail::hasState())(sol[0]))::value, void>
 {
     VTKReader vtu(fileName);
     using PrimaryVariables = typename SolutionVector::block_type;
     using Scalar = typename PrimaryVariables::field_type;
+    constexpr auto dim = FVGridGeometry::GridView::dimension;
 
     for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
     {
-        const auto pvName = pvNamesFunc(pvIdx);
-        const auto vec = vtu.readData<std::vector<Scalar>>(pvName, dataType);
+        const auto pvName = pvNameFunc(pvIdx);
+        auto vec = vtu.readData<std::vector<Scalar>>(pvName, dataType);
         if (vec.size() != sol.size())
             DUNE_THROW(Dune::IOError, "Size mismatch between solution vector and read data (" << sol.size() << " != " << vec.size() << ")");
 
-        for (std::size_t i = 0; i < sol.size(); ++i)
-            sol[i][pvIdx] = vec[i];
-    }
-}
-
-/*!
- * \ingroup InputOutput
- * \brief read from a parallel file into a solution vector without state
- */
-template <class SolutionVector, class PvNamesFunc, class FVGridGeometry>
-auto loadSolutionFromParallelVtkFile(const std::string fileName,
-                                     const VTKReader::DataType& dataType,
-                                     PvNamesFunc&& pvNamesFunc,
-                                     SolutionVector& sol,
-                                     const FVGridGeometry& fvGridGeometry)
--> typename std::enable_if_t<!decltype(isValid(Detail::hasState())(sol[0]))::value, void>
-{
-    VTKReader vtu(fileName);
-    using PrimaryVariables = typename SolutionVector::block_type;
-    using Scalar = typename PrimaryVariables::field_type;
-    using GridView = typename FVGridGeometry::GridView;
-    static constexpr auto dim = GridView::dimension;
-
-    for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
-    {
-        const auto pvName = pvNamesFunc(pvIdx);
-        auto vec = vtu.readData<std::vector<Scalar>>(pvName, dataType);
-
-        std::vector<bool> visited;
-        if (dataType == VTKReader::DataType::pointData)
+        if (dataType == VTKReader::DataType::cellData)
         {
-            visited.resize(fvGridGeometry.gridView().size(dim), false);
+            std::size_t i = 0;
+            for (const auto& element : elements(fvGridGeometry.gridView(), Dune::Partitions::interior))
+            {
+                const auto eIdx = fvGridGeometry.elementMapper().index(element);
+                sol[eIdx][pvIdx] = vec[i++];
+            }
         }
+        // for staggered face data (which is written out as VTK point data) we just read in the vector
+        else if (dataType == VTKReader::DataType::pointData && FVGridGeometry::discMethod == DiscretizationMethod::staggered)
+        {
+            if (sol.size() != vec.size())
+                DUNE_THROW(Dune::InvalidStateException, "Solution size (" << sol.size() << ") does not match input size (" << vec.size() << ")!");
 
-        std::size_t i = 0;
-        for (const auto& element : elements(fvGridGeometry.gridView(), Dune::Partitions::interior))
+            for (std::size_t i = 0; i < sol.size(); ++i)
+                sol[i][pvIdx] = vec[i];
+        }
+        else
         {
-            if (dataType == VTKReader::DataType::cellData)
-            {
-                auto eIdx = fvGridGeometry.elementMapper().index(element);
-                sol[eIdx][pvIdx] = vec[i];
-                ++i;
-            }
-            else
+            std::size_t i = 0;
+            std::vector<bool> visited(fvGridGeometry.gridView().size(dim), false);
+            for (const auto& element : elements(fvGridGeometry.gridView(), Dune::Partitions::interior))
             {
                 for (int vIdxLocal = 0; vIdxLocal < element.subEntities(dim); ++vIdxLocal)
                 {
-                    auto vIdxGlobal = fvGridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
+                    const auto vIdxGlobal = fvGridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
                     if (!visited[vIdxGlobal])
                     {
-                        sol[vIdxGlobal][pvIdx] = vec[i];
-                        ++i;
+                        sol[vIdxGlobal][pvIdx] = vec[i++];
                         visited[vIdxGlobal] = true;
                     }
                 }
@@ -181,108 +166,77 @@ auto loadSolutionFromParallelVtkFile(const std::string fileName,
 
 /*!
  * \ingroup InputOutput
- * \brief read from a sequential file into a solution vector with state
- */
-template <class SolutionVector, class PvNamesFunc>
-auto loadSolutionFromSequentialVtkFile(const std::string fileName,
-                                       const VTKReader::DataType& dataType,
-                                       PvNamesFunc&& pvNamesFunc,
-                                       SolutionVector& sol)
--> typename std::enable_if_t<decltype(isValid(Detail::hasState())(sol[0]))::value, void>
-{
-    VTKReader vtu(fileName);
-    const auto vec = vtu.readData<std::vector<int>>("phase presence", dataType);
-    std::unordered_set<int> states;
-    for (size_t i = 0; i < sol.size(); ++i)
-    {
-        const int state = vec[i];
-        sol[i].setState(state);
-        states.insert(state);
-    }
-
-    using PrimaryVariables = typename SolutionVector::block_type;
-    using Scalar = typename PrimaryVariables::field_type;
-    for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
-    {
-        if (pvNamesFunc(pvIdx, 1) == pvNamesFunc(pvIdx, 2))
-        {
-            const auto vec = vtu.readData<std::vector<Scalar>>(pvNamesFunc(pvIdx, 1), dataType);
-            for (size_t i = 0; i < sol.size(); ++i)
-                sol[i][pvIdx] = vec[i];
-        }
-        else
-        {
-            std::unordered_map<int, std::vector<Scalar>> switchedPvsSol;
-            for (const auto& state : states)
-                switchedPvsSol[state] = vtu.readData<std::vector<Scalar>>(pvNamesFunc(pvIdx, state), dataType);
-
-            for (size_t i = 0; i < sol.size(); ++i)
-                sol[i][pvIdx] = switchedPvsSol[sol[i].state()][i];
-        }
-    }
-}
-
-/*!
- * \ingroup InputOutput
- * \brief read from a parallel file into a solution vector with state
+ * \brief read from a sequential file into a solution vector with primary variables with state
  */
-template <class SolutionVector, class PvNamesFunc, class FVGridGeometry>
-auto loadSolutionFromParallelVtkFile(const std::string fileName,
-                                     const VTKReader::DataType& dataType,
-                                     PvNamesFunc&& pvNamesFunc,
-                                     SolutionVector& sol,
-                                     const FVGridGeometry& fvGridGeometry)
+template <class SolutionVector, class PvNameFunc, class FVGridGeometry>
+auto loadSolutionFromVtkFile(SolutionVector& sol,
+                             const std::string fileName,
+                             PvNameFunc&& pvNameFunc,
+                             const FVGridGeometry& fvGridGeometry,
+                             const VTKReader::DataType& dataType)
 -> typename std::enable_if_t<decltype(isValid(Detail::hasState())(sol[0]))::value, void>
 {
     VTKReader vtu(fileName);
 
-    // get states at dof locations
+    // get states at each dof location
     const auto stateAtDof = vtu.readData<std::vector<int>>("phase presence", dataType);
 
-    // determine which states are present
+    // determine all states that are present
     std::unordered_set<int> states;
-    for (size_t i = 0; i < sol.size(); ++i)
+    for (size_t i = 0; i < stateAtDof.size(); ++i)
         states.insert(stateAtDof[i]);
 
-    // declare a vector to determine which vertices have been visited already
-    std::vector<bool> visited;
-    using GridView = typename FVGridGeometry::GridView;
-    static constexpr auto dim = GridView::dimension;
-    if (dataType == VTKReader::DataType::pointData)
-        visited.resize(fvGridGeometry.gridView().size(dim));
-
     using PrimaryVariables = typename SolutionVector::block_type;
     using Scalar = typename PrimaryVariables::field_type;
     for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
     {
-        std::unordered_map<int, std::vector<Scalar>> switchedPvsSol;
+        // check if the primary variable is state invariant
+        bool isStateInvariant = true;
         for (const auto& state : states)
-            switchedPvsSol[state] = vtu.readData<std::vector<Scalar>>(pvNamesFunc(pvIdx, state), dataType);
+            isStateInvariant = isStateInvariant && pvNameFunc(pvIdx, state) == pvNameFunc(pvIdx, *states.begin());
+
+        std::unordered_map<int, std::vector<Scalar>> data;
+        if (isStateInvariant)
+            data[0] = vtu.readData<std::vector<Scalar>>(pvNameFunc(pvIdx, *states.begin()), dataType);
+
+        // the primary variable changes with the state
+        // read the data for all occuring states from the file
+        else
+            for (const auto& state : states)
+                data[state] = vtu.readData<std::vector<Scalar>>(pvNameFunc(pvIdx, state), dataType);
 
-        std::fill(visited.begin(), visited.end(), false);
-        std::size_t i = 0;
-        for (const auto& element : elements(fvGridGeometry.gridView(), Dune::Partitions::interior))
+        // sanity check
+        for (const auto& d : data)
+            if (d.second.size() != sol.size())
+                DUNE_THROW(Dune::IOError, "Size mismatch between solution vector and read data (" << sol.size() << " != " << d.second.size() << ")");
+
+        if (dataType == VTKReader::DataType::cellData)
         {
-            if (dataType == VTKReader::DataType::cellData)
+            std::size_t i = 0;
+            for (const auto& element : elements(fvGridGeometry.gridView(), Dune::Partitions::interior))
             {
-                auto eIdx = fvGridGeometry.elementMapper().index(element);
-                auto state = stateAtDof[i];
-                sol[eIdx].setState(state);
-                sol[eIdx][pvIdx] = switchedPvsSol[state][i];
-                ++i;
+                const auto eIdx = fvGridGeometry.elementMapper().index(element);
+                const auto state = isStateInvariant ? 0 : stateAtDof[i];
+                sol[eIdx][pvIdx] = data[state][i++];
+                sol[eIdx][pvIdx].setState(state);
             }
-            else
+        }
+        else
+        {
+            std::size_t i = 0;
+            constexpr int dim = FVGridGeometry::GridView::dimension;
+            std::vector<bool> visited(fvGridGeometry.gridView().size(dim), false);
+            for (const auto& element : elements(fvGridGeometry.gridView(), Dune::Partitions::interior))
             {
                 for (int vIdxLocal = 0; vIdxLocal < element.subEntities(dim); ++vIdxLocal)
                 {
-                    auto vIdxGlobal = fvGridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
+                    const auto vIdxGlobal = fvGridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
                     if (!visited[vIdxGlobal])
                     {
-                        auto state = stateAtDof[i];
-                        sol[vIdxGlobal].setState(state);
-                        sol[vIdxGlobal][pvIdx] = switchedPvsSol[state][i];
+                        const auto state = isStateInvariant ? 0 : stateAtDof[i];
+                        sol[vIdxGlobal][pvIdx] = data[state][i++];
+                        sol[vIdxGlobal][pvIdx].setState(state);
                         visited[vIdxGlobal] = true;
-                        ++i;
                     }
                 }
             }
@@ -293,150 +247,98 @@ auto loadSolutionFromParallelVtkFile(const std::string fileName,
 /*!
  * \ingroup InputOutput
  * \brief helper function to determine the primary variable names of a model with privar state
+ * \note use this as input for the load solution function
  */
 template<class ModelTraits, class FluidSystem>
-std::string primaryVariableName(int pvIdx, int state)
+std::string pvNameWithState(int pvIdx, int state, const std::string& paramGroup = "")
 {
     static auto numStates = (1 << ModelTraits::numPhases()) - 1;
     const auto paramNameWithState = "LoadSolution.PriVarNamesState" + std::to_string(state);
 
-    if (hasParam("LoadSolution.PriVarNames") && !hasParam(paramNameWithState))
-    {
+    if (hasParamInGroup(paramGroup, "LoadSolution.PriVarNames") && !hasParamInGroup(paramGroup, paramNameWithState))
         DUNE_THROW(Dune::NotImplemented, "please provide LoadSolution.PriVarNamesState1..." << numStates
-                   << " or remove LoadSolution.PriVarNames to use default names");
-    }
-    else if (hasParam(paramNameWithState))
-    {
-        const auto pvNames = getParam<std::vector<std::string>>(paramNameWithState);
-        return pvNames[pvIdx];
-    }
-    else
-        return ModelTraits::template primaryVariableName<FluidSystem>(pvIdx, state);
-}
-
-/*!
- * \ingroup InputOutput
- * \brief helper function to determine the cell primary variable names of a staggered model with privar state
- */
-template<class ModelTraits, class FluidSystem>
-std::string primaryVariableNameCell(int pvIdx, int state)
-{
-    static auto numStates = (1 << ModelTraits::numPhases()) - 1;
-    const auto paramNameWithState = "LoadSolution.PriVarNamesCellState" + std::to_string(state);
+                   << " or remove LoadSolution.PriVarNames to use the model's default primary variable names");
 
-    if (hasParam("LoadSolution.PriVarNamesCell") && !hasParam(paramNameWithState))
-    {
-        DUNE_THROW(Dune::NotImplemented, "please provide LoadSolution.PriVarNamesCellState1..." << numStates
-                   << " or remove LoadSolution.PriVarNamesCell to use default names");
-    }
-    else if (hasParam(paramNameWithState))
+    else if (hasParamInGroup(paramGroup, paramNameWithState))
     {
-        const auto pvNames = getParam<std::vector<std::string>>(paramNameWithState);
+        static const auto pvNames = getParamFromGroup<std::vector<std::string>>(paramGroup, paramNameWithState);
         return pvNames[pvIdx];
     }
-    else
-        return ModelTraits::template primaryVariableNameCell<FluidSystem>(pvIdx, state);
-}
 
-/*!
- * \ingroup InputOutput
- * \brief helper function to determine the primary variable names of a model
- */
-template<class ModelTraits>
-std::string primaryVariableName(int pvIdx)
-{
-    if (hasParam("LoadSolution.PriVarNames"))
-    {
-        static auto pvNames = getParam<std::vector<std::string>>("LoadSolution.PriVarNames");
-        return pvNames[pvIdx];
-    }
     else
-        return ModelTraits::primaryVariableName(pvIdx);
-}
-
-/*!
- * \ingroup InputOutput
- * \brief helper function to determine the cell primary variable names of a staggered model
- */
-template<class ModelTraits>
-std::string primaryVariableNameCell(int pvIdx)
-{
-    if (hasParam("LoadSolution.PriVarNamesCell"))
-    {
-        static auto pvNames = getParam<std::vector<std::string>>("LoadSolution.PriVarNamesCell");
-        return pvNames[pvIdx];
-    }
-    else
-        return ModelTraits::primaryVariableNameCell(pvIdx);
+        return ModelTraits::template primaryVariableName<FluidSystem>(pvIdx, state);
 }
 
 /*!
  * \ingroup InputOutput
- * \brief helper function to determine the face primary variable names of a staggered model
+ * \brief helper function to determine the primary variable names of a model without state
+ * \note use this as input for the load solution function
  */
 template<class ModelTraits>
-std::string primaryVariableNameFace(int pvIdx)
+std::string pvName(int pvIdx, int state, const std::string& paramGroup = "")
 {
-    if (hasParam("LoadSolution.PriVarNamesFace"))
+    if (hasParamInGroup(paramGroup, "LoadSolution.PriVarNames"))
     {
-        static auto pvNames = getParam<std::vector<std::string>>("LoadSolution.PriVarNamesFace");
+        static const auto pvNames = getParamFromGroup<std::vector<std::string>>(paramGroup, "LoadSolution.PriVarNames");
         return pvNames[pvIdx];
     }
     else
-        return ModelTraits::primaryVariableNameFace(pvIdx);
+        return ModelTraits::primaryVariableName(pvIdx, 0);
 }
 
-
 /*!
  * \ingroup InputOutput
  * \brief load a solution vector from file
- * \note Supports the following file extensions: *.vtu *.vtp
+ * \note Supports the following file extensions: *.vtu *.vtp *.pvtu, *.pvtp
+ * \param sol the solution vector to read from file
+ * \param fileName the file name of the file to read from
+ * \param pvNameFunc a function with the signature std::string(int pvIdx)
+ *        in case the primary variables have a state the signature is std::string(int pvIdx, int state)
+ * \param fvGridGeometry the grid geometry of the discretization method used
  */
-template <class SolutionVector, class PvNamesFunc, class FVGridGeometry>
-void loadSolution(const std::string& fileName,
-                  DiscretizationMethod discMethod,
-                  PvNamesFunc&& pvNamesFunc,
-                  SolutionVector& sol,
+template <class SolutionVector, class PvNameFunc, class FVGridGeometry>
+void loadSolution(SolutionVector& sol,
+                  const std::string& fileName,
+                  PvNameFunc&& pvNameFunc,
                   const FVGridGeometry& fvGridGeometry)
 {
     const auto extension = fileName.substr(fileName.find_last_of(".") + 1);
-    auto dataType = discMethod == DiscretizationMethod::box ?
+    auto dataType = FVGridGeometry::discMethod == DiscretizationMethod::box ?
                     VTKReader::DataType::pointData : VTKReader::DataType::cellData;
 
     if (extension == "vtu" || extension == "vtp")
     {
-        if (discMethod == DiscretizationMethod::staggered && extension == "vtp")
+        if (FVGridGeometry::discMethod == DiscretizationMethod::staggered && extension == "vtp")
             dataType = VTKReader::DataType::pointData;
 
-        loadSolutionFromSequentialVtkFile(fileName, dataType, pvNamesFunc, sol);
+        loadSolutionFromVtkFile(sol, fileName, pvNameFunc, fvGridGeometry, dataType);
     }
     else if (extension == "pvtu" || extension == "pvtp")
     {
-        if (discMethod == DiscretizationMethod::staggered)
-            DUNE_THROW(Dune::NotImplemented, "reading staggered solution from a parallel Vtk file");
+        if (FVGridGeometry::discMethod == DiscretizationMethod::staggered)
+            DUNE_THROW(Dune::NotImplemented, "reading staggered solution from a parallel vtk file");
 
-        loadSolutionFromParallelVtkFile(fileName, dataType, pvNamesFunc, sol, fvGridGeometry);
+        loadSolutionFromVtkFile(sol, fileName, pvNameFunc, fvGridGeometry, dataType);
     }
     else
-        DUNE_THROW(Dune::NotImplemented, "loadSolution for extension " << extension);
+        DUNE_THROW(Dune::NotImplemented, "loadSolution for file with extension " << extension);
 
-    // synchronize values on ghost and overlap dofs
+    // communicate solution on ghost and overlap dofs
     if (fvGridGeometry.gridView().comm().size() > 1)
     {
-        using PrimaryVariables = typename SolutionVector::block_type;
         using GridView = typename FVGridGeometry::GridView;
         if (dataType == VTKReader::DataType::cellData)
         {
-            LoadSolutionDataHandle<PrimaryVariables, SolutionVector, GridView, 0>
-              dataHandle(sol, fvGridGeometry.gridView());
+            LoadSolutionDataHandle<SolutionVector, typename FVGridGeometry::ElementMapper, 0>
+                dataHandle(sol, fvGridGeometry.elementMapper());
             fvGridGeometry.gridView().communicate(dataHandle,
                                                   Dune::InteriorBorder_All_Interface,
                                                   Dune::ForwardCommunication);
         }
         else
         {
-            LoadSolutionDataHandle<PrimaryVariables, SolutionVector, GridView, GridView::dimension>
-              dataHandle(sol, fvGridGeometry.gridView());
+            LoadSolutionDataHandle<SolutionVector, typename FVGridGeometry::VertexMapper, GridView::dimension>
+                dataHandle(sol, fvGridGeometry.vertexMapper());
             fvGridGeometry.gridView().communicate(dataHandle,
                                                   Dune::InteriorBorder_All_Interface,
                                                   Dune::ForwardCommunication);
@@ -444,6 +346,6 @@ void loadSolution(const std::string& fileName,
     }
 }
 
-} // namespace Dumux
+} // end namespace Dumux
 
 #endif
diff --git a/test/freeflow/navierstokes/CMakeLists.txt b/test/freeflow/navierstokes/CMakeLists.txt
index ebdd77f2bc083cdc70a6139adc19d998a3fe603f..96ec34ce148c7f3398f0d99aa844b84f2d835d80 100644
--- a/test/freeflow/navierstokes/CMakeLists.txt
+++ b/test/freeflow/navierstokes/CMakeLists.txt
@@ -73,14 +73,13 @@ dune_add_test(NAME test_channel_navierstokes
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes -Vtk.WriteFaceData 1")
 
 dune_add_test(NAME test_channel_navierstokes_restart
-              SOURCES test_channel.cc
-              COMPILE_DEFINITIONS ENABLE_NAVIERSTOKES=1
+              TARGET test_channel_navierstokes
               CMAKE_GUARD HAVE_UMFPACK
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/channel-navierstokes-reference.vtu
                                      ${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes_restart-00001.vtu
-                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes -Vtk.WriteFaceData 1 -TimeLoop.DtInitial 1 -Restart.Time 1 -Restart.FileCell test_channel_navierstokes-00001.vtu -Restart.FileFace test_channel_navierstokes-face-00001.vtp  -Problem.Name test_channel_navierstokes_restart")
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_navierstokes -Vtk.WriteFaceData 1 -TimeLoop.DtInitial 1 -Restart.Time 1 -CellCenter.Restart.File test_channel_navierstokes-00001.vtu -Face.Restart.File test_channel_navierstokes-face-00001.vtp  -Problem.Name test_channel_navierstokes_restart")
 
 # the restart test has to run after the test that produces the corresponding vtu file
 set_tests_properties(test_channel_navierstokes_restart PROPERTIES DEPENDS test_channel_navierstokes)
diff --git a/test/freeflow/navierstokes/test_channel.cc b/test/freeflow/navierstokes/test_channel.cc
index 9b85acf622edabb1276458812d033dcef013b38e..6776a9595a2e563c282b03f88ca91567fee1470d 100644
--- a/test/freeflow/navierstokes/test_channel.cc
+++ b/test/freeflow/navierstokes/test_channel.cc
@@ -143,13 +143,14 @@ int main(int argc, char** argv) try
     {
         using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
-        auto fileNameCell = getParam<std::string>("Restart.FileCell");
-        loadSolution(fileNameCell, FVGridGeometry::discMethod,
-                     primaryVariableNameCell<ModelTraits>, x[Dune::index_constant<0>{}], *fvGridGeometry);
-
-        auto fileNameFace = getParam<std::string>("Restart.FileFace");
-        loadSolution(fileNameFace, FVGridGeometry::discMethod,
-                     primaryVariableNameFace<ModelTraits>, x[Dune::index_constant<1>{}], *fvGridGeometry);
+        auto fileNameCell = getParamFromGroup<std::string>("CellCenter", "Restart.File");
+        loadSolution(x[FVGridGeometry::cellCenterIdx()], fileNameCell,
+                     [](int pvIdx){ return "p"; }, // test option with lambda
+                     *fvGridGeometry);
+
+        auto fileNameFace = getParamFromGroup<std::string>("Face", "Restart.File");
+        loadSolution(x[FVGridGeometry::faceIdx()], fileNameFace,
+                     ModelTraits::primaryVariableNameFace, *fvGridGeometry);
     }
     else
         problem->applyInitialSolution(x);
diff --git a/test/porousmediumflow/1p/implicit/test_1pnifv.cc b/test/porousmediumflow/1p/implicit/test_1pnifv.cc
index 9e2abdf602f27b3f35b96b7d5bc190a23911e9d3..22116b726c60e61f5a63e768bed664c1d33e1f7d 100644
--- a/test/porousmediumflow/1p/implicit/test_1pnifv.cc
+++ b/test/porousmediumflow/1p/implicit/test_1pnifv.cc
@@ -132,8 +132,8 @@ int main(int argc, char** argv) try
     if (restartTime > 0)
     {
         using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
-        auto fileName = getParam<std::string>("Restart.File");
-        loadSolution(fileName, FVGridGeometry::discMethod, primaryVariableName<ModelTraits>, x, *fvGridGeometry);
+        const auto fileName = getParam<std::string>("Restart.File");
+        loadSolution(x, fileName, primaryVariableName<ModelTraits>, *fvGridGeometry);
     }
     else
         problem->applyInitialSolution(x);
diff --git a/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc b/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
index 06f13ebd405266964f95baecff8ea8424803614a..adfbcfc2c38992857b26c7cef90911e84846b3ec 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
+++ b/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
@@ -137,8 +137,8 @@ int main(int argc, char** argv) try
     if (restartTime > 0)
     {
         using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
-        auto fileName = getParam<std::string>("Restart.File");
-        loadSolution(fileName, FVGridGeometry::discMethod, primaryVariableName<ModelTraits>, x, *fvGridGeometry);
+        const auto fileName = getParam<std::string>("Restart.File");
+        loadSolution(x, fileName, primaryVariableName<ModelTraits>, *fvGridGeometry);
     }
     else
         problem->applyInitialSolution(x);
diff --git a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
index a1190fff940078988b8c24118dd2e4a11f4f1b9e..134ab6cb94f100971cd696bd9649dbc95387c062 100644
--- a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
+++ b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
@@ -106,8 +106,8 @@ int main(int argc, char** argv) try
     {
         using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
         using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-        auto fileName = getParam<std::string>("Restart.File");
-        loadSolution(fileName, FVGridGeometry::discMethod, primaryVariableName<ModelTraits, FluidSystem>, x, *fvGridGeometry);
+        const auto fileName = getParam<std::string>("Restart.File");
+        loadSolution(x, fileName, primaryVariableName<ModelTraits, FluidSystem>, *fvGridGeometry);
     }
     else
         problem->applyInitialSolution(x);