Add restart functionality

include/base/restart.h

@@ -0,0 +1,153 @@
+#pragma once
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/solution_transfer.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/tria.h>
+
+// for implementing a check-pointing and restart mechanism.
+#include <boost/archive/binary_iarchive.hpp>
+#include <boost/archive/binary_oarchive.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Utilities
+{
+  /**
+   * @brief Create files required for a restart consisting of the triangulation, vectors, and some metadata
+   *
+   * @tparam PartitionerPtr shred_ptr holding a fully distributed  Partitioner with ghost elements
+   *
+   * @param triangulation The triangulation used to create the checkpoint
+   * @param dof_handler DoFHandeler associated to the global vectors
+   * @param partitioner A valid shared_ptr with partitioner used to create tmp vectors if necessary
+   * @param vectors A list of vectors, which are stored in the checkpoint
+   * @param name File name(s) for the restart files
+   * @param t Absolute time associated to the data vectors
+   * @param t Timestep associated to the data vectors
+   */
+  template <int dim, typename VectorType, typename PartitionerPtr>
+  void
+  create_restart_snapshot(
+    const dealii::parallel::distributed::Triangulation<dim> &triangulation,
+    const dealii::DoFHandler<dim>                           &dof_handler,
+    PartitionerPtr                                           partitioner,
+    const std::vector<const VectorType *>                   &vectors,
+    const std::string                                       &name,
+    const double                                             t,
+    unsigned int                                             timestep)
+  {
+    // To load a restart, we need read access to ghost entries.
+    // Thus, we might need temporary vectors to comply with the required
+    // partitioner layout (fully distributed vector with read access to ghost
+    // entries). We use lazy allocation, if required.
+    std::vector<VectorType>         tmp_vectors(vectors.size());
+    std::vector<const VectorType *> tmp_vectors_ptr;
+
+    for (std::size_t i = 0; i < vectors.size(); ++i)
+      {
+        if (!vectors[i]->partitioners_are_globally_compatible(
+              *partitioner.get()))
+          {
+            // in case our vector carries anyway only local data, i.e., no ghost
+            // values at all, we need to create a temporary vector
+            tmp_vectors[i].reinit(partitioner);
+            tmp_vectors[i] = *(vectors[i]);
+            tmp_vectors_ptr.emplace_back(&tmp_vectors[i]);
+          }
+        else
+          {
+            // if the partitioner is compatible, we can just use the existing
+            // vector
+            tmp_vectors_ptr.emplace_back(vectors[i]);
+          }
+      }
+
+    dealii::parallel::distributed::SolutionTransfer<dim, VectorType>
+      solution_transfer(dof_handler);
+
+    solution_transfer.prepare_for_serialization(tmp_vectors_ptr);
+    triangulation.save(name + ".mesh");
+
+    if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+      {
+        std::ofstream file(name + ".metadata", std::ios::binary);
+        boost::archive::binary_oarchive oa(file);
+        oa << t;
+        oa << timestep;
+      }
+  }
+
+
+  /**
+   * @brief Load a restart snapshot previously stored using create_restart_snapshot
+   *
+   * @param dof_handler DoFHandler associated to the vectors
+   * @param vectors Vectors container to load data into
+   * @param name Name of the file bundle
+   * @param t Absolute time stored in the metadata
+   * @param timestep Timestep stored in the metadata
+   */
+  template <int dim, typename VectorType>
+  void
+  load_restart_snapshot(const dealii::DoFHandler<dim> &dof_handler,
+                        std::vector<VectorType *>     &vectors,
+                        const std::string             &name,
+                        double                        &t,
+                        unsigned int                  &timestep)
+  {
+    auto partitioner = std::make_shared<Utilities::MPI::Partitioner>(
+      dof_handler.locally_owned_dofs(),
+      DoFTools::extract_locally_relevant_dofs(dof_handler),
+      MPI_COMM_WORLD);
+
+    // To load a restart, we need write access to ghost entries
+    // Thus, we might need temporary vectors to comply with the required
+    // partitioner layout (fully distributed vector with write access to ghost
+    // entries). We use lazy allocation, if required.
+    std::vector<VectorType>   tmp_vectors(vectors.size());
+    std::vector<VectorType *> tmp_vectors_ptr;
+
+    for (std::size_t i = 0; i < vectors.size(); ++i)
+      {
+        // With LA:d:V, we can grant write access using zero_out_ghost_values
+        // (this wouldn't work with distributed PETSc vectors or similar)
+        if (vectors[i]->has_ghost_elements())
+          {
+            vectors[i]->zero_out_ghost_values();
+            tmp_vectors_ptr.emplace_back(vectors[i]);
+          }
+        else
+          {
+            // in case our vector carries anyway only local data, i.e., no ghost
+            // values at all, we need to create a temporary vector
+            tmp_vectors[i].reinit(partitioner);
+            tmp_vectors_ptr.emplace_back(&tmp_vectors[i]);
+          }
+      }
+
+    // triangulation.load(name + "-checkpoint.mesh");
+    dealii::parallel::distributed::SolutionTransfer<dim, VectorType>
+      solution_transfer(dof_handler);
+    solution_transfer.deserialize(tmp_vectors_ptr);
+
+    for (unsigned int i = 0; i < tmp_vectors_ptr.size(); ++i)
+      {
+        // Copy over data from the temporary vectors, if necessary
+        if (!vectors[i]->has_ghost_elements())
+          vectors[i]->copy_locally_owned_data_from(*(tmp_vectors_ptr[i]));
+        // ... and update all ghost values
+        vectors[i]->update_ghost_values();
+      }
+
+    // Last but not least, retrieve the time-stamp data
+    std::ifstream                   file(name + ".metadata", std::ios::binary);
+    boost::archive::binary_iarchive ia(file);
+    ia >> t;
+    ia >> timestep;
+  }
+} // namespace Utilities
+
+DEAL_II_NAMESPACE_CLOSE

