From 5549b58e1fbebf0f584341d039fede2a44ae0889 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Tue, 16 Jan 2024 10:35:57 -0600 Subject: [PATCH] Geometron (#2744) --- include/openmc/cell.h | 14 +- include/openmc/dagmc.h | 8 +- include/openmc/geometry.h | 16 +- include/openmc/particle.h | 14 +- include/openmc/particle_data.h | 440 ++++++++++++++++++++------------- include/openmc/plot.h | 8 +- include/openmc/surface.h | 7 +- include/openmc/universe.h | 3 +- src/cell.cpp | 7 +- src/dagmc.cpp | 12 +- src/geometry.cpp | 67 ++--- src/particle.cpp | 18 +- src/particle_data.cpp | 23 +- src/plot.cpp | 43 ++-- src/surface.cpp | 4 +- src/universe.cpp | 3 +- 16 files changed, 403 insertions(+), 284 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 39acaf06740..c405f536b5d 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -40,6 +40,7 @@ constexpr int32_t OP_UNION {std::numeric_limits::max() - 4}; //============================================================================== class Cell; +class GeometryState; class ParentCell; class CellInstance; class Universe; @@ -82,7 +83,7 @@ class Region { //! Find the oncoming boundary of this cell. std::pair distance( - Position r, Direction u, int32_t on_surface, Particle* p) const; + Position r, Direction u, int32_t on_surface) const; //! Get the BoundingBox for this cell. BoundingBox bounding_box(int32_t cell_id) const; @@ -183,7 +184,7 @@ class Cell { //! Find the oncoming boundary of this cell. virtual std::pair distance( - Position r, Direction u, int32_t on_surface, Particle* p) const = 0; + Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0; //! Write all information needed to reconstruct the cell to an HDF5 group. //! \param group_id An HDF5 group id. @@ -260,7 +261,8 @@ class Cell { //! \param[in] instance of the cell to find parent cells for //! \param[in] p particle used to do a fast search for parent cells //! \return parent cells - vector find_parent_cells(int32_t instance, Particle& p) const; + vector find_parent_cells( + int32_t instance, GeometryState& p) const; //! Determine the path to this cell instance in the geometry hierarchy //! \param[in] instance of the cell to find parent cells for @@ -332,10 +334,10 @@ class CSGCell : public Cell { // Methods vector surfaces() const override { return region_.surfaces(); } - std::pair distance( - Position r, Direction u, int32_t on_surface, Particle* p) const override + std::pair distance(Position r, Direction u, + int32_t on_surface, GeometryState* p) const override { - return region_.distance(r, u, on_surface, p); + return region_.distance(r, u, on_surface); } bool contains(Position r, Direction u, int32_t on_surface) const override diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 5ed426e5afb..0b23e567abc 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -42,7 +42,7 @@ class DAGSurface : public Surface { double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; - Direction reflect(Position r, Direction u, Particle* p) const override; + Direction reflect(Position r, Direction u, GeometryState* p) const override; inline void to_hdf5_inner(hid_t group_id) const override {}; @@ -63,8 +63,8 @@ class DAGCell : public Cell { bool contains(Position r, Direction u, int32_t on_surface) const override; - std::pair distance( - Position r, Direction u, int32_t on_surface, Particle* p) const override; + std::pair distance(Position r, Direction u, + int32_t on_surface, GeometryState* p) const override; BoundingBox bounding_box() const override; @@ -143,7 +143,7 @@ class DAGUniverse : public Universe { //! string of the ID ranges for entities of dimension \p dim std::string dagmc_ids_for_dim(int dim) const; - bool find_cell(Particle& p) const override; + bool find_cell(GeometryState& p) const override; void to_hdf5(hid_t universes_group) const override; diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 001e58c4cc3..107cc7d1f3e 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -11,7 +11,7 @@ namespace openmc { class BoundaryInfo; -class Particle; +class GeometryState; //============================================================================== // Global variables @@ -39,7 +39,7 @@ inline bool coincident(double d1, double d2) //! Check for overlapping cells at a particle's position. //============================================================================== -bool check_cell_overlap(Particle& p, bool error = true); +bool check_cell_overlap(GeometryState& p, bool error = true); //============================================================================== //! Get the cell instance for a particle at the specified universe level @@ -50,7 +50,7 @@ bool check_cell_overlap(Particle& p, bool error = true); //! should be computed. \return The instance of the cell at the specified level. //============================================================================== -int cell_instance_at_level(const Particle& p, int level); +int cell_instance_at_level(const GeometryState& p, int level); //============================================================================== //! Locate a particle in the geometry tree and set its geometry data fields. @@ -60,20 +60,22 @@ int cell_instance_at_level(const Particle& p, int level); //! \return True if the particle's location could be found and ascribed to a //! valid geometry coordinate stack. //============================================================================== -bool exhaustive_find_cell(Particle& p); -bool neighbor_list_find_cell(Particle& p); // Only usable on surface crossings +bool exhaustive_find_cell(GeometryState& p, bool verbose = false); +bool neighbor_list_find_cell( + GeometryState& p, bool verbose = false); // Only usable on surface crossings //============================================================================== //! Move a particle into a new lattice tile. //============================================================================== -void cross_lattice(Particle& p, const BoundaryInfo& boundary); +void cross_lattice( + GeometryState& p, const BoundaryInfo& boundary, bool verbose = false); //============================================================================== //! Find the next boundary a particle will intersect. //============================================================================== -BoundaryInfo distance_to_boundary(Particle& p); +BoundaryInfo distance_to_boundary(GeometryState& p); } // namespace openmc diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 80a561f4d05..e423880f87c 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -5,7 +5,6 @@ //! \brief Particle type #include -#include #include #include "openmc/constants.h" @@ -103,17 +102,8 @@ class Particle : public ParticleData { //! mark a particle as lost and create a particle restart file //! \param message A warning message to display - void mark_as_lost(const char* message); - - void mark_as_lost(const std::string& message) - { - mark_as_lost(message.c_str()); - } - - void mark_as_lost(const std::stringstream& message) - { - mark_as_lost(message.str()); - } + virtual void mark_as_lost(const char* message) override; + using GeometryState::mark_as_lost; //! create a particle restart HDF5 file void write_restart() const; diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 37eb05bdb04..cb68fc45740 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -194,6 +194,157 @@ struct BoundaryInfo { lattice_translation {}; //!< which way lattice indices will change }; +/* + * Contains all geometry state information for a particle. + */ +class GeometryState { +public: + GeometryState(); + + /* + * GeometryState does not store any ID info, so give some reasonable behavior + * here. The Particle class redefines this. This is only here for the error + * reporting behavior that occurs in geometry.cpp. The explanation for + * mark_as_lost is the same. + */ + virtual void mark_as_lost(const char* message); + void mark_as_lost(const std::string& message); + void mark_as_lost(const std::stringstream& message); + + // resets all coordinate levels for the particle + void clear() + { + for (auto& level : coord_) + level.reset(); + n_coord_ = 1; + } + + // Initialize all internal state from position and direction + void init_from_r_u(Position r_a, Direction u_a) + { + clear(); + surface() = 0; + material() = C_NONE; + r() = r_a; + u() = u_a; + r_last_current() = r_a; + r_last() = r_a; + u_last() = u_a; + } + + // Unique ID. This is not geometric info, but the + // error reporting in geometry.cpp requires this. + // We could save this to implement it in Particle, + // but that would require virtuals. + int64_t& id() { return id_; } + const int64_t& id() const { return id_; } + + // Number of current coordinate levels + int& n_coord() { return n_coord_; } + const int& n_coord() const { return n_coord_; } + + // Offset for distributed properties + int& cell_instance() { return cell_instance_; } + const int& cell_instance() const { return cell_instance_; } + + // Coordinates for all nesting levels + LocalCoord& coord(int i) { return coord_[i]; } + const LocalCoord& coord(int i) const { return coord_[i]; } + const vector& coord() const { return coord_; } + + // Innermost universe nesting coordinates + LocalCoord& lowest_coord() { return coord_[n_coord_ - 1]; } + const LocalCoord& lowest_coord() const { return coord_[n_coord_ - 1]; } + + // Last coordinates on all nesting levels, before crossing a surface + int& n_coord_last() { return n_coord_last_; } + const int& n_coord_last() const { return n_coord_last_; } + int& cell_last(int i) { return cell_last_[i]; } + const int& cell_last(int i) const { return cell_last_[i]; } + + // Coordinates of last collision or reflective/periodic surface + // crossing for current tallies + Position& r_last_current() { return r_last_current_; } + const Position& r_last_current() const { return r_last_current_; } + + // Previous direction and spatial coordinates before a collision + Position& r_last() { return r_last_; } + const Position& r_last() const { return r_last_; } + Position& u_last() { return u_last_; } + const Position& u_last() const { return u_last_; } + + // Accessors for position in global coordinates + Position& r() { return coord_[0].r; } + const Position& r() const { return coord_[0].r; } + + // Accessors for position in local coordinates + Position& r_local() { return coord_[n_coord_ - 1].r; } + const Position& r_local() const { return coord_[n_coord_ - 1].r; } + + // Accessors for direction in global coordinates + Direction& u() { return coord_[0].u; } + const Direction& u() const { return coord_[0].u; } + + // Accessors for direction in local coordinates + Direction& u_local() { return coord_[n_coord_ - 1].u; } + const Direction& u_local() const { return coord_[n_coord_ - 1].u; } + + // Surface that the particle is on + int& surface() { return surface_; } + const int& surface() const { return surface_; } + + // Boundary information + BoundaryInfo& boundary() { return boundary_; } + +#ifdef DAGMC + // DagMC state variables + moab::DagMC::RayHistory& history() { return history_; } + Direction& last_dir() { return last_dir_; } +#endif + + // material of current and last cell + int& material() { return material_; } + const int& material() const { return material_; } + int& material_last() { return material_last_; } + const int& material_last() const { return material_last_; } + + // temperature of current and last cell + double& sqrtkT() { return sqrtkT_; } + const double& sqrtkT() const { return sqrtkT_; } + double& sqrtkT_last() { return sqrtkT_last_; } + +private: + int64_t id_ {-1}; //!< Unique ID + + int n_coord_ {1}; //!< number of current coordinate levels + int cell_instance_; //!< offset for distributed properties + vector coord_; //!< coordinates for all levels + + int n_coord_last_ {1}; //!< number of current coordinates + vector cell_last_; //!< coordinates for all levels + + Position r_last_current_; //!< coordinates of the last collision or + //!< reflective/periodic surface crossing for + //!< current tallies + Position r_last_; //!< previous coordinates + Direction u_last_; //!< previous direction coordinates + + int surface_ {0}; //!< index for surface particle is on + + BoundaryInfo boundary_; //!< Info about the next intersection + + int material_ {-1}; //!< index for current material + int material_last_ {-1}; //!< index for last material + + double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV + double sqrtkT_last_ {0.0}; //!< last temperature + +#ifdef DAGMC + moab::DagMC::RayHistory history_; + Direction last_dir_; +#endif +}; + //============================================================================ //! Defines how particle data is laid out in memory //============================================================================ @@ -229,163 +380,112 @@ struct BoundaryInfo { * Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18. * https://doi.org/10.1016/j.anucene.2017.11.032. */ -class ParticleData { -public: - //---------------------------------------------------------------------------- - // Constructors - ParticleData(); - +class ParticleData : public GeometryState { private: //========================================================================== - // Data members (accessor methods are below) + // Data members -- see public: below for descriptions - // Cross section caches - vector neutron_xs_; //!< Microscopic neutron cross sections - vector photon_xs_; //!< Microscopic photon cross sections - MacroXS macro_xs_; //!< Macroscopic cross sections - CacheDataMG mg_xs_cache_; //!< Multigroup XS cache + vector neutron_xs_; + vector photon_xs_; + MacroXS macro_xs_; + CacheDataMG mg_xs_cache_; - int64_t id_; //!< Unique ID - ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.) + ParticleType type_ {ParticleType::neutron}; - int n_coord_ {1}; //!< number of current coordinate levels - int cell_instance_; //!< offset for distributed properties - vector coord_; //!< coordinates for all levels + double E_; + double E_last_; + int g_ {0}; + int g_last_; - // Particle coordinates before crossing a surface - int n_coord_last_ {1}; //!< number of current coordinates - vector cell_last_; //!< coordinates for all levels + double wgt_ {1.0}; + double mu_; + double time_ {0.0}; + double time_last_ {0.0}; + double wgt_last_ {1.0}; - // Energy data - double E_; //!< post-collision energy in eV - double E_last_; //!< pre-collision energy in eV - int g_ {C_NONE}; //!< post-collision energy group (MG only) - int g_last_; //!< pre-collision energy group (MG only) + bool fission_ {false}; + TallyEvent event_; + int event_nuclide_; + int event_mt_; + int delayed_group_ {0}; - // Other physical data - double wgt_ {1.0}; //!< particle weight - double mu_; //!< angle of scatter - double time_ {0.0}; //!< time in [s] - double time_last_ {0.0}; //!< previous time in [s] + int n_bank_ {0}; + int n_bank_second_ {0}; + double wgt_bank_ {0.0}; + int n_delayed_bank_[MAX_DELAYED_GROUPS]; - // Other physical data - Position r_last_current_; //!< coordinates of the last collision or - //!< reflective/periodic surface crossing for - //!< current tallies - Position r_last_; //!< previous coordinates - Direction u_last_; //!< previous direction coordinates - double wgt_last_ {1.0}; //!< pre-collision particle weight - - // What event took place - bool fission_ {false}; //!< did particle cause implicit fission - TallyEvent event_; //!< scatter, absorption - int event_nuclide_; //!< index in nuclides array - int event_mt_; //!< reaction MT - int delayed_group_ {0}; //!< delayed group - - // Post-collision physical data - int n_bank_ {0}; //!< number of fission sites banked - int n_bank_second_ {0}; //!< number of secondary particles banked - double wgt_bank_ {0.0}; //!< weight of fission sites banked - int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission - //!< sites banked - - // Indices for various arrays - int surface_ {0}; //!< index for surface particle is on - int cell_born_ {-1}; //!< index for cell particle was born in - int material_ {-1}; //!< index for current material - int material_last_ {-1}; //!< index for last material + int cell_born_ {-1}; - // Boundary information - BoundaryInfo boundary_; + int n_collision_ {0}; - // Temperature of current cell - double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV - double sqrtkT_last_ {0.0}; //!< last temperature - - // Statistical data - int n_collision_ {0}; //!< number of collisions - - // Track output bool write_track_ {false}; - // Current PRNG state - uint64_t seeds_[N_STREAMS]; // current seeds - int stream_; // current RNG stream + uint64_t seeds_[N_STREAMS]; + int stream_; - // Secondary particle bank vector secondary_bank_; - int64_t current_work_; // current work index + int64_t current_work_; - vector flux_derivs_; // for derivatives for this particle + vector flux_derivs_; - vector filter_matches_; // tally filter matches + vector filter_matches_; - vector tracks_; // tracks for outputting to file + vector tracks_; - vector nu_bank_; // bank of most recently fissioned particles + vector nu_bank_; - vector pht_storage_; // interim pulse-height results + vector pht_storage_; - // Global tally accumulators double keff_tally_absorption_ {0.0}; double keff_tally_collision_ {0.0}; double keff_tally_tracklength_ {0.0}; double keff_tally_leakage_ {0.0}; - bool trace_ {false}; //!< flag to show debug information + bool trace_ {false}; - double collision_distance_; // distance to particle's next closest collision + double collision_distance_; - int n_event_ {0}; // number of events executed in this particle's history + int n_event_ {0}; - // Weight window information - int n_split_ {0}; // Number of times this particle has been split - double ww_factor_ { - 0.0}; // Particle-specific factor for on-the-fly weight window adjustment + int n_split_ {0}; + double ww_factor_ {0.0}; -// DagMC state variables -#ifdef DAGMC - moab::DagMC::RayHistory history_; - Direction last_dir_; -#endif - - int64_t n_progeny_ {0}; // Number of progeny produced by this particle + int64_t n_progeny_ {0}; public: + //---------------------------------------------------------------------------- + // Constructors + ParticleData(); + //========================================================================== // Methods and accessors - NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; } + // Cross section caches + NuclideMicroXS& neutron_xs(int i) + { + return neutron_xs_[i]; + } // Microscopic neutron cross sections const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; } + + // Microscopic photon cross sections ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; } + + // Macroscopic cross sections MacroXS& macro_xs() { return macro_xs_; } const MacroXS& macro_xs() const { return macro_xs_; } + + // Multigroup macroscopic cross sections CacheDataMG& mg_xs_cache() { return mg_xs_cache_; } const CacheDataMG& mg_xs_cache() const { return mg_xs_cache_; } - int64_t& id() { return id_; } - const int64_t& id() const { return id_; } + // Particle type (n, p, e, gamma, etc) ParticleType& type() { return type_; } const ParticleType& type() const { return type_; } - int& n_coord() { return n_coord_; } - const int& n_coord() const { return n_coord_; } - int& cell_instance() { return cell_instance_; } - const int& cell_instance() const { return cell_instance_; } - LocalCoord& coord(int i) { return coord_[i]; } - const LocalCoord& coord(int i) const { return coord_[i]; } - const vector& coord() const { return coord_; } - - LocalCoord& lowest_coord() { return coord_[n_coord_ - 1]; } - const LocalCoord& lowest_coord() const { return coord_[n_coord_ - 1]; } - - int& n_coord_last() { return n_coord_last_; } - const int& n_coord_last() const { return n_coord_last_; } - int& cell_last(int i) { return cell_last_[i]; } - const int& cell_last(int i) const { return cell_last_[i]; } - + // Current particle energy, energy before collision, + // and corresponding multigroup group indices. Energy + // units are eV. double& E() { return E_; } const double& E() const { return E_; } double& E_last() { return E_last_; } @@ -395,113 +495,119 @@ class ParticleData { int& g_last() { return g_last_; } const int& g_last() const { return g_last_; } + // Statistic weight of particle. Setting to zero + // indicates that the particle is dead. double& wgt() { return wgt_; } double wgt() const { return wgt_; } + double& wgt_last() { return wgt_last_; } + const double& wgt_last() const { return wgt_last_; } + bool alive() const { return wgt_ != 0.0; } + + // Polar scattering angle after a collision double& mu() { return mu_; } const double& mu() const { return mu_; } + + // Tracks the time of a particle as it traverses the problem. + // Units are seconds. double& time() { return time_; } const double& time() const { return time_; } double& time_last() { return time_last_; } const double& time_last() const { return time_last_; } - bool alive() const { return wgt_ != 0.0; } - Position& r_last_current() { return r_last_current_; } - const Position& r_last_current() const { return r_last_current_; } - Position& r_last() { return r_last_; } - const Position& r_last() const { return r_last_; } - Position& u_last() { return u_last_; } - const Position& u_last() const { return u_last_; } - double& wgt_last() { return wgt_last_; } - const double& wgt_last() const { return wgt_last_; } - - bool& fission() { return fission_; } + // What event took place, described in greater detail below TallyEvent& event() { return event_; } const TallyEvent& event() const { return event_; } - int& event_nuclide() { return event_nuclide_; } + bool& fission() { return fission_; } // true if implicit fission + int& event_nuclide() { return event_nuclide_; } // index of collision nuclide const int& event_nuclide() const { return event_nuclide_; } - int& event_mt() { return event_mt_; } - int& delayed_group() { return delayed_group_; } + int& event_mt() { return event_mt_; } // MT number of collision + int& delayed_group() { return delayed_group_; } // delayed group - int& n_bank() { return n_bank_; } - int& n_bank_second() { return n_bank_second_; } - double& wgt_bank() { return wgt_bank_; } - int* n_delayed_bank() { return n_delayed_bank_; } - int& n_delayed_bank(int i) { return n_delayed_bank_[i]; } + // Post-collision data + int& n_bank() { return n_bank_; } // number of banked fission sites + int& n_bank_second() + { + return n_bank_second_; + } // number of secondaries banked + double& wgt_bank() { return wgt_bank_; } // weight of banked fission sites + int* n_delayed_bank() + { + return n_delayed_bank_; + } // number of delayed fission sites + int& n_delayed_bank(int i) + { + return n_delayed_bank_[i]; + } // number of delayed fission sites - int& surface() { return surface_; } - const int& surface() const { return surface_; } + // Index of cell particle is born in int& cell_born() { return cell_born_; } const int& cell_born() const { return cell_born_; } - int& material() { return material_; } - const int& material() const { return material_; } - int& material_last() { return material_last_; } - const int& material_last() const { return material_last_; } - - BoundaryInfo& boundary() { return boundary_; } - - double& sqrtkT() { return sqrtkT_; } - const double& sqrtkT() const { return sqrtkT_; } - double& sqrtkT_last() { return sqrtkT_last_; } + // index of the current and last material + // Total number of collisions suffered by particle int& n_collision() { return n_collision_; } const int& n_collision() const { return n_collision_; } + // whether this track is to be written bool& write_track() { return write_track_; } + + // RNG state uint64_t& seeds(int i) { return seeds_[i]; } uint64_t* seeds() { return seeds_; } int& stream() { return stream_; } + // secondary particle bank SourceSite& secondary_bank(int i) { return secondary_bank_[i]; } decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; } + + // Current simulation work index int64_t& current_work() { return current_work_; } const int64_t& current_work() const { return current_work_; } + + // Used in tally derivatives double& flux_derivs(int i) { return flux_derivs_[i]; } const double& flux_derivs(int i) const { return flux_derivs_[i]; } + + // Matches of tallies decltype(filter_matches_)& filter_matches() { return filter_matches_; } FilterMatch& filter_matches(int i) { return filter_matches_[i]; } + + // Tracks to output to file decltype(tracks_)& tracks() { return tracks_; } + + // Bank of recently fissioned particles decltype(nu_bank_)& nu_bank() { return nu_bank_; } NuBank& nu_bank(int i) { return nu_bank_[i]; } + + // Interim pulse height tally storage vector& pht_storage() { return pht_storage_; } + // Global tally accumulators double& keff_tally_absorption() { return keff_tally_absorption_; } double& keff_tally_collision() { return keff_tally_collision_; } double& keff_tally_tracklength() { return keff_tally_tracklength_; } double& keff_tally_leakage() { return keff_tally_leakage_; } + // Shows debug info bool& trace() { return trace_; } + + // Distance to the next collision double& collision_distance() { return collision_distance_; } + + // Number of events particle has undergone int& n_event() { return n_event_; } + // Number of times variance reduction has caused a particle split int n_split() const { return n_split_; } int& n_split() { return n_split_; } + // Particle-specific factor for on-the-fly weight window adjustment double ww_factor() const { return ww_factor_; } double& ww_factor() { return ww_factor_; } -#ifdef DAGMC - moab::DagMC::RayHistory& history() { return history_; } - Direction& last_dir() { return last_dir_; } -#endif - + // Number of progeny produced by this particle int64_t& n_progeny() { return n_progeny_; } - // Accessors for position in global coordinates - Position& r() { return coord_[0].r; } - const Position& r() const { return coord_[0].r; } - - // Accessors for position in local coordinates - Position& r_local() { return coord_[n_coord_ - 1].r; } - const Position& r_local() const { return coord_[n_coord_ - 1].r; } - - // Accessors for direction in global coordinates - Direction& u() { return coord_[0].u; } - const Direction& u() const { return coord_[0].u; } - - // Accessors for direction in local coordinates - Direction& u_local() { return coord_[n_coord_ - 1].u; } - const Direction& u_local() const { return coord_[n_coord_ - 1].u; } - //! Gets the pointer to the particle's current PRN seed uint64_t* current_seed() { return seeds_ + stream_; } const uint64_t* current_seed() const { return seeds_ + stream_; } @@ -513,14 +619,6 @@ class ParticleData { micro.last_E = 0.0; } - //! resets all coordinate levels for the particle - void clear() - { - for (auto& level : coord_) - level.reset(); - n_coord_ = 1; - } - //! Get track information based on particle's current state TrackState get_track_state() const; diff --git a/include/openmc/plot.h b/include/openmc/plot.h index ba82f43e27a..7b66036fb71 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -121,7 +121,7 @@ struct IdData { IdData(size_t h_res, size_t v_res); // Methods - void set_value(size_t y, size_t x, const Particle& p, int level); + void set_value(size_t y, size_t x, const GeometryState& p, int level); void set_overlap(size_t y, size_t x); // Members @@ -133,7 +133,7 @@ struct PropertyData { PropertyData(size_t h_res, size_t v_res); // Methods - void set_value(size_t y, size_t x, const Particle& p, int level); + void set_value(size_t y, size_t x, const GeometryState& p, int level); void set_overlap(size_t y, size_t x); // Members @@ -205,7 +205,7 @@ T SlicePlotBase::get_map() const #pragma omp parallel { - Particle p; + GeometryState p; p.r() = xyz; p.u() = dir; p.coord(0).universe = model::root_universe; @@ -291,7 +291,7 @@ class ProjectionPlot : public PlottableInterface { * find a distance to the boundary in a non-standard surface intersection * check. It's an exhaustive search over surfaces in the top-level universe. */ - static int advance_to_boundary_from_void(Particle& p); + static int advance_to_boundary_from_void(GeometryState& p); /* Checks if a vector of two TrackSegments is equivalent. We define this * to mean not having matching intersection lengths, but rather having diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 3bb84c6e00c..350775123c1 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -105,12 +105,13 @@ class Surface { //! Determine the direction of a ray reflected from the surface. //! \param[in] r The point at which the ray is incident. //! \param[in] u Incident direction of the ray - //! \param[inout] p Pointer to the particle + //! \param[inout] p Pointer to the particle. Only DAGMC uses this. //! \return Outgoing direction of the ray - virtual Direction reflect(Position r, Direction u, Particle* p) const; + virtual Direction reflect( + Position r, Direction u, GeometryState* p = nullptr) const; virtual Direction diffuse_reflect( - Position r, Direction u, uint64_t* seed) const; + Position r, Direction u, uint64_t* seed, GeometryState* p = nullptr) const; //! Evaluate the equation describing the surface. //! diff --git a/include/openmc/universe.h b/include/openmc/universe.h index b98d2d57f64..26f33cb383d 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -9,6 +9,7 @@ namespace openmc { class DAGUniverse; #endif +class GeometryState; class Universe; class UniversePartitioner; @@ -32,7 +33,7 @@ class Universe { //! \param group_id An HDF5 group id. virtual void to_hdf5(hid_t group_id) const; - virtual bool find_cell(Particle& p) const; + virtual bool find_cell(GeometryState& p) const; BoundingBox bounding_box() const; diff --git a/src/cell.cpp b/src/cell.cpp index cba52f1b55f..00083d967b3 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -784,7 +784,7 @@ std::string Region::str() const //============================================================================== std::pair Region::distance( - Position r, Direction u, int32_t on_surface, Particle* p) const + Position r, Direction u, int32_t on_surface) const { double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; @@ -1301,14 +1301,15 @@ vector Cell::find_parent_cells( { // create a temporary particle - Particle dummy_particle {}; + GeometryState dummy_particle {}; dummy_particle.r() = r; dummy_particle.u() = {0., 0., 1.}; return find_parent_cells(instance, dummy_particle); } -vector Cell::find_parent_cells(int32_t instance, Particle& p) const +vector Cell::find_parent_cells( + int32_t instance, GeometryState& p) const { // look up the particle's location exhaustive_find_cell(p); diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 4ffe4e6526a..561fdca81ba 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -404,7 +404,7 @@ int32_t DAGUniverse::implicit_complement_idx() const return cell_idx_offset_ + dagmc_instance_->index_by_handle(ic) - 1; } -bool DAGUniverse::find_cell(Particle& p) const +bool DAGUniverse::find_cell(GeometryState& p) const { // if the particle isn't in any of the other DagMC // cells, place it in the implicit complement @@ -583,9 +583,8 @@ DAGCell::DAGCell(std::shared_ptr dag_ptr, int32_t dag_idx) }; std::pair DAGCell::distance( - Position r, Direction u, int32_t on_surface, Particle* p) const + Position r, Direction u, int32_t on_surface, GeometryState* p) const { - Expects(p); // if we've changed direction or we're not on a surface, // reset the history and update last direction if (u != p->last_dir()) { @@ -635,10 +634,9 @@ std::pair DAGCell::distance( p->material() == MATERIAL_VOID ? "-1 (VOID)" : std::to_string(model::materials[p->material()]->id()); - auto lost_particle_msg = fmt::format( + p->mark_as_lost(fmt::format( "No intersection found with DAGMC cell {}, filled with material {}", id_, - material_id); - p->mark_as_lost(lost_particle_msg); + material_id)); } return {dist, surf_idx}; @@ -723,7 +721,7 @@ Direction DAGSurface::normal(Position r) const return dir; } -Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const +Direction DAGSurface::reflect(Position r, Direction u, GeometryState* p) const { Expects(p); p->history().reset_to_last_intersection(); diff --git a/src/geometry.cpp b/src/geometry.cpp index 47152e17295..a16addd4848 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -32,7 +32,7 @@ vector overlap_check_count; // Non-member functions //============================================================================== -bool check_cell_overlap(Particle& p, bool error) +bool check_cell_overlap(GeometryState& p, bool error) { int n_coord = p.n_coord(); @@ -63,7 +63,7 @@ bool check_cell_overlap(Particle& p, bool error) //============================================================================== -int cell_instance_at_level(const Particle& p, int level) +int cell_instance_at_level(const GeometryState& p, int level) { // throw error if the requested level is too deep for the geometry if (level > model::n_coord_levels) { @@ -99,7 +99,8 @@ int cell_instance_at_level(const Particle& p, int level) //============================================================================== -bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) +bool find_cell_inner( + GeometryState& p, const NeighborList* neighbor_list, bool verbose) { // Find which cell of this universe the particle is in. Use the neighbor list // to shorten the search if one was provided. @@ -156,7 +157,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) i_cell = p.lowest_coord().cell; // Announce the cell that the particle is entering. - if (found && (settings::verbosity >= 10 || p.trace())) { + if (found && verbose) { auto msg = fmt::format(" Entering cell {}", model::cells[i_cell]->id_); write_message(msg, 1); } @@ -243,10 +244,9 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) if (lat.outer_ != NO_OUTER_UNIVERSE) { coord.universe = lat.outer_; } else { - warning(fmt::format("Particle {} is outside lattice {} but the " - "lattice has no defined outer universe.", + p.mark_as_lost(fmt::format( + "Particle {} left lattice {}, but it has no outer definition.", p.id(), lat.id_)); - return false; } } } @@ -259,7 +259,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) //============================================================================== -bool neighbor_list_find_cell(Particle& p) +bool neighbor_list_find_cell(GeometryState& p, bool verbose) { // Reset all the deeper coordinate levels. @@ -274,20 +274,20 @@ bool neighbor_list_find_cell(Particle& p) // Search for the particle in that cell's neighbor list. Return if we // found the particle. - bool found = find_cell_inner(p, &c.neighbors_); + bool found = find_cell_inner(p, &c.neighbors_, verbose); if (found) return found; // The particle could not be found in the neighbor list. Try searching all // cells in this universe, and update the neighbor list if we find a new // neighboring cell. - found = find_cell_inner(p, nullptr); + found = find_cell_inner(p, nullptr, verbose); if (found) c.neighbors_.push_back(p.coord(coord_lvl).cell); return found; } -bool exhaustive_find_cell(Particle& p) +bool exhaustive_find_cell(GeometryState& p, bool verbose) { int i_universe = p.lowest_coord().universe; if (i_universe == C_NONE) { @@ -299,17 +299,17 @@ bool exhaustive_find_cell(Particle& p) for (int i = p.n_coord(); i < model::n_coord_levels; i++) { p.coord(i).reset(); } - return find_cell_inner(p, nullptr); + return find_cell_inner(p, nullptr, verbose); } //============================================================================== -void cross_lattice(Particle& p, const BoundaryInfo& boundary) +void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose) { auto& coord {p.lowest_coord()}; auto& lat {*model::lattices[coord.lattice]}; - if (settings::verbosity >= 10 || p.trace()) { + if (verbose) { write_message( fmt::format(" Crossing lattice {}. Current position ({},{},{}). r={}", lat.id_, coord.lattice_i[0], coord.lattice_i[1], coord.lattice_i[2], @@ -336,10 +336,11 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary) // The particle is outside the lattice. Search for it from the base coords. p.n_coord() = 1; bool found = exhaustive_find_cell(p); - if (!found && p.alive()) { - p.mark_as_lost(fmt::format("Could not locate particle {} after " - "crossing a lattice boundary", - p.id())); + + if (!found) { + p.mark_as_lost(fmt::format("Particle {} could not be located after " + "crossing a boundary of lattice {}", + p.id(), lat.id_)); } } else { @@ -352,10 +353,10 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary) // this case, search for it from the base coords. p.n_coord() = 1; bool found = exhaustive_find_cell(p); - if (!found && p.alive()) { - p.mark_as_lost(fmt::format("Could not locate particle {} after " - "crossing a lattice boundary", - p.id())); + if (!found) { + p.mark_as_lost(fmt::format("Particle {} could not be located after " + "crossing a boundary of lattice {}", + p.id(), lat.id_)); } } } @@ -363,7 +364,7 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary) //============================================================================== -BoundaryInfo distance_to_boundary(Particle& p) +BoundaryInfo distance_to_boundary(GeometryState& p) { BoundaryInfo info; double d_lat = INFINITY; @@ -408,8 +409,9 @@ BoundaryInfo distance_to_boundary(Particle& p) level_lat_trans = lattice_distance.second; if (d_lat < 0) { - p.mark_as_lost(fmt::format( - "Particle {} had a negative distance to a lattice boundary", p.id())); + p.mark_as_lost(fmt::format("Particle {} had a negative distance " + "to a lattice boundary.", + p.id())); } } @@ -463,18 +465,19 @@ BoundaryInfo distance_to_boundary(Particle& p) extern "C" int openmc_find_cell( const double* xyz, int32_t* index, int32_t* instance) { - Particle p; + GeometryState geom_state; - p.r() = Position {xyz}; - p.u() = {0.0, 0.0, 1.0}; + geom_state.r() = Position {xyz}; + geom_state.u() = {0.0, 0.0, 1.0}; - if (!exhaustive_find_cell(p)) { - set_errmsg(fmt::format("Could not find cell at position {}.", p.r())); + if (!exhaustive_find_cell(geom_state)) { + set_errmsg( + fmt::format("Could not find cell at position {}.", geom_state.r())); return OPENMC_E_GEOMETRY; } - *index = p.lowest_coord().cell; - *instance = p.cell_instance(); + *index = geom_state.lowest_coord().cell; + *instance = geom_state.cell_instance(); return 0; } diff --git a/src/particle.cpp b/src/particle.cpp index 8992bc8683b..549de5cff65 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -277,7 +277,9 @@ void Particle::event_cross_surface() boundary().lattice_translation[1] != 0 || boundary().lattice_translation[2] != 0) { // Particle crosses lattice boundary - cross_lattice(*this, boundary()); + + bool verbose = settings::verbosity >= 10 || trace(); + cross_lattice(*this, boundary(), verbose); event() = TallyEvent::LATTICE; } else { // Particle crosses surface @@ -406,7 +408,8 @@ void Particle::event_revive_from_secondary() // have to determine it before the energy of the secondary particle can be // removed from the pulse-height of this cell. if (lowest_coord().cell == C_NONE) { - if (!exhaustive_find_cell(*this)) { + bool verbose = settings::verbosity >= 10 || trace(); + if (!exhaustive_find_cell(*this, verbose)) { mark_as_lost("Could not find the cell containing particle " + std::to_string(id())); return; @@ -556,7 +559,8 @@ void Particle::cross_surface() } #endif - if (neighbor_list_find_cell(*this)) + bool verbose = settings::verbosity >= 10 || trace(); + if (neighbor_list_find_cell(*this, verbose)) return; // ========================================================================== @@ -564,7 +568,7 @@ void Particle::cross_surface() // Remove lower coordinate levels n_coord() = 1; - bool found = exhaustive_find_cell(*this); + bool found = exhaustive_find_cell(*this, verbose); if (settings::run_mode != RunMode::PLOTTING && (!found)) { // If a cell is still not found, there are two possible causes: 1) there is @@ -579,7 +583,7 @@ void Particle::cross_surface() // Couldn't find next cell anywhere! This probably means there is an actual // undefined region in the geometry. - if (!exhaustive_find_cell(*this)) { + if (!exhaustive_find_cell(*this, verbose)) { mark_as_lost("After particle " + std::to_string(id()) + " crossed surface " + std::to_string(surf->id_) + " it could not be located in any cell and it did not leak."); @@ -654,8 +658,8 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u) // (unless we're using a dagmc model, which has exactly one universe) n_coord() = 1; if (surf.geom_type_ != GeometryType::DAG && !neighbor_list_find_cell(*this)) { - this->mark_as_lost("Couldn't find particle after reflecting from surface " + - std::to_string(surf.id_) + "."); + mark_as_lost("Couldn't find particle after reflecting from surface " + + std::to_string(surf.id_) + "."); return; } diff --git a/src/particle_data.cpp b/src/particle_data.cpp index 0fc7202c5bf..d8a5bb9d99c 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -1,6 +1,9 @@ #include "openmc/particle_data.h" +#include + #include "openmc/cell.h" +#include "openmc/error.h" #include "openmc/geometry.h" #include "openmc/material.h" #include "openmc/nuclide.h" @@ -12,6 +15,21 @@ namespace openmc { +void GeometryState::mark_as_lost(const std::string& message) +{ + mark_as_lost(message.c_str()); +} + +void GeometryState::mark_as_lost(const std::stringstream& message) +{ + mark_as_lost(message.str()); +} + +void GeometryState::mark_as_lost(const char* message) +{ + fatal_error(message); +} + void LocalCoord::rotate(const vector& rotation) { r = r.rotate(rotation); @@ -30,13 +48,16 @@ void LocalCoord::reset() rotated = false; } -ParticleData::ParticleData() +GeometryState::GeometryState() { // Create and clear coordinate levels coord_.resize(model::n_coord_levels); cell_last_.resize(model::n_coord_levels); clear(); +} +ParticleData::ParticleData() +{ zero_delayed_bank(); // Every particle starts with no accumulated flux derivative. Note that in diff --git a/src/plot.cpp b/src/plot.cpp index 8d995b0a9a2..bf733cff48f 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -45,7 +45,7 @@ constexpr int32_t OVERLAP {-3}; IdData::IdData(size_t h_res, size_t v_res) : data_({v_res, h_res, 3}, NOT_FOUND) {} -void IdData::set_value(size_t y, size_t x, const Particle& p, int level) +void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level) { // set cell data if (p.n_coord() <= level) { @@ -78,7 +78,8 @@ PropertyData::PropertyData(size_t h_res, size_t v_res) : data_({v_res, h_res, 2}, NOT_FOUND) {} -void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level) +void PropertyData::set_value( + size_t y, size_t x, const GeometryState& p, int level) { Cell* c = model::cells.at(p.lowest_coord().cell).get(); data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN; @@ -1084,7 +1085,7 @@ void ProjectionPlot::set_output_path(pugi::xml_node node) // Advances to the next boundary from outside the geometry // Returns -1 if no intersection found, and the surface index // if an intersection was found. -int ProjectionPlot::advance_to_boundary_from_void(Particle& p) +int ProjectionPlot::advance_to_boundary_from_void(GeometryState& p) { constexpr double scoot = 1e-5; double min_dist = {INFINITY}; @@ -1234,20 +1235,8 @@ void ProjectionPlot::create_output() const const int n_threads = num_threads(); const int tid = thread_num(); - SourceSite s; // Where particle starts from (camera) - s.E = 1; - s.wgt = 1; - s.delayed_group = 0; - s.particle = ParticleType::photon; // just has to be something reasonable - s.parent_id = 1; - s.progeny_id = 2; - s.r = camera_position_; - - Particle p; - s.u.x = 1.0; - s.u.y = 0.0; - s.u.z = 0.0; - p.from_source(&s); + GeometryState p; + p.u() = {1.0, 0.0, 0.0}; int vert = tid; for (int iter = 0; iter <= pixels_[1] / n_threads; iter++) { @@ -1263,6 +1252,10 @@ void ProjectionPlot::create_output() const for (int horiz = 0; horiz < pixels_[0]; ++horiz) { + // Projection mode below decides ray starting conditions + Position init_r; + Direction init_u; + // Generate the starting position/direction of the ray if (orthographic_width_ == 0.0) { // perspective projection double this_phi = @@ -1273,18 +1266,22 @@ void ProjectionPlot::create_output() const camera_local_vec.x = std::cos(this_phi) * std::sin(this_mu); camera_local_vec.y = std::sin(this_phi) * std::sin(this_mu); camera_local_vec.z = std::cos(this_mu); - s.u = camera_local_vec.rotate(camera_to_model); + init_u = camera_local_vec.rotate(camera_to_model); + init_r = camera_position_; } else { // orthographic projection - s.u = looking_direction; + init_u = looking_direction; double x_pix_coord = (static_cast(horiz) - p0 / 2.0) / p0; double y_pix_coord = (static_cast(vert) - p1 / 2.0) / p0; - s.r = camera_position_ + - cam_yaxis * x_pix_coord * orthographic_width_ + - cam_zaxis * y_pix_coord * orthographic_width_; + + init_r = camera_position_; + init_r += cam_yaxis * x_pix_coord * orthographic_width_; + init_r += cam_zaxis * y_pix_coord * orthographic_width_; } - p.from_source(&s); // put particle at camera + // Resets internal geometry state of particle + p.init_from_r_u(init_r, init_u); + bool hitsomething = false; bool intersection_found = true; int loop_counter = 0; diff --git a/src/surface.cpp b/src/surface.cpp index 394a4b36630..12fef070e18 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -129,7 +129,7 @@ bool Surface::sense(Position r, Direction u) const return f > 0.0; } -Direction Surface::reflect(Position r, Direction u, Particle* p) const +Direction Surface::reflect(Position r, Direction u, GeometryState* p) const { // Determine projection of direction onto normal and squared magnitude of // normal. @@ -140,7 +140,7 @@ Direction Surface::reflect(Position r, Direction u, Particle* p) const } Direction Surface::diffuse_reflect( - Position r, Direction u, uint64_t* seed) const + Position r, Direction u, uint64_t* seed, GeometryState* p) const { // Diffuse reflect direction according to the normal. // cosine distribution diff --git a/src/universe.cpp b/src/universe.cpp index 67e1d1354c3..b4ef6b4b265 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -3,6 +3,7 @@ #include #include "openmc/hdf5_interface.h" +#include "openmc/particle.h" namespace openmc { @@ -36,7 +37,7 @@ void Universe::to_hdf5(hid_t universes_group) const close_group(group); } -bool Universe::find_cell(Particle& p) const +bool Universe::find_cell(GeometryState& p) const { const auto& cells { !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};