diff --git a/dumux/io/loadsolution.hh b/dumux/io/loadsolution.hh
index 6c95e18946386b749eeeb9728c93b43cfb530fe6..d267af335d3d48cec48d0628fc95bf7ec429a3ab 100644
--- a/dumux/io/loadsolution.hh
+++ b/dumux/io/loadsolution.hh
@@ -174,7 +174,7 @@ void loadSolution(const std::string& fileName,
 {
     const auto extension = fileName.substr(fileName.find_last_of(".") + 1);
 
-    if (extension == "vtu")
+    if (extension == "vtu" || extension == "vtp")
     {
         const auto dataType = discMethod == DiscretizationMethod::box
                               ? VTUReader::DataType::pointData : VTUReader::DataType::cellData;
diff --git a/dumux/io/vtk/vtureader.hh b/dumux/io/vtk/vtureader.hh
index 5ce9f7e6a6442f1c730fee06bcd3cafa4ff8368d..a4e35815bb2dc796925df834fcf04ce40cae0427 100644
--- a/dumux/io/vtk/vtureader.hh
+++ b/dumux/io/vtk/vtureader.hh
@@ -67,20 +67,30 @@ public:
     Container readData(const std::string& name, const DataType& type)
     {
         using namespace tinyxml2;
-        XMLElement *pieceNode = doc_.FirstChildElement("VTKFile")->FirstChildElement("UnstructuredGrid")->FirstChildElement("Piece");
-        if (pieceNode == nullptr)
-            DUNE_THROW(Dune::IOError, "Couldn't get Piece node in " << fileName_ << ".");
 
-        XMLElement *dataNode = nullptr;
+        XMLElement *dataNode = doc_.FirstChildElement("VTKFile");
+        if (dataNode == nullptr)
+            DUNE_THROW(Dune::IOError, "Couldn't get 'VTKFile' node in " << fileName_ << ".");
+
+        dataNode = dataNode->FirstChildElement("UnstructuredGrid");
+        if (dataNode == nullptr)
+            dataNode = doc_.FirstChildElement("VTKFile")->FirstChildElement("PolyData");
+        if (dataNode == nullptr)
+            DUNE_THROW(Dune::IOError, "Couldn't get 'UnstructuredGrid' or 'PolyData' node in " << fileName_ << ".");
+
+        dataNode = dataNode->FirstChildElement("Piece");
+        if (dataNode == nullptr)
+            DUNE_THROW(Dune::IOError, "Couldn't get 'Piece' node in " << fileName_ << ".");
+
         if (type == DataType::pointData)
-            dataNode = pieceNode->FirstChildElement("PointData");
+            dataNode = dataNode->FirstChildElement("PointData");
         else if (type == DataType::cellData)
-            dataNode = pieceNode->FirstChildElement("CellData");
+            dataNode = dataNode->FirstChildElement("CellData");
         else
             DUNE_THROW(Dune::IOError, "Only cell and point data are supported.");
 
         if (dataNode == nullptr)
-            DUNE_THROW(Dune::IOError, "Couldn't get data node in " << fileName_ << ".");
+            DUNE_THROW(Dune::IOError, "Couldn't get 'PointData' or 'CellData' node in " << fileName_ << ".");
 
         // loop over XML node siblings to find the correct data array
         XMLElement *dataArray = dataNode->FirstChildElement("DataArray");
diff --git a/dumux/porousmediumflow/1p/model.hh b/dumux/porousmediumflow/1p/model.hh
index 581c01acf6c4cac9d94b6f96db1cbb7d4cd76b31..0bcc9c71bf605490ee58f86af564ec9fc282907b 100644
--- a/dumux/porousmediumflow/1p/model.hh
+++ b/dumux/porousmediumflow/1p/model.hh
@@ -74,6 +74,11 @@ struct OnePModelTraits
     static constexpr bool enableAdvection() { return true; }
     static constexpr bool enableMolecularDiffusion() { return false; }
     static constexpr bool enableEnergyBalance() { return false; }
+
+    static std::string primaryVariableName(int pvIdx = 0, int state = 0)
+    {
+        return "pressure";
+    }
 };
 
 /*!
diff --git a/dumux/porousmediumflow/2p/model.hh b/dumux/porousmediumflow/2p/model.hh
index ed4a77e9166fe7bf2ba087db36e8d122e25894d3..6f0d7a893d7c2da284c03281af8605f2f22d3702 100644
--- a/dumux/porousmediumflow/2p/model.hh
+++ b/dumux/porousmediumflow/2p/model.hh
@@ -99,7 +99,7 @@ struct TwoPModelTraits
     static constexpr bool enableMolecularDiffusion() { return false; }
     static constexpr bool enableEnergyBalance() { return false; }
 
-    static std::string primaryVariableName(int pvIdx)
+    static std::string primaryVariableName(int pvIdx, int state = 0)
     {
         if (priVarFormulation() == TwoPFormulation::p0s1)
             return pvIdx == 0 ? "pw" : "Sn";
diff --git a/dumux/porousmediumflow/nonisothermal/model.hh b/dumux/porousmediumflow/nonisothermal/model.hh
index ed6a92516366e3ef27675735f479602b8ea1ceff..e734cba2d760bbc8cd72529b8022d4ef0a6ba8eb 100644
--- a/dumux/porousmediumflow/nonisothermal/model.hh
+++ b/dumux/porousmediumflow/nonisothermal/model.hh
@@ -74,6 +74,14 @@ struct PorousMediumFlowNIModelTraits : public IsothermalTraits
     static constexpr bool enableEnergyBalance() { return true; }
     //! The indices related to the non-isothermal model
     using Indices = EnergyIndices< typename IsothermalTraits::Indices, numEq()>;
+
+    static std::string primaryVariableName(int pvIdx, int state = 0)
+    {
+        if (pvIdx < numEq() - 1)
+            return IsothermalTraits::primaryVariableName(pvIdx, state);
+        else
+            return "temperature";
+    }
 };
 
 } // end namespace Dumux
diff --git a/test/porousmediumflow/1p/implicit/CMakeLists.txt b/test/porousmediumflow/1p/implicit/CMakeLists.txt
index d021037df340eede379d9c43e0ec25fb9212f3da..e93c0df4ac424731ae3ba92eeab87d69065cd50d 100644
--- a/test/porousmediumflow/1p/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/CMakeLists.txt
@@ -63,6 +63,19 @@ dune_add_test(NAME test_1pnibox_convection
                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_1pnibox_convection test_1pnifv_convection.input -Problem.Name 1pnibox_convection"
                        --zeroThreshold {"velocity":1e-15})
 
+dune_add_test(NAME test_1pnibox_convection_restart
+              SOURCES test_1pnifv.cc
+              COMPILE_DEFINITIONS TYPETAG=OnePNIConvectionBoxTypeTag
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS --script fuzzy
+                       --files ${CMAKE_SOURCE_DIR}/test/references/1pniboxconvection-reference.vtp
+                               ${CMAKE_CURRENT_BINARY_DIR}/1pnibox_convection_restart-00005.vtp
+                       --command "${CMAKE_CURRENT_BINARY_DIR}/test_1pnibox_convection test_1pnifv_convection.input -Problem.Name 1pnibox_convection_restart -TimeLoop.DtInitial 1000 -Restart.Time 11761.8 -Restart.File 1pnibox_convection-00006.vtp"
+                       --zeroThreshold {"velocity":1e-15})
+
+# the restart test has to run after the test that produces the corresponding vtu file
+set_tests_properties(test_1pnibox_convection_restart PROPERTIES DEPENDS test_1pnibox_convection)
+
 dune_add_test(NAME test_1pnicctpfa_conduction
               SOURCES test_1pnifv.cc
               COMPILE_DEFINITIONS TYPETAG=OnePNIConductionCCTpfaTypeTag
diff --git a/test/porousmediumflow/1p/implicit/test_1pnifv.cc b/test/porousmediumflow/1p/implicit/test_1pnifv.cc
index 5ca82d837bbef3385eef014b6e47b2fd8ab708fc..56ff9df872669fc0837b8357564c619cf6074a33 100644
--- a/test/porousmediumflow/1p/implicit/test_1pnifv.cc
+++ b/test/porousmediumflow/1p/implicit/test_1pnifv.cc
@@ -51,6 +51,7 @@
 
 #include <dumux/io/vtkoutputmodule.hh>
 #include <dumux/io/grid/gridmanager.hh>
+#include <dumux/io/loadsolution.hh>
 
 /*!
  * \brief Provides an interface for customizing error messages associated with
@@ -122,10 +123,20 @@ int main(int argc, char** argv) try
     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 = getParam<Scalar>("Restart.Time", 0);
+
     // the solution vector
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     SolutionVector x(fvGridGeometry->numDofs());
-    problem->applyInitialSolution(x);
+    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);
+    }
+    else
+        problem->applyInitialSolution(x);
     auto xOld = x;
 
     // the grid variables
@@ -140,7 +151,7 @@ int main(int argc, char** argv) try
     vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
     VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
     vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
-    vtkWriter.write(0.0);
+    vtkWriter.write(restartTime);
 
     // output every vtkOutputInterval time step
     const int vtkOutputInterval = getParam<int>("Problem.OutputInterval");