[PDI] Allow restarting simulations from a previously saved hdf5 file

gyselax · Nov 15, 2023 · 90b99db · 90b99db
1 parent 4b44e07
commit 90b99db
Show file tree

Hide file tree

Showing 22 changed files with 503 additions and 61 deletions.
diff --git a/simulations/geometryXVx/bump_on_tail/bumpontail_fem_uniform.cpp b/simulations/geometryXVx/bump_on_tail/bumpontail_fem_uniform.cpp
@@ -30,6 +30,7 @@
 #include "params.yaml.hpp"
 #include "pdi_out.yml.hpp"
 #include "predcorr.hpp"
+#include "restartinitialization.hpp"
 #include "singlemodeperturbinitialization.hpp"
 #include "species_info.hpp"
 #include "spline_interpolator.hpp"
@@ -54,6 +55,7 @@ int main(int argc, char** argv)
 {
     ddc::ScopeGuard scope(argc, argv);
 
+    long int iter_start(0);
     PC_tree_t conf_voicexx;
     if (argc == 2) {
         conf_voicexx = PC_parse_path(fs::path(argv[1]).c_str());
@@ -63,8 +65,15 @@ int main(int argc, char** argv)
             file << params_yaml;
             return EXIT_SUCCESS;
         }
+    } else if (argc == 4) {
+        if (argv[1] == std::string_view("--iter-restart")) {
+            iter_start = std::strtol(argv[2], NULL, 10);
+            conf_voicexx = PC_parse_path(fs::path(argv[3]).c_str());
+        }
     } else {
         cerr << "usage: " << argv[0] << " [--dump-config] <config_file.yml>" << endl;
+        cerr << "or to perform a restart" << argv[0] << " [--iter-restart] <iter> <config_file.yml>"
+             << endl;
         return EXIT_FAILURE;
     }
     PC_errhandler(PC_NULL_HANDLER);
@@ -155,11 +164,23 @@ int main(int argc, char** argv)
             std::move(mean_velocity_bot));
     init_fequilibrium(allfequilibrium);
     DFieldSpXVx allfdistribu(meshSpXVx);
-    SingleModePerturbInitialization const
-            init(allfequilibrium,
-                 ddc::discrete_space<IDimSp>().perturb_modes(),
-                 ddc::discrete_space<IDimSp>().perturb_amplitudes());
-    init(allfdistribu);
+
+    PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
+    PDI_init(conf_pdi);
+
+    ddc::expose_to_pdi("iter_start", iter_start);
+
+    double time_start(0);
+    if (iter_start == 0) {
+        SingleModePerturbInitialization const
+                init(allfequilibrium,
+                     ddc::discrete_space<IDimSp>().perturb_modes(),
+                     ddc::discrete_space<IDimSp>().perturb_amplitudes());
+        init(allfdistribu);
+    } else {
+        RestartInitialization const restart(iter_start, time_start);
+        restart(allfdistribu);
+    }
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -169,10 +190,6 @@ int main(int argc, char** argv)
     double const time_diag = PCpp_double(conf_voicexx, ".Output.time_diag");
     int const nbstep_diag = int(time_diag / deltat);
 
-    PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
-
-    PDI_init(conf_pdi);
-
     ConstantExtrapolationBoundaryValue<BSplinesX> bv_x_min(x_min);
     ConstantExtrapolationBoundaryValue<BSplinesX> bv_x_max(x_max);
 
@@ -230,7 +247,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, deltat, nbiter);
+    predcorr(allfdistribu, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

diff --git a/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp b/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp
@@ -26,6 +26,7 @@
 #include "params.yaml.hpp"
 #include "pdi_out.yml.hpp"
 #include "predcorr.hpp"
+#include "restartinitialization.hpp"
 #include "singlemodeperturbinitialization.hpp"
 #include "species_info.hpp"
 #include "spline_interpolator.hpp"
@@ -50,6 +51,7 @@ int main(int argc, char** argv)
 {
     ddc::ScopeGuard scope(argc, argv);
 
+    long int iter_start(0);
     PC_tree_t conf_voicexx;
     if (argc == 2) {
         conf_voicexx = PC_parse_path(fs::path(argv[1]).c_str());
@@ -59,8 +61,15 @@ int main(int argc, char** argv)
             file << params_yaml;
             return EXIT_SUCCESS;
         }
+    } else if (argc == 4) {
+        if (argv[1] == std::string_view("--iter-restart")) {
+            iter_start = std::strtol(argv[2], NULL, 10);
+            conf_voicexx = PC_parse_path(fs::path(argv[3]).c_str());
+        }
     } else {
         cerr << "usage: " << argv[0] << " [--dump-config] <config_file.yml>" << endl;
+        cerr << "or to perform a restart" << argv[0] << " [--iter-restart] <iter> <config_file.yml>"
+             << endl;
         return EXIT_FAILURE;
     }
     PC_errhandler(PC_NULL_HANDLER);
@@ -147,11 +156,23 @@ int main(int argc, char** argv)
             std::move(mean_velocity_bot));
     init_fequilibrium(allfequilibrium);
     DFieldSpXVx allfdistribu(meshSpXVx);
