[Output] Overhaul writer logic, rename writeVTK

- The configuration parameter io.writeVTK is now io.vtk.write
- The configuration parameters restarts.first and restarts.spacing
  were moved into the io namespace
- Restarts are stored in a separate file
- VTK data, restarts (in HDF5), and the other HDF5 data are now
  each controlled through a parameter:
    io.data.write, io.restarts.write, io.vtk.write

In particular, it is thus now possible to write VTK data while
opening the data and restart files in read-only mode (or not at all);
This makes it possible to, by keeping the restarts file (and the binary!)
around, quickly (re-)compute a a small segment from the history.

The likeliest scenario: VTK files are huge and for a computation with
tens of thousands of timesteps (e.g. 30000), we know already (from the
data, or just expect it) that an event occurred during the last 500
timesteps. So we recompute them from a restart in the following mode:

  -io.restarts.first 29500
  -io.restarts.write false
  -io.data.write false
  -io.vtk.write true

Note that care needs to be taken to use a number for io.restarts.first
which is compatible with the old io.restarts.spacing

src/hdf5-writer.hh

@@ -8,7 +8,6 @@
 #include "hdf5/frictionalboundary-writer.hh"
 #include "hdf5/iteration-writer.hh"
 #include "hdf5/patchinfo-writer.hh"
-#include "hdf5/restart-io.hh"
 #include "hdf5/surface-writer.hh"
 #include "hdf5/time-writer.hh"
 
