changed: redo restarting support

store restart data separate from visualization data.
this has several advantages:
- no unusual striding in visualization files due to multistep methods
- we can store data which are not on a per-control-point basis
- simplify logic

Apps/Common/AppCommon.h

@@ -37,53 +37,21 @@ namespace SIM
     int dim; //!< Dimensionality of simulation
   };
 
-  //! \brief Handles application restarts.
-  //! \param[in] simulator The top SIMbase instance of your application
-  //! \param[in] solver The SIMSolver instance of your application
-  //! \param[in] restartfile The file to read from
-  //! \param[in] interval The stride in the input file
-  //! \param[in] steps The number of time steps to load
-  template<class Simulator, class Solver>
-  void handleRestart(Simulator& simulator, Solver& solver,
-                     const std::string& restartfile,
-                     int interval = 1, int steps = 1)
-  {
-    DataExporter reader(true,interval,steps);
-    XMLWriter* xml = new XMLWriter(restartfile,solver.getProcessAdm());
-    HDF5Writer* hdf = new HDF5Writer(restartfile,solver.getProcessAdm(),true);
-    reader.registerWriter(xml);
-    reader.registerWriter(hdf);
-    simulator.registerFields(reader);
-    int max = reader.getTimeLevel();
-    // correct loaded level if we stopped in the middle of a "stride" level
-    if ((max+1) % (steps+1))
-      max -= (max % (steps+1))+steps;
-    double time;
-    hdf->openFile(max-steps);
-    hdf->readDouble(max-steps,"timeinfo","SIMbase-1",time);
-    solver.fastForward(time/solver.getTimePrm().time.dt);
-    for (int i=steps;i>=0;--i) {
-      reader.loadTimeLevel(max-i,xml,hdf);
-      solver.postSolve(solver.getTimePrm(),true);
-      if (i > 0) solver.advanceStep();
-    }
-    xml->writeTimeInfo(0, interval, steps, solver.getTimePrm());
-  }
-
   //! \brief Handles application data output.
   //! \param[in] simulator The top SIMbase instance of your application
   //! \param[in] solver The SIMSolver instance of your application
   //! \param[in] hdf5file The file to save to
   //! \param[in] append Whether or not to append to file
   //! \param[in] interval The stride in the output file
-  //! \param[in] steps The number of time steps to dump in onw row
+  //! \param[in] restartInterval The stride in the restart file
   template<class Simulator, class Solver>
   DataExporter* handleDataOutput(Simulator& simulator, Solver& solver,
                                  const std::string& hdf5file,
                                  bool append = false,
-                                 int interval = 1, int steps = 1)
+                                 int interval = 1,
+                                 int restartInterval = 0)
   {
-    DataExporter* writer = new DataExporter(true,interval,steps);
+    DataExporter* writer = new DataExporter(true,interval,restartInterval);
     XMLWriter* xml = new XMLWriter(hdf5file,solver.getProcessAdm());
     HDF5Writer* hdf = new HDF5Writer(hdf5file,solver.getProcessAdm(),append);
     writer->registerWriter(xml);

Apps/Common/SIMCoupled.h

@@ -158,6 +158,20 @@ public:
     return result;
   }
 
+  //! \brief Serialize internal state for restarting purposes.
+  //! \param data Container for serialized data
+  bool serialize(DataExporter::SerializeData& data)
+  {
+    return S1.serialize(data) && S2.serialize(data);
+  }
+
+  //! \brief Set internal state from a serialized state.
+  //! \param[in] data Container for serialized data
+  bool deSerialize(const DataExporter::SerializeData& data)
+  {
+    return S1.deSerialize(data) && S2.deSerialize(data);
+  }
+
 protected:
   T1& S1; //!< First substep
   T2& S2; //!< Second substep

Apps/Common/SIMExplicitRK.h

@@ -13,6 +13,7 @@
 #ifndef SIM_EXPLICIT_RK_H_
 #define SIM_EXPLICIT_RK_H_
 
+#include "DataExporter.h"
 #include "TimeIntUtils.h"
 #include "TimeStep.h"
 
@@ -141,6 +142,20 @@ public:
     return solver.saveStep(tp, nBlock);
   }
 
+  //! \brief Serialize internal state for restarting purposes.
+  //! \param data Container for serialized data
+  bool serialize(DataExporter::SerializeData& data)
+  {
+    return solver.serialize(data);
+  }
+
+  //! \brief Set internal state from a serialized state.
+  //! \param[in] data Container for serialized data
+  bool deSerialize(const DataExporter::SerializeData& data)
+  {
+    return solver.deSerialize(data);
+  }
+
 protected:
   Solver& solver; //!< Reference to simulator
   RKTableaux RK;  //!< Tableaux of Runge-Kutta coefficients

Apps/Common/SIMSemi3D.h

