diff --git a/dumux/io/loadsolution.hh b/dumux/io/loadsolution.hh
new file mode 100644
index 0000000000000000000000000000000000000000..25132ba57d6cd6656acfa245ed3792747efa567a
--- /dev/null
+++ b/dumux/io/loadsolution.hh
@@ -0,0 +1,216 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup InputOutput
+ * \brief read from a file into a solution vector
+ */
+#ifndef DUMUX_LOADSOLUTION_HH
+#define DUMUX_LOADSOLUTION_HH
+
+#include <string>
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/parameters.hh>
+#include <dumux/common/typetraits/isvalid.hh>
+#include <dumux/io/tinyxml2/tinyxml2.h>
+
+namespace Dumux {
+
+//! helper struct detecting if a PrimaryVariables object has a state() function
+struct hasState
+{
+    template<class PrimaryVariables>
+    auto operator()(PrimaryVariables&& priVars)
+    -> decltype(priVars.state())
+    {}
+};
+
+/*!
+ * \ingroup InputOutput
+ * \brief helper class to read from a file into a solution vector
+ */
+class LoadSolutionHelper {
+public:
+    template <class FVGridGeometry, class SolutionVector, class PvNamesFunc>
+    static void loadFromVtkFile(const std::string fileName,
+                                const FVGridGeometry& fvGridGeometry,
+                                PvNamesFunc pvNamesFunc,
+                                SolutionVector& sol)
+    {
+        using namespace tinyxml2;
+
+        XMLDocument xmlDoc;
+        auto eResult = xmlDoc.LoadFile(fileName.c_str());
+        if (eResult != XML_SUCCESS)
+            DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName << ".");
+
+        XMLElement *pieceNode = xmlDoc.FirstChildElement("VTKFile")->FirstChildElement("UnstructuredGrid")->FirstChildElement("Piece");
+        if (pieceNode == nullptr)
+            DUNE_THROW(Dune::IOError, "Couldn't get Piece node in " << fileName << ".");
+
+        XMLElement *dataNode = nullptr;
+        if (FVGridGeometry::discMethod == DiscretizationMethod::box)
+            dataNode = pieceNode->FirstChildElement("PointData");
+        else
+            dataNode = pieceNode->FirstChildElement("CellData");
+
+        if (dataNode == nullptr)
+            DUNE_THROW(Dune::IOError, "Couldn't get data node in " << fileName << ".");
+
+        setPrimaryVariables_(dataNode, fvGridGeometry, pvNamesFunc, sol);
+    }
+
+private:
+
+    template <class Scalar, class FVGridGeometry>
+    static std::vector<Scalar> extractDataToVector_(tinyxml2::XMLElement *xmlNode,
+                                                    const std::string& name,
+                                                    const FVGridGeometry& fvGridGeometry)
+    {
+        // loop over XML node siblings to find the correct data array
+        tinyxml2::XMLElement *dataArray = xmlNode->FirstChildElement("DataArray");
+        for (; dataArray != nullptr; dataArray = dataArray->NextSiblingElement("DataArray"))
+        {
+            const char *attributeText = dataArray->Attribute("Name");
+
+            if (attributeText == nullptr)
+                DUNE_THROW(Dune::IOError, "Couldn't get Name attribute.");
+
+            if (std::string(attributeText) == name)
+                break;
+        }
+        if (dataArray == nullptr)
+            DUNE_THROW(Dune::IOError, "Couldn't find the data array " << name << ".");
+
+        std::stringstream dataStream(dataArray->GetText());
+        std::vector<Scalar> vec(fvGridGeometry.numDofs());
+        for (auto& val : vec)
+        {
+            dataStream >> val;
+        }
+
+        return vec;
+    }
+
+    template<class SolutionVector, class FVGridGeometry, class PvNamesFunc>
+    static auto setPrimaryVariables_(tinyxml2::XMLElement *node,
+                                     const FVGridGeometry& fvGridGeometry,
+                                     PvNamesFunc pvNamesFunc,
+                                     SolutionVector& sol)
+    -> typename std::enable_if_t<!decltype(isValid(hasState())(sol[0]))::value, void>
+    {
+        using PrimaryVariables = typename SolutionVector::block_type;
+        using Scalar = typename PrimaryVariables::field_type;
+
+        for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
+        {
+            auto pvName = pvNamesFunc(pvIdx);
+
+            auto vec = extractDataToVector_<Scalar>(node, pvName, fvGridGeometry);
+
+            for (size_t i = 0; i < sol.size(); ++i)
+                sol[i][pvIdx] = vec[i];
+        }
+    }
+
+    template<class SolutionVector, class FVGridGeometry, class PvNamesFunc>
+    static auto setPrimaryVariables_(tinyxml2::XMLElement *node,
+                                     const FVGridGeometry& fvGridGeometry,
+                                     PvNamesFunc pvNamesFunc,
+                                     SolutionVector& sol)
+    -> typename std::enable_if_t<decltype(isValid(hasState())(sol[0]))::value, void>
+    {
+        auto vec = extractDataToVector_<int>(node, "phase presence", fvGridGeometry);
+
+        int minState = std::numeric_limits<int>::max();
+        int maxState = std::numeric_limits<int>::min();
+        for (size_t i = 0; i < sol.size(); ++i)
+        {
+            int state = vec[i];
+            sol[i].setState(state);
+            using std::min;
+            minState = min(minState, state);
+            using std::max;
+            maxState = max(maxState, 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))
+            {
+                auto vec = extractDataToVector_<Scalar>(node, pvNamesFunc(pvIdx, 1), fvGridGeometry);
+
+                for (size_t i = 0; i < sol.size(); ++i)
+                    sol[i][pvIdx] = vec[i];
+            }
+            else
+            {
+                std::vector<std::vector<Scalar>> switchedPvsSol;
+                for (int state = minState; state <= maxState; ++state)
+                    switchedPvsSol.push_back(extractDataToVector_<Scalar>(node,
+                                                                          pvNamesFunc(pvIdx, state),
+                                                                          fvGridGeometry));
+
+                for (size_t i = 0; i < sol.size(); ++i)
+                    sol[i][pvIdx] = switchedPvsSol[sol[i].state() - minState][i];
+            }
+        }
+    }
+};
+
+template<class ModelTraits, class FluidSystem>
+std::string primaryVariableName(int pvIdx, int state)
+{
+    if (haveParam("LoadSolution.PrimaryVariableNames"))
+        DUNE_THROW(Dune::NotImplemented, "reading PrimaryVariableNames from input file");
+    else
+        return ModelTraits::template primaryVariableName<FluidSystem>(pvIdx, state);
+}
+
+template<class ModelTraits>
+std::string primaryVariableName(int pvIdx)
+{
+    if (haveParam("LoadSolution.PrimaryVariableNames"))
+        DUNE_THROW(Dune::NotImplemented, "reading PrimaryVariableNames from input file");
+    else
+        return ModelTraits::primaryVariableName(pvIdx);
+}
+
+template <class FVGridGeometry, class SolutionVector, class PvNamesFunc>
+void loadSolution(const std::string& fileName,
+                  const FVGridGeometry& fvGridGeometry,
+                  PvNamesFunc pvNamesFunc,
+                  SolutionVector& sol)
+{
+    auto extension = fileName.substr(fileName.find_last_of(".") + 1);
+
+    if (extension == "vtu")
+        LoadSolutionHelper::loadFromVtkFile(fileName, fvGridGeometry, pvNamesFunc, sol);
+    else
+        DUNE_THROW(Dune::NotImplemented, "loadSolution for extension " << extension);
+}
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/io/restart.hh b/dumux/io/restart.hh
index 124dc4bdb7c211cd75bc7725b95b8e02cede177b..b78ceed772b61d0fc18d3c2af0317ded770e908b 100644
--- a/dumux/io/restart.hh
+++ b/dumux/io/restart.hh
@@ -35,20 +35,8 @@
 #include <fstream>
 #include <sstream>
 