@@ -32,8 +31,8 @@ class HDF5Writer {
         patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
 #endif
         surfaceWriter_(file_, vertexCoordinates, surface),
-        frictionalBoundaryWriter_(file_, vertexCoordinates, frictionalBoundary),
-        restartWriter_(file_, vertexCoordinates.size()) {
+        frictionalBoundaryWriter_(file_, vertexCoordinates,
+                                  frictionalBoundary) {
   }
 
   template <class Friction>
@@ -53,10 +52,6 @@ class HDF5Writer {
     iterationWriter_.write(programState.timeStep, iterationCount);
   }
 
-  void writeRestart(ProgramState const &programState) {
-    restartWriter_.write(programState);
-  }
-
 private:
   HDF5::Grouplike &file_;
 
@@ -67,7 +62,5 @@ class HDF5Writer {
 #endif
   SurfaceWriter<ProgramState, GridView> surfaceWriter_;
   FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
-
-  RestartIO<ProgramState> restartWriter_;
 };
 #endif

src/hdf5/restart-io.cc

@@ -5,18 +5,17 @@
 #include "restart-io.hh"
 
 template <class ProgramState>
-RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &file, size_t vertexCount)
-    : group_(file, "restarts"),
-      displacementWriter_(group_, "displacement", vertexCount,
+RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
+    : displacementWriter_(group, "displacement", vertexCount,
                           Vector::block_type::dimension),
-      velocityWriter_(group_, "velocity", vertexCount,
+      velocityWriter_(group, "velocity", vertexCount,
                       Vector::block_type::dimension),
-      accelerationWriter_(group_, "acceleration", vertexCount,
+      accelerationWriter_(group, "acceleration", vertexCount,
                           Vector::block_type::dimension),
-      stateWriter_(group_, "state", vertexCount),
-      weightedNormalStressWriter_(group_, "weightedNormalStress", vertexCount),
-      relativeTimeWriter_(group_, "relativeTime"),
-      relativeTimeIncrementWriter_(group_, "relativeTimeIncrement") {}
+      stateWriter_(group, "state", vertexCount),
+      weightedNormalStressWriter_(group, "weightedNormalStress", vertexCount),
+      relativeTimeWriter_(group, "relativeTime"),
+      relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {}
 
 template <class ProgramState>
 void RestartIO<ProgramState>::write(ProgramState const &programState) {

src/hdf5/restart-io.hh

@@ -9,15 +9,13 @@ template <class ProgramState> class RestartIO {
   using Vector = typename ProgramState::Vector;
 
 public:
-  RestartIO(HDF5::Grouplike &file, size_t vertexCount);
+  RestartIO(HDF5::Grouplike &group, size_t vertexCount);
 
   void write(ProgramState const &programState);
 
   void read(size_t timeStep, ProgramState &programState);
 
 private:
-  HDF5::Group group_;
-
   HDF5::SequenceIO<2> displacementWriter_;
   HDF5::SequenceIO<2> velocityWriter_;
   HDF5::SequenceIO<2> accelerationWriter_;

src/one-body-problem.cc

@@ -176,15 +176,29 @@ int main(int argc, char *argv[]) {
           _ell += gravityFunctional;
         };
 
-    ProgramState<Vector, ScalarVector> programState(leafVertexCount);
-    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");
-
-    HDF5::File ioFile("output.h5");
-    if (firstRestart != 0)
-      RestartIO<ProgramState<Vector, ScalarVector>>(ioFile, leafVertexCount)
-          .read(firstRestart, programState);
+    using MyProgramState = ProgramState<Vector, ScalarVector>;
+    MyProgramState programState(leafVertexCount);
+    auto const firstRestart = parset.get<size_t>("io.restarts.first");
+    auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
+    auto const writeRestarts = parset.get<bool>("io.restarts.write");
+    auto const writeData = parset.get<bool>("io.data.write");
+    bool const handleRestarts = writeRestarts or firstRestart > 0;
+
+    auto dataFile =
+        writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+    auto restartFile = handleRestarts
+                           ? std::make_unique<HDF5::File>(
+                                 "restarts.h5",
+                                 writeRestarts ? HDF5::Access::READWRITE
+                                               : HDF5::Access::READONLY)
+                           : nullptr;
+    auto restartIO = handleRestarts
+                         ? std::make_unique<RestartIO<MyProgramState>>(
+                               *restartFile, leafVertexCount)
+                         : nullptr;
+
+    if (firstRestart > 0) // automatically adjusts the time and timestep
+      restartIO->read(firstRestart, programState);
     else
       programState.setupInitialConditions(parset, computeExternalForces,
                                           matrices, myAssembler, dirichletNodes,
@@ -205,24 +219,30 @@ int main(int argc, char *argv[]) {
         vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
     }
 
-    HDF5Writer<ProgramState<Vector, ScalarVector>,
-               typename MyAssembler::VertexBasis, GridView>
-        hdf5Writer(ioFile, vertexCoordinates, myAssembler.vertexBasis, surface,
-                   frictionalBoundary, weakPatch);
-    MyVTKWriter<typename MyAssembler::VertexBasis,
-                typename MyAssembler::CellBasis> const
-        vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+    using MyVertexBasis = typename MyAssembler::VertexBasis;
+    auto dataWriter =
+        writeData ? std::make_unique<
+                        HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
+                        *dataFile, vertexCoordinates, myAssembler.vertexBasis,
+                        surface, frictionalBoundary, weakPatch)
+                  : nullptr;
+    MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
+        myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
 
     IterationRegister iterationCount;
     auto const report = [&](bool initial = false) {
-      hdf5Writer.reportSolution(programState, *myGlobalFriction);
-      if (!initial)
-        hdf5Writer.reportIterations(programState, iterationCount);
-
-      if (!initial and programState.timeStep % restartSpacing == 0)
-        hdf5Writer.writeRestart(programState);
+      if (writeData) {
+        dataWriter->reportSolution(programState, *myGlobalFriction);
+        if (!initial)
+          dataWriter->reportIterations(programState, iterationCount);
+        dataFile->flush();
+      }
 
-      ioFile.flush();
+      if (writeRestarts and !initial and
+          programState.timeStep % restartSpacing == 0) {
+        restartIO->write(programState);
+        restartFile->flush();
+      }
 
       if (parset.get<bool>("io.printProgress"))
         std::cout << "timeStep = " << std::setw(6) << programState.timeStep
@@ -230,7 +250,7 @@ int main(int argc, char *argv[]) {
                   << ", tau = " << std::setw(12) << programState.relativeTau
                   << std::endl;
 
-      if (parset.get<bool>("io.writeVTK")) {
+      if (parset.get<bool>("io.vtk.write")) {
         ScalarVector stress;
         myAssembler.assembleVonMisesStress(body.getYoungModulus(),
                                            body.getPoissonRatio(),

src/one-body-problem.cfg

@@ -2,8 +2,12 @@
 gravity         = 9.81  # [m/s^2]
 
 [io]
+data.write      = true
 printProgress   = false
-writeVTK        = false
+restarts.first  = 0
+restarts.spacing= 20
+restarts.write  = true
+vtk.write       = false
 
 [problem]
 finalTime       = 1000  # [s]
@@ -35,11 +39,6 @@ b               = 0.017 # [ ]
 a               = 0.020 # [ ]
 b               = 0.005 # [ ]
 
-[restarts]
-template = restart_%06d
-first = 0
-spacing = 20
-
 [timeSteps]
 scheme = newmark