diff --git a/include/openmc/particle_restart.h b/include/openmc/particle_restart.h index 24ea237a42b..39d7ced0fa0 100644 --- a/include/openmc/particle_restart.h +++ b/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 diff --git a/include/openmc/track_output.h b/include/openmc/track_output.h index 2380fe44055..f5be1bd9a8f 100644 --- a/include/openmc/track_output.h +++ b/include/openmc/track_output.h @@ -2,6 +2,7 @@ #define OPENMC_TRACK_OUTPUT_H #include "openmc/particle.h" +#include 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 diff --git a/src/particle.cpp b/src/particle.cpp index 779ae18da9e..e63c32eb660 100644 --- a/src/particle.cpp +++ b/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."); + } } } diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index c226d51ec2a..0e2e3e577be 100644 --- a/src/particle_restart.cpp +++ b/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 diff --git a/src/simulation.cpp b/src/simulation.cpp index 89aa9ca0ef3..16c2d20149d 100644 --- a/src/simulation.cpp +++ b/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; diff --git a/src/track_output.cpp b/src/track_output.cpp index 02e09223645..fdb9129d7df 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -10,7 +10,7 @@ #include "openmc/tensor.h" #include -#include +// #include #include // for size_t #include @@ -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(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(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(); }