include/base/time_handler.h

@@ -45,6 +45,12 @@ class Time
     time_current += delta_t;
     ++timestep;
   }
+  void
+  set_time(double absolute_time, unsigned int nth_timestep)
+  {
+    time_current = absolute_time;
+    timestep     = nth_timestep;
+  }
 
 private:
   unsigned int timestep;

include/base/utilities.h

@@ -104,7 +104,42 @@ namespace Utilities
       << std::endl;
   }
 
-
+  /**
+   * @brief Helper function to convert different representations of the same number
+   * (in double rpecision) to a common format by truncating trailing zeros.
+   *
+   * Examples:
+   * number = 0.0100 -> 0.01
+   * number = 42.0100 -> 42.01
+   * number = 42.00 -> 42
+   *
+   *
+   * @param number the number to format
+   * @return std::string the string representation
+   */
+  std::string
+  format_time_stamp_to_string(double number)
+  {
+    // We need to apply some cosmetics to t to make it usable
+    std::string str = std::to_string(number);
+    for (int i = str.size() - 1; i >= 1; i--)
+      {
+        if (str.at(i) == '0')
+          {
+            str.pop_back(); // Remove if last digit is '0'.
+          }
+        else if (str.at(i) == '.')
+          {
+            str.pop_back(); // Remove dot.
+            break;          // Break after '.' is removed.
+          }
+        else
+          {
+            break; // Or break before a digit is removed.
+          }
+      }
+    return str;
+  }
 
   /**
    * @brief Rounds a given number up to a defined precision.

include/parameter/parameter_handling.h

@@ -311,8 +311,9 @@ namespace Parameters
   // Set the timestep size $ \varDelta t $ and the simulation end-time.
   struct Time
   {
-    double delta_t  = 0.1;
-    double end_time = 1.;
+    double delta_t    = 0.1;
+    double end_time   = 1.;
+    double start_time = 0;
 
     void
     add_parameters(ParameterHandler &prm);
@@ -323,6 +324,10 @@ namespace Parameters
   {
     prm.enter_subsection("Time");
     {
+      prm.add_parameter("Start time",
+                        start_time,
+                        "Restart simulation or start from the beginning",
+                        Patterns::Double(0.));
       prm.add_parameter("End time", end_time, "End time", Patterns::Double());
 
       prm.add_parameter("Time step size",

include/solid_mechanics/mf_elasticity.h

@@ -62,6 +62,7 @@ static const unsigned int debug_level = 0;
 #include <adapter/precice_adapter.h>
 #include <base/fe_integrator.h>
 #include <base/q_equidistant.h>
+#include <base/restart.h>
 #include <base/time_handler.h>
 #include <base/utilities.h>
 #include <cases/case_base.h>
@@ -79,8 +80,6 @@ using namespace dealii;
 // into it:
 namespace FSI
 {
-  using namespace dealii;
-
   // The Solid class is the central class.
   template <int dim, typename Number>
   class Solid
@@ -146,6 +145,13 @@ namespace FSI
     void
     output_results(const unsigned int result_number) const;
 
+    void
+    write_restart();
+
+
+    void
+    load_restart();
+
     // Set up an Additional data object
     template <typename AdditionalData>
     void
@@ -558,7 +564,8 @@ namespace FSI
     testcase = testcase_;
     make_grid();
     system_setup();
-    output_results(0);
+    if (parameters.start_time != 0)
+      output_results(0);
     time.increment();
 
     // We then declare the incremental solution update $\varDelta
@@ -573,6 +580,9 @@ namespace FSI
       mf_data_reference);
     precice_adapter->initialize(total_displacement);
 
+    if (parameters.start_time != 0)
+      load_restart();
+
     // At the beginning, we reset the solution update for this time step...
     while (precice_adapter->is_coupling_ongoing())
       {
@@ -614,6 +624,11 @@ namespace FSI
             time.increment();
           }
       }
+    // write the restart after we finished
+    // we decrement the time in the write_restart function, as we call the
+    // function after we incremented the time for the next iteration
+    write_restart();
+
 
     // for post-processing, print average CG iterations over the whole run:
     timer_out << std::endl
@@ -635,7 +650,18 @@ namespace FSI
   {
     Assert(testcase.get() != nullptr, ExcInternalError());
     testcase->make_coarse_grid_and_bcs(triangulation);
-    triangulation.refine_global(parameters.n_global_refinement);
+
+    if (parameters.start_time != 0)
+      {
+        std::string restart_file(
+          parameters.output_folder + "restart-t-" +
+          Utilities::format_time_stamp_to_string(parameters.start_time) +
+          ".mesh");
+        pcout << "--     Looking for restart file \"" << restart_file << "\"\n";
+        triangulation.load(restart_file);
+      }
+    else
+      triangulation.refine_global(parameters.n_global_refinement);
 
     if (testcase->body_force.get() == nullptr)
       testcase->body_force =
@@ -2041,6 +2067,79 @@ namespace FSI
     return std::make_tuple(lin_it, lin_res, cond_number);
   }
 
+
+
+  template <int dim, typename Number>
+  void
+  Solid<dim, Number>::write_restart()
+  {
+    // first, decrement the time, as we call the function after an (unnecessary)
+    // increment
+    double       t_decrement        = time.current() - time.get_delta_t();
+    unsigned int timestep_decrement = time.get_timestep() - 1;
+
+    auto t_string = Utilities::format_time_stamp_to_string(t_decrement);
+
+    std::string restart_file(parameters.output_folder + "restart-t-" +
+                             t_string);
+    pcout << "--     Creating restart files \"" + restart_file +
+               "\" for t = " + t_string
+          << " ( timestep " << timestep_decrement << " ) "
+          << "\n";
+    std::vector<const VectorType *> in_vectors({&total_displacement,
+                                                &velocity,
+                                                &acceleration,
+                                                &old_displacement,
+                                                &velocity_old,
+                                                &acceleration_old});
+
+    // Make sure to pass updated vectors into the function
+    for (auto &in : in_vectors)
+      in->update_ghost_values();
+
+    Utilities::create_restart_snapshot<dim, VectorType>(
+      triangulation,
+      dof_handler,
+      total_displacement.get_partitioner(),
+      in_vectors,
+      restart_file,
+      t_decrement,
+      timestep_decrement);
+  }
+
+
+
+  template <int dim, typename Number>
+  void
+  Solid<dim, Number>::load_restart()
+  {
+    std::vector<VectorType *> in_vectors({&total_displacement,
+                                          &velocity,
+                                          &acceleration,
+                                          &old_displacement,
+                                          &velocity_old,
+                                          &acceleration_old});
+
+    std::string restart_file(
+      parameters.output_folder + "restart-t-" +
+      Utilities::format_time_stamp_to_string(parameters.start_time));
+
+    double       loaded_time     = 0;
+    unsigned int loaded_timestep = 0;
+    Utilities::load_restart_snapshot<dim, VectorType>(
+      dof_handler, in_vectors, restart_file, loaded_time, loaded_timestep);
+
+    pcout << "--     Restarting computation from files \"" + restart_file +
+               "\" at t = "
+          << loaded_time << " ( timestep " << loaded_timestep << " ) "
+          << "\n";
+    time.set_time(loaded_time, loaded_timestep);
+    // we loaded the time we already computed, so let's move on to the next step
+    time.increment();
+  }
+
+
+
   template <int dim, typename Number>
   void
   Solid<dim, Number>::output_results(const unsigned int result_number) const