From 10d32af3feb55d4b65b6bf39b41548f1b796c863 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 17 Aug 2023 16:04:48 -0400 Subject: [PATCH 01/17] factor out geometry info from particle to geometron --- include/openmc/particle_data.h | 378 +++++++++++++++++++-------------- src/particle_data.cpp | 5 +- 2 files changed, 227 insertions(+), 156 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 37eb05bdb04..27edd606549 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -194,6 +194,96 @@ struct BoundaryInfo { lattice_translation {}; //!< which way lattice indices will change }; +/* + * Contains all geometry state information for a particle. + */ +class Geometron { +public: + Geometron(); + +public: + // 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_; } + + // resets all coordinate levels for the particle + void clear() + { + for (auto& level : coord_) + level.reset(); + n_coord_ = 1; + } + +private: + int n_coord_ {1}; + int cell_instance_; + vector coord_; + + int n_coord_last_ {1}; + vector cell_last_; + + Position r_last_current_; + Position r_last_; + Direction u_last_; + + int surface_ {0}; + + BoundaryInfo boundary_; +}; + //============================================================================ //! Defines how particle data is laid out in memory //============================================================================ @@ -229,163 +319,127 @@ 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 Geometron { 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 - - int64_t id_; //!< Unique ID - ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.) - - int n_coord_ {1}; //!< number of current coordinate levels - int cell_instance_; //!< offset for distributed properties - vector coord_; //!< coordinates for all levels - - // Particle coordinates before crossing a surface - int n_coord_last_ {1}; //!< number of current coordinates - vector cell_last_; //!< coordinates for all levels - - // 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) - - // 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] - - // 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 + vector neutron_xs_; + vector photon_xs_; + MacroXS macro_xs_; + CacheDataMG mg_xs_cache_; - // Boundary information - BoundaryInfo boundary_; + int64_t id_; + ParticleType type_ {ParticleType::neutron}; + + double E_; + double E_last_; + int g_ {0}; + int g_last_; - // Temperature of current cell - double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV - double sqrtkT_last_ {0.0}; //!< last temperature + double wgt_ {1.0}; + double mu_; + double time_ {0.0}; + double time_last_ {0.0}; + double wgt_last_ {1.0}; - // Statistical data - int n_collision_ {0}; //!< number of collisions + bool fission_ {false}; + TallyEvent event_; + int event_nuclide_; + int event_mt_; + int delayed_group_ {0}; + + int n_bank_ {0}; + int n_bank_second_ {0}; + double wgt_bank_ {0.0}; + int n_delayed_bank_[MAX_DELAYED_GROUPS]; + + int cell_born_ {-1}; + int material_ {-1}; + int material_last_ {-1}; + + double sqrtkT_ {-1.0}; + double sqrtkT_last_ {0.0}; + + int n_collision_ {0}; - // 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_; } + // Unique ID 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 +449,135 @@ 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_; } + + // index of the current and last material 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_; } - + // temperature of current and last cell double& sqrtkT() { return sqrtkT_; } const double& sqrtkT() const { return sqrtkT_; } double& sqrtkT_last() { return sqrtkT_last_; } + // 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 + // DagMC state variables 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 +589,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/src/particle_data.cpp b/src/particle_data.cpp index 0fc7202c5bf..b80b09e8599 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -30,13 +30,16 @@ void LocalCoord::reset() rotated = false; } -ParticleData::ParticleData() +Geometron::Geometron() { // 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 From f3560ea6b3fe5aefe39eb203ea59079c12b48d2f Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Fri, 18 Aug 2023 14:03:15 -0400 Subject: [PATCH 02/17] factor particle geometry state out into Geometron --- include/openmc/cell.h | 11 +++---- include/openmc/dagmc.h | 6 ++-- include/openmc/geometry.h | 53 +++++++++++++++++++++++++++++----- include/openmc/particle_data.h | 53 +++++++++++++++++----------------- include/openmc/surface.h | 2 +- include/openmc/universe.h | 3 +- src/boundary_condition.cpp | 2 +- src/cell.cpp | 6 ++-- src/dagmc.cpp | 14 ++------- src/geometry.cpp | 50 +++++++++++++++----------------- src/particle.cpp | 47 ++++++++++++++++++++++++++---- src/surface.cpp | 2 +- src/universe.cpp | 3 +- 13 files changed, 159 insertions(+), 93 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 39acaf06740..00eae4d520a 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 Geometron; 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, Geometron* 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,7 @@ 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, Geometron& p) const; //! Determine the path to this cell instance in the geometry hierarchy //! \param[in] instance of the cell to find parent cells for @@ -333,9 +334,9 @@ class CSGCell : public Cell { vector surfaces() const override { return region_.surfaces(); } std::pair distance( - Position r, Direction u, int32_t on_surface, Particle* p) const override + Position r, Direction u, int32_t on_surface, Geometron* 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..7aee45648f4 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) const override; inline void to_hdf5_inner(hid_t group_id) const override {}; @@ -64,7 +64,7 @@ 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; + Position r, Direction u, int32_t on_surface, Geometron* 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(Geometron& 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..f44e5c5865a 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -11,7 +11,7 @@ namespace openmc { class BoundaryInfo; -class Particle; +class Geometron; //============================================================================== // 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(Geometron& 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 Geometron& p, int level); //============================================================================== //! Locate a particle in the geometry tree and set its geometry data fields. @@ -60,20 +60,59 @@ 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(Geometron& p, bool verbose = false); +bool neighbor_list_find_cell( + Geometron& 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( + Geometron& 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(Geometron& p); + +/* Geometry routines can potentially throw this type of exception + */ +class ParticleLost : public std::exception { +public: + enum class Reason { + negative_lattice_distance, // should not have negative distance to lattice + bad_boundary_crossing, // could not locate after crossing a boundary + no_universe_outside_lattice, // undefined behavior. All space should be + // defined. + no_dagmc_intersection + } reason; + + // Extra data to be passed up for exception handling or for detailed + // error reporting. For instance no_dagmc_intersection will give the + // DAGMC cell ID. + int id; + + ParticleLost(Reason reason_a, int id_a = 0) : reason(reason_a), id(id_a) {} + + // Handles uncaught exception messages, e.g. Geometron with no ID. + // These will not show up in the usual particle tracking, as more + // detailed error reporting is provided in the Particle class. + const char* what() + { + switch (reason) { + case Reason::negative_lattice_distance: + return "Negative distance to a lattice"; + case Reason::bad_boundary_crossing: + return "Crossed a boundary and lost particle"; + case Reason::no_universe_outside_lattice: + return "Outside lattice but no outer region defined"; + case Reason::no_dagmc_intersection: + return "No intersection found with DAGMC cell"; + } + } +}; } // namespace openmc diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 27edd606549..28a0be91a36 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -267,6 +267,22 @@ class Geometron { n_coord_ = 1; } +#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_; } + + // temperature of current and last cell + double& sqrtkT() { return sqrtkT_; } + const double& sqrtkT() const { return sqrtkT_; } + double& sqrtkT_last() { return sqrtkT_last_; } + private: int n_coord_ {1}; int cell_instance_; @@ -282,6 +298,17 @@ class Geometron { int surface_ {0}; BoundaryInfo boundary_; + + int material_ {-1}; + int material_last_ {-1}; + + double sqrtkT_ {-1.0}; + double sqrtkT_last_ {0.0}; + +#ifdef DAGMC + moab::DagMC::RayHistory history_; + Direction last_dir_; +#endif }; //============================================================================ @@ -355,11 +382,6 @@ class ParticleData : public Geometron { int n_delayed_bank_[MAX_DELAYED_GROUPS]; int cell_born_ {-1}; - int material_ {-1}; - int material_last_ {-1}; - - double sqrtkT_ {-1.0}; - double sqrtkT_last_ {0.0}; int n_collision_ {0}; @@ -396,11 +418,6 @@ class ParticleData : public Geometron { int n_split_ {0}; double ww_factor_ {0.0}; -#ifdef DAGMC - moab::DagMC::RayHistory history_; - Direction last_dir_; -#endif - int64_t n_progeny_ {0}; public: @@ -498,16 +515,6 @@ class ParticleData : public Geometron { const int& cell_born() const { return cell_born_; } // index of the current and last material - 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_; } - // Total number of collisions suffered by particle int& n_collision() { return n_collision_; } const int& n_collision() const { return n_collision_; } @@ -569,12 +576,6 @@ class ParticleData : public Geometron { double ww_factor() const { return ww_factor_; } double& ww_factor() { return ww_factor_; } -#ifdef DAGMC - // DagMC state variables - 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_; } diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 3bb84c6e00c..6620303edd4 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -107,7 +107,7 @@ class Surface { //! \param[in] u Incident direction of the ray //! \param[inout] p Pointer to the particle //! \return Outgoing direction of the ray - virtual Direction reflect(Position r, Direction u, Particle* p) const; + virtual Direction reflect(Position r, Direction u) const; virtual Direction diffuse_reflect( Position r, Direction u, uint64_t* seed) const; diff --git a/include/openmc/universe.h b/include/openmc/universe.h index b98d2d57f64..bd872eb0535 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -9,6 +9,7 @@ namespace openmc { class DAGUniverse; #endif +class Geometron; 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(Geometron& p) const; BoundingBox bounding_box() const; diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 42a8e27a164..193c842128a 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -25,7 +25,7 @@ void VacuumBC::handle_particle(Particle& p, const Surface& surf) const void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const { - Direction u = surf.reflect(p.r(), p.u(), &p); + Direction u = surf.reflect(p.r(), p.u()); u /= u.norm(); // Handle the effects of the surface albedo on the particle's weight. diff --git a/src/cell.cpp b/src/cell.cpp index c96abf39bef..c0999214677 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()}; @@ -1279,14 +1279,14 @@ vector Cell::find_parent_cells( { // create a temporary particle - Particle dummy_particle {}; + Geometron 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, Geometron& p) const { // look up the particle's location exhaustive_find_cell(p); diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 4ffe4e6526a..002a6b3c83c 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(Geometron& 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, Geometron* 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()) { @@ -631,14 +630,7 @@ std::pair DAGCell::distance( // isn't found in a volume that is not the implicit complement. In the case // that the DAGMC model is the root universe of the geometry, even a missing // intersection in the implicit complement should trigger this condition. - std::string material_id = - p->material() == MATERIAL_VOID - ? "-1 (VOID)" - : std::to_string(model::materials[p->material()]->id()); - auto lost_particle_msg = fmt::format( - "No intersection found with DAGMC cell {}, filled with material {}", id_, - material_id); - p->mark_as_lost(lost_particle_msg); + throw ParticleLost(ParticleLost::Reason::no_dagmc_intersection, id_); } return {dist, surf_idx}; diff --git a/src/geometry.cpp b/src/geometry.cpp index 47152e17295..2c67b0b3350 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(Geometron& 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 Geometron& 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( + Geometron& 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,8 @@ 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.id(), lat.id_)); - return false; + throw ParticleLost( + ParticleLost::Reason::no_universe_outside_lattice, lat.id_); } } } @@ -259,7 +258,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) //============================================================================== -bool neighbor_list_find_cell(Particle& p) +bool neighbor_list_find_cell(Geometron& p, bool verbose) { // Reset all the deeper coordinate levels. @@ -274,20 +273,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(Geometron& p, bool verbose) { int i_universe = p.lowest_coord().universe; if (i_universe == C_NONE) { @@ -299,17 +298,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(Geometron& 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 +335,9 @@ 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) { + throw ParticleLost(ParticleLost::Reason::bad_boundary_crossing, lat.id_); } } else { @@ -352,10 +350,9 @@ 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) { + throw ParticleLost( + ParticleLost::Reason::bad_boundary_crossing, lat.id_); } } } @@ -363,7 +360,7 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary) //============================================================================== -BoundaryInfo distance_to_boundary(Particle& p) +BoundaryInfo distance_to_boundary(Geometron& p) { BoundaryInfo info; double d_lat = INFINITY; @@ -408,8 +405,7 @@ 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())); + throw ParticleLost(ParticleLost::Reason::negative_lattice_distance); } } @@ -463,7 +459,7 @@ BoundaryInfo distance_to_boundary(Particle& p) extern "C" int openmc_find_cell( const double* xyz, int32_t* index, int32_t* instance) { - Particle p; + Geometron p; p.r() = Position {xyz}; p.u() = {0.0, 0.0, 1.0}; diff --git a/src/particle.cpp b/src/particle.cpp index 8992bc8683b..c92e7b31d19 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -205,7 +205,38 @@ void Particle::event_calculate_xs() void Particle::event_advance() { // Find the distance to the nearest boundary - boundary() = distance_to_boundary(*this); + try { + boundary() = distance_to_boundary(*this); + } catch (ParticleLost msg) { + // Lost particles do not mean we should terminate the code. + // Instead we accept that they're lost and only allow a certain + // amount to geometrically fail. + switch (msg.reason) { + case ParticleLost::Reason::negative_lattice_distance: + mark_as_lost(fmt::format( + "Particle {} had a negative distance to a lattice boundary", id())); + break; + case ParticleLost::Reason::bad_boundary_crossing: + mark_as_lost(fmt::format("Particle {} could not be located after " + "crossing a boundary of lattice {}", + id(), msg.id)); + break; + case ParticleLost::Reason::no_universe_outside_lattice: + mark_as_lost(fmt::format( + "Particle {} left lattice {}, but it has no outer definition.", id(), + msg.id)); + break; + case ParticleLost::Reason::no_dagmc_intersection: + std::string material_id = + material() == MATERIAL_VOID + ? "-1 (VOID)" + : std::to_string(model::materials[material()]->id()); + mark_as_lost(fmt::format( + "No intersection found with DAGMC cell {}, filled with material {}", + msg.id, material_id)); + break; + } + } // Sample a distance to collision if (type() == ParticleType::electron || type() == ParticleType::positron) { @@ -277,7 +308,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 +439,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 +590,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 +599,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 +614,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."); diff --git a/src/surface.cpp b/src/surface.cpp index 394a4b36630..e41eca4ffdb 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) const { // Determine projection of direction onto normal and squared magnitude of // normal. diff --git a/src/universe.cpp b/src/universe.cpp index 67e1d1354c3..bf89d4069d0 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(Geometron& p) const { const auto& cells { !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())}; From b2688cf7ec95accbb58184e0137deda12fc32076 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Mon, 21 Aug 2023 11:03:51 -0400 Subject: [PATCH 03/17] use geometron in plotting --- include/openmc/particle_data.h | 30 +++++++++++++++++++++--------- include/openmc/plot.h | 6 +++--- src/plot.cpp | 28 ++++++++-------------------- 3 files changed, 32 insertions(+), 32 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 28a0be91a36..8ae934bdf0d 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -201,7 +201,27 @@ class Geometron { public: Geometron(); -public: + // 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; + } + // Number of current coordinate levels int& n_coord() { return n_coord_; } const int& n_coord() const { return n_coord_; } @@ -259,14 +279,6 @@ class Geometron { // Boundary information BoundaryInfo& boundary() { return boundary_; } - // resets all coordinate levels for the particle - void clear() - { - for (auto& level : coord_) - level.reset(); - n_coord_ = 1; - } - #ifdef DAGMC // DagMC state variables moab::DagMC::RayHistory& history() { return history_; } diff --git a/include/openmc/plot.h b/include/openmc/plot.h index ba82f43e27a..8cee79f9b7c 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 Geometron& 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 Geometron& 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; + Geometron p; p.r() = xyz; p.u() = dir; p.coord(0).universe = model::root_universe; diff --git a/src/plot.cpp b/src/plot.cpp index 8d995b0a9a2..1381848206e 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 Geometron& p, int level) { // set cell data if (p.n_coord() <= level) { @@ -78,7 +78,7 @@ 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 Geometron& p, int level) { Cell* c = model::cells.at(p.lowest_coord().cell).get(); data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN; @@ -1234,20 +1234,9 @@ 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); + Geometron p; + p.r() = camera_position_; + p.u() = {1.0, 0.0, 0.0}; int vert = tid; for (int iter = 0; iter <= pixels_[1] / n_threads; iter++) { @@ -1273,18 +1262,17 @@ 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); + p.u() = camera_local_vec.rotate(camera_to_model); } else { // orthographic projection - s.u = looking_direction; + p.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_ + + p.u() = camera_position_ + cam_yaxis * x_pix_coord * orthographic_width_ + cam_zaxis * y_pix_coord * orthographic_width_; } - p.from_source(&s); // put particle at camera bool hitsomething = false; bool intersection_found = true; int loop_counter = 0; From d4cfb9aa4a6766c75ead6792608de711ba16f6ef Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Sat, 9 Sep 2023 18:00:58 -0400 Subject: [PATCH 04/17] pass geometron pointer for dagmc surface reflections --- include/openmc/dagmc.h | 2 +- include/openmc/surface.h | 7 ++++--- src/boundary_condition.cpp | 2 +- src/dagmc.cpp | 2 +- src/surface.cpp | 4 ++-- 5 files changed, 9 insertions(+), 8 deletions(-) diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 7aee45648f4..645c1a160db 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) const override; + Direction reflect(Position r, Direction u, Geometron* p) const override; inline void to_hdf5_inner(hid_t group_id) const override {}; diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 6620303edd4..14ebc9e0f32 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) const; + virtual Direction reflect( + Position r, Direction u, Geometron* p = nullptr) const; virtual Direction diffuse_reflect( - Position r, Direction u, uint64_t* seed) const; + Position r, Direction u, uint64_t* seed, Geometron* p = nullptr) const; //! Evaluate the equation describing the surface. //! diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 193c842128a..42a8e27a164 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -25,7 +25,7 @@ void VacuumBC::handle_particle(Particle& p, const Surface& surf) const void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const { - Direction u = surf.reflect(p.r(), p.u()); + Direction u = surf.reflect(p.r(), p.u(), &p); u /= u.norm(); // Handle the effects of the surface albedo on the particle's weight. diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 002a6b3c83c..0478d25179f 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -715,7 +715,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, Geometron* p) const { Expects(p); p->history().reset_to_last_intersection(); diff --git a/src/surface.cpp b/src/surface.cpp index e41eca4ffdb..a6b403be995 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) const +Direction Surface::reflect(Position r, Direction u, Geometron* p) const { // Determine projection of direction onto normal and squared magnitude of // normal. @@ -140,7 +140,7 @@ Direction Surface::reflect(Position r, Direction u) const } Direction Surface::diffuse_reflect( - Position r, Direction u, uint64_t* seed) const + Position r, Direction u, uint64_t* seed, Geometron* p) const { // Diffuse reflect direction according to the normal. // cosine distribution From bd1c2630442cdb6e24f20db66694366dbe7cd33a Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 5 Oct 2023 11:17:05 -0400 Subject: [PATCH 05/17] make mark_as_lost virtual --- include/openmc/particle.h | 6 +++--- include/openmc/particle_data.h | 22 +++++++++++++++++----- src/particle_data.cpp | 17 +++++++++++++++++ 3 files changed, 37 insertions(+), 8 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 80a561f4d05..7ec4556a5b0 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -103,14 +103,14 @@ 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); + virtual void mark_as_lost(const char* message) override; - void mark_as_lost(const std::string& message) + virtual void mark_as_lost(const std::string& message) override { mark_as_lost(message.c_str()); } - void mark_as_lost(const std::stringstream& message) + virtual void mark_as_lost(const std::stringstream& message) override { mark_as_lost(message.str()); } diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 8ae934bdf0d..3358ad8923c 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -201,6 +201,14 @@ class Geometron { public: Geometron(); + // Geometron 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); + virtual void mark_as_lost(const std::string& message); + virtual void mark_as_lost(const std::stringstream& message); + // resets all coordinate levels for the particle void clear() { @@ -222,6 +230,13 @@ class Geometron { 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_; } @@ -296,6 +311,8 @@ class Geometron { double& sqrtkT_last() { return sqrtkT_last_; } private: + int64_t id_ {-1}; + int n_coord_ {1}; int cell_instance_; vector coord_; @@ -368,7 +385,6 @@ class ParticleData : public Geometron { MacroXS macro_xs_; CacheDataMG mg_xs_cache_; - int64_t id_; ParticleType type_ {ParticleType::neutron}; double E_; @@ -458,10 +474,6 @@ class ParticleData : public Geometron { CacheDataMG& mg_xs_cache() { return mg_xs_cache_; } const CacheDataMG& mg_xs_cache() const { return mg_xs_cache_; } - // Unique ID - 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_; } diff --git a/src/particle_data.cpp b/src/particle_data.cpp index b80b09e8599..30e84fd0be5 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -1,5 +1,7 @@ #include "openmc/particle_data.h" +#include + #include "openmc/cell.h" #include "openmc/geometry.h" #include "openmc/material.h" @@ -12,6 +14,21 @@ namespace openmc { +void Geometron::mark_as_lost(const char* message) +{ + fatal_error(message); +} + +void Geometron::mark_as_lost(const std::string& message) +{ + mark_as_lost(message.c_str()); +} + +void Geometron::mark_as_lost(const std::stringstream& message) +{ + mark_as_lost(message.str()); +} + void LocalCoord::rotate(const vector& rotation) { r = r.rotate(rotation); From c393ead7835fabf7ae0b54c6d05b82f2144b2947 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Fri, 6 Oct 2023 11:01:31 -0400 Subject: [PATCH 06/17] fix dagmc build maybe --- src/dagmc.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 0478d25179f..9cbea5b3b14 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -630,7 +630,13 @@ std::pair DAGCell::distance( // isn't found in a volume that is not the implicit complement. In the case // that the DAGMC model is the root universe of the geometry, even a missing // intersection in the implicit complement should trigger this condition. - throw ParticleLost(ParticleLost::Reason::no_dagmc_intersection, id_); + std::string material_id = + p->material() == MATERIAL_VOID + ? "-1 (VOID)" + : std::to_string(model::materials[p->material()]->id()); + p->mark_as_lost(fmt::format( + "No intersection found with DAGMC cell {}, filled with material {}", id_, + material_id)); } return {dist, surf_idx}; From 4dc0b764cadff9117cd2ef45de2d91eab071f0dd Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Sun, 22 Oct 2023 15:58:02 -0400 Subject: [PATCH 07/17] finish cherry picking --- include/openmc/plot.h | 2 +- src/plot.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 8cee79f9b7c..70544f97853 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -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(Geometron& 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/src/plot.cpp b/src/plot.cpp index 1381848206e..eb67d2d9535 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1084,7 +1084,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(Geometron& p) { constexpr double scoot = 1e-5; double min_dist = {INFINITY}; From d045070d3ee457a19fb32c86f72553b974937700 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 5 Oct 2023 10:45:05 -0400 Subject: [PATCH 08/17] more palatable error handling --- include/openmc/geometry.h | 37 ---------------------------------- include/openmc/particle.h | 10 --------- include/openmc/particle_data.h | 10 +++++++-- src/geometry.cpp | 18 +++++++++++------ src/particle.cpp | 33 +----------------------------- src/particle_data.cpp | 11 +--------- 6 files changed, 22 insertions(+), 97 deletions(-) diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index f44e5c5865a..c4bb19100bb 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -77,43 +77,6 @@ void cross_lattice( BoundaryInfo distance_to_boundary(Geometron& p); -/* Geometry routines can potentially throw this type of exception - */ -class ParticleLost : public std::exception { -public: - enum class Reason { - negative_lattice_distance, // should not have negative distance to lattice - bad_boundary_crossing, // could not locate after crossing a boundary - no_universe_outside_lattice, // undefined behavior. All space should be - // defined. - no_dagmc_intersection - } reason; - - // Extra data to be passed up for exception handling or for detailed - // error reporting. For instance no_dagmc_intersection will give the - // DAGMC cell ID. - int id; - - ParticleLost(Reason reason_a, int id_a = 0) : reason(reason_a), id(id_a) {} - - // Handles uncaught exception messages, e.g. Geometron with no ID. - // These will not show up in the usual particle tracking, as more - // detailed error reporting is provided in the Particle class. - const char* what() - { - switch (reason) { - case Reason::negative_lattice_distance: - return "Negative distance to a lattice"; - case Reason::bad_boundary_crossing: - return "Crossed a boundary and lost particle"; - case Reason::no_universe_outside_lattice: - return "Outside lattice but no outer region defined"; - case Reason::no_dagmc_intersection: - return "No intersection found with DAGMC cell"; - } - } -}; - } // namespace openmc #endif // OPENMC_GEOMETRY_H diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 7ec4556a5b0..52e017493a8 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -105,16 +105,6 @@ class Particle : public ParticleData { //! \param message A warning message to display virtual void mark_as_lost(const char* message) override; - virtual void mark_as_lost(const std::string& message) override - { - mark_as_lost(message.c_str()); - } - - virtual void mark_as_lost(const std::stringstream& message) override - { - mark_as_lost(message.str()); - } - //! 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 3358ad8923c..c81d11f53c5 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -206,8 +206,14 @@ class Geometron { // reporting behavior that occurs in geometry.cpp. The explanation for // mark_as_lost is the same. virtual void mark_as_lost(const char* message); - virtual void mark_as_lost(const std::string& message); - virtual void mark_as_lost(const std::stringstream& 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.c_str()); + } // resets all coordinate levels for the particle void clear() diff --git a/src/geometry.cpp b/src/geometry.cpp index 2c67b0b3350..2c2ebcb591d 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -244,8 +244,9 @@ bool find_cell_inner( if (lat.outer_ != NO_OUTER_UNIVERSE) { coord.universe = lat.outer_; } else { - throw ParticleLost( - ParticleLost::Reason::no_universe_outside_lattice, lat.id_); + p.mark_as_lost(fmt::format( + "Particle {} left lattice {}, but it has no outer definition.", + p.id(), lat.id_)); } } } @@ -337,7 +338,9 @@ void cross_lattice(Geometron& p, const BoundaryInfo& boundary, bool verbose) bool found = exhaustive_find_cell(p); if (!found) { - throw ParticleLost(ParticleLost::Reason::bad_boundary_crossing, lat.id_); + p.mark_as_lost(fmt::format("Particle {} could not be located after " + "crossing a boundary of lattice {}", + p.id(), lat.id_)); } } else { @@ -351,8 +354,9 @@ void cross_lattice(Geometron& p, const BoundaryInfo& boundary, bool verbose) p.n_coord() = 1; bool found = exhaustive_find_cell(p); if (!found) { - throw ParticleLost( - ParticleLost::Reason::bad_boundary_crossing, lat.id_); + p.mark_as_lost(fmt::format("Particle {} could not be located after " + "crossing a boundary of lattice {}", + p.id(), lat.id_)); } } } @@ -405,7 +409,9 @@ BoundaryInfo distance_to_boundary(Geometron& p) level_lat_trans = lattice_distance.second; if (d_lat < 0) { - throw ParticleLost(ParticleLost::Reason::negative_lattice_distance); + p.mark_as_lost(fmt::format("Particle {} had a negative distance " + "to a lattice boundary.", + p.id())); } } diff --git a/src/particle.cpp b/src/particle.cpp index c92e7b31d19..50dcf664dcf 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -205,38 +205,7 @@ void Particle::event_calculate_xs() void Particle::event_advance() { // Find the distance to the nearest boundary - try { - boundary() = distance_to_boundary(*this); - } catch (ParticleLost msg) { - // Lost particles do not mean we should terminate the code. - // Instead we accept that they're lost and only allow a certain - // amount to geometrically fail. - switch (msg.reason) { - case ParticleLost::Reason::negative_lattice_distance: - mark_as_lost(fmt::format( - "Particle {} had a negative distance to a lattice boundary", id())); - break; - case ParticleLost::Reason::bad_boundary_crossing: - mark_as_lost(fmt::format("Particle {} could not be located after " - "crossing a boundary of lattice {}", - id(), msg.id)); - break; - case ParticleLost::Reason::no_universe_outside_lattice: - mark_as_lost(fmt::format( - "Particle {} left lattice {}, but it has no outer definition.", id(), - msg.id)); - break; - case ParticleLost::Reason::no_dagmc_intersection: - std::string material_id = - material() == MATERIAL_VOID - ? "-1 (VOID)" - : std::to_string(model::materials[material()]->id()); - mark_as_lost(fmt::format( - "No intersection found with DAGMC cell {}, filled with material {}", - msg.id, material_id)); - break; - } - } + boundary() = distance_to_boundary(*this); // Sample a distance to collision if (type() == ParticleType::electron || type() == ParticleType::positron) { diff --git a/src/particle_data.cpp b/src/particle_data.cpp index 30e84fd0be5..6b6c0bc61e6 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -3,6 +3,7 @@ #include #include "openmc/cell.h" +#include "openmc/error.h" #include "openmc/geometry.h" #include "openmc/material.h" #include "openmc/nuclide.h" @@ -19,16 +20,6 @@ void Geometron::mark_as_lost(const char* message) fatal_error(message); } -void Geometron::mark_as_lost(const std::string& message) -{ - mark_as_lost(message.c_str()); -} - -void Geometron::mark_as_lost(const std::stringstream& message) -{ - mark_as_lost(message.str()); -} - void LocalCoord::rotate(const vector& rotation) { r = r.rotate(rotation); From 29816baef0420ebf256808402bb0a61f5a3d668e Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 4 Jan 2024 18:02:50 -0500 Subject: [PATCH 09/17] fully remove ParticleLostException --- include/openmc/particle.h | 1 - include/openmc/particle_data.h | 11 +++------- src/particle.cpp | 39 ++++++++++++++++++++-------------- src/particle_data.cpp | 10 +++++++++ 4 files changed, 36 insertions(+), 25 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 52e017493a8..63e502fb724 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" diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index c81d11f53c5..a494ada37a6 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -206,14 +206,8 @@ class Geometron { // 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) - { - mark_as_lost(message.c_str()); - } - void mark_as_lost(const std::stringstream& message) - { - mark_as_lost(message.c_str()); - } + 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() @@ -310,6 +304,7 @@ class Geometron { 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_; } diff --git a/src/particle.cpp b/src/particle.cpp index 50dcf664dcf..f24c0da0563 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -159,7 +159,7 @@ void Particle::event_calculate_xs() // beginning of the history and again for any secondary particles if (lowest_coord().cell == C_NONE) { if (!exhaustive_find_cell(*this)) { - mark_as_lost( + Geometron::mark_as_lost( "Could not find the cell containing particle " + std::to_string(id())); return; } @@ -410,8 +410,13 @@ void Particle::event_revive_from_secondary() if (lowest_coord().cell == C_NONE) { 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())); + // Note: Geometron:: prefix is necessary here because it's ambiguous + // whether Particle::mark_as_lost(const char*) should be called with + // implicit argument conversion or Geometron::mark_as_lost(const + // std::string&) should be called. + Geometron::mark_as_lost( + "Could not find the cell containing particle " + + std::to_string(id())); return; } // Set birth cell attribute @@ -584,9 +589,10 @@ void Particle::cross_surface() // undefined region in the geometry. 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."); + Geometron::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."); return; } } @@ -622,8 +628,8 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u) { // Do not handle reflective boundary conditions on lower universes if (n_coord() != 1) { - mark_as_lost("Cannot reflect particle " + std::to_string(id()) + - " off surface in a lower universe."); + Geometron::mark_as_lost("Cannot reflect particle " + std::to_string(id()) + + " off surface in a lower universe."); return; } @@ -658,8 +664,9 @@ 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_) + "."); + Geometron::mark_as_lost( + "Couldn't find particle after reflecting from surface " + + std::to_string(surf.id_) + "."); return; } @@ -677,7 +684,7 @@ void Particle::cross_periodic_bc( { // Do not handle periodic boundary conditions on lower universes if (n_coord() != 1) { - mark_as_lost( + Geometron::mark_as_lost( "Cannot transfer particle " + std::to_string(id()) + " across surface in a lower universe. Boundary conditions must be " "applied to root universe."); @@ -705,11 +712,11 @@ void Particle::cross_periodic_bc( n_coord() = 1; if (!neighbor_list_find_cell(*this)) { - mark_as_lost("Couldn't find particle after hitting periodic " - "boundary on surface " + - std::to_string(surf.id_) + - ". The normal vector " - "of one periodic surface may need to be reversed."); + Geometron::mark_as_lost("Couldn't find particle after hitting periodic " + "boundary on surface " + + std::to_string(surf.id_) + + ". The normal vector " + "of one periodic surface may need to be reversed."); return; } diff --git a/src/particle_data.cpp b/src/particle_data.cpp index 6b6c0bc61e6..1a8ebbf99fe 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -15,6 +15,16 @@ namespace openmc { +void Geometron::mark_as_lost(const std::string& message) +{ + mark_as_lost(message.c_str()); +} + +void Geometron::mark_as_lost(const std::stringstream& message) +{ + mark_as_lost(message.str()); +} + void Geometron::mark_as_lost(const char* message) { fatal_error(message); From 294fa139647c59bc9137c3f7e7106d1308fa4dd9 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 4 Jan 2024 18:13:27 -0500 Subject: [PATCH 10/17] appease disparate cpplinter version --- src/plot.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/plot.cpp b/src/plot.cpp index eb67d2d9535..be0b1b25aa8 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1268,9 +1268,10 @@ void ProjectionPlot::create_output() const double x_pix_coord = (static_cast(horiz) - p0 / 2.0) / p0; double y_pix_coord = (static_cast(vert) - p1 / 2.0) / p0; - p.u() = camera_position_ + - cam_yaxis * x_pix_coord * orthographic_width_ + - cam_zaxis * y_pix_coord * orthographic_width_; + + p.u() = camera_position_; + p.u() += cam_yaxis * x_pix_coord * orthographic_width_; + p.u() += cam_zaxis * y_pix_coord * orthographic_width_; } bool hitsomething = false; From 2381d0e71a029dd882531563b9846795ea18a334 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Sat, 6 Jan 2024 14:25:53 -0600 Subject: [PATCH 11/17] fix particle initialization in projection plots --- src/plot.cpp | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/plot.cpp b/src/plot.cpp index be0b1b25aa8..0f1682e18b2 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1235,7 +1235,6 @@ void ProjectionPlot::create_output() const const int tid = thread_num(); Geometron p; - p.r() = camera_position_; p.u() = {1.0, 0.0, 0.0}; int vert = tid; @@ -1252,6 +1251,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 = @@ -1262,18 +1265,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); - p.u() = camera_local_vec.rotate(camera_to_model); + init_u = camera_local_vec.rotate(camera_to_model); + init_r = camera_position_; } else { // orthographic projection - p.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; - p.u() = camera_position_; - p.u() += cam_yaxis * x_pix_coord * orthographic_width_; - p.u() += 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_; } + // 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; From 046294631eb4429f9d2516d33244316d0c2754d0 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Sat, 6 Jan 2024 14:29:54 -0600 Subject: [PATCH 12/17] regold proj plot test for tiny floating point error --- include/openmc/particle.h | 1 + include/openmc/particle_data.h | 10 +++-- src/particle.cpp | 39 ++++++++----------- .../plot_projections/results_true.dat | 2 +- 4 files changed, 24 insertions(+), 28 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 63e502fb724..09d8a10f1ce 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -103,6 +103,7 @@ class Particle : public ParticleData { //! mark a particle as lost and create a particle restart file //! \param message A warning message to display virtual void mark_as_lost(const char* message) override; + using Geometron::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 a494ada37a6..d874a75d8b7 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -201,10 +201,12 @@ class Geometron { public: Geometron(); - // Geometron 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. + /* + * Geometron 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); diff --git a/src/particle.cpp b/src/particle.cpp index f24c0da0563..549de5cff65 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -159,7 +159,7 @@ void Particle::event_calculate_xs() // beginning of the history and again for any secondary particles if (lowest_coord().cell == C_NONE) { if (!exhaustive_find_cell(*this)) { - Geometron::mark_as_lost( + mark_as_lost( "Could not find the cell containing particle " + std::to_string(id())); return; } @@ -410,13 +410,8 @@ void Particle::event_revive_from_secondary() if (lowest_coord().cell == C_NONE) { bool verbose = settings::verbosity >= 10 || trace(); if (!exhaustive_find_cell(*this, verbose)) { - // Note: Geometron:: prefix is necessary here because it's ambiguous - // whether Particle::mark_as_lost(const char*) should be called with - // implicit argument conversion or Geometron::mark_as_lost(const - // std::string&) should be called. - Geometron::mark_as_lost( - "Could not find the cell containing particle " + - std::to_string(id())); + mark_as_lost("Could not find the cell containing particle " + + std::to_string(id())); return; } // Set birth cell attribute @@ -589,10 +584,9 @@ void Particle::cross_surface() // undefined region in the geometry. if (!exhaustive_find_cell(*this, verbose)) { - Geometron::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."); + 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."); return; } } @@ -628,8 +622,8 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u) { // Do not handle reflective boundary conditions on lower universes if (n_coord() != 1) { - Geometron::mark_as_lost("Cannot reflect particle " + std::to_string(id()) + - " off surface in a lower universe."); + mark_as_lost("Cannot reflect particle " + std::to_string(id()) + + " off surface in a lower universe."); return; } @@ -664,9 +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)) { - Geometron::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; } @@ -684,7 +677,7 @@ void Particle::cross_periodic_bc( { // Do not handle periodic boundary conditions on lower universes if (n_coord() != 1) { - Geometron::mark_as_lost( + mark_as_lost( "Cannot transfer particle " + std::to_string(id()) + " across surface in a lower universe. Boundary conditions must be " "applied to root universe."); @@ -712,11 +705,11 @@ void Particle::cross_periodic_bc( n_coord() = 1; if (!neighbor_list_find_cell(*this)) { - Geometron::mark_as_lost("Couldn't find particle after hitting periodic " - "boundary on surface " + - std::to_string(surf.id_) + - ". The normal vector " - "of one periodic surface may need to be reversed."); + mark_as_lost("Couldn't find particle after hitting periodic " + "boundary on surface " + + std::to_string(surf.id_) + + ". The normal vector " + "of one periodic surface may need to be reversed."); return; } diff --git a/tests/regression_tests/plot_projections/results_true.dat b/tests/regression_tests/plot_projections/results_true.dat index 1a67da3002a..a7cd235d3b5 100644 --- a/tests/regression_tests/plot_projections/results_true.dat +++ b/tests/regression_tests/plot_projections/results_true.dat @@ -1 +1 @@ -24fb0f41ee018ea086962dbd6bcd0b536d11d4b34644bfef4f0e74f8b462fe41a84af39c7ff79046d5d7cfe209084eac54712fa0ec01038e97eb43df1abd0334 \ No newline at end of file +cc00d9d011dc7349a165d347901ca6662456f2433271e56984d471f2399aac52a375c7488ab2902596df2a876f5e3abf5e9aebd625d17e3d9e4a5f935819d45b \ No newline at end of file From b66468ae17d3737fd7558a9ebd1ae58ebd175a90 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Fri, 12 Jan 2024 19:58:11 -0600 Subject: [PATCH 13/17] reintroduce end-of-line comments --- include/openmc/particle_data.h | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index d874a75d8b7..784d0b01763 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -314,28 +314,30 @@ class Geometron { double& sqrtkT_last() { return sqrtkT_last_; } private: - int64_t id_ {-1}; + int64_t id_ {-1}; //!< Unique ID - int n_coord_ {1}; - int cell_instance_; - vector coord_; + 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}; - vector cell_last_; + int n_coord_last_ {1}; //!< number of current coordinates + vector cell_last_; //!< coordinates for all levels - Position r_last_current_; - Position r_last_; - Direction u_last_; + 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}; + int surface_ {0}; //!< index for surface particle is on - BoundaryInfo boundary_; + BoundaryInfo boundary_; //!< Info about the next intersection - int material_ {-1}; - int material_last_ {-1}; + int material_ {-1}; //!< index for current material + int material_last_ {-1}; //!< index for last material - double sqrtkT_ {-1.0}; - double sqrtkT_last_ {0.0}; + double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV + double sqrtkT_last_ {0.0}; //!< last temperature #ifdef DAGMC moab::DagMC::RayHistory history_; From 926f2ba733e318faed7a8c26a5ced831d6582d8b Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Fri, 12 Jan 2024 19:59:46 -0600 Subject: [PATCH 14/17] address comments about p --- src/geometry.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/geometry.cpp b/src/geometry.cpp index 2c2ebcb591d..6bb545693c2 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -465,18 +465,19 @@ BoundaryInfo distance_to_boundary(Geometron& p) extern "C" int openmc_find_cell( const double* xyz, int32_t* index, int32_t* instance) { - Geometron p; + Geometron 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; } From b9ff5d5fb3d093a4a167d6d992ee3b181fafe0dc Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Fri, 12 Jan 2024 20:04:45 -0600 Subject: [PATCH 15/17] Geometron -> GeometryState --- include/openmc/cell.h | 11 ++++++----- include/openmc/dagmc.h | 8 ++++---- include/openmc/geometry.h | 14 +++++++------- include/openmc/particle.h | 2 +- include/openmc/particle_data.h | 8 ++++---- include/openmc/plot.h | 8 ++++---- include/openmc/surface.h | 4 ++-- include/openmc/universe.h | 4 ++-- src/cell.cpp | 5 +++-- src/dagmc.cpp | 6 +++--- src/geometry.cpp | 16 ++++++++-------- src/particle_data.cpp | 8 ++++---- src/plot.cpp | 9 +++++---- src/surface.cpp | 4 ++-- src/universe.cpp | 2 +- 15 files changed, 56 insertions(+), 53 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 00eae4d520a..c405f536b5d 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -40,7 +40,7 @@ constexpr int32_t OP_UNION {std::numeric_limits::max() - 4}; //============================================================================== class Cell; -class Geometron; +class GeometryState; class ParentCell; class CellInstance; class Universe; @@ -184,7 +184,7 @@ class Cell { //! Find the oncoming boundary of this cell. virtual std::pair distance( - Position r, Direction u, int32_t on_surface, Geometron* 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. @@ -261,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, Geometron& 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 @@ -333,8 +334,8 @@ class CSGCell : public Cell { // Methods vector surfaces() const override { return region_.surfaces(); } - std::pair distance( - Position r, Direction u, int32_t on_surface, Geometron* 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); } diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 645c1a160db..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, Geometron* 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, Geometron* 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(Geometron& 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 c4bb19100bb..107cc7d1f3e 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -11,7 +11,7 @@ namespace openmc { class BoundaryInfo; -class Geometron; +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(Geometron& 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(Geometron& p, bool error = true); //! should be computed. \return The instance of the cell at the specified level. //============================================================================== -int cell_instance_at_level(const Geometron& 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,22 +60,22 @@ int cell_instance_at_level(const Geometron& 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(Geometron& p, bool verbose = false); +bool exhaustive_find_cell(GeometryState& p, bool verbose = false); bool neighbor_list_find_cell( - Geometron& p, bool verbose = false); // Only usable on surface crossings + GeometryState& p, bool verbose = false); // Only usable on surface crossings //============================================================================== //! Move a particle into a new lattice tile. //============================================================================== void cross_lattice( - Geometron& p, const BoundaryInfo& boundary, bool verbose = false); + GeometryState& p, const BoundaryInfo& boundary, bool verbose = false); //============================================================================== //! Find the next boundary a particle will intersect. //============================================================================== -BoundaryInfo distance_to_boundary(Geometron& p); +BoundaryInfo distance_to_boundary(GeometryState& p); } // namespace openmc diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 09d8a10f1ce..e423880f87c 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -103,7 +103,7 @@ class Particle : public ParticleData { //! mark a particle as lost and create a particle restart file //! \param message A warning message to display virtual void mark_as_lost(const char* message) override; - using Geometron::mark_as_lost; + 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 784d0b01763..cb68fc45740 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -197,12 +197,12 @@ struct BoundaryInfo { /* * Contains all geometry state information for a particle. */ -class Geometron { +class GeometryState { public: - Geometron(); + GeometryState(); /* - * Geometron does not store any ID info, so give some reasonable behavior + * 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. @@ -380,7 +380,7 @@ class Geometron { * Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18. * https://doi.org/10.1016/j.anucene.2017.11.032. */ -class ParticleData : public Geometron { +class ParticleData : public GeometryState { private: //========================================================================== // Data members -- see public: below for descriptions diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 70544f97853..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 Geometron& 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 Geometron& 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 { - Geometron 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(Geometron& 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 14ebc9e0f32..350775123c1 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -108,10 +108,10 @@ class Surface { //! \param[inout] p Pointer to the particle. Only DAGMC uses this. //! \return Outgoing direction of the ray virtual Direction reflect( - Position r, Direction u, Geometron* p = nullptr) const; + Position r, Direction u, GeometryState* p = nullptr) const; virtual Direction diffuse_reflect( - Position r, Direction u, uint64_t* seed, Geometron* p = nullptr) 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 bd872eb0535..26f33cb383d 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -9,7 +9,7 @@ namespace openmc { class DAGUniverse; #endif -class Geometron; +class GeometryState; class Universe; class UniversePartitioner; @@ -33,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(Geometron& p) const; + virtual bool find_cell(GeometryState& p) const; BoundingBox bounding_box() const; diff --git a/src/cell.cpp b/src/cell.cpp index c0999214677..9d4aea7b891 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1279,14 +1279,15 @@ vector Cell::find_parent_cells( { // create a temporary particle - Geometron 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, Geometron& 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 9cbea5b3b14..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(Geometron& 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,7 +583,7 @@ DAGCell::DAGCell(std::shared_ptr dag_ptr, int32_t dag_idx) }; std::pair DAGCell::distance( - Position r, Direction u, int32_t on_surface, Geometron* p) const + Position r, Direction u, int32_t on_surface, GeometryState* p) const { // if we've changed direction or we're not on a surface, // reset the history and update last direction @@ -721,7 +721,7 @@ Direction DAGSurface::normal(Position r) const return dir; } -Direction DAGSurface::reflect(Position r, Direction u, Geometron* 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 6bb545693c2..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(Geometron& 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(Geometron& p, bool error) //============================================================================== -int cell_instance_at_level(const Geometron& 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) { @@ -100,7 +100,7 @@ int cell_instance_at_level(const Geometron& p, int level) //============================================================================== bool find_cell_inner( - Geometron& p, const NeighborList* neighbor_list, bool verbose) + 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. @@ -259,7 +259,7 @@ bool find_cell_inner( //============================================================================== -bool neighbor_list_find_cell(Geometron& p, bool verbose) +bool neighbor_list_find_cell(GeometryState& p, bool verbose) { // Reset all the deeper coordinate levels. @@ -287,7 +287,7 @@ bool neighbor_list_find_cell(Geometron& p, bool verbose) return found; } -bool exhaustive_find_cell(Geometron& p, bool verbose) +bool exhaustive_find_cell(GeometryState& p, bool verbose) { int i_universe = p.lowest_coord().universe; if (i_universe == C_NONE) { @@ -304,7 +304,7 @@ bool exhaustive_find_cell(Geometron& p, bool verbose) //============================================================================== -void cross_lattice(Geometron& p, const BoundaryInfo& boundary, bool verbose) +void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose) { auto& coord {p.lowest_coord()}; auto& lat {*model::lattices[coord.lattice]}; @@ -364,7 +364,7 @@ void cross_lattice(Geometron& p, const BoundaryInfo& boundary, bool verbose) //============================================================================== -BoundaryInfo distance_to_boundary(Geometron& p) +BoundaryInfo distance_to_boundary(GeometryState& p) { BoundaryInfo info; double d_lat = INFINITY; @@ -465,7 +465,7 @@ BoundaryInfo distance_to_boundary(Geometron& p) extern "C" int openmc_find_cell( const double* xyz, int32_t* index, int32_t* instance) { - Geometron geom_state; + GeometryState geom_state; geom_state.r() = Position {xyz}; geom_state.u() = {0.0, 0.0, 1.0}; diff --git a/src/particle_data.cpp b/src/particle_data.cpp index 1a8ebbf99fe..d8a5bb9d99c 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -15,17 +15,17 @@ namespace openmc { -void Geometron::mark_as_lost(const std::string& message) +void GeometryState::mark_as_lost(const std::string& message) { mark_as_lost(message.c_str()); } -void Geometron::mark_as_lost(const std::stringstream& message) +void GeometryState::mark_as_lost(const std::stringstream& message) { mark_as_lost(message.str()); } -void Geometron::mark_as_lost(const char* message) +void GeometryState::mark_as_lost(const char* message) { fatal_error(message); } @@ -48,7 +48,7 @@ void LocalCoord::reset() rotated = false; } -Geometron::Geometron() +GeometryState::GeometryState() { // Create and clear coordinate levels coord_.resize(model::n_coord_levels); diff --git a/src/plot.cpp b/src/plot.cpp index 0f1682e18b2..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 Geometron& 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 Geometron& 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(Geometron& p) +int ProjectionPlot::advance_to_boundary_from_void(GeometryState& p) { constexpr double scoot = 1e-5; double min_dist = {INFINITY}; @@ -1234,7 +1235,7 @@ void ProjectionPlot::create_output() const const int n_threads = num_threads(); const int tid = thread_num(); - Geometron p; + GeometryState p; p.u() = {1.0, 0.0, 0.0}; int vert = tid; diff --git a/src/surface.cpp b/src/surface.cpp index a6b403be995..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, Geometron* 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, Geometron* p) const } Direction Surface::diffuse_reflect( - Position r, Direction u, uint64_t* seed, Geometron* p) 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 bf89d4069d0..b4ef6b4b265 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -37,7 +37,7 @@ void Universe::to_hdf5(hid_t universes_group) const close_group(group); } -bool Universe::find_cell(Geometron& p) const +bool Universe::find_cell(GeometryState& p) const { const auto& cells { !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())}; From 0fc15e7943f82bf1f33499f1eff23d91380e5f79 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Sat, 13 Jan 2024 22:05:04 -0600 Subject: [PATCH 16/17] git screwup oops --- tests/regression_tests/plot_projections/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/plot_projections/results_true.dat b/tests/regression_tests/plot_projections/results_true.dat index a7cd235d3b5..d4e1ba896e2 100644 --- a/tests/regression_tests/plot_projections/results_true.dat +++ b/tests/regression_tests/plot_projections/results_true.dat @@ -1 +1 @@ -cc00d9d011dc7349a165d347901ca6662456f2433271e56984d471f2399aac52a375c7488ab2902596df2a876f5e3abf5e9aebd625d17e3d9e4a5f935819d45b \ No newline at end of file +24fb0f41ee018ea086962dbd6bcd0b536d11d4b34644bfef4f0e74f8b462fe41a84af39c7ff79046d5d7cfe209084eac54712fa0ec01038e97eb43df1abd0334 From f9e338f26b6bc62296280e3ba0de33bafb9794a5 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Sat, 13 Jan 2024 23:07:42 -0600 Subject: [PATCH 17/17] no newline on file... ah --- tests/regression_tests/plot_projections/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/plot_projections/results_true.dat b/tests/regression_tests/plot_projections/results_true.dat index d4e1ba896e2..1a67da3002a 100644 --- a/tests/regression_tests/plot_projections/results_true.dat +++ b/tests/regression_tests/plot_projections/results_true.dat @@ -1 +1 @@ -24fb0f41ee018ea086962dbd6bcd0b536d11d4b34644bfef4f0e74f8b462fe41a84af39c7ff79046d5d7cfe209084eac54712fa0ec01038e97eb43df1abd0334 +24fb0f41ee018ea086962dbd6bcd0b536d11d4b34644bfef4f0e74f8b462fe41a84af39c7ff79046d5d7cfe209084eac54712fa0ec01038e97eb43df1abd0334 \ No newline at end of file