Implemented checkpoint-restart

include/lac/lac_initializer.h

@@ -3,6 +3,7 @@
 
 // This includes all types we know of.
 #include "lac/lac_type.h"
+#include <deal.II/base/mpi.h>
 #include <mpi.h>
 
 /**
@@ -15,20 +16,20 @@ class ScopedLACInitializer
   ScopedLACInitializer(const std::vector<types::global_dof_index> &dofs_per_block,
                        const std::vector<IndexSet> &owned,
                        const std::vector<IndexSet> &relevant,
-                       const MPI_Comm &comm = MPI_COMM_WORLD):
+                       const MPI_Comm extcomm= MPI_COMM_WORLD):
     dofs_per_block(dofs_per_block),
     owned(owned),
     relevant(relevant),
-    comm(comm)
-  {};
+    comm(Utilities::MPI::duplicate_communicator(extcomm))
+  {}
 
   /**
    * Initialize a non ghosted TrilinosWrappers::MPI::BlockVector.
    */
   void operator() (TrilinosWrappers::MPI::BlockVector &v, bool fast=false)
   {
     v.reinit(owned, comm, fast);
-  };
+  }
 
 
   /**
@@ -37,15 +38,15 @@ class ScopedLACInitializer
   void ghosted(TrilinosWrappers::MPI::BlockVector &v, bool fast=false)
   {
     v.reinit(owned, relevant, comm, fast);
-  };
+  }
 
   /**
    * Initialize a serial BlockVector<double>.
    */
   void operator() (BlockVector<double> &v, bool fast=false)
   {
     v.reinit(dofs_per_block, fast);
-  };
+  }
 
 
   /**
@@ -56,7 +57,7 @@ class ScopedLACInitializer
   {
     Assert(false, ExcInternalError("You tried to create a ghosted vector in a serial run."));
     (void)fast;
-  };
+  }
 
   /**
    * Initialize a Trilinos Sparsity Pattern.
@@ -112,7 +113,7 @@ class ScopedLACInitializer
   /**
    * MPI Communicator.
    */
-  const MPI_Comm &comm;
+  MPI_Comm comm;
 };
 
 

include/pidomus.h

@@ -17,23 +17,17 @@
 
 
 #include <deal.II/base/timer.h>
-// #include <deal.II/base/parameter_handler.h>
 
 #include <deal.II/lac/trilinos_precondition.h>
 #include <deal.II/lac/trilinos_block_sparse_matrix.h>
 #include <deal.II/lac/linear_operator.h>
 
-
-
 #include <mpi.h>
 
-
-// #include <deal.II/lac/precondition.h>
-
-
 #include "base_interface.h"
 #include "simulator_access.h"
 #include "pidomus_signals.h"
+#include "copy_data.h"
 
 #include <deal2lkit/parsed_grid_generator.h>
 #include <deal2lkit/parsed_finite_element.h>
@@ -72,7 +66,7 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
 
   piDoMUS (const std::string &name,
            const BaseInterface<dim, spacedim, LAC> &energy,
-           const MPI_Comm &comm = MPI_COMM_WORLD);
+           const MPI_Comm comm = MPI_COMM_WORLD);
 
   virtual void declare_parameters(ParameterHandler &prm);
   virtual void parse_parameters_call_back();
@@ -307,7 +301,7 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
 
   void set_constrained_dofs_to_zero(typename LAC::VectorType &v) const;
 
-  const MPI_Comm &comm;
+  MPI_Comm comm;
   const BaseInterface<dim, spacedim, LAC>    &interface;
 
   unsigned int n_cycles;
@@ -522,6 +516,86 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
                   const typename LAC::VectorType &solution_dot);
 
 