@@ -131,8 +131,7 @@ public:
     int plane = 1 + startCtx;
     for (size_t i = 0; i < m_planes.size(); i++, plane++)
       if (plane_exporters.size() <= i) {
-        DataExporter* exp = new DataExporter(true,exporter.getStride(),
+        DataExporter* exp = new DataExporter(true,exporter.getStride());
-                                             exporter.getOrder());
         std::stringstream str;
         str << "_plane" << plane;
         XMLWriter* xml = new XMLWriter(name+str.str(),
@@ -172,6 +171,12 @@ public:
     return true;
   }
 
+  //! \brief No serialization support.
+  bool serialize(DataExporter::SerializeData& data) { return false; }
+
+  //! \brief No deserialization support.
+  bool deSerialize(const DataExporter::SerializeData& data) { return false; }
+
   //! \brief Solves the nonlinear equations by Newton-Raphson iterations.
   bool solveStep(TimeStep& tp)
   {

Apps/Common/SIMSolver.h

@@ -15,8 +15,8 @@
 #define _SIM_SOLVER_H_
 
 #include "SIMadmin.h"
-#include "DataExporter.h"
 #include "TimeStep.h"
+#include "HDF5Writer.h"
 #include "IFEM.h"
 #include "tinyxml.h"
 
@@ -100,6 +100,20 @@ public:
     return 0;
   }
 
+  //! \brief Serialize internal state for restarting purposes.
+  //! \param data Container for serialized data
+  bool serialize(DataExporter::SerializeData& data)
+  {
+    return tp.serialize(data) && S1.serialize(data);
+  }
+
+  //! \brief Set internal state from a serialized state.
+  //! \param[in] data Container for serialized data
+  bool deSerialize(const DataExporter::SerializeData& data)
+  {
+    return tp.deSerialize(data) && S1.deSerialize(data);
+  }
+
 protected:
   //! \brief Parses a data section from an input stream.
   virtual bool parse(char* keyw, std::istream& is) { return tp.parse(keyw,is); }
@@ -143,12 +157,44 @@ protected:
     if (saveRes && !S1.saveStep(tp,nBlock))
       return false;
 
-    if (saveRes && exporter)
+    if (saveRes && exporter) {
-      exporter->dumpTimeLevel(&tp,newMesh);
+      DataExporter::SerializeData data;
+      if (exporter->dumpForRestart(&tp) && this->serialize(data))
+        return exporter->dumpTimeLevel(&tp,newMesh,&data);
+      else // no restart dump, or serialization failure
+        return exporter->dumpTimeLevel(&tp,newMesh);
+    }
 
     return true;
   }
 
+public:
+  //! \brief Handles application restarts by reading a serialized solver state.
+  //! \param[in] restartFile File to read restart state from
+  //! \param[in] restartStep Index of the time step to read restart state for
+  //! \return One-based time step index of the restart state read.
+  //! If zero, no restart specified. If negative, read failure.
+  int restart(const std::string& restartFile, int restartStep)
+  {
+    if (restartFile.empty()) return 0;
+
+    DataExporter::SerializeData data;
+    HDF5Writer hdf(restartFile,adm,true);
+    if ((restartStep = hdf.readRestartData(data,restartStep)) >= 0)
+    {
+      IFEM::cout <<"\n === Restarting from a serialized state ==="
+                 <<"\n     file = "<< restartFile
+                 <<"\n     step = "<< restartStep << std::endl;
+      if (this->deSerialize(data))
+        return restartStep+1;
+      else
+        restartStep = -2;
+    }
+
+    std::cerr <<" *** SIMSolver: Failed to read restart data."<< std::endl;
+    return restartStep;
+  }
+
 private:
   bool saveDivergedSol; //!< If \e true, save also the diverged solution to VTF
 

cmake/Modules/FindIFEMDeps.cmake

@@ -222,6 +222,17 @@ if(IFEM_USE_ISTL)
   endif()
 endif()
 
+# Cereal
+if(IFEM_USE_CEREAL)
+  find_path(CEREAL_INCLUDE_DIRS NAMES cereal/cereal.hpp)
+  if(CEREAL_INCLUDE_DIRS)
+    list(APPEND IFEM_DEPINCLUDES ${CEREAL_INCLUDE_DIRS})
+    list(APPEND IFEM_DEFINITIONS -DHAS_CEREAL=1)
+    message(STATUS "Cereal serialization support enabled")
+    set(CEREAL_FOUND 1)
+  endif()
+endif()
+
 # Portability issues
 include(CheckFunctionExists)
 set(CMAKE_REQUIRED_DEFINITIONS)

cmake/Modules/IFEMOptions.cmake

@@ -11,6 +11,7 @@ OPTION(IFEM_USE_SAMG           "Compile with SAMG support?"           OFF)
 OPTION(IFEM_USE_HDF5           "Compile with HDF5 support?"           ON)
 OPTION(IFEM_USE_VTFWRITER      "Compile with VTFWriter support?"      ON)
 OPTION(IFEM_USE_UMFPACK        "Compile with UMFPACK support?"        ON)
+OPTION(IFEM_USE_CEREAL         "Compile with cereal support?"         ON)
 OPTION(IFEM_AS_SUBMODULE       "Compile IFEM as a submodule of apps?" OFF)
 OPTION(IFEM_WHOLE_PROG_OPTIM   "Compile IFEM with link-time optimizations?" OFF)
 OPTION(IFEM_TEST_MEMCHECK      "Run tests through valgrind?"          OFF)

cmake/Templates/IFEMConfig.cmake.in

@@ -1,7 +1,7 @@
 if(NOT IFEM_FOUND)
   set(IFEM_INCLUDE_DIRS @IFEM_INCLUDE_DIRS@)
   set(IFEM_LIBRARIES -L@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_LIBDIR@ -lIFEM @IFEM_DEPLIBS@)
-  set(IFEM_DEFINITIONS "@IFEM_DEFINITIONS@")
+  set(IFEM_DEFINITIONS @IFEM_DEFINITIONS@)
   set(IFEM_CXX_FLAGS "@IFEM_CXX_FLAGS@")
 
   list(APPEND CMAKE_MODULE_PATH @CMAKE_INSTALL_PREFIX@/lib/IFEM)
@@ -12,6 +12,7 @@ if(NOT IFEM_FOUND)
   set(ISTL_FOUND @ISTL_FOUND@)
   set(MPI_FOUND @MPI_FOUND@)
   set(HDF5_FOUND @HDF5_FOUND@)
+  set(CEREAL_FOUND @CEREAL_FOUND@)
 
   set(IFEM_REGTEST_SCRIPT ${IFEM_PATH}/scripts/regtest.sh.in)
   set(IFEM_CHECKCOMMITS_SCRIPT ${IFEM_PATH}/CheckCommits.cmake)

src/SIM/SIMinput.C

@@ -464,6 +464,8 @@ bool SIMinput::parse (const TiXmlElement* elem)
   bool result = true;
   if (!strcasecmp(elem->Value(),"discretization"))
     result = opt.parseDiscretizationTag(elem);
+  else if (!strcasecmp(elem->Value(),"restart"))
+    result = opt.parseRestartTag(elem);
   else if (!strcasecmp(elem->Value(),"initialcondition"))
     result = this->parseICTag(elem);
   else if (!strcasecmp(elem->Value(),"linearsolver"))
