added restart I/O feature

include/etc.hpp

@@ -27,4 +27,15 @@ inline void DeletePointers(Array2D<T*> &ptr_array)
 
 bool FileExists(const std::string& name);
 
+template<typename ... Args>
+std::string string_format( const std::string& format, Args ... args )
+{
+    int size_s = std::snprintf( nullptr, 0, format.c_str(), args ... ) + 1; // Extra space for '\0'
+    if( size_s <= 0 ){ throw std::runtime_error( "Error during formatting." ); }
+    auto size = static_cast<size_t>( size_s );
+    std::unique_ptr<char[]> buf( new char[ size ] );
+    std::snprintf( buf.get(), size, format.c_str(), args ... );
+    return std::string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside
+}
+
 #endif

include/multiblock_solver.hpp

@@ -237,7 +237,9 @@ friend class ParameterizedProblem;
    virtual void SaveVisualization(const int step, const double time);
 
    void SaveSolution(std::string filename = "");
+   void SaveSolutionWithTime(const std::string filename, const int step, const double time);
    void LoadSolution(const std::string &filename);
+   void LoadSolutionWithTime(const std::string &filename, int &step, double &time);
    void CopySolution(BlockVector *input_sol);
 
    virtual void AllocateROMTensorElems()

include/unsteady_ns_solver.hpp

@@ -26,6 +26,8 @@ friend class SteadyNSOperator;
    int time_order = -1;
    // report interval for time integration
    int report_interval = 0;
+   // restart save interval
+   int restart_interval = 0;
 
    // BDFk/EXTk coefficients.
    /* use first order for now. */

src/multiblock_solver.cpp

@@ -807,6 +807,25 @@ void MultiBlockSolver::SaveSolution(std::string filename)
    printf("Done!\n");
 }
 
+void MultiBlockSolver::SaveSolutionWithTime(std::string filename, const int step, const double time)
+{
+   SaveSolution(filename);
+   printf("Saving time/time step ...");
+
+   hid_t file_id;
+   herr_t errf = 0;
+   file_id = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
+   assert(file_id >= 0);
+
+   // TODO: currently we only need solution vector. But we can add more data as we need.
+   hdf5_utils::WriteAttribute(file_id, "timestep", step);
+   hdf5_utils::WriteAttribute(file_id, "time", time);
+
+   errf = H5Fclose(file_id);
+   assert(errf >= 0);
+   printf("Done!\n");
+}
+
 void MultiBlockSolver::LoadSolution(const std::string &filename)
 {
    // solution vector U must be instantiated beforehand.
@@ -839,6 +858,25 @@ void MultiBlockSolver::LoadSolution(const std::string &filename)
       mfem_error("solution size does not match!\n");
 }
 
+void MultiBlockSolver::LoadSolutionWithTime(const std::string &filename, int &step, double &time)
+{
+   LoadSolution(filename);
+   printf("Loading time/time step ...");
+
+   hid_t file_id;
+   herr_t errf = 0;
+   file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+   assert(file_id >= 0);
+
+   // TODO: currently we only need solution vector. But we can add more data as we need.
+   hdf5_utils::ReadAttribute(file_id, "timestep", step);
+   hdf5_utils::ReadAttribute(file_id, "time", time);
+
+   errf = H5Fclose(file_id);
+   assert(errf >= 0);
+   printf("Done!\n");
+}
+
 void MultiBlockSolver::CopySolution(BlockVector *input_sol)
 {
    assert(input_sol->NumBlocks() == U->NumBlocks());

src/parameter.cpp

@@ -11,17 +11,6 @@
 using namespace mfem;
 using namespace std;
 
-template<typename ... Args>
-std::string string_format( const std::string& format, Args ... args )
-{
-    int size_s = std::snprintf( nullptr, 0, format.c_str(), args ... ) + 1; // Extra space for '\0'
-    if( size_s <= 0 ){ throw std::runtime_error( "Error during formatting." ); }
-    auto size = static_cast<size_t>( size_s );
-    std::unique_ptr<char[]> buf( new char[ size ] );
-    std::snprintf( buf.get(), size, format.c_str(), args ... );
-    return std::string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside
-}
-
 DoubleParam::DoubleParam(const std::string &input_key, YAML::Node option)
    : Parameter(input_key),
      log_scale(config.GetOptionFromDict<bool>("log_scale", false, option)),

src/unsteady_ns_solver.cpp

@@ -26,6 +26,9 @@ UnsteadyNSSolver::UnsteadyNSSolver()
    time_order = config.GetOption<int>("time-integration/bdf_order", 1);
    report_interval = config.GetOption<int>("time-integration/report_interval", 0);
 
+   if (save_sol)
+      restart_interval = config.GetOption<int>("save_solution/restart_interval", 0);
+
    if (time_order != 1)
       mfem_error("UnsteadyNSSolver supports only first-order time integration for now!\n");
 }
@@ -77,12 +80,22 @@ void UnsteadyNSSolver::AssembleOperator()
 bool UnsteadyNSSolver::Solve()
 {
    bool converged = true;
+   int initial_step = 0;
+   double time = 0.0;
+
+   bool use_restart = config.GetOption<bool>("solver/use_restart", false);
+   std::string restart_file, file_fmt;
+   file_fmt = "%s/%s_%08d.h5";
+   if (use_restart)
+   {
+      restart_file = config.GetRequiredOption<std::string>("solver/restart_file");
+      LoadSolutionWithTime(restart_file, initial_step, time);
+   }
 
    InitializeTimeIntegration();
 
    SortByVariables(*U, *U_step);
 
-   double time = 0.0;
    SaveVisualization(0, time);
 
    double cfl = 0.0;
@@ -97,6 +110,12 @@ bool UnsteadyNSSolver::Solve()
 
       if (((step+1) % visual.time_interval) == 0)
          SaveVisualization(step+1, time);
+
+      if ((save_sol) && (((step+1) % restart_interval) == 0))
+      {
+         restart_file = string_format(file_fmt, sol_dir.c_str(), sol_prefix.c_str(), step+1);
+         SaveSolutionWithTime(restart_file, step+1, time);
+      }
    }
 
    SortBySubdomains(*U_step, *U);