[io] add capability to restart from a VTK file

Add the TinyXML-2 parser. Add functions in the `Restart` class
that read from a given VTK file into a `SolutionVector`. The
names of the primary variables in the VTK file are handed over by
the problem. Adapt `applyInitialSolution` of the FVProblem such
that it invokes the restart functionality if a start time larger
than 0 is handed over.

dumux/common/fvproblem.hh

@@ -476,10 +476,18 @@ public:
      * \brief Applies the initial solution for all degrees of freedom of the grid.
      * \param sol the initial solution vector
      */
-    void applyInitialSolution(SolutionVector& sol) const
+    void applyInitialSolution(SolutionVector& sol, Scalar restartTime = 0) const
     {
-        // set the initial values by forwarding to a specialized method
-        applyInitialSolutionImpl_(sol, std::integral_constant<bool, isBox>());
+        if (restartTime > 0)
+        {
+            // set initial values for a restarted simulation
+            asImp_().applyRestartSolution(sol);
+        }
+        else
+        {
+            // set the initial values by forwarding to a specialized method
+            applyInitialSolutionImpl_(sol, std::integral_constant<bool, isBox>());
+        }
     }
 
     /*!

dumux/io/restart.hh

@@ -35,6 +35,8 @@
 #include <fstream>
 #include <sstream>
 
+#include <dumux/io/tinyxml2/tinyxml2.h>
+
 namespace Dumux {
 /*!
  * \ingroup InputOutput
@@ -287,6 +289,67 @@ public:
                    "Restart::restartFileList()");
     }
 
+    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 FVGridGeometry, class SolutionVector>
+    static void loadSolutionFromVtkFile(const FVGridGeometry& fvGridGeometry,
+                                        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.");
+
+        XMLElement *pieceNode = xmlDoc.FirstChildElement("VTKFile")->FirstChildElement("UnstructuredGrid")->FirstChildElement("Piece");
+        if (pieceNode == nullptr)
+            DUNE_THROW(Dune::IOError, "Couldn't get Piece node.");
+
+        XMLElement *cellDataNode = pieceNode->FirstChildElement("CellData");
+        if (cellDataNode == nullptr)
+            DUNE_THROW(Dune::IOError, "Couldn't get CellData node.");
+
+        using PrimaryVariables = typename SolutionVector::block_type;
+        using Scalar = typename PrimaryVariables::field_type;
+        for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
+        {
+            auto vec = extractDataToVector<Scalar>(cellDataNode, pvNames[pvIdx], fvGridGeometry);
+
+            for (size_t i = 0; i < sol.size(); ++i)
+                sol[i][pvIdx] = vec[i];
+        }
+    }
 
 private:
     std::string fileName_;

test/porousmediumflow/2p/implicit/incompressible/problem.hh

@@ -29,6 +29,8 @@
 #include <dumux/discretization/cellcentered/tpfa/properties.hh>
 #include <dumux/discretization/cellcentered/mpfa/properties.hh>
 
+#include <dumux/io/restart.hh>
+
 #include <dumux/material/components/trichloroethene.hh>
 #include <dumux/material/components/simpleh2o.hh>
 #include <dumux/material/fluidsystems/1pliquid.hh>
@@ -230,6 +232,12 @@ public:
         return 293.15; // 10°C
     }
 
+    template <class SolutionVector>
+    void applyRestartSolution(SolutionVector& sol) const
+    {
+        Restart::loadSolutionFromVtkFile(this->fvGridGeometry(), {"pw", "Sn"}, sol);
+    }
+
 private:
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
     {

test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc

@@ -25,6 +25,7 @@
 
 #include <ctime>
 #include <iostream>
+#include <sstream>
 
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/timer.hh>
@@ -120,10 +121,19 @@ int main(int argc, char** argv) try
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     auto problem = std::make_shared<Problem>(fvGridGeometry);
 
+    // get some time loop parameters
+    using Scalar = typename GET_PROP_TYPE(TypeTag, 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 = getParam<Scalar>("Restart.Time", 0);
+
     // the solution vector
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     SolutionVector x(fvGridGeometry->numDofs());
-    problem->applyInitialSolution(x);
+    problem->applyInitialSolution(x, restartTime);
     auto xOld = x;
 
     // maybe update the interface parameters
@@ -135,17 +145,6 @@ int main(int argc, char** argv) try
     auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
     gridVariables->init(x, xOld);
 
-    // get some time loop parameters
-    using Scalar = typename GET_PROP_TYPE(TypeTag, 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");
-
     // intialize the vtk output module
     using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
 
@@ -156,7 +155,7 @@ int main(int argc, char** argv) try
     using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput);
     vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
     VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
-    vtkWriter.write(0.0);
+    vtkWriter.write(restartTime);
 
     // instantiate time loop
     auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);