-    SingleModePerturbInitialization const
-            init(allfequilibrium,
-                 ddc::discrete_space<IDimSp>().perturb_modes(),
-                 ddc::discrete_space<IDimSp>().perturb_amplitudes());
-    init(allfdistribu);
+
+    PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
+    PDI_init(conf_pdi);
+
+    ddc::expose_to_pdi("iter_start", iter_start);
+
+    double time_start(0);
+    if (iter_start == 0) {
+        SingleModePerturbInitialization const
+                init(allfequilibrium,
+                     ddc::discrete_space<IDimSp>().perturb_modes(),
+                     ddc::discrete_space<IDimSp>().perturb_amplitudes());
+        init(allfdistribu);
+    } else {
+        RestartInitialization const restart(iter_start, time_start);
+        restart(allfdistribu);
+    }
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -161,10 +182,6 @@ int main(int argc, char** argv)
     double const time_diag = PCpp_double(conf_voicexx, ".Output.time_diag");
     int const nbstep_diag = int(time_diag / deltat);
 
-    PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
-
-    PDI_init(conf_pdi);
-
     ConstantExtrapolationBoundaryValue<BSplinesX> bv_x_min(x_min);
     ConstantExtrapolationBoundaryValue<BSplinesX> bv_x_max(x_max);
 
@@ -219,7 +236,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, deltat, nbiter);
+    predcorr(allfdistribu, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

diff --git a/simulations/geometryXVx/bump_on_tail/pdi_out.yml.hpp b/simulations/geometryXVx/bump_on_tail/pdi_out.yml.hpp
@@ -5,6 +5,7 @@ constexpr char const* const PDI_CFG = R"PDI_CFG(
   Nx : int
   Nvx : int
   iter : int
+  iter_start : int
   time_saved : double
   nbstep_diag: int
   iter_saved : int
@@ -55,7 +56,7 @@ constexpr char const* const PDI_CFG = R"PDI_CFG(
     on_data:
       iter:
         - set:
-          - iter_saved: '${iter}/${nbstep_diag}'
+          - iter_saved: '${iter_start} + ${iter}/${nbstep_diag}'
     on_finalize:
       - release: [iter_saved]
   decl_hdf5:
@@ -68,5 +69,8 @@ constexpr char const* const PDI_CFG = R"PDI_CFG(
       when: '${iter} % ${nbstep_diag} = 0'
       collision_policy: replace_and_warn
       write: [time_saved, fdistribu, electrostatic_potential]
+    - file: 'VOICEXX_${iter_start:05}.h5'
+      on_event: restart
+      read: [time_saved, fdistribu]
   #trace: ~
 )PDI_CFG";
diff --git a/simulations/geometryXVx/landau/landau_fem_uniform.cpp b/simulations/geometryXVx/landau/landau_fem_uniform.cpp
@@ -29,6 +29,7 @@
 #include "params.yaml.hpp"
 #include "pdi_out.yml.hpp"
 #include "predcorr.hpp"
+#include "restartinitialization.hpp"
 #include "singlemodeperturbinitialization.hpp"
 #include "species_info.hpp"
 #include "spline_interpolator.hpp"
@@ -53,6 +54,7 @@ int main(int argc, char** argv)
 {
     ddc::ScopeGuard scope(argc, argv);
 
+    long int iter_start(0);
     PC_tree_t conf_voicexx;
     if (argc == 2) {
         conf_voicexx = PC_parse_path(fs::path(argv[1]).c_str());
@@ -62,8 +64,15 @@ int main(int argc, char** argv)
             file << params_yaml;
             return EXIT_SUCCESS;
         }
+    } else if (argc == 4) {
+        if (argv[1] == std::string_view("--iter-restart")) {
+            iter_start = std::strtol(argv[2], NULL, 10);
+            conf_voicexx = PC_parse_path(fs::path(argv[3]).c_str());
+        }
     } else {
         cerr << "usage: " << argv[0] << " [--dump-config] <config_file.yml>" << endl;
+        cerr << "or to perform a restart" << argv[0] << " [--iter-restart] <iter> <config_file.yml>"
+             << endl;
         return EXIT_FAILURE;
     }
     PC_errhandler(PC_NULL_HANDLER);
@@ -153,11 +162,23 @@ int main(int argc, char** argv)
             std::move(mean_velocity_eq));
     init_fequilibrium(allfequilibrium);
     DFieldSpXVx allfdistribu(meshSpXVx);
-    SingleModePerturbInitialization const
-            init(allfequilibrium,
-                 ddc::discrete_space<IDimSp>().perturb_modes(),
-                 ddc::discrete_space<IDimSp>().perturb_amplitudes());
-    init(allfdistribu);
+
+    PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
+    PDI_init(conf_pdi);
+
+    ddc::expose_to_pdi("iter_start", iter_start);
+
+    double time_start(0);
+    if (iter_start == 0) {
+        SingleModePerturbInitialization const
+                init(allfequilibrium,
+                     ddc::discrete_space<IDimSp>().perturb_modes(),
+                     ddc::discrete_space<IDimSp>().perturb_amplitudes());
+        init(allfdistribu);
+    } else {
+        RestartInitialization const restart(iter_start, time_start);
+        restart(allfdistribu);
+    }
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -167,10 +188,6 @@ int main(int argc, char** argv)
     double const time_diag = PCpp_double(conf_voicexx, ".Output.time_diag");
     int const nbstep_diag = int(time_diag / deltat);
 