+  // checkpoint restart
+
+  // used in output solution
+  unsigned int step_number;
+
+  unsigned int old_step;
+
+  /** resume compuatation
+   * from a previously created snapshot
+   */
+  bool resume_computation;
+
+  /**
+   * create snapshot throughout the simulation
+   */
+  bool save_snapshot;
+
+  /**
+   * this prefix is put in front of the
+   * common names used for checkpoint- restart
+   */
+  std::string snap_prefix;
+
+  /**
+    * Save the state of this program to a set of files in the output
+    * directory. In reality, however, only some variables are stored (in
+    * particular the mesh, the solution vectors, etc) whereas others can
+    * either be re-generated (matrices, DoFHandler objects, etc) or are
+    * read from the input parameter file.
+    *
+    * This function is implemented in
+    * <code>source/checkpoint_restart.cc</code>.
+    */
+  void create_snapshot() const;
+
+  void save_solutions_and_triangulation(const LADealII::VectorType &y,
+                                        const LADealII::VectorType &y_dot,
+                                        const LADealII::VectorType &locally_relevant_y_expl,
+                                        const LADealII::VectorType &,
+                                        const LADealII::VectorType &) const;
+
+  void save_solutions_and_triangulation(const LATrilinos::VectorType &y,
+                                        const LATrilinos::VectorType &y_dot,
+                                        const LATrilinos::VectorType &locally_relevant_y_expl,
+                                        const LATrilinos::VectorType &,
+                                        const LATrilinos::VectorType &) const;
+
+
+  /**
+   * Restore the state of this program from a set of files in the output
+   * directory.
+   *
+   * This function is implemented in
+   * <code>source/checkpoint_restart.cc</code>.
+   */
+  void resume_from_snapshot();
+
+  void load_solutions(LADealII::VectorType &y,
+                      LADealII::VectorType &y_expl,
+                      LADealII::VectorType &y_dot);
+
+  void load_solutions(LATrilinos::VectorType &y,
+                      LATrilinos::VectorType &y_expl,
+                      LATrilinos::VectorType &y_dot);
+
+public:
+  /**
+   * Save a number of variables using BOOST serialization mechanism.
+   *
+   * This function is implemented in
+   * <code>source/checkpoint_restart.cc</code>.
+   */
+  template <class Archive>
+  void serialize (Archive &ar, const unsigned int /*version*/);
+
+private:
+
+
+
+
   /**
    * Struct containing the signals
    */

source/pidomus.cc

@@ -44,11 +44,11 @@ using namespace deal2lkit;
 template <int dim, int spacedim, typename LAC>
 piDoMUS<dim, spacedim, LAC>::piDoMUS (const std::string &name,
                                       const BaseInterface<dim, spacedim, LAC> &interface,
-                                      const MPI_Comm &communicator)
+                                      const MPI_Comm communicator)
   :
   ParameterAcceptor(name),
   SundialsInterface<typename LAC::VectorType>(communicator),
-  comm(communicator),
+  comm(Utilities::MPI::duplicate_communicator(communicator)),
   interface(interface),
   pcout (std::cout,
          (Utilities::MPI::this_mpi_process(comm)
@@ -100,7 +100,9 @@ piDoMUS<dim, spacedim, LAC>::piDoMUS (const std::string &name,
 
   ida(*this),
   euler(*this),
-  we_are_parallel(Utilities::MPI::n_mpi_processes(comm) > 1)
+  we_are_parallel(Utilities::MPI::n_mpi_processes(comm) > 1),
+  step_number(0),
+  old_step(0)
 {
 
   interface.initialize_simulator (*this);
@@ -123,15 +125,31 @@ void piDoMUS<dim, spacedim, LAC>::run ()
 
   interface.connect_to_signals();
 
-  for (current_cycle = 0; current_cycle < n_cycles; ++current_cycle)
+  if (resume_computation)
     {
-      if (current_cycle == 0)
+      resume_from_snapshot();
+
+      syncronize(current_time,solution,solution_dot);
+      ida.set_initial_time(current_time);
+      euler.set_initial_time(current_time+current_dt);
+    }
+
+  current_cycle = (resume_computation? current_cycle : 0);
+  old_step = (resume_computation? step_number : 0);
+
+  for (; current_cycle < n_cycles; ++current_cycle)
+    {
+      if (!resume_computation)
         {
-          make_grid_fe();
-          setup_dofs(true);
+          if (current_cycle == 0)
+            {
+              make_grid_fe();
+              setup_dofs(true);
+            }
+          else
+            refine_mesh();
         }
-      else
-        refine_mesh();
+      resume_computation = false;
 
       constraints.distribute(solution);
       constraints_dot.distribute(solution_dot);