-#include <dumux/common/typetraits/isvalid.hh>
-#include <dumux/io/tinyxml2/tinyxml2.h>
-
 namespace Dumux {
 
-//! helper struct detecting if a PrimaryVariables object has a state() function
-struct hasState
-{
-    template<class PrimaryVariables>
-    auto operator()(PrimaryVariables&& priVars)
-    -> decltype(priVars.state())
-    {}
-};
-
 /*!
  * \ingroup InputOutput
  * \brief Load or save a state of a model to/from the harddisk.
@@ -300,150 +288,6 @@ public:
                    "Restart::restartFileList()");
     }
 
-    template <class FVGridGeometry, class SolutionVector>
-    static void loadSolutionFromVtkFile(const FVGridGeometry& fvGridGeometry,
-                                        const std::vector<std::string>& pvNames,
-                                        SolutionVector& sol)
-    {
-        using namespace tinyxml2;
-        auto fileName = getParam<std::string>("Restart.File");
-
-        XMLDocument xmlDoc;
-        auto eResult = xmlDoc.LoadFile(fileName.c_str());
-        if (eResult != XML_SUCCESS)
-            DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName << ".");
-
-        XMLElement *pieceNode = xmlDoc.FirstChildElement("VTKFile")->FirstChildElement("UnstructuredGrid")->FirstChildElement("Piece");
-        if (pieceNode == nullptr)
-            DUNE_THROW(Dune::IOError, "Couldn't get Piece node in " << fileName << ".");
-
-        XMLElement *dataNode = nullptr;
-        if (FVGridGeometry::discMethod == DiscretizationMethod::box)
-            dataNode = pieceNode->FirstChildElement("PointData");
-        else
-            dataNode = pieceNode->FirstChildElement("CellData");
-
-        if (dataNode == nullptr)
-            DUNE_THROW(Dune::IOError, "Couldn't get data node in " << fileName << ".");
-
-        setPrimaryVariables_(dataNode, fvGridGeometry, pvNames, sol);
-    }
-
-private:
-
-    template <class Scalar, class FVGridGeometry>
-    static std::vector<Scalar> extractDataToVector_(tinyxml2::XMLElement *xmlNode,
-                                                    const std::string& name,
-                                                    const FVGridGeometry& fvGridGeometry)
-    {
-        // loop over XML node siblings to find the correct data array
-        tinyxml2::XMLElement *dataArray = xmlNode->FirstChildElement("DataArray");
-        for (; dataArray != nullptr; dataArray = dataArray->NextSiblingElement("DataArray"))
-        {
-            const char *attributeText = dataArray->Attribute("Name");
-
-            if (attributeText == nullptr)
-                DUNE_THROW(Dune::IOError, "Couldn't get Name attribute.");
-
-            if (std::string(attributeText) == name)
-                break;
-        }
-        if (dataArray == nullptr)
-            DUNE_THROW(Dune::IOError, "Couldn't find the data array " << name << ".");
-
-        std::stringstream dataStream(dataArray->GetText());
-        std::vector<Scalar> vec(fvGridGeometry.numDofs());
-        for (auto& val : vec)
-        {
-            dataStream >> val;
-        }
-
-        return vec;
-    }
-
-    template<class SolutionVector, class FVGridGeometry>
-    static auto setPrimaryVariables_(tinyxml2::XMLElement *node,
-                                     const FVGridGeometry& fvGridGeometry,
-                                     std::vector<std::string> pvNames,
-                                     SolutionVector& sol)
-    -> typename std::enable_if_t<!decltype(isValid(hasState())(sol[0]))::value, void>
-    {
-        using PrimaryVariables = typename SolutionVector::block_type;
-        using Scalar = typename PrimaryVariables::field_type;
-
-        for (size_t pvIdx = 0; pvIdx < pvNames.size(); ++pvIdx)
-        {
-            auto pvName = pvNames[pvIdx];
-
-            auto vec = extractDataToVector_<Scalar>(node, pvName, fvGridGeometry);
-
-            for (size_t i = 0; i < sol.size(); ++i)
-                sol[i][pvIdx] = vec[i];
-        }
-    }
-
-    template<class SolutionVector, class FVGridGeometry>
-    static auto setPrimaryVariables_(tinyxml2::XMLElement *node,
-                                     const FVGridGeometry& fvGridGeometry,
-                                     std::vector<std::string> pvNames,
-                                     SolutionVector& sol)
-    -> typename std::enable_if_t<decltype(isValid(hasState())(sol[0]))::value, void>
-    {
-        auto phasePresenceIt = std::find(pvNames.begin(), pvNames.end(), "phase presence");
-        if (phasePresenceIt != pvNames.end())
-        {
-            auto vec = extractDataToVector_<int>(node, "phase presence", fvGridGeometry);
-
-            for (size_t i = 0; i < sol.size(); ++i)
-                sol[i].setState(vec[i]);
-
-            pvNames.erase(phasePresenceIt);
-        }
-
-        using PrimaryVariables = typename SolutionVector::block_type;
-        using Scalar = typename PrimaryVariables::field_type;
-        for (size_t pvIdx = 0; pvIdx < pvNames.size(); ++pvIdx)
-        {
-            auto pvName = pvNames[pvIdx];
-
-            auto switchedPvNames = getSwitchedPvNames_(pvName);
-
-            if (switchedPvNames[0] == pvName)
-            {
-                auto vec = extractDataToVector_<Scalar>(node, pvName, fvGridGeometry);
-
-                for (size_t i = 0; i < sol.size(); ++i)
-                    sol[i][pvIdx] = vec[i];
-            }
-            else
-            {
-                std::vector<std::vector<Scalar>> switchedPvsSol;
-                for (size_t switchedPvIdx = 0; switchedPvIdx < switchedPvNames.size(); ++switchedPvIdx)
-                    switchedPvsSol.push_back(extractDataToVector_<Scalar>(node,
-                                                                          switchedPvNames[switchedPvIdx],
-                                                                          fvGridGeometry));
-
-                for (size_t i = 0; i < sol.size(); ++i)
-                    sol[i][pvIdx] = switchedPvsSol[sol[i].state()-1][i];
-            }
-        }
-    }
-
-    static std::vector<std::string> getSwitchedPvNames_(std::string pvName)
-    {
-        std::vector<std::string> switchedPvNames;
-        int lastPos = -1;
-        size_t pos = 0;
-        do
-        {
-            pos = pvName.find('/', lastPos + 1);
-            switchedPvNames.push_back(pvName.substr(lastPos+1, pos - lastPos - 1));
-            lastPos = pos;
-        } while (pos != std::string::npos);
-
-        return switchedPvNames;
-    }
-
     std::string fileName_;
     std::ifstream inStream_;
     std::ofstream outStream_;
diff --git a/dumux/porousmediumflow/2p/model.hh b/dumux/porousmediumflow/2p/model.hh
index 73516e6f1c42e78e33d8a7017d03ed73b18f6253..ed4a77e9166fe7bf2ba087db36e8d122e25894d3 100644
--- a/dumux/porousmediumflow/2p/model.hh
+++ b/dumux/porousmediumflow/2p/model.hh
@@ -98,20 +98,13 @@ struct TwoPModelTraits
     static constexpr bool enableAdvection() { return true; }
     static constexpr bool enableMolecularDiffusion() { return false; }
     static constexpr bool enableEnergyBalance() { return false; }
-};
 
-template<TwoPFormulation formulation>
-struct TwoPPrimaryVariableNames
-{
-    static std::vector<std::string> get()
+    static std::string primaryVariableName(int pvIdx)
     {
-        switch (formulation)
-        {
-            case TwoPFormulation::p0s1:
-                return {"pw", "Sn"};
-            case TwoPFormulation::p1s0:
-                return {"pn", "Sw"};
-        }
+        if (priVarFormulation() == TwoPFormulation::p0s1)
+            return pvIdx == 0 ? "pw" : "Sn";
+        else
+            return pvIdx == 0 ? "pn" : "Sw";
     }
 };
 
@@ -164,9 +157,6 @@ NEW_TYPE_TAG(TwoPNI, INHERITS_FROM(TwoP));
 SET_PROP(TwoP, Formulation)
 { static constexpr auto value = TwoPFormulation::p0s1; };
 
-NEW_PROP_TAG(PrimaryVariableNames);
-SET_TYPE_PROP(TwoP, PrimaryVariableNames, TwoPPrimaryVariableNames<GET_PROP_VALUE(TypeTag, Formulation)>);
-
 SET_TYPE_PROP(TwoP, LocalResidual, ImmiscibleLocalResidual<TypeTag>);         //!< Use the immiscible local residual operator for the 2p model
 
 //! The model traits class
diff --git a/dumux/porousmediumflow/2p2c/model.hh b/dumux/porousmediumflow/2p2c/model.hh
index 5ed84ac163bef2819716b60d9ae44aa256e85ab7..2ec91b72e7e0a9168fd3a924c8fea12e430744ab 100644
--- a/dumux/porousmediumflow/2p2c/model.hh
+++ b/dumux/porousmediumflow/2p2c/model.hh
@@ -77,48 +77,51 @@
 #ifndef DUMUX_2P2C_MODEL_HH
 #define DUMUX_2P2C_MODEL_HH
 
-#include <dune/common/fvector.hh>
+#include <array>
 
 // property forward declarations
 #include <dumux/common/properties.hh>
-#include <dumux/porousmediumflow/properties.hh>
+
+#include <dumux/porousmediumflow/2pnc/model.hh>
 #include <dumux/porousmediumflow/2p/formulation.hh>
 #include <dumux/porousmediumflow/nonisothermal/model.hh>
-#include <dumux/porousmediumflow/nonisothermal/indices.hh>
 #include <dumux/porousmediumflow/nonisothermal/vtkoutputfields.hh>
-
-#include <dumux/material/spatialparams/fv.hh>
 #include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
-#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
-
-#include <dumux/porousmediumflow/compositional/localresidual.hh>
-#include <dumux/porousmediumflow/compositional/switchableprimaryvariables.hh>
-
-#include <dumux/porousmediumflow/2pnc/model.hh>
 
 #include "volumevariables.hh"
 
 namespace Dumux {
 
-template<TwoPFormulation formulation, class FluidSystem>
-struct TwoPTwoCPrimaryVariableNames
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief Specifies a number properties of two-phase two-component models.
+ *
+ * \tparam f The two-phase formulation used
+ * \tparam useM Boolean to specify if moles or masses are balanced
+ * \tparam replCompEqIdx the equation which is replaced by the total mass balance (none if replCompEqIdx >= numComponents)
+ */
+template<TwoPFormulation formulation, bool useMol, int replCompEqIdx = 2>
+struct TwoPTwoCModelTraits : public TwoPNCModelTraits</*numComps=*/2, useMol, /*setMFracForFirstPhase=*/true, formulation, replCompEqIdx>
 {
-    static std::vector<std::string> get()
+    template <class FluidSystem>
+    static std::string primaryVariableName(int pvIdx, int state)
     {
-        const auto p0s1SwitchedPvNames =
-          "x_" + FluidSystem::phaseName(0) + "^" + FluidSystem::componentName(1)
-          + "/x_" + FluidSystem::phaseName(1) + "^" + FluidSystem::componentName(0)
-          + "/Sn";
-        const auto p1s0SwitchedPvNames =
-          "x_" + FluidSystem::phaseName(0) + "^" + FluidSystem::componentName(1)
-          + "/x_" + FluidSystem::phaseName(1) + "^" + FluidSystem::componentName(0)
-          + "/Sw";
+        static const std::string xString = useMol ? "x" : "X";
+        static const std::array<std::string, 3> p0s1SwitchedPvNames = {{
+            xString + "_" + FluidSystem::phaseName(0) + "^" + FluidSystem::componentName(1),
+            xString + "_" + FluidSystem::phaseName(1) + "^" + FluidSystem::componentName(0),
+            "S_n"}};
+        static const std::array<std::string, 3> p1s0SwitchedPvNames = {{
+            xString + "_" + FluidSystem::phaseName(0) + "^" + FluidSystem::componentName(1),
+            xString + "_" + FluidSystem::phaseName(1) + "^" + FluidSystem::componentName(0),
+            "S_w"}};
+
         switch (formulation)
         {
-            case TwoPFormulation::p0s1:
-                return {"pw", p0s1SwitchedPvNames, "phase presence"};
-            case TwoPFormulation::p1s0:
-                return {"pn", p1s0SwitchedPvNames, "phase presence"};
+        case TwoPFormulation::p0s1:
+            return pvIdx == 0 ? "p_w" : p0s1SwitchedPvNames[state-1];
+        case TwoPFormulation::p1s0:
+            return pvIdx == 0 ? "p_n" : p1s0SwitchedPvNames[state-1];
         }
     }
 };
@@ -135,14 +138,20 @@ NEW_TYPE_TAG(TwoPTwoCNI, INHERITS_FROM(TwoPTwoC));
 // Property values
 //////////////////////////////////////////////////////////////////
 
-NEW_PROP_TAG(PrimaryVariableNames);
-SET_PROP(TwoPTwoC, PrimaryVariableNames)
+/*!
+ * \brief Set the model traits property.
+ */
+SET_PROP(TwoPTwoC, ModelTraits)
 {
 private:
-    static constexpr auto formulation = GET_PROP_VALUE(TypeTag, Formulation);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    static_assert(FluidSystem::numComponents == 2, "Only fluid systems with 2 components are supported by the 2p-2c model!");
+    static_assert(FluidSystem::numPhases == 2, "Only fluid systems with 2 phases are supported by the 2p-2c model!");
+
 public:
-    using type = TwoPTwoCPrimaryVariableNames<formulation, FluidSystem>;
+    using type = TwoPTwoCModelTraits< GET_PROP_VALUE(TypeTag, Formulation),
+                                      GET_PROP_VALUE(TypeTag, UseMoles),
+                                      GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx) >;
 };
 
 //! Use the 2p2c VolumeVariables