@@ -1176,3 +1178,11 @@ bool SIMinput::hasIC (const std::string& name) const
 
   return false;
 }
+
+
+bool SIMinput::deSerialize (const std::map<std::string,std::string>&)
+{
+  std::cerr <<" *** SIMinput::deSerialize: Must be implemented in sub-class.\n"
+            <<"     Restart not supported for "<< this->getName() << std::endl;
+  return false;
+}

src/SIM/SIMinput.h

@@ -222,6 +222,9 @@ public:
   //! nullptr means this.
   bool setInitialConditions(SIMdependency* fieldHolder = nullptr);
 
+  //! \brief Deserialization support (for simulation restart).
+  virtual bool deSerialize(const std::map<std::string,std::string>&);
+
 private:
   //! \brief Sets initial conditions from a file.
   //! \param fieldHolder The SIM-object to inject the initial conditions into

src/SIM/SIMoptions.C

@@ -49,6 +49,7 @@ SIMoptions::SIMoptions ()
   pSolOnly = false;
   enableController = false;
   restartInc = 0;
+  restartStep = -1;
 
   nGauss[0] = nGauss[1] = 4;
   nViz[0] = nViz[1] = nViz[2] = 2;
@@ -172,6 +173,17 @@ bool SIMoptions::parseOutputTag (const TiXmlElement* elem)
 }
 
 
+bool SIMoptions::parseRestartTag (const TiXmlElement* elem)
+{
+  if (!strcasecmp(elem->Value(),"restart")) {
+    utl::getAttribute(elem,"file",restartFile);
+    utl::getAttribute(elem,"step",restartStep);
+  }
+
+  return true;
+}
+
+
 bool SIMoptions::parseConsoleTag (const TiXmlElement* elem)
 {
   if (!strcasecmp(elem->Value(),"logging")) {
@@ -269,6 +281,14 @@ bool SIMoptions::parseOldOptions (int argc, char** argv, int& i)
   }
   else if (!strcmp(argv[i],"-saveInc") && i < argc-1)
     dtSave = atof(argv[++i]);
+  else if (!strcmp(argv[i],"-restart") && i < argc-1)
+  {
+    restartFile = strtok(argv[++i],".");
+    if (i+1 < argc && argv[i+1][0] != '-')
+      restartStep = atoi(argv[++i]);
+  }
+  else if (!strcmp(argv[i],"-restartInc") && i < argc-1)
+    restartInc = atoi(argv[++i]);
   else if (!strcmp(argv[i],"-eig") && i < argc-1)
     eig = atoi(argv[++i]);
   else if (!strcmp(argv[i],"-nev") && i < argc-1)

src/SIM/SIMoptions.h

@@ -51,6 +51,8 @@ public:
   bool parseEigSolTag(const TiXmlElement* elem);
   //! \brief Parses a subelement of the \a resultoutput XML-tag.
   bool parseOutputTag(const TiXmlElement* elem);
+  //! \brief Parses the \a restart XML-tag.
+  bool parseRestartTag(const TiXmlElement* elem);
   //! \brief Parses a projection method XML-tag.
   bool parseProjectionMethod(const char* ptype);
 
@@ -83,11 +85,16 @@ public:
   int format;    //!< VTF-file format (-1=NONE, 0=ASCII, 1=BINARY)
   int nViz[3];   //!< Number of visualization points over each knot-span
   int saveInc;   //!< Number of load/time increments between each result output
-  int restartInc;//!< Number of load/time increments between each restart output
   double dtSave; //!< Time interval between each result output
   bool pSolOnly; //!< If \e true, don't save secondary solution variables
 
   std::string hdf5; //!< Prefix for HDF5-file
+
+  // Restart options
+  int         restartInc;  //!< Number of increments between each restart output
+  int         restartStep; //!< Index to the actual state to restart from
+  std::string restartFile; //!< File to read restart state data from
+
   bool enableController; //!< Whether or not to enable external program control
 
   int printPid; //!< PID to print info to screen for

src/SIM/SIMoutput.h

@@ -269,6 +269,9 @@ public:
   //! \brief Checks whether point result files have been defined or not.
   bool hasPointResultFile() const;
 
+  //! \brief Serialization support.
+  virtual bool serialize(std::map<std::string,std::string>&) { return false; }
+
 protected:
   //! \brief Preprocesses the result sampling points.
   virtual void preprocessResultPoints();

src/SIM/TimeStep.C

@@ -16,6 +16,10 @@
 #include "IFEM.h"
 #include "tinyxml.h"
 #include <sstream>
+#ifdef HAS_CEREAL
+#include <cereal/cereal.hpp>
+#include <cereal/archives/binary.hpp>
+#endif
 
 
 TimeStep::TimeStep () : step(0), iter(time.it), lstep(0)
@@ -310,3 +314,52 @@ bool TimeStep::cutback ()
              << time.dt << std::endl;
   return true;
 }
+
+
+#ifdef HAS_CEREAL
+/*!
+  \brief Serialize to/from archive.
+  \param ar Input or output archive
+*/
+
+template<class T> void doSerializeOps(T& ar, TimeStep& tp)
+{
+  ar(tp.step);
+  ar(tp.starTime);
+  ar(tp.maxCFL);
+  ar(tp.time.t);
+  ar(tp.time.dt);
+  ar(tp.time.dtn);
+  ar(tp.time.CFL);
+  ar(tp.time.first);
+}
+#endif
+
+
+bool TimeStep::serialize(std::map<std::string,std::string>& data)
+{
+#ifdef HAS_CEREAL
+  std::ostringstream str;
+  cereal::BinaryOutputArchive ar(str);
+  doSerializeOps(ar,*this);
+  data.insert(std::make_pair("TimeStep", str.str()));
+  return true;
+#endif
+  return false;
+}
+
+
+bool TimeStep::deSerialize(const std::map<std::string,std::string>& data)
+{
+#ifdef HAS_CEREAL
+  std::stringstream str;
+  auto it = data.find("TimeStep");
+  if (it != data.end()) {
+    str << it->second;
+    cereal::BinaryInputArchive ar(str);
+    doSerializeOps(ar,*this);
+    return true;
+  }
+#endif
+  return false;
+}

