[Extend] Implement restarts

src/CMakeLists.txt

@@ -20,9 +20,12 @@ file(MAKE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/sand-wedge-data")
 dune_symlink_to_source_files("sand-wedge-data/boundaryconditions.py")
 dune_symlink_to_source_files("sand-wedge-data/parset.cfg")
 
-find_package(Boost REQUIRED system filesystem)
+find_package(Boost REQUIRED system filesystem serialization)
 include_directories(${Boost_INCLUDE_DIR})
 
+# dataio.hh expects this to be set
+set_property(DIRECTORY APPEND PROPERTY COMPILE_DEFINITIONS "HAVE_BOOST_SERIALIZATION")
+
 foreach(_dim 2 3)
   set(_target sand-wedge-${_dim}D)
   add_executable(${_target} ${SOURCE_FILES})
@@ -30,6 +33,7 @@ foreach(_dim 2 3)
   add_dune_ug_flags(${_target})
 
   target_link_libraries(${_target} ${Boost_FILESYSTEM_LIBRARY})
+  target_link_libraries(${_target} ${Boost_SERIALIZATION_LIBRARY})
   target_link_libraries(${_target} ${Boost_SYSTEM_LIBRARY})
 
   set_property(TARGET ${_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")

src/program_state.hh

@@ -4,6 +4,7 @@
 #include <dune/common/parametertree.hh>
 
 #include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/dunedataio.hh>
 
 #include <dune/tectonic/body.hh>
 
@@ -110,4 +111,22 @@ template <class Vector, class ScalarVector> class ProgramState {
   double relativeTau;
   size_t timeStep;
 };
+
+namespace boost {
+namespace serialization {
+  template <class Archive, class Vector, class ScalarVector>
+  void serialize(Archive &ar, ProgramState<Vector, ScalarVector> &ps,
+                 const unsigned int version) {
+    ar &ps.u;
+    ar &ps.v;
+    ar &ps.a;
+    ar &ps.alpha;
+    ar &ps.normalStress;
+    ar &ps.relativeTime;
+    ar &ps.relativeTau;
+    ar &ps.timeStep;
+  }
+}
+}
+
 #endif

src/sand-wedge-data/parset.cfg

@@ -41,6 +41,11 @@ refinementTolerance = 1e-5
 number = 100000
 scheme = newmark
 
+[restarts]
+template = restart_%06d
+first = 0
+spacing = 20
+
 [u0.solver]
 tolerance         = 1e-8
 maximumIterations = 100000

src/sand-wedge.cc

@@ -39,6 +39,7 @@
 #include <dune/fufem/boundarypatch.hh>
 #pragma clang diagnostic pop
 
+#include <dune/fufem/dataio.hh>
 #include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/sharedpointermap.hh>
@@ -210,9 +211,17 @@ int main(int argc, char *argv[]) {
         };
 
     ProgramState<Vector, ScalarVector> programState(leafVertexCount);
-    programState.setupInitialConditions(parset, computeExternalForces, matrices,
-                                        myAssembler, dirichletNodes, noNodes,
-                                        frictionalBoundary, body);
+    auto const firstRestart = parset.get<size_t>("restarts.first");
+    auto const restartSpacing = parset.get<size_t>("restarts.spacing");
+    auto const restartTemplate = parset.get<std::string>("restarts.template");
+    if (firstRestart != 0)
+      DataIO::loadData(programState,
+                       str(boost::format(restartTemplate) % firstRestart));
+    else
+      programState.setupInitialConditions(parset, computeExternalForces,
+                                          matrices, myAssembler, dirichletNodes,
+                                          noNodes, frictionalBoundary, body);
+    assert(programState.timeStep == firstRestart);
 
     MyGlobalFrictionData<LocalVector> frictionInfo(
         parset.sub("boundary.friction"), weakPatch);
@@ -274,6 +283,10 @@ int main(int argc, char *argv[]) {
       specialVelocityWriter.write(velocity);
       specialDisplacementWriter.write(displacement);
 
+      if (programState.timeStep % restartSpacing == 0)
+        DataIO::writeData(programState, str(boost::format(restartTemplate) %
+                                            programState.timeStep));
+
       if (parset.get<bool>("io.writeVTK")) {
         ScalarVector stress;
         myAssembler.assembleVonMisesStress(body.getYoungModulus(),