Lost particle track files

include/openmc/particle_restart.h

@@ -1,9 +1,12 @@
 #ifndef OPENMC_PARTICLE_RESTART_H
 #define OPENMC_PARTICLE_RESTART_H
 
+#include "openmc/particle.h"
+
 namespace openmc {
 
 void run_particle_restart();
+void run_lost_particle_track(Particle& lost);
 
 } // namespace openmc
 

include/openmc/track_output.h

@@ -2,6 +2,7 @@
 #define OPENMC_TRACK_OUTPUT_H
 
 #include "openmc/particle.h"
+#include <hdf5.h>
 
 namespace openmc {
 
@@ -36,6 +37,12 @@ void write_particle_track(Particle& p);
 //! \param[in] p  Current particle
 void finalize_particle_track(Particle& p);
 
+// Check whether --track is calling or lost particles are calling automatically
+extern bool lost_particle_track_file_open;
+extern hid_t track_file;
+extern hid_t track_dtype;
+extern thread_local bool in_lost_track;
+
 } // namespace openmc
 
 #endif // OPENMC_TRACK_OUTPUT_H

src/particle.cpp

@@ -20,6 +20,7 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/particle_data.h"
+#include "openmc/particle_restart.h"
 #include "openmc/photon.h"
 #include "openmc/physics.h"
 #include "openmc/physics_mg.h"
@@ -837,11 +838,24 @@ void Particle::cross_periodic_bc(
 
 void Particle::mark_as_lost(const char* message)
 {
+  // Skip if we are already replaying a lost particle
+  if (in_lost_track) {
+    wgt() = 0.0;
+    return;
+  }
+
   // Print warning and write lost particle file
-  warning(message);
+#pragma omp critical(PrintErrorMessage)
+  {
+    warning(message);
+  }
   if (settings::max_write_lost_particles < 0 ||
       simulation::n_lost_particles < settings::max_write_lost_particles) {
     write_restart();
+    if (settings::run_mode != RunMode::PARTICLE &&
+        !settings::write_all_tracks) {
+      run_lost_particle_track(*this);
+    }
   }
   // Increment number of lost particles
   wgt() = 0.0;
@@ -852,11 +866,21 @@ void Particle::mark_as_lost(const char* message)
   auto n = simulation::current_batch * settings::gen_per_batch *
            simulation::work_per_rank;
 
-  // Abort the simulation if the maximum number of lost particles has been
-  // reached
-  if (simulation::n_lost_particles >= settings::max_lost_particles &&
-      simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {
-    fatal_error("Maximum number of lost particles has been reached.");
+// Abort the simulation if the maximum number of lost particles reached
+#pragma omp critical(TrackFile)
+  {
+    if (simulation::n_lost_particles >= settings::max_lost_particles &&
+        simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {
+      // close track file manually
+      if (track_file >= 0) {
+        H5Tclose(track_dtype);
+        file_close(track_file);
+        track_file = -1;
+        track_dtype = -1;
+      }
+      lost_particle_track_file_open = false;
+      fatal_error("Maximum number of lost particles has been reached.");
+    }
   }
 }
 

src/particle_restart.cpp

@@ -7,11 +7,11 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
-#include "openmc/particle.h"
 #include "openmc/photon.h"
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
+#include "openmc/source.h"
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/track_output.h"
@@ -139,4 +139,58 @@ void run_particle_restart()
   }
 }
 
+// Function for automatically creating tracks file for lost particles
+void run_lost_particle_track(Particle& lost)
+{
+  if (in_lost_track)
+    return;
+  in_lost_track = true;
+
+#pragma omp critical(TrackFile)
+  {
+    if (track_file < 0 && !lost_particle_track_file_open) {
+      open_track_file();
+      lost_particle_track_file_open = true;
+    }
+  }
+
+  Particle p;
+  p.id() = lost.id();
+
+  int64_t i = lost.current_work();
+  SourceSite site;
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    site = simulation::source_bank[i];
+  } else if (settings::run_mode == RunMode::FIXED_SOURCE &&
+             settings::use_shared_secondary_bank &&
+             i < simulation::shared_secondary_bank_read.size()) {
+    site = simulation::shared_secondary_bank_read[i];
+  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
+    int64_t id = compute_transport_seed(compute_particle_id(i + 1));
+    uint64_t seed = init_seed(id, STREAM_SOURCE);
+    site = sample_external_source(&seed);
+  }
+
+  p.wgt() = site.wgt;
+  p.E() = site.E;
+  p.r() = site.r;
+  p.u() = site.u;
+  p.time() = site.time;
+  p.type() = lost.type();
+
+  int64_t particle_seed = compute_transport_seed(p.id());
+  init_particle_seeds(particle_seed, p.seeds());
+
+  if (settings::run_CE) {
+    p.invalidate_neutron_xs();
+  }
+
+  p.write_track() = true;
+  add_particle_track(p);
+  transport_history_based_single_particle(p);
+  // finalize_particle_track called internally by transport, don't call again
+
+  in_lost_track = false;
+}
+
 } // namespace openmc

src/simulation.cpp