src/SIM/TimeStep.h

@@ -18,6 +18,8 @@
 #include <vector>
 #include <cstddef>
 #include <iostream>
+#include <map>
+#include <string>
 
 class TiXmlElement;
 
@@ -59,6 +61,13 @@ public:
   //! \return \e false Cannot do further cut-back, time step size too small
   bool cutback();
 
+  //! \brief Serialize internal state for restarting purposes.
+  //! \param data Container for serialized data
+  bool serialize(std::map<std::string,std::string>& data);
+  //! \brief Set internal state from a serialized state.
+  //! \param[in] data Container for serialized data
+  bool deSerialize(const std::map<std::string,std::string>& data);
+
   int        step; //!< Time step counter
   int&       iter; //!< Iteration counter
   TimeDomain time; //!< Time domain data

src/Utility/DataExporter.C

@@ -95,13 +95,24 @@ bool DataExporter::setFieldValue (const std::string& name,
 }
 
 
-bool DataExporter::dumpTimeLevel (const TimeStep* tp, bool geometryUpdated)
+bool DataExporter::dumpForRestart (const TimeStep* tp) const
+{
+  return m_nrestart > 0 && tp && tp->step % m_nrestart == 0;
+}
+
+
+bool DataExporter::dumpTimeLevel (const TimeStep* tp, bool geometryUpdated,
+                                  SerializeData* serializeData)
 {
   // ignore multiple calls for the same time step
   if (tp && tp->step == m_last_step)
     return true;
 
-  if (tp && tp->step % m_ndump && tp->step % m_ndump > m_order)
+  bool writeData = !tp || tp->step % m_ndump == 0;
+  bool writeRestart = serializeData && this->dumpForRestart(tp);
+  int restartLevel = writeRestart ? tp->step / m_nrestart : 0;
+
+  if (!writeRestart && !writeData)
     return true;
 
   if (tp)
@@ -123,7 +134,8 @@ bool DataExporter::dumpTimeLevel (const TimeStep* tp, bool geometryUpdated)
           (*it2)->writeVector(m_level,*it);
           break;
         case SIM:
-          (*it2)->writeSIM(m_level,*it,geometryUpdated,it->second.prefix);
+          if (writeData)
+            (*it2)->writeSIM(m_level,*it,geometryUpdated,it->second.prefix);
           break;
         case NODALFORCES:
           (*it2)->writeNodalForces(m_level,*it);
@@ -141,7 +153,9 @@ bool DataExporter::dumpTimeLevel (const TimeStep* tp, bool geometryUpdated)
       }
     }
     if (tp)
-      (*it2)->writeTimeInfo(m_level,m_order,m_ndump,*tp);
+      (*it2)->writeTimeInfo(m_level,m_ndump,*tp);
+    if (writeRestart)
+      (*it2)->writeRestartData(restartLevel, *serializeData);
 
     (*it2)->closeFile(m_level);
   }
@@ -156,57 +170,6 @@ bool DataExporter::dumpTimeLevel (const TimeStep* tp, bool geometryUpdated)
 }
 
 
