Nitsche checkpoint particles (#1192)

Description
In the previous version of the code, the lethe-fluid-nitsche solver would crash when trying to checkpoint a simulation that was already restarted. This was because the checkpointing mechanism was poorly implemented. Indeed, instead of checkpointing the particle_handler itself, the solver would actually re-load the triangulation of the solid object and generate new particles. This is not a good approach and is not maintanable.
Solution
I fought agaisnt my laziness and re-implemented the thing to actually checkpoint and load the particles. This now enables things to restart multiple times.
Testing
To test this feature, I updated the existing restart test so that it actually tries to keep on checkpointing. This would throw an assert in Debug. This does not happen anymore and things are all fine and dandy
Documentation
Nothing here needs to be documented, but some issues can be closed because of this.

Co-authored-by: Amishga Alphonius <[email protected]>

CHANGELOG.md

@@ -3,6 +3,12 @@
 All notable changes to the Lethe project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
+## [Master] - 2024-07-09
+
+### Fixed
+
+- MINOR The lethe-fluid-nitsche solver could not checkpoint once it had been restarted. This has been fixed by ensuring that the particles are properly checkpointed and read from the restart files. [#1192](https://github.com/chaos-polymtl/lethe/pull/1192)
+
 ## [Master] - 2024-07-08
 
 ### Changed

applications/lethe-fluid-nitsche/nitsche_navier_stokes.cc

@@ -43,9 +43,6 @@ main(int argc, char *argv[])
           prm.parse_input(argv[1]);
           NSparam.parse(prm);
 
-          AssertThrow(NSparam.nitsche->number_solids > 0,
-                      NoSolidWarning("gls_nitsche_navier_stokes_22"));
-
           GLSNitscheNavierStokesSolver<2> problem_22(NSparam);
           problem_22.solve();
         }
@@ -59,9 +56,6 @@ main(int argc, char *argv[])
           prm.parse_input(argv[1]);
           NSparam.parse(prm);
 
-          AssertThrow(NSparam.nitsche->number_solids > 0,
-                      NoSolidWarning("gls_nitsche_navier_stokes_22"));
-
           GLSNitscheNavierStokesSolver<3> problem_33(NSparam);
           problem_33.solve();
         }

applications_tests/lethe-fluid-nitsche/CMakeLists.txt