-    PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
-
-    PDI_init(conf_pdi);
-
     ConstantExtrapolationBoundaryValue<BSplinesX> bv_x_min(x_min);
     ConstantExtrapolationBoundaryValue<BSplinesX> bv_x_max(x_max);
 
@@ -229,7 +246,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, deltat, nbiter);
+    predcorr(allfdistribu, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

diff --git a/simulations/geometryXVx/landau/landau_fft.cpp b/simulations/geometryXVx/landau/landau_fft.cpp
@@ -26,6 +26,7 @@
 #include "params.yaml.hpp"
 #include "pdi_out.yml.hpp"
 #include "predcorr.hpp"
+#include "restartinitialization.hpp"
 #include "singlemodeperturbinitialization.hpp"
 #include "species_info.hpp"
 #include "spline_interpolator.hpp"
@@ -50,6 +51,7 @@ int main(int argc, char** argv)
 {
     ddc::ScopeGuard scope(argc, argv);
 
+    long int iter_start(0);
     PC_tree_t conf_voicexx;
     if (argc == 2) {
         conf_voicexx = PC_parse_path(fs::path(argv[1]).c_str());
@@ -59,8 +61,15 @@ int main(int argc, char** argv)
             file << params_yaml;
             return EXIT_SUCCESS;
         }
+    } else if (argc == 4) {
+        if (argv[1] == std::string_view("--iter-restart")) {
+            iter_start = std::strtol(argv[2], NULL, 10);
+            conf_voicexx = PC_parse_path(fs::path(argv[3]).c_str());
+        }
     } else {
         cerr << "usage: " << argv[0] << " [--dump-config] <config_file.yml>" << endl;
+        cerr << "or to perform a restart" << argv[0] << " [--iter-restart] <iter> <config_file.yml>"
+             << endl;
         return EXIT_FAILURE;
     }
     PC_errhandler(PC_NULL_HANDLER);
@@ -150,11 +159,23 @@ int main(int argc, char** argv)
             std::move(mean_velocity_eq));
     init_fequilibrium(allfequilibrium);
     DFieldSpXVx allfdistribu(meshSpXVx);
-    SingleModePerturbInitialization const
-            init(allfequilibrium,
-                 ddc::discrete_space<IDimSp>().perturb_modes(),
-                 ddc::discrete_space<IDimSp>().perturb_amplitudes());
-    init(allfdistribu);
+
+    PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
+    PDI_init(conf_pdi);
+
+    ddc::expose_to_pdi("iter_start", iter_start);
+
+    double time_start(0);
+    if (iter_start == 0) {
+        SingleModePerturbInitialization const
+                init(allfequilibrium,
+                     ddc::discrete_space<IDimSp>().perturb_modes(),
+                     ddc::discrete_space<IDimSp>().perturb_amplitudes());
+        init(allfdistribu);
+    } else {
+        RestartInitialization const restart(iter_start, time_start);
+        restart(allfdistribu);
+    }
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -164,10 +185,6 @@ int main(int argc, char** argv)
     double const time_diag = PCpp_double(conf_voicexx, ".Output.time_diag");
     int const nbstep_diag = int(time_diag / deltat);
 
-    PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
-
-    PDI_init(conf_pdi);
-
     ConstantExtrapolationBoundaryValue<BSplinesX> bv_x_min(x_min);
     ConstantExtrapolationBoundaryValue<BSplinesX> bv_x_max(x_max);
 
@@ -222,7 +239,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, deltat, nbiter);
+    predcorr(allfdistribu, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

diff --git a/simulations/geometryXVx/landau/pdi_out.yml.hpp b/simulations/geometryXVx/landau/pdi_out.yml.hpp
@@ -5,6 +5,7 @@ constexpr char const* const PDI_CFG = R"PDI_CFG(
   Nx : int
   Nvx : int
   iter : int
+  iter_start : int
   time_saved : double
   nbstep_diag: int
   iter_saved : int
@@ -55,7 +56,7 @@ constexpr char const* const PDI_CFG = R"PDI_CFG(
     on_data:
       iter:
         - set:
-          - iter_saved: '${iter}/${nbstep_diag}'
+          - iter_saved: '${iter_start} + ${iter}/${nbstep_diag}'
     on_finalize:
       - release: [iter_saved]
   decl_hdf5:
@@ -68,5 +69,8 @@ constexpr char const* const PDI_CFG = R"PDI_CFG(
       when: '${iter} % ${nbstep_diag} = 0'
       collision_policy: replace_and_warn
       write: [time_saved, fdistribu, electrostatic_potential]
+    - file: 'VOICEXX_${iter_start:05}.h5'
+      on_event: restart
+      read: [time_saved, fdistribu]
   #trace: ~
 )PDI_CFG";