-bool DataExporter::loadTimeLevel (int level, DataWriter* info,
-                                  DataWriter* input)
-{
-  if (!input)
-    if (m_writers.empty())
-      return false;
-    else if (m_dataReader)
-      input = m_dataReader;
-    else
-      input = m_writers.front();
-
-  if (!info)
-    if (m_infoReader)
-      info = m_infoReader;
-    else if (m_writers.size() < 2)
-      return false;
-    else
-      info = m_writers[1];
-
-  int level2=level;
-  if (level == -1)
-    if ((m_level = info->getLastTimeLevel()) < 0)
-      return false;
-    else
-      level2 = m_level;
-
-  bool ok = true;
-  input->openFile(level2);
-  std::map<std::string,FileEntry>::iterator it;
-  for (it = m_entry.begin(); it != m_entry.end() && ok; ++it) {
-    if (!it->second.data)
-      ok = false;
-    else switch (it->second.field)
-    {
-      case SIM:
-        ok = input->readSIM(level2,*it);
-        break;
-      default:
-        break;
-    }
-  }
-  input->closeFile(level2);
-  // if we load the last time level, we want to advance
-  // if we load a specified time level, we do not want to advance
-  if (level == -1)
-    m_level++;
-
-  return ok;
-}
-
-
 int DataExporter::getTimeLevel ()
 {
   if (m_level == -1)
@@ -236,13 +199,7 @@ void DataExporter::setNormPrefixes(const char** prefix)
 
 int DataExporter::realTimeLevel(int filelevel) const
 {
-  return realTimeLevel(filelevel,m_order,m_ndump);
+  return filelevel*m_ndump;
-}
-
-
-int DataExporter::realTimeLevel(int filelevel, int order, int interval) const
-{
-  return filelevel/order*interval;
 }
 
 

src/Utility/DataExporter.h

@@ -15,8 +15,9 @@
 #define _DATA_EXPORTER_H
 
 #include "ControlFIFO.h"
-#include <vector>
 #include <map>
+#include <string>
+#include <vector>
 
 class DataWriter;
 class ProcessAdm;
@@ -33,6 +34,8 @@ class TimeStep;
 class DataExporter : public ControlCallback
 {
 public:
+  typedef std::map<std::string, std::string> SerializeData; //!< Convenience typedef
+
    //! \brief Supported field types
   enum FieldType {
     VECTOR,
@@ -51,7 +54,6 @@ public:
     NORMS        = 8,   //!< Storage of norms
     EIGENMODES   = 16,  //!< Storage of eigenmodes
     ONCE         = 32,  //!< Only write field once
-    RESTART      = 64,  //!< Write restart info
     GRID         = 128, //!< Always store an updated grid
     REDUNDANT    = 256  //!< Field is redundantly calculated on all processes
   };
@@ -74,11 +76,10 @@ public:
   //! \brief Default constructor.
   //! \param[in] dynWriters If \e true, delete the writers on destruction
   //! \param[in] ndump Interval between dumps
-  //! \param[in] order The temporal order of simulations
+  //! \param[in] nrestart Restart stride. 0 to disable
-  //! (always dumps order solutions in a row)
+  DataExporter(bool dynWriters = false, int ndump=1, int nrestart=0) :
-  DataExporter(bool dynWriters = false, int ndump=1, int order=1) :
+    m_delete(dynWriters), m_level(-1), m_ndump(ndump),
-    m_delete(dynWriters), m_level(-1), m_ndump(ndump), m_order(order),
+    m_last_step(-1), m_nrestart(nrestart), m_infoReader(0), m_dataReader(0) {}
-    m_last_step(-1), m_infoReader(0), m_dataReader(0) {}
 
   //! \brief The destructor deletes the writers if \a dynWriters was \e true.
   virtual ~DataExporter();
@@ -112,7 +113,9 @@ public:
   //! \brief Dumps all registered fields using the registered writers.
   //! \param[in] tp Current time stepping info
   //! \param[in] geometryUpdated Whether or not geometries are updated
-  bool dumpTimeLevel(const TimeStep* tp=nullptr, bool geometryUpdated=false);
+  //! \param[in] serializeData Serialized data from simulators for restart files
+  bool dumpTimeLevel(const TimeStep* tp=nullptr, bool geometryUpdated=false,
+                     SerializeData* serializeData = nullptr);
 
   //! \brief Loads last time level with first registered writer by default.
   //! \param[in] level Time level to load, defaults to last time level
@@ -124,10 +127,8 @@ public:
   //! \brief Returns the current time level of the exporter.
   int getTimeLevel();
 
-  //! \brief Calculates the real time level taking order and ndump into account.
+  //! \brief Calculates the real time level taking ndump into account.
   int realTimeLevel(int filelevel) const;
-  //! \brief Calculates the real time level taking order and ndump into account.
-  int realTimeLevel(int filelevel, int order, int interval) const;
 
   //! \brief Sets the prefices used for norm output.
   void setNormPrefixes(const char** prefix);
@@ -140,8 +141,11 @@ public:
   //! \brief Return name from data writer
   std::string getName() const;
 
+  //! \brief Returns visualization data stride
   int getStride() const { return m_ndump; }
-  int getOrder() const { return m_order; }
+
+  //! \brief Returns whether current step should be saved for restart or not.
+  bool dumpForRestart(const TimeStep* tp) const;
 
 protected:
   //! \brief Internal helper function.
@@ -155,8 +159,8 @@ protected:
   bool m_delete;    //!< If true, we are in charge of freeing up datawriters
   int  m_level;     //!< Current time level
   int  m_ndump;     //!< Time level stride for dumping
-  int  m_order;     //!< The temporal order used (needed to facilitate restart)
+  int  m_last_step; //!< Last time step we dumped for
-  int  m_last_step; //!< Last time step we dumped for.
+  int  m_nrestart;  //!< Stride for restart data dumping
 
   DataWriter* m_infoReader; //!< DataWriter to read data information from
   DataWriter* m_dataReader; //!< DataWriter to read numerical data from
@@ -234,17 +238,11 @@ public:
   virtual void writeBasis(int level, const DataEntry& entry,
                           const std::string& prefix) = 0;
 
-  //! \brief Reads data from a file into s SIM object.
-  //! \param[in] level The time level to read the data at
-  //! \param[in] entry The DataEntry describing the SIM
-  virtual bool readSIM(int level, const DataEntry& entry) = 0;
-
   //! \brief Writes time stepping info to file.
   //! \param[in] level The time level to write the info at
-  //! \param[in] order The temporal order
   //! \param[in] interval The number of time steps between each data dump
   //! \param[in] tp The current time stepping info
-  virtual bool writeTimeInfo(int level, int order, int interval,
+  virtual bool writeTimeInfo(int level, int interval,
                              const TimeStep& tp) = 0;
 
   //! \brief Sets the prefices used for norm output.
@@ -253,6 +251,11 @@ public:
   //! \brief Returns the name of the file
   const std::string& getName() const { return m_name; }
 
+  //! \brief Write restart data.
+  //! \param[in] step Level to write restart data at
+  //! \param[in] data Data to write
+  virtual bool writeRestartData(int level, const DataExporter::SerializeData& data) = 0;
+
 protected:
   std::string  m_name;   //!< File name
   const char** m_prefix; //!< The norm prefixes

src/Utility/HDF5Writer.C

@@ -38,7 +38,8 @@
 
 HDF5Writer::HDF5Writer (const std::string& name, const ProcessAdm& adm,
                         bool append, bool keepOpen)
-  : DataWriter(name,adm,".hdf5"), m_file(0), m_keepOpen(keepOpen)
+  : DataWriter(name,adm,".hdf5"), m_file(0), m_restart_file(0),
+    m_keepOpen(keepOpen)
 #ifdef HAVE_MPI
   , m_adm(adm)
 #endif
@@ -52,6 +53,11 @@ HDF5Writer::HDF5Writer (const std::string& name, const ProcessAdm& adm,
     m_flag = H5F_ACC_RDWR;
   else
     m_flag = H5F_ACC_TRUNC;
+  m_restart_name = name + "_restart.hdf5";
+  if (stat(m_restart_name.c_str(),&temp) == 0)
+    m_restart_flag = H5F_ACC_RDWR;
+  else
+    m_restart_flag = H5F_ACC_TRUNC;
 #endif
 }
 
@@ -89,10 +95,11 @@ int HDF5Writer::getLastTimeLevel ()
 }
 
 
-void HDF5Writer::openFile(int level)
+void HDF5Writer::openFile(int level, bool restart)
 {
-  if (m_file)
+  if ((m_file && !restart) || (m_restart_file && restart))
     return;
+  int file = 0;
 #ifdef HAS_HDF5
   hid_t acc_tpl = H5P_DEFAULT;
 #ifdef HAVE_MPI
@@ -100,13 +107,15 @@ void HDF5Writer::openFile(int level)
   acc_tpl = H5Pcreate(H5P_FILE_ACCESS);
   H5Pset_fapl_mpio(acc_tpl, *m_adm.getCommunicator(), info);
 #endif
+  unsigned int flag = restart ? m_restart_flag : m_flag;
+  std::string fname = restart ? m_restart_name : m_name;
 
-  if (m_flag == H5F_ACC_TRUNC)
+  if (flag == H5F_ACC_TRUNC)
-    m_file = H5Fcreate(m_name.c_str(),m_flag,H5P_DEFAULT,acc_tpl);
+    file = H5Fcreate(fname.c_str(),m_flag,H5P_DEFAULT,acc_tpl);
   else {
     // check free disk space - to protect against corrupting files
     // due to out of space condition
-    if (m_flag == H5F_ACC_RDWR) {
+    if (flag == H5F_ACC_RDWR) {
 #ifdef HAVE_GET_CURRENT_DIR_NAME
       char* cwd = get_current_dir_name();
       struct statvfs vfs;
@@ -118,9 +127,9 @@ void HDF5Writer::openFile(int level)
       free(cwd);
 #endif
     }
-    m_file = H5Fopen(m_name.c_str(),m_flag,acc_tpl);
+    file = H5Fopen(fname.c_str(),flag,acc_tpl);
   }
-  if (m_file <= 0)
+  if (file <= 0)
   {
     std::cerr <<" *** HDF5Writer: Failed to open "<< m_name << std::endl;
     return;
@@ -128,12 +137,16 @@ void HDF5Writer::openFile(int level)
 
   std::stringstream str;
   str << '/' << level;
-  if (!checkGroupExistence(m_file,str.str().c_str()))
+  if (!checkGroupExistence(file,str.str().c_str()))
-    H5Gclose(H5Gcreate2(m_file,str.str().c_str(),0,H5P_DEFAULT,H5P_DEFAULT));
+    H5Gclose(H5Gcreate2(file,str.str().c_str(),0,H5P_DEFAULT,H5P_DEFAULT));
 #ifdef HAVE_MPI
   H5Pclose(acc_tpl);
 #endif
 #endif
+  if (restart)
+    m_restart_file = file;
+  else
+   m_file = file;
 }
 
 
@@ -190,6 +203,23 @@ void HDF5Writer::readArray(int group, const std::string& name,
 }
 
 
+void HDF5Writer::readArray(int group, const std::string& name,
+                           int& len,  char*& data)
+{
+#ifdef HAS_HDF5
+  hid_t set = H5Dopen2(group,name.c_str(),H5P_DEFAULT);
+  hsize_t siz = H5Dget_storage_size(set);
+  len = siz;
+  data = new char[siz];
+  H5Dread(set,H5T_NATIVE_CHAR,H5S_ALL,H5S_ALL,H5P_DEFAULT,data);
+  H5Dclose(set);
+#else
+  len = 0;
+  std::cout << "HDF5Writer: compiled without HDF5 support, no data read" << std::endl;
+#endif
+}
+
+
 void HDF5Writer::readString(const std::string& name, std::string& out, bool close)
 {
 #ifdef HAS_HDF5
@@ -289,49 +319,6 @@ void HDF5Writer::writeVector(int level, const DataEntry& entry)
 }
 
 
-bool HDF5Writer::readSIM (int level, const DataEntry& entry)
-{
-  SIMbase* sim = static_cast<SIMbase*>(const_cast<void*>(entry.second.data));
-  Vector* sol = static_cast<Vector*>(const_cast<void*>(entry.second.data2));
-  if (!sim || !sol) return false;
-  if (!(abs(entry.second.results) & DataExporter::RESTART))
-    return true;
-
-  bool ok = true;
-#ifdef HAS_HDF5
-  std::string name;
-  if (sim->mixedProblem())
-    name = entry.first;
-  else
-    name = entry.second.description + " restart";
-
-  for (int i = 0; i < sim->getNoPatches() && ok; ++i) {
-    std::stringstream str;
-    str << '/' << level << '/' << i+1;
-    hid_t group2 = H5Gopen2(m_file,str.str().c_str(),H5P_DEFAULT);
-    int loc = sim->getLocalPatchIndex(i+1);
-    if (loc > 0) {
-      double* tmp = nullptr; int siz = 0;
-      readArray(group2,name,siz,tmp);
-      ok = sim->injectPatchSolution(*sol,Vector(tmp,siz),loc-1);
-      if (hasGeometries(level, sim->getName()+"-1")) {
-        std::string out;
-        std::stringstream geom;
-        geom << '/' << level << "/basis/" << sim->getName() << "-1" << "/" << i+1;
-        readString(geom.str(), out, false);
-        std::stringstream str;
-        str << out;
-        sim->getPatch(loc)->read(str);
-      }
-      delete[] tmp;
-    }
-    H5Gclose(group2);
-  }
-#endif
-  return ok;
-}
-
-
 bool HDF5Writer::readVector(int level, const std::string& name,
                             int patch, std::vector<double>& vec)
 {
@@ -435,13 +422,6 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
     if (loc > 0 && (!(abs(results) & DataExporter::REDUNDANT) ||
                     sim->getGlobalProcessID() == 0)) // we own the patch
     {
-      if (abs(results) & DataExporter::RESTART) {
-        Vector psol;
-        int ncmps = entry.second.ncmps;
-        sim->extractPatchSolution(*sol,psol,loc-1,ncmps);
-        writeArray(group2, entry.second.description+" restart",
-                           psol.size(), psol.ptr(), H5T_NATIVE_DOUBLE);
-      }
       if (abs(results) & DataExporter::PRIMARY) {
         Vector psol;
         int ncmps = entry.second.ncmps;
@@ -534,10 +514,6 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
     else // must write empty dummy records for the other patches
     {
       double dummy;
-      if (abs(results) & DataExporter::RESTART) {
-        writeArray(group2, entry.second.description+" restart",
-                   0, &dummy, H5T_NATIVE_DOUBLE);
-      }
       if (abs(results) & DataExporter::PRIMARY) {
         if (entry.second.results < 0) {
           writeArray(group2, entry.second.description,
@@ -682,8 +658,7 @@ bool HDF5Writer::checkGroupExistence(int parent, const char* path)
 
 // TODO: implement for variable time steps.
 // (named time series to allow different timelevels for different fields)
-bool HDF5Writer::writeTimeInfo (int level, int order, int interval,
+bool HDF5Writer::writeTimeInfo (int level, int interval, const TimeStep& tp)
-                                const TimeStep& tp)
 {
 #ifdef HAS_HDF5
   std::stringstream str;
@@ -797,3 +772,105 @@ bool HDF5Writer::readVector(int level, const std::string& name,
   closeFile(level);
   return ok;
 }
+
+
+bool HDF5Writer::writeRestartData(int level, const DataExporter::SerializeData& data)
+{
+#ifdef HAS_HDF5
+  if (level > 0)
+    m_restart_flag = H5F_ACC_RDWR;
+  openFile(level, true);
+  int pid = 0;
+  int ptot = 1;
+#ifdef HAVE_MPI
+  pid = m_adm.getProcId();
+  ptot = m_adm.getNoProcs();
+#endif
+  for (int p = 0; p < ptot; ++p) {
+    for (auto& it : data) {
+      std::stringstream str;
+      str << level << '/' << p;
+
+      int group;
+      if (checkGroupExistence(m_restart_file,str.str().c_str()))
+        group = H5Gopen2(m_restart_file,str.str().c_str(),H5P_DEFAULT);
+      else
+        group = H5Gcreate2(m_restart_file,str.str().c_str(),0,H5P_DEFAULT,H5P_DEFAULT);
+
+      if (!H5Lexists(group, it.first.c_str(), 0)) {
+        if (pid == p)
+          writeArray(group, it.first, it.second.size(), it.second.data(), H5T_NATIVE_CHAR);
+        else
+          writeArray(group, it.first, 0, nullptr, H5T_NATIVE_CHAR);
+      }
+
+      H5Gclose(group);
+    }
+  }
+  H5Fclose(m_restart_file);
+  m_restart_file = 0;
+#else
+  std::cout << "HDF5Writer: compiled without HDF5 support, no data written" << std::endl;
+#endif
+
+  return true;
+}
+
+
+struct read_restart_ctx {
+  HDF5Writer* w;
+  DataExporter::SerializeData* data;
+};
+
+
+#ifdef HAS_HDF5
+static herr_t read_restart_data(hid_t group_id, const char* member_name, void* data)
+{
+  read_restart_ctx* ctx = static_cast<read_restart_ctx*>(data);
+
+  char* c;
+  int len;
+  ctx->w->readArray(group_id, member_name, len, c);
+  std::string tmp;
+  tmp.resize(len);
+  tmp.assign(c, len);
+  ctx->data->insert(std::make_pair(std::string(member_name),tmp));
+  return 0;
+}
+#endif
+
+
+int HDF5Writer::readRestartData(DataExporter::SerializeData& data, int level)
+{
+#ifdef HAS_HDF5
+  openFile(0);
+
+  if (level == -1) {
+    while (true) {
+      std::stringstream str;
+      str << '/' << level+1;
+      if (!checkGroupExistence(m_file, str.str().c_str()))
+        break;
+      ++level;
+    }
+  }
+
+  std::stringstream str;
+  str << '/' << level << '/'
+#ifdef HAVE_MPI
+   << m_adm.getProcId();
+#else
+  << 0;
+#endif
+  int idx = 0;
+  read_restart_ctx ctx;
+  ctx.w = this;
+  ctx.data = &data;
+  int it = H5Giterate(m_file, str.str().c_str(), &idx, read_restart_data, &ctx);
+  closeFile(0);
+  return it < 0 ? it : level;
+#else
+  std::cout << "HDF5Writer: compiled without HDF5 support, no data read" << std::endl;
+  return -1;
+#endif
+}

src/Utility/HDF5Writer.h

@@ -45,7 +45,12 @@ public:
 
   //! \brief Opens the file at a given time level.
   //! \param[in] level The requested time level
-  virtual void openFile(int level);
+  virtual void openFile(int level) { openFile(level, false); }
+
+  //! \brief Opens the file at a given time level.
+  //! \param[in] level The requested time level
+  //! \param[in] restart If true, open restart file
+  void openFile(int level, bool restart);
 
   //! \brief Closes the file.
   //! \param[in] level Level we just wrote to the file
@@ -74,11 +79,6 @@ public:
   virtual void writeSIM(int level, const DataEntry& entry,
                         bool geometryUpdated, const std::string& prefix);
 
-  //! \brief Reads data from a file into a SIM.
-  //! \param[in] level The time level to read the data at
-  //! \param[in] entry The DataEntry describing the SIM
-  virtual bool readSIM(int level, const DataEntry& entry);
-
   //! \brief Writes nodal forces to file.
   //! \param[in] level The time level to write the data at
   //! \param[in] entry The DataEntry describing the vector
@@ -100,11 +100,9 @@ public:
 
   //! \brief Writes time stepping info to file.
   //! \param[in] level The time level to write the info at
-  //! \param[in] order The temporal order
   //! \param[in] interval The number of time steps between each data dump
   //! \param[in] tp The current time stepping info
-  virtual bool writeTimeInfo(int level, int order, int interval,
+  virtual bool writeTimeInfo(int level, int interval, const TimeStep& tp);
-                             const TimeStep& tp);
 
   //! \brief Reads a text string.
   //! \param[in] name The name (path in HDF5 file) to the string
@@ -141,6 +139,24 @@ public:
   //! \param[in] basisName Check for a particular basis
   bool hasGeometries(int level, const std::string& basisName = "");
 
+  //! \brief Write restart data.
+  //! \param level Level to write data at
+  //! \param data Data to write
+  bool writeRestartData(int level, const DataExporter::SerializeData& data);
+
+  //! \brief Read restart data from file.
+  //! \param data The map to store data in
+  //! \param level Level to read (-1 to read last level in file)
+  //! \returns Negative value on error, else restart level loaded
+  int readRestartData(DataExporter::SerializeData& data, int level = -1);
+
+  //! \brief Internal helper function. Reads an array into an array of chars.
+  //! \param[in] group The HDF5 group to read data from
+  //! \param[in] name The name of the array
+  //! \param[in] len The length of the data to read
+  //! \param[out] data The array to read data into
+  void readArray(int group, const std::string& name, int& len, char*& data);
+
 protected:
   //! \brief Internal helper function. Writes a data array to HDF5 file.
   //! \param[in] group The HDF5 group to write data into
@@ -182,8 +198,11 @@ protected:
 
 private:
   int          m_file; //!< The HDF5 handle for our file
+  int  m_restart_file; //!< The HDF5 handle for our restart file
   unsigned int m_flag; //!< The file flags to open HDF5 file with
+  unsigned int m_restart_flag; //!< The file flags to open the restart file with
   bool     m_keepOpen; //!< If \e true, we always keep the file open
+  std::string m_restart_name; //!< The restart file to use
 #ifdef HAVE_MPI
   const ProcessAdm& m_adm;   //!< Pointer to process adm in use
 #endif

src/Utility/XMLWriter.C

@@ -28,7 +28,7 @@ XMLWriter::XMLWriter (const std::string& name, const ProcessAdm& adm) :
   m_doc = nullptr;
   m_node = nullptr;
   m_dt = 0;
-  m_order = m_interval = 1;
+  m_interval = 1;
 }
 
 
@@ -73,7 +73,6 @@ void XMLWriter::closeFile(int level, bool force)
   TiXmlElement element3("timestep");
   sprintf(temp,"%f",m_dt);
   element3.SetAttribute("constant","1");
-  element3.SetAttribute("order",m_order);
   element3.SetAttribute("interval",m_interval);
   pNewNode = m_node->InsertEndChild(element3);
   TiXmlText value2(temp);
@@ -185,12 +184,6 @@ void XMLWriter::writeKnotspan(int level, const DataEntry& entry,
 }
 
 
-bool XMLWriter::readSIM (int level, const DataEntry& entry)
-{
-  return true;
-}
-
-
 void XMLWriter::writeSIM (int level, const DataEntry& entry, bool,
                           const std::string& prefix)
 {
@@ -215,22 +208,12 @@ void XMLWriter::writeSIM (int level, const DataEntry& entry, bool,
   else
     basisname = prefix+sim->getName()+"-1";
 
-  // restart vector
-  if (results & DataExporter::RESTART) {
-    addField(prefix+"restart",entry.second.description,basisname,
-             cmps,sim->getNoPatches(),"restart");
-  }
-
   if (results & DataExporter::PRIMARY) {
     if (entry.second.results < 0)
       addField(entry.second.description, entry.second.description,
                basisname, cmps, sim->getNoPatches(), "field");
     else if (sim->mixedProblem())
     {
-      // primary solution vector
-      addField(prefix+entry.first,entry.second.description,basisname,
-               0,sim->getNoPatches(),"restart");
-
       // primary solution fields
       for (int b=1; b <= sim->getNoBasis(); ++b) {
         std::stringstream str;
@@ -297,11 +280,9 @@ void XMLWriter::writeSIM (int level, const DataEntry& entry, bool,
 }
 
 
-bool XMLWriter::writeTimeInfo (int level, int order, int interval,
+bool XMLWriter::writeTimeInfo (int level, int interval, const TimeStep& tp)
-                               const TimeStep& tp)
 {
   m_dt = tp.time.dt;
-  m_order = order;
   m_interval = interval;
   return true;
 }

src/Utility/XMLWriter.h

@@ -40,7 +40,6 @@ public:
     int    patches;    //!< Number of patches in field
     int    components; //!< Number of components in field
     double timestep;   //!< The time step associated with the field
-    int    order;      //!< The temporal order associated with the field
     int    interval;   //!< The dumping interval for the field
     bool   once;       //!< If true, field is only stored at first time level
   };
@@ -86,11 +85,6 @@ public:
   virtual void writeSIM(int level, const DataEntry& entry,
                         bool, const std::string& prefix);
 
-  //! \brief Reads data from a file into SIM.
-  //! \param[in] level The time level to read the data at
-  //! \param[in] entry The DataEntry describing the SIM
-  virtual bool readSIM(int level, const DataEntry& entry);
-
   //! \brief Writes nodal forces to file.
   //! \param[in] level The time level to write the data at
   //! \param[in] entry The DataEntry describing the vector
@@ -105,11 +99,9 @@ public:
 
   //! \brief Writes time stepping info to file.
   //! \param[in] level The time level to write the info at
-  //! \param[in] order The temporal order
   //! \param[in] interval The number of time steps between each data dump
   //! \param[in] tp The current time stepping info
-  virtual bool writeTimeInfo(int level, int order, int interval,
+  virtual bool writeTimeInfo(int level, int interval, const TimeStep& tp);
-                             const TimeStep& tp);
 
   //! \brief Write a basis to file
   //! \param[in] level The time level to write the basis at
@@ -117,6 +109,13 @@ public:
   //! \param[in] prefix Prefix for basis
   virtual void writeBasis(int level, const DataEntry& entry,
                           const std::string& prefix) {}
+
+  //! \brief No XML for restart data
+  //! \param level Level to write data at
+  //! \param data Data to write
+  virtual bool writeRestartData(int level,
+                                const DataExporter::SerializeData& data)
+  { return true; }
 protected:
   //! \brief Internal helper function adding an XML-element to the file
   //! \param[in] name The name of the field to add
@@ -137,7 +136,6 @@ private:
   TiXmlNode*     m_node; //!< The document's root node
 
   double m_dt;       //!< The current (constant) time step size
-  int    m_order;    //!< The temporal order
   int    m_interval; //!< Stride for dumping (dump at every m_interval'th level)
 };