@@ -80,6 +80,11 @@ int openmc_simulation_init()
   if (simulation::initialized)
     return 0;
 
+  // Reset lost particle track file state for new simulation
+  lost_particle_track_file_open = false;
+  in_lost_track = false;
+  simulation::n_lost_particles = 0;
+
   // Initialize nuclear data (energy limits, log grid)
   if (settings::run_CE) {
     initialize_data();
@@ -188,6 +193,10 @@ int openmc_simulation_finalize()
   if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
     close_track_file();
   }
+  if (lost_particle_track_file_open) {
+    close_track_file();
+    lost_particle_track_file_open = false;
+  }
 
   // Increment total number of generations
   simulation::total_gen += simulation::current_batch * settings::gen_per_batch;

src/track_output.cpp

@@ -10,7 +10,7 @@
 
 #include "openmc/tensor.h"
 #include <fmt/core.h>
-#include <hdf5.h>
+// #include <hdf5.h>
 
 #include <cstddef> // for size_t
 #include <string>
@@ -21,9 +21,11 @@ namespace openmc {
 // Global variables
 //==============================================================================
 
-hid_t track_file;     //! HDF5 identifier for track file
-hid_t track_dtype;    //! HDF5 identifier for track datatype
-int n_tracks_written; //! Number of tracks written
+hid_t track_file {-1};  //! HDF5 identifier for track file
+hid_t track_dtype {-1}; //! HDF5 identifier for track datatype
+int n_tracks_written;   //! Number of tracks written
+bool lost_particle_track_file_open {false};
+thread_local bool in_lost_track {false};
 
 //==============================================================================
 // Non-member functions
@@ -42,6 +44,14 @@ void write_particle_track(Particle& p)
 
 void open_track_file()
 {
+  // Close existing file handle if still open from a previous run
+  if (track_file >= 0) {
+    H5Tclose(track_dtype);
+    file_close(track_file);
+    track_file = -1;
+    track_dtype = -1;
+  }
+
   // Open file and write filetype/version -- when MPI is enabled and there is
   // more than one rank, each rank writes its own file
 #ifdef OPENMC_MPI
@@ -82,11 +92,16 @@ void open_track_file()
 
 void close_track_file()
 {
-  H5Tclose(track_dtype);
-  file_close(track_file);
-
-  // Reset number of tracks written
-  n_tracks_written = 0;
+#pragma omp critical(TrackFile)
+  {
+    if (track_file >= 0) {
+      H5Tclose(track_dtype);
+      file_close(track_file);
+      track_file = -1;
+      track_dtype = -1;
+      n_tracks_written = 0;
+    }
+  }
 }
 
 bool check_track_criteria(const Particle& p)
@@ -128,29 +143,40 @@ void finalize_particle_track(Particle& p)
   }
   offsets.push_back(offset);
 
-#pragma omp critical(FinalizeParticleTrack)
+  // Thread safe by ensuring opening and closing of HDF5 file is not
+  // simultaneous
+#pragma omp critical(TrackFile)
   {
-    // Create name for dataset
-    std::string dset_name = fmt::format("track_{}_{}_{}",
-      simulation::current_batch, simulation::current_gen, p.id());
-
-    // Write array of TrackState to file
-    hsize_t dims[] {static_cast<hsize_t>(tracks.size())};
-    hid_t dspace = H5Screate_simple(1, dims, nullptr);
-    hid_t dset = H5Dcreate(track_file, dset_name.c_str(), track_dtype, dspace,
-      H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-    H5Dwrite(dset, track_dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, tracks.data());
-
-    // Write attributes
-    write_attribute(dset, "n_particles", p.tracks().size());
-    write_attribute(dset, "offsets", offsets);
-    write_attribute(dset, "particles", particles);
-
-    // Free resources
-    H5Dclose(dset);
-    H5Sclose(dspace);
-  }
+    // Guard against writing to a closed file
+    if (track_file >= 0) {
+      static int track_index = 0;
+
+      // Create name for dataset
+      std::string dset_name = fmt::format("track_{}_{}_{}",
+        simulation::current_batch, simulation::current_gen, p.id());
+
+      // In case ID already exists (no duplicates in HDF5)
+      if (H5Lexists(track_file, dset_name.c_str(), H5P_DEFAULT) > 0) {
+        dset_name = fmt::format("{}_{}", dset_name, track_index++);
+      }
 
+      // Write array of TrackState to file
+      hsize_t dims[] {static_cast<hsize_t>(tracks.size())};
+      hid_t dspace = H5Screate_simple(1, dims, nullptr);
+      hid_t dset = H5Dcreate(track_file, dset_name.c_str(), track_dtype, dspace,
+        H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+      H5Dwrite(dset, track_dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, tracks.data());
+
+      // Write attributes
+      write_attribute(dset, "n_particles", p.tracks().size());
+      write_attribute(dset, "offsets", offsets);
+      write_attribute(dset, "particles", particles);
+
+      // Free resources
+      H5Dclose(dset);
+      H5Sclose(dspace);
+    }
+  }
   // Clear particle tracks
   p.tracks().clear();
 }