diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index debb516bb7b190057c222e4488f448070e3219e9..16f0fe5ce202b0d240eb6bce35ee8f6dcdef580d 100644
--- a/src/CMakeLists.txt
+++ b/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}")
diff --git a/src/program_state.hh b/src/program_state.hh
index 92aeb766048326a2e5f0632b25207cf4d650dd1d..6ce0eba0d6943af81f2bbcc50ed27441d2474c37 100644
--- a/src/program_state.hh
+++ b/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
diff --git a/src/sand-wedge-data/parset.cfg b/src/sand-wedge-data/parset.cfg
index 52a760fea86955be74ea78e2aa18bbe1354c8ad6..3b38ac859fc7fd63006f3aab433a90871d9b18df 100644
--- a/src/sand-wedge-data/parset.cfg
+++ b/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
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 675110ac0258b096e635dc7c883aeb10f942df85..b3107785ccbd792c9eecfc8b4a47ed2bb354ea81 100644
--- a/src/sand-wedge.cc
+++ b/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(),