@@ -13,6 +13,8 @@ file(COPY two-bar-mixer-restart/restart.triangulation_variable.data DESTINATION
 file(COPY two-bar-mixer-restart/restart.triangulation.info DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/two-bar-mixer-restart.${_build_type}/mpirun=1/")
 file(COPY two-bar-mixer-restart/restart_solid_particles_00.pvdhandler DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/two-bar-mixer-restart.${_build_type}/mpirun=1/")
 file(COPY two-bar-mixer-restart/restart_solid_particles_01.pvdhandler DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/two-bar-mixer-restart.${_build_type}/mpirun=1/")
+file(COPY two-bar-mixer-restart/restart_solid_00.particles DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/two-bar-mixer-restart.${_build_type}/mpirun=1/")
+file(COPY two-bar-mixer-restart/restart_solid_01.particles DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/two-bar-mixer-restart.${_build_type}/mpirun=1/")
 file(COPY two-bar-mixer-restart/restart_solid_00.triangulation DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/two-bar-mixer-restart.${_build_type}/mpirun=1/")
 file(COPY two-bar-mixer-restart/restart_solid_00.triangulation.info DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/two-bar-mixer-restart.${_build_type}/mpirun=1/")
 file(COPY two-bar-mixer-restart/restart_solid_00.triangulation_fixed.data DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/two-bar-mixer-restart.${_build_type}/mpirun=1/")

applications_tests/lethe-fluid-nitsche/two-bar-mixer-start.prm → applications_tests/lethe-fluid-nitsche/generators/two-bar-mixer-start.prm

@@ -72,6 +72,14 @@ subsection initial conditions
   end
 end
 
+#---------------------------------------------------
+# Boundary Conditions
+#---------------------------------------------------
+
+subsection boundary conditions
+  set number = 0
+end
+
 #---------------------------------------------------
 # Physical Properties
 #---------------------------------------------------

applications_tests/lethe-fluid-nitsche/two-bar-mixer-restart.prm

@@ -33,7 +33,7 @@ end
 
 subsection restart
   # Enable checkpointing
-  set checkpoint = false
+  set checkpoint = true
 
   # Enable restart
   set restart = true

applications_tests/lethe-fluid-nitsche/two-bar-mixer-restart/restart_solid_00.particles

@@ -0,0 +1 @@
+0 0 16 4 16

applications_tests/lethe-fluid-nitsche/two-bar-mixer-restart/restart_solid_01.particles

@@ -0,0 +1 @@
+0 0 16 4 16

source/core/solid_base.cc

@@ -722,14 +722,28 @@ template <int dim, int spacedim>
 void
 SolidBase<dim, spacedim>::write_checkpoint(std::string prefix)
 {
-#if (DEAL_II_VERSION_MAJOR < 10 && DEAL_II_VERSION_MINOR < 4)
-  parallel::distributed::
-    SolutionTransfer<dim, GlobalVectorType, DoFHandler<dim, spacedim>>
-      system_trans_vectors(this->displacement_dh);
-#else
+  // First, manage the serialization of the particle information.
+  // The particles are attached to the background triangulation and not
+  // the solid triangulation. Consequently, they only need to be prepared
+  // for serialization and will be serialized when the background triangulation
+  // itself is serialized by the main solver.
+  std::string                   prefix_particles = prefix + "_particles";
+  std::ostringstream            oss;
+  boost::archive::text_oarchive oa(oss, boost::archive::no_header);
+  oa << *(this->solid_particle_handler);
+
+  // Write additional particle information for deserialization
+  std::string   particle_filename = prefix + ".particles";
+  std::ofstream output(particle_filename.c_str());
+  output << oss.str() << std::endl;
+
+  //   Prepare particle handler for serialization
+  this->solid_particle_handler->prepare_for_serialization();
+
+  // Prepare serialization of the solid triangulation which is
+  // only used for visualization purposes
   parallel::distributed::SolutionTransfer<dim, GlobalVectorType, spacedim>
     system_trans_vectors(this->displacement_dh);
-#endif
 
   std::vector<const GlobalVectorType *> sol_set_transfer;
   displacement_relevant = displacement;
@@ -773,14 +787,8 @@ SolidBase<dim, spacedim>::read_checkpoint(std::string prefix)
   std::vector<GlobalVectorType *> x_system(1);
   x_system[0] = &(displacement);
 
-#if (DEAL_II_VERSION_MAJOR < 10 && DEAL_II_VERSION_MINOR < 4)
-  parallel::distributed::
-    SolutionTransfer<dim, GlobalVectorType, DoFHandler<dim, spacedim>>
-      system_trans_vectors(this->displacement_dh);
-#else
   parallel::distributed::SolutionTransfer<dim, GlobalVectorType, spacedim>
     system_trans_vectors(this->displacement_dh);
-#endif
 
   system_trans_vectors.deserialize(x_system);
   displacement_relevant = displacement;
@@ -791,8 +799,28 @@ SolidBase<dim, spacedim>::read_checkpoint(std::string prefix)
   if (!param->solid_mesh.simplex)
     move_solid_triangulation_with_displacement();
 
-  // We did not checkpoint particles, we re-create them from scratch
-  setup_particles();
+  // Gather particle serialization information and deserialize them
+  std::string   particle_filename = prefix + ".particles";
+  std::ifstream input(particle_filename.c_str());
+  AssertThrow(input, ExcFileNotOpen(particle_filename));
+
+  std::string buffer;
+  std::getline(input, buffer);
+  std::istringstream            iss(buffer);
+  boost::archive::text_iarchive ia(iss, boost::archive::no_header);
+
+
+  // At this stage, the pointer for the particle handler has not been allocated
+  // it thus needs to be done by calling the setup
+  // TODO, this function should be called in the constructor
+  this->setup_particles_handler();
+  ia >> *(this->solid_particle_handler);
+
+  this->solid_particle_handler->deserialize();
+
+  // Number of particles is reset to the number of particles in the
+  // solid_particle_handler after deserialization
+  initial_number_of_particles = solid_particle_handler->n_global_particles();
 }
 
 

source/solvers/gls_nitsche_navier_stokes.cc

@@ -114,13 +114,8 @@ GLSNitscheNavierStokesSolver<dim, spacedim>::assemble_nitsche_restriction()
           local_rhs    = 0;
 
 
-#if (DEAL_II_VERSION_MAJOR < 10 && DEAL_II_VERSION_MINOR < 4)
-          const auto &cell =
-            particle->get_surrounding_cell(*this->triangulation);
-#else
-          const auto &cell = particle->get_surrounding_cell();
-#endif
-          double h_cell = 0;
+          const auto &cell   = particle->get_surrounding_cell();
+          double      h_cell = 0;
           if (dim == 2)
             h_cell = std::sqrt(4. * cell->measure() / M_PI) /
                      this->velocity_fem_degree;
@@ -977,9 +972,16 @@ template <int dim, int spacedim>
 void
 GLSNitscheNavierStokesSolver<dim, spacedim>::write_checkpoint()
 {
-  // Prepare solid particles for serialization
-  for (auto solid : solids)
-    solid->get_solid_particle_handler()->prepare_for_serialization();
+  // Write solid base checkpoint
+  for (unsigned int i_solid = 0; i_solid < solids.size(); ++i_solid)
+    {
+      std::string filename =
+        this->simulation_parameters.simulation_control.output_folder +
+        this->simulation_parameters.restart_parameters.filename + "_solid_" +
+        Utilities::int_to_string(i_solid, 2);
+
+      this->solids[i_solid]->write_checkpoint(filename);
+    }
 
   // Call regular checkpointing routine
   this->GLSNavierStokesSolver<spacedim>::write_checkpoint();
@@ -1005,16 +1007,6 @@ GLSNitscheNavierStokesSolver<dim, spacedim>::write_checkpoint()
             Utilities::int_to_string(i_solid, 2));
         }
     }
-
-  // Write solid base checkpoint
-  for (unsigned int i_solid = 0; i_solid < solids.size(); ++i_solid)
-    {
-      std::string filename =
-        this->simulation_parameters.simulation_control.output_folder +
-        this->simulation_parameters.restart_parameters.filename + "_solid_" +
-        Utilities::int_to_string(i_solid, 2);
-      this->solids[i_solid]->write_checkpoint(filename);
-    }
 }
 
 template <int dim, int spacedim>
@@ -1024,10 +1016,12 @@ GLSNitscheNavierStokesSolver<dim, spacedim>::read_checkpoint()
   this->GLSNavierStokesSolver<spacedim>::read_checkpoint();
 
   TimerOutput::Scope t(this->computing_timer,
-                       "Reset Nitsche solid mesh and particles");
-
+                       "Reload Nitsche solid mesh and particles");
 
   // Reload initial configurations
+  // This must be done after the background triangulation is read
+  // because the solid particle_handler is deserialized physically
+  // after the triangulation is read.
   for (unsigned int i_solid = 0; i_solid < solids.size(); ++i_solid)
     {
       std::string prefix =
@@ -1038,6 +1032,7 @@ GLSNitscheNavierStokesSolver<dim, spacedim>::read_checkpoint()
       solids[i_solid]->read_checkpoint(prefix);
     }
 
+
   // Reload particle and solid pvd handlers
   std::string prefix =
     this->simulation_parameters.simulation_control.output_folder +