@@ -182,11 +191,9 @@ private:
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static_assert(FluidSystem::numComponents == 2, "Only fluid systems with 2 components are supported by the 2p2c model!");
     static_assert(FluidSystem::numPhases == 2, "Only fluid systems with 2 phases are supported by the 2p2c model!");
-    using IsothermalTraits = TwoPNCModelTraits<FluidSystem::numComponents,
-                                               GET_PROP_VALUE(TypeTag, UseMoles),
-                                               /*setMoleFractionsForFirstPhase=*/true,
-                                               GET_PROP_VALUE(TypeTag, Formulation),
-                                               GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx)>;
+    using IsothermalTraits = TwoPTwoCModelTraits<GET_PROP_VALUE(TypeTag, Formulation),
+                                                 GET_PROP_VALUE(TypeTag, UseMoles),
+                                                 GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx)>;
 public:
     using type = PorousMediumFlowNIModelTraits<IsothermalTraits>;
 };
diff --git a/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc b/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
index 6419e8b9be5efec6ccb6adee55c2f85ef9ec33a4..634dbe3b821f511570a116d146617edec2d06c1d 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
+++ b/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
@@ -51,7 +51,7 @@
 
 #include <dumux/io/vtkoutputmodule.hh>
 #include <dumux/io/grid/gridmanager.hh>
-#include <dumux/io/restart.hh>
+#include <dumux/io/loadsolution.hh>
 
 /*!
  * \brief Provides an interface for customizing error messages associated with
@@ -136,8 +136,9 @@ int main(int argc, char** argv) try
     SolutionVector x(fvGridGeometry->numDofs());
     if (restartTime > 0)
     {
-        using PvNames = typename GET_PROP_TYPE(TypeTag, PrimaryVariableNames);
-        Restart::loadSolutionFromVtkFile(*fvGridGeometry, PvNames::get(), x);
+        using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+        auto fileName = getParam<std::string>("Restart.File");
+        loadSolution(fileName, *fvGridGeometry, primaryVariableName<ModelTraits>, x);
     }
     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 fa0fe82b67d3bb52ddf8742baab0d9938a11ca26..f55102ebfd39a11cbbb74e96151e9323c4aef41a 100644
--- a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
+++ b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
@@ -47,7 +47,7 @@
 
 #include <dumux/io/vtkoutputmodule.hh>
 #include <dumux/io/grid/gridmanager.hh>
-#include <dumux/io/restart.hh>
+#include <dumux/io/loadsolution.hh>
 
 // the problem definitions
 #include "injectionproblem.hh"
@@ -104,8 +104,10 @@ int main(int argc, char** argv) try
     SolutionVector x(fvGridGeometry->numDofs());
     if (restartTime > 0)
     {
-        using PvNames = typename GET_PROP_TYPE(TypeTag, PrimaryVariableNames);
-        Restart::loadSolutionFromVtkFile(*fvGridGeometry, PvNames::get(), x);
+        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, primaryVariableName<ModelTraits, FluidSystem>, x);
     }
     else
         problem->applyInitialSolution(x);