From 2af4c9cd9219e30a85de6262bf3ec7eb4a08652f Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Fri, 16 Apr 2021 15:35:33 -0400 Subject: [PATCH] use base class to handle layout of particle data --- CMakeLists.txt | 1 + include/openmc/bank.h | 6 +- include/openmc/event.h | 2 +- include/openmc/geometry.h | 6 +- include/openmc/mesh.h | 9 +- include/openmc/particle.h | 424 +--------------- include/openmc/particle_data.h | 451 ++++++++++++++++++ include/openmc/physics.h | 2 +- include/openmc/reaction_product.h | 2 +- include/openmc/simulation.h | 3 - include/openmc/source.h | 16 +- include/openmc/state_point.h | 3 +- include/openmc/tallies/filter_particle.h | 6 +- src/bank.cpp | 10 +- src/bremsstrahlung.cpp | 6 +- src/eigenvalue.cpp | 2 +- src/initialize.cpp | 2 +- src/material.cpp | 6 +- src/mesh.cpp | 12 +- src/mgxs_interface.cpp | 2 +- src/nuclide.cpp | 7 +- src/output.cpp | 8 +- src/particle.cpp | 423 ++++++++-------- src/particle_data.cpp | 35 ++ src/particle_restart.cpp | 12 +- src/photon.cpp | 8 +- src/physics.cpp | 48 +- src/physics_mg.cpp | 7 +- src/reaction.cpp | 2 +- src/reaction_product.cpp | 2 +- src/simulation.cpp | 25 +- src/source.cpp | 14 +- src/state_point.cpp | 29 +- src/tallies/derivative.cpp | 6 +- src/tallies/filter_particle.cpp | 7 +- src/tallies/tally.cpp | 2 +- src/tallies/tally_scoring.cpp | 7 +- .../source_dlopen/source_sampling.cpp | 8 +- .../parameterized_source_sampling.cpp | 8 +- 39 files changed, 834 insertions(+), 795 deletions(-) create mode 100644 include/openmc/particle_data.h create mode 100644 src/particle_data.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f3cafcc5c..1b5c1e7e1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -295,6 +295,7 @@ list(APPEND libopenmc_SOURCES src/nuclide.cpp src/output.cpp src/particle.cpp + src/particle_data.cpp src/particle_restart.cpp src/photon.cpp src/physics.cpp diff --git a/include/openmc/bank.h b/include/openmc/bank.h index 4f6f95c19..1eae1774b 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -16,11 +16,11 @@ namespace openmc { namespace simulation { -extern std::vector source_bank; +extern std::vector source_bank; -extern SharedArray surf_source_bank; +extern SharedArray surf_source_bank; -extern SharedArray fission_bank; +extern SharedArray fission_bank; extern std::vector progeny_per_particle; diff --git a/include/openmc/event.h b/include/openmc/event.h index 7428eedd5..f45b2a18a 100644 --- a/include/openmc/event.h +++ b/include/openmc/event.h @@ -25,7 +25,7 @@ namespace openmc { // consistent locality improvements. struct EventQueueItem{ int64_t idx; //!< particle index in event-based particle buffer - Particle::Type type; //!< particle type + ParticleType type; //!< particle type int64_t material; //!< material that particle is in double E; //!< particle energy diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 0681952ad..9ca0a9e9e 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -6,11 +6,13 @@ #include #include -#include "openmc/particle.h" - +#include "openmc/constants.h" namespace openmc { +class BoundaryInfo; +class Particle; + //============================================================================== // Global variables //============================================================================== diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 2af2e87df..5904abb6e 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -158,8 +158,8 @@ public: //! \param[in] Pointer to bank sites //! \param[in] Number of bank sites //! \param[out] Whether any bank sites are outside the mesh - xt::xtensor count_sites(const Particle::Bank* bank, - int64_t length, bool* outside) const; + xt::xtensor count_sites( + const ParticleBank* bank, int64_t length, bool* outside) const; //! Get bin given mesh indices // @@ -256,9 +256,8 @@ public: //! \param[in] bank Array of bank sites //! \param[out] Whether any bank sites are outside the mesh //! \return Array indicating number of sites in each mesh/energy bin - xt::xtensor count_sites(const Particle::Bank* bank, - int64_t length, - bool* outside) const; + xt::xtensor count_sites( + const ParticleBank* bank, int64_t length, bool* outside) const; // Data members double volume_frac_; //!< Volume fraction of each mesh element diff --git a/include/openmc/particle.h b/include/openmc/particle.h index f0e146743..77dec18b6 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -11,6 +11,7 @@ #include #include "openmc/constants.h" +#include "openmc/particle_data.h" #include "openmc/position.h" #include "openmc/random_lcg.h" #include "openmc/tallies/filter_match.h" @@ -22,173 +23,22 @@ namespace openmc { -//============================================================================== -// Constants -//============================================================================== - -// Since cross section libraries come with different numbers of delayed groups -// (e.g. ENDF/B-VII.1 has 6 and JEFF 3.1.1 has 8 delayed groups) and we don't -// yet know what cross section library is being used when the tallies.xml file -// is read in, we want to have an upper bound on the size of the array we -// use to store the bins for delayed group tallies. -constexpr int MAX_DELAYED_GROUPS {8}; - -constexpr double CACHE_INVALID {-1.0}; - -//============================================================================== -// Class declarations -//============================================================================== - -// Forward declare the Surface class for use in function arguments. +// Forward declare the Surface class for use in Particle::cross_vacuum_bc, etc. class Surface; -class LocalCoord { -public: - void rotate(const std::vector& rotation); - - //! clear data from a single coordinate level - void reset(); - - Position r; //!< particle position - Direction u; //!< particle direction - int cell {-1}; - int universe {-1}; - int lattice {-1}; - array lattice_i {-1, -1, -1}; - bool rotated {false}; //!< Is the level rotated? -}; - -//============================================================================== -//! Cached microscopic cross sections for a particular nuclide at the current -//! energy -//============================================================================== - -struct NuclideMicroXS { - // Microscopic cross sections in barns - double total; //!< total cross section - double absorption; //!< absorption (disappearance) - double fission; //!< fission - double nu_fission; //!< neutron production from fission - - double elastic; //!< If sab_frac is not 1 or 0, then this value is - //!< averaged over bound and non-bound nuclei - double thermal; //!< Bound thermal elastic & inelastic scattering - double thermal_elastic; //!< Bound thermal elastic scattering - double photon_prod; //!< microscopic photon production xs - - // Cross sections for depletion reactions (note that these are not stored in - // macroscopic cache) - double reaction[DEPLETION_RX.size()]; - - // Indicies and factors needed to compute cross sections from the data tables - int index_grid; //!< Index on nuclide energy grid - int index_temp; //!< Temperature index for nuclide - double interp_factor; //!< Interpolation factor on nuc. energy grid - int index_sab {-1}; //!< Index in sab_tables - int index_temp_sab; //!< Temperature index for sab_tables - double sab_frac; //!< Fraction of atoms affected by S(a,b) - bool use_ptable; //!< In URR range with probability tables? - - // Energy and temperature last used to evaluate these cross sections. If - // these values have changed, then the cross sections must be re-evaluated. - double last_E {0.0}; //!< Last evaluated energy - double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant - //!< * temperature (eV)) -}; - -//============================================================================== -//! Cached microscopic photon cross sections for a particular element at the -//! current energy -//============================================================================== - -struct ElementMicroXS { - int index_grid; //!< index on element energy grid - double last_E {0.0}; //!< last evaluated energy in [eV] - double interp_factor; //!< interpolation factor on energy grid - double total; //!< microscopic total photon xs - double coherent; //!< microscopic coherent xs - double incoherent; //!< microscopic incoherent xs - double photoelectric; //!< microscopic photoelectric xs - double pair_production; //!< microscopic pair production xs -}; - -//============================================================================== -// MACROXS contains cached macroscopic cross sections for the material a -// particle is traveling through -//============================================================================== - -struct MacroXS { - double total; //!< macroscopic total xs - double absorption; //!< macroscopic absorption xs - double fission; //!< macroscopic fission xs - double nu_fission; //!< macroscopic production xs - double photon_prod; //!< macroscopic photon production xs - - // Photon cross sections - double coherent; //!< macroscopic coherent xs - double incoherent; //!< macroscopic incoherent xs - double photoelectric; //!< macroscopic photoelectric xs - double pair_production; //!< macroscopic pair production xs -}; - -//============================================================================== -// Information about nearest boundary crossing -//============================================================================== - -struct BoundaryInfo { - double distance {INFINITY}; //!< distance to nearest boundary - int surface_index {0}; //!< if boundary is surface, index in surfaces vector - int coord_level; //!< coordinate level after crossing boundary - std::array lattice_translation {}; //!< which way lattice indices will change -}; - //============================================================================ //! State of a particle being transported through geometry +//! This class defines actions particles can take. Its base +//! class defines particle data layout in memory. //============================================================================ -class Particle { +class Particle : public ParticleData { public: - //========================================================================== - // Aliases and type definitions - - //! Particle types - enum class Type { - neutron, photon, electron, positron - }; - - //! Saved ("banked") state of a particle - //! NOTE: This structure's MPI type is built in initialize_mpi() of - //! initialize.cpp. Any changes made to the struct here must also be - //! made when building the Bank MPI type in initialize_mpi(). - //! NOTE: This structure is also used on the python side, and is defined - //! in lib/core.py. Changes made to the type here must also be made to the - //! python defintion. - struct Bank { - Position r; - Direction u; - double E; - double wgt; - int delayed_group; - int surf_id; - Type particle; - int64_t parent_id; - int64_t progeny_id; - }; - - //! Saved ("banked") state of a particle, for nu-fission tallying - struct NuBank { - double E; //!< particle energy - double wgt; //!< particle weight - int delayed_group; //!< particle delayed group - }; //========================================================================== // Constructors - Particle(); - - //! resets all coordinate levels for the particle - void clear(); + Particle() = default; //! create a secondary particle // @@ -198,7 +48,7 @@ public: //! \param u Direction of the secondary particle //! \param E Energy of the secondary particle in [eV] //! \param type Particle type - void create_secondary(double wgt, Direction u, double E, Type type); + void create_secondary(double wgt, Direction u, double E, ParticleType type); //! initialize from a source site // @@ -206,7 +56,7 @@ public: //! site may have been produced from an external source, from fission, or //! simply as a secondary particle. //! \param src Source site data - void from_source(const Bank* src); + void from_source(const ParticleBank* src); // Coarse-grained particle events void event_calculate_xs(); @@ -255,269 +105,15 @@ public: //! create a particle restart HDF5 file void write_restart() const; - - //! 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_;} - - //! Force recalculation of neutron xs by setting last energy to zero - void invalidate_neutron_xs() - { - for (auto& micro : neutron_xs_) - micro.last_E = 0.0; - } - -private: - //========================================================================== - // Data members (accessor methods are below) - - // Cross section caches - std::vector neutron_xs_; //!< Microscopic neutron cross sections - std::vector photon_xs_; //!< Microscopic photon cross sections - MacroXS macro_xs_; //!< Macroscopic cross sections - - int64_t id_; //!< Unique ID - Type type_ {Type::neutron}; //!< Particle type (n, p, e, etc.) - - int n_coord_ {1}; //!< number of current coordinate levels - int cell_instance_; //!< offset for distributed properties - std::vector coord_; //!< coordinates for all levels - - // Particle coordinates before crossing a surface - int n_coord_last_ {1}; //!< number of current coordinates - std::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_ {0}; //!< 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 - bool alive_ {true}; //!< is particle alive? - - // 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 - double wgt_absorb_ {0.0}; //!< weight absorbed for survival biasing - - // 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 - - // Boundary information - BoundaryInfo boundary_; - - // Temperature of current cell - double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV - double sqrtkT_last_ {0.0}; //!< last temperature - - // Statistical data - int n_collision_ {0}; //!< number of collisions - - // Track output - bool write_track_ {false}; - - // Current PRNG state - uint64_t seeds_[N_STREAMS]; // current seeds - int stream_; // current RNG stream - - // Secondary particle bank - std::vector secondary_bank_; - - int64_t current_work_; // current work index - - std::vector flux_derivs_; // for derivatives for this particle - - std::vector filter_matches_; // tally filter matches - - std::vector> tracks_; // tracks for outputting to file - - std::vector nu_bank_; // bank of most recently fissioned particles - - // 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 - - double collision_distance_; // distance to particle's next closest collision - - int n_event_ {0}; // number of events executed in this particle's history - - // 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 - -public: - //========================================================================== - // Methods and accessors - - NuclideMicroXS& neutron_xs(const int& i) { return neutron_xs_[i]; } - const NuclideMicroXS& neutron_xs(const int& i) const - { - return neutron_xs_[i]; - } - ElementMicroXS& photon_xs(const int& i) { return photon_xs_[i]; } - MacroXS& macro_xs() { return macro_xs_; } - const MacroXS& macro_xs() const { return macro_xs_; } - - int64_t& id() { return id_; } - Type& type() { return type_; } - const Type& 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(const int& i) { return coord_[i]; } - const LocalCoord& coord(const int& i) const { return coord_[i]; } - - int& n_coord_last() { return n_coord_last_; } - const int& n_coord_last() const { return n_coord_last_; } - int& cell_last(const int& i) { return cell_last_[i]; } - const int& cell_last(const int& i) const { return cell_last_[i]; } - - double& E() { return E_; } - const double& E() const { return E_; } - double& E_last() { return E_last_; } - const double& E_last() const { return E_last_; } - int& g() { return g_; } - const int& g() const { return g_; } - int& g_last() { return g_last_; } - const int& g_last() const { return g_last_; } - - double& wgt() { return wgt_; } - double& mu() { return mu_; } - const double& mu() const { return mu_; } - bool& alive() { return alive_; } - - 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_; } - double& wgt_absorb() { return wgt_absorb_; } - const double& wgt_absorb() const { return wgt_absorb_; } - - bool& fission() { return fission_; } - TallyEvent& event() { return event_; } - const TallyEvent& event() const { return event_; } - int& event_nuclide() { return event_nuclide_; } - const int& event_nuclide() const { return event_nuclide_; } - int& event_mt() { return event_mt_; } - int& delayed_group() { return 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(const int& i) { return n_delayed_bank_[i]; } - - int& surface() { return surface_; } - const int& surface() const { return surface_; } - int& cell_born() { return cell_born_; } - const int& cell_born() const { return cell_born_; } - int& material() { return material_; } - const int& material() const { return material_; } - int& material_last() { return material_last_; } - - BoundaryInfo& boundary() { return boundary_; } - - double& sqrtkT() { return sqrtkT_; } - const double& sqrtkT() const { return sqrtkT_; } - double& sqrtkT_last() { return sqrtkT_last_; } - - int& n_collision() { return n_collision_; } - const int& n_collision() const { return n_collision_; } - - bool& write_track() { return write_track_; } - uint64_t& seeds(const int& i) { return seeds_[i]; } - uint64_t* seeds() { return seeds_; } - int& stream() { return stream_; } - - Particle::Bank& secondary_bank(const int& i) { return secondary_bank_[i]; } - decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; } - int64_t& current_work() { return current_work_; } - decltype(flux_derivs_)& flux_derivs() { return flux_derivs_; } - const decltype(flux_derivs_)& flux_derivs() const { return flux_derivs_; } - decltype(filter_matches_)& filter_matches() { return filter_matches_; } - FilterMatch& filter_matches(const int& i) { return filter_matches_[i]; } - decltype(tracks_)& tracks() { return tracks_; } - decltype(nu_bank_)& nu_bank() { return nu_bank_; } - NuBank& nu_bank(const int& i) { return nu_bank_[i]; } - - 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_; } - - bool& trace() { return trace_; } - double& collision_distance() { return collision_distance_; } - int& n_event() { return n_event_; } - -#ifdef DAGMC - moab::DagMC::RayHistory& rayhistory() { return history_; } - Direction& last_dir() { return last_dir_; } -#endif - - 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; } }; //============================================================================ //! Functions //============================================================================ -std::string particle_type_to_str(Particle::Type type); +std::string particle_type_to_str(ParticleType type); -Particle::Type str_to_particle_type(std::string str); +ParticleType str_to_particle_type(std::string str); } // namespace openmc diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h new file mode 100644 index 000000000..8da54dde4 --- /dev/null +++ b/include/openmc/particle_data.h @@ -0,0 +1,451 @@ +#ifndef OPENMC_PARTICLE_REPRESENTATION_H +#define OPENMC_PARTICLE_REPRESENTATION_H + +#include "openmc/constants.h" +#include "openmc/position.h" +#include "openmc/random_lcg.h" +#include "openmc/tallies/filter_match.h" +#include + +namespace openmc { + +//============================================================================== +// Constants +//============================================================================== + +// Since cross section libraries come with different numbers of delayed groups +// (e.g. ENDF/B-VII.1 has 6 and JEFF 3.1.1 has 8 delayed groups) and we don't +// yet know what cross section library is being used when the tallies.xml file +// is read in, we want to have an upper bound on the size of the array we +// use to store the bins for delayed group tallies. +constexpr int MAX_DELAYED_GROUPS {8}; + +constexpr double CACHE_INVALID {-1.0}; + +//========================================================================== +// Aliases and type definitions + +//! Particle types +enum class ParticleType { neutron, photon, electron, positron }; + +//! Saved ("banked") state of a particle +//! NOTE: This structure's MPI type is built in initialize_mpi() of +//! initialize.cpp. Any changes made to the struct here must also be +//! made when building the Bank MPI type in initialize_mpi(). +//! NOTE: This structure is also used on the python side, and is defined +//! in lib/core.py. Changes made to the type here must also be made to the +//! python defintion. +struct ParticleBank { + Position r; + Direction u; + double E; + double wgt; + int delayed_group; + int surf_id; + ParticleType particle; + int64_t parent_id; + int64_t progeny_id; +}; + +//! Saved ("banked") state of a particle, for nu-fission tallying +struct NuBank { + double E; //!< particle energy + double wgt; //!< particle weight + int delayed_group; //!< particle delayed group +}; + +class LocalCoord { +public: + void rotate(const std::vector& rotation); + + //! clear data from a single coordinate level + void reset(); + + Position r; //!< particle position + Direction u; //!< particle direction + int cell {-1}; + int universe {-1}; + int lattice {-1}; + int lattice_x {-1}; + int lattice_y {-1}; + int lattice_z {-1}; + bool rotated {false}; //!< Is the level rotated? +}; + +//============================================================================== +//! Cached microscopic cross sections for a particular nuclide at the current +//! energy +//============================================================================== + +struct NuclideMicroXS { + // Microscopic cross sections in barns + double total; //!< total cross section + double absorption; //!< absorption (disappearance) + double fission; //!< fission + double nu_fission; //!< neutron production from fission + + double elastic; //!< If sab_frac is not 1 or 0, then this value is + //!< averaged over bound and non-bound nuclei + double thermal; //!< Bound thermal elastic & inelastic scattering + double thermal_elastic; //!< Bound thermal elastic scattering + double photon_prod; //!< microscopic photon production xs + + // Cross sections for depletion reactions (note that these are not stored in + // macroscopic cache) + double reaction[DEPLETION_RX.size()]; + + // Indicies and factors needed to compute cross sections from the data tables + int index_grid; //!< Index on nuclide energy grid + int index_temp; //!< Temperature index for nuclide + double interp_factor; //!< Interpolation factor on nuc. energy grid + int index_sab {-1}; //!< Index in sab_tables + int index_temp_sab; //!< Temperature index for sab_tables + double sab_frac; //!< Fraction of atoms affected by S(a,b) + bool use_ptable; //!< In URR range with probability tables? + + // Energy and temperature last used to evaluate these cross sections. If + // these values have changed, then the cross sections must be re-evaluated. + double last_E {0.0}; //!< Last evaluated energy + double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant + //!< * temperature (eV)) +}; + +//============================================================================== +//! Cached microscopic photon cross sections for a particular element at the +//! current energy +//============================================================================== + +struct ElementMicroXS { + int index_grid; //!< index on element energy grid + double last_E {0.0}; //!< last evaluated energy in [eV] + double interp_factor; //!< interpolation factor on energy grid + double total; //!< microscopic total photon xs + double coherent; //!< microscopic coherent xs + double incoherent; //!< microscopic incoherent xs + double photoelectric; //!< microscopic photoelectric xs + double pair_production; //!< microscopic pair production xs +}; + +//============================================================================== +// MacroXS contains cached macroscopic cross sections for the material a +// particle is traveling through +//============================================================================== + +struct MacroXS { + double total; //!< macroscopic total xs + double absorption; //!< macroscopic absorption xs + double fission; //!< macroscopic fission xs + double nu_fission; //!< macroscopic production xs + double photon_prod; //!< macroscopic photon production xs + + // Photon cross sections + double coherent; //!< macroscopic coherent xs + double incoherent; //!< macroscopic incoherent xs + double photoelectric; //!< macroscopic photoelectric xs + double pair_production; //!< macroscopic pair production xs +}; + +//============================================================================== +// Information about nearest boundary crossing +//============================================================================== + +struct BoundaryInfo { + double distance {INFINITY}; //!< distance to nearest boundary + int surface_index {0}; //!< if boundary is surface, index in surfaces vector + int coord_level; //!< coordinate level after crossing boundary + std::array + lattice_translation {}; //!< which way lattice indices will change +}; + +//============================================================================ +//! Defines how particle data is laid out in memory +//============================================================================ + +class ParticleData { + +public: + ParticleData(); + +private: + //========================================================================== + // Data members (accessor methods are below) + + // Cross section caches + std::vector + neutron_xs_; //!< Microscopic neutron cross sections + std::vector photon_xs_; //!< Microscopic photon cross sections + MacroXS macro_xs_; //!< Macroscopic cross sections + + 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 + std::vector coord_; //!< coordinates for all levels + + // Particle coordinates before crossing a surface + int n_coord_last_ {1}; //!< number of current coordinates + std::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_ {0}; //!< 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 + bool alive_ {true}; //!< is particle alive? + + // 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 + double wgt_absorb_ {0.0}; //!< weight absorbed for survival biasing + + // 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 + + // Boundary information + BoundaryInfo boundary_; + + // Temperature of current cell + double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV + double sqrtkT_last_ {0.0}; //!< last temperature + + // Statistical data + int n_collision_ {0}; //!< number of collisions + + // Track output + bool write_track_ {false}; + + // Current PRNG state + uint64_t seeds_[N_STREAMS]; // current seeds + int stream_; // current RNG stream + + // Secondary particle bank + std::vector secondary_bank_; + + int64_t current_work_; // current work index + + std::vector flux_derivs_; // for derivatives for this particle + + std::vector filter_matches_; // tally filter matches + + std::vector> tracks_; // tracks for outputting to file + + std::vector nu_bank_; // bank of most recently fissioned particles + + // 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 + + double collision_distance_; // distance to particle's next closest collision + + int n_event_ {0}; // number of events executed in this particle's history + +// 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 + +public: + //========================================================================== + // Methods and accessors + + NuclideMicroXS& neutron_xs(const int& i) { return neutron_xs_[i]; } + const NuclideMicroXS& neutron_xs(const int& i) const + { + return neutron_xs_[i]; + } + ElementMicroXS& photon_xs(const int& i) { return photon_xs_[i]; } + MacroXS& macro_xs() { return macro_xs_; } + const MacroXS& macro_xs() const { return macro_xs_; } + + int64_t& id() { return id_; } + const int64_t& id() const { return id_; } + 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(const int& i) { return coord_[i]; } + const LocalCoord& coord(const int& i) const { return coord_[i]; } + + int& n_coord_last() { return n_coord_last_; } + const int& n_coord_last() const { return n_coord_last_; } + int& cell_last(const int& i) { return cell_last_[i]; } + const int& cell_last(const int& i) const { return cell_last_[i]; } + + double& E() { return E_; } + const double& E() const { return E_; } + double& E_last() { return E_last_; } + const double& E_last() const { return E_last_; } + int& g() { return g_; } + const int& g() const { return g_; } + int& g_last() { return g_last_; } + const int& g_last() const { return g_last_; } + + double& wgt() { return wgt_; } + double& mu() { return mu_; } + const double& mu() const { return mu_; } + bool& alive() { return alive_; } + + 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_; } + double& wgt_absorb() { return wgt_absorb_; } + const double& wgt_absorb() const { return wgt_absorb_; } + + bool& fission() { return fission_; } + TallyEvent& event() { return event_; } + const TallyEvent& event() const { return event_; } + int& event_nuclide() { return event_nuclide_; } + const int& event_nuclide() const { return event_nuclide_; } + int& event_mt() { return event_mt_; } + int& delayed_group() { return 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(const int& i) { return n_delayed_bank_[i]; } + + int& surface() { return surface_; } + const int& surface() const { return surface_; } + int& cell_born() { return cell_born_; } + const int& cell_born() const { return cell_born_; } + int& material() { return material_; } + const int& material() const { return material_; } + int& material_last() { return material_last_; } + + BoundaryInfo& boundary() { return boundary_; } + + double& sqrtkT() { return sqrtkT_; } + const double& sqrtkT() const { return sqrtkT_; } + double& sqrtkT_last() { return sqrtkT_last_; } + + int& n_collision() { return n_collision_; } + const int& n_collision() const { return n_collision_; } + + bool& write_track() { return write_track_; } + uint64_t& seeds(const int& i) { return seeds_[i]; } + uint64_t* seeds() { return seeds_; } + int& stream() { return stream_; } + + ParticleBank& secondary_bank(const int& i) { return secondary_bank_[i]; } + decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; } + int64_t& current_work() { return current_work_; } + const int64_t& current_work() const { return current_work_; } + double& flux_derivs(const int& i) { return flux_derivs_[i]; } + const double& flux_derivs(const int& i) const { return flux_derivs_[i]; } + decltype(filter_matches_)& filter_matches() { return filter_matches_; } + FilterMatch& filter_matches(const int& i) { return filter_matches_[i]; } + decltype(tracks_)& tracks() { return tracks_; } + decltype(nu_bank_)& nu_bank() { return nu_bank_; } + NuBank& nu_bank(const int& i) { return nu_bank_[i]; } + + 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_; } + + bool& trace() { return trace_; } + double& collision_distance() { return collision_distance_; } + int& n_event() { return n_event_; } + +#ifdef DAGMC + moab::DagMC::RayHistory& rayhistory() { return history_; } + Direction& last_dir() { return last_dir_; } +#endif + + 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_; } + + //! Force recalculation of neutron xs by setting last energy to zero + void invalidate_neutron_xs() + { + for (auto& micro : neutron_xs_) + micro.last_E = 0.0; + } + + //! resets all coordinate levels for the particle + void clear() + { + for (auto& level : coord_) + level.reset(); + n_coord_ = 1; + } + + void zero_delayed_bank() + { + for (int& n : n_delayed_bank_) { + n = 0; + } + } + + void zero_flux_derivs() + { + for (double& d : flux_derivs_) { + d = 0; + } + } +}; + +} // namespace openmc + +#endif // OPENMC_PARTICLE_REPRESENTATION_H diff --git a/include/openmc/physics.h b/include/openmc/physics.h index b6b2fc204..efafa9550 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -81,7 +81,7 @@ Direction sample_cxs_target_velocity(double awr, double E, Direction u, double k uint64_t* seed); void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, - Particle::Bank* site, uint64_t* seed); + ParticleBank* site, uint64_t* seed); //! handles all reactions with a single secondary neutron (other than fission), //! i.e. level scattering, (n,np), (n,na), etc. diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index c842b6f96..b072ae485 100644 --- a/include/openmc/reaction_product.h +++ b/include/openmc/reaction_product.h @@ -45,7 +45,7 @@ public: //! \param[inout] seed Pseudorandom seed pointer void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const; - Particle::Type particle_; //!< Particle type + ParticleType particle_; //!< Particle type EmissionMode emission_mode_; //!< Emission mode double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s] std::unique_ptr yield_; //!< Yield as a function of energy diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 32e2c6730..cd604dcb6 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -69,9 +69,6 @@ void initialize_generation(); //! Full initialization of a particle history void initialize_history(Particle& p, int64_t index_source); -//! Helper function for initialize_history() that is called independently elsewhere -void initialize_history_partial(Particle& p); - //! Finalize a batch //! //! Handles synchronization and accumulation of tallies, calculation of Shannon diff --git a/include/openmc/source.h b/include/openmc/source.h index 12e815699..24c5332e5 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -36,7 +36,7 @@ public: virtual ~Source() = default; // Methods that must be implemented - virtual Particle::Bank sample(uint64_t* seed) const = 0; + virtual ParticleBank sample(uint64_t* seed) const = 0; // Methods that can be overridden virtual double strength() const { return 1.0; } @@ -55,10 +55,10 @@ public: //! Sample from the external source distribution //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - Particle::Bank sample(uint64_t* seed) const override; + ParticleBank sample(uint64_t* seed) const override; // Properties - Particle::Type particle_type() const { return particle_; } + ParticleType particle_type() const { return particle_; } double strength() const override { return strength_; } // Make observing pointers available @@ -67,7 +67,7 @@ public: Distribution* energy() const { return energy_.get(); } private: - Particle::Type particle_ {Particle::Type::neutron}; //!< Type of particle emitted + ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted double strength_ {1.0}; //!< Source strength UPtrSpace space_; //!< Spatial distribution UPtrAngle angle_; //!< Angular distribution @@ -84,10 +84,10 @@ public: explicit FileSource(std::string path); // Methods - Particle::Bank sample(uint64_t* seed) const override; + ParticleBank sample(uint64_t* seed) const override; private: - std::vector sites_; //!< Source sites from a file + std::vector sites_; //!< Source sites from a file }; //============================================================================== @@ -101,7 +101,7 @@ public: ~CustomSourceWrapper(); // Defer implementation to custom source library - Particle::Bank sample(uint64_t* seed) const override + ParticleBank sample(uint64_t* seed) const override { return custom_source_->sample(seed); } @@ -125,7 +125,7 @@ extern "C" void initialize_source(); //! source strength //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled source site -Particle::Bank sample_external_source(uint64_t* seed); +ParticleBank sample_external_source(uint64_t* seed); void free_memory_source(); diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index f93fa00aa..7256192d4 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -15,7 +15,8 @@ void load_state_point(); std::vector calculate_surf_source_size(); void write_source_point(const char* filename, bool surf_source_bank = false); void write_source_bank(hid_t group_id, bool surf_source_bank); -void read_source_bank(hid_t group_id, std::vector& sites, bool distribute); +void read_source_bank( + hid_t group_id, std::vector& sites, bool distribute); void write_tally_results_nr(hid_t file_id); void restart_set_keff(); void write_unstructured_mesh_results(); diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h index c7eb3ce41..e41ff3846 100644 --- a/include/openmc/tallies/filter_particle.h +++ b/include/openmc/tallies/filter_particle.h @@ -37,15 +37,15 @@ public: //---------------------------------------------------------------------------- // Accessors - const std::vector& particles() const { return particles_; } + const std::vector& particles() const { return particles_; } - void set_particles(gsl::span particles); + void set_particles(gsl::span particles); private: //---------------------------------------------------------------------------- // Data members - std::vector particles_; + std::vector particles_; }; } // namespace openmc diff --git a/src/bank.cpp b/src/bank.cpp index 958676d6e..9b77133a1 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -15,15 +15,15 @@ namespace openmc { namespace simulation { -std::vector source_bank; +std::vector source_bank; -SharedArray surf_source_bank; +SharedArray surf_source_bank; // The fission bank is allocated as a SharedArray, rather than a vector, as it will // be shared by all threads in the simulation. It will be allocated to a fixed // maximum capacity in the init_fission_bank() function. Then, Elements will be // added to it by using SharedArray's special thread_safe_append() function. -SharedArray fission_bank; +SharedArray fission_bank; // Each entry in this vector corresponds to the number of progeny produced // this generation for the particle located at that index. This vector is @@ -79,8 +79,8 @@ void sort_fission_bank() // We need a scratch vector to make permutation of the fission bank into // sorted order easy. Under normal usage conditions, the fission bank is // over provisioned, so we can use that as scratch space. - Particle::Bank* sorted_bank; - std::vector sorted_bank_holder; + ParticleBank* sorted_bank; + std::vector sorted_bank_holder; // If there is not enough space, allocate a temporary vector and point to it if (simulation::fission_bank.size() > simulation::fission_bank.capacity() / 2) { diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index f96b60c1e..8fe8ea713 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -31,13 +31,13 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) if (p.material() == MATERIAL_VOID) return; - int photon = static_cast(Particle::Type::photon); + int photon = static_cast(ParticleType::photon); if (p.E() < settings::energy_cutoff[photon]) return; // Get bremsstrahlung data for this material and particle type BremsstrahlungData* mat; - if (p.type() == Particle::Type::positron) { + if (p.type() == ParticleType::positron) { mat = &model::materials[p.material()]->ttb_->positron; } else { mat = &model::materials[p.material()]->ttb_->electron; @@ -110,7 +110,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) if (w > settings::energy_cutoff[photon]) { // Create secondary photon - p.create_secondary(p.wgt(), p.u(), w, Particle::Type::photon); + p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon); *E_lost += w; } } diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index f5d8a9c80..fd1f9bf43 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -137,7 +137,7 @@ void synchronize_bank() // Allocate temporary source bank -- we don't really know how many fission // sites were created, so overallocate by a factor of 3 int64_t index_temp = 0; - std::vector temp_sites(3*simulation::work_per_rank); + std::vector temp_sites(3 * simulation::work_per_rank); for (int64_t i = 0; i < simulation::fission_bank.size(); i++ ) { const auto& site = simulation::fission_bank[i]; diff --git a/src/initialize.cpp b/src/initialize.cpp index e84678feb..1a1e7b2ab 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -132,7 +132,7 @@ void initialize_mpi(MPI_Comm intracomm) mpi::master = (mpi::rank == 0); // Create bank datatype - Particle::Bank b; + ParticleBank b; MPI_Aint disp[9]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); diff --git a/src/material.cpp b/src/material.cpp index 4cc2f7056..89103e4c2 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -753,9 +753,9 @@ void Material::calculate_xs(Particle& p) const p.macro_xs().fission = 0.0; p.macro_xs().nu_fission = 0.0; - if (p.type() == Particle::Type::neutron) { + if (p.type() == ParticleType::neutron) { this->calculate_neutron_xs(p); - } else if (p.type() == Particle::Type::photon) { + } else if (p.type() == ParticleType::photon) { this->calculate_photon_xs(p); } } @@ -763,7 +763,7 @@ void Material::calculate_xs(Particle& p) const void Material::calculate_neutron_xs(Particle& p) const { // Find energy index on energy grid - int neutron = static_cast(Particle::Type::neutron); + int neutron = static_cast(ParticleType::neutron); int i_grid = std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing; diff --git a/src/mesh.cpp b/src/mesh.cpp index 9e974477c..d6cadf12f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -283,10 +283,8 @@ int StructuredMesh::n_surface_bins() const return 4 * n_dimension_ * n_bins(); } -xt::xtensor -StructuredMesh::count_sites(const Particle::Bank* bank, - int64_t length, - bool* outside) const +xt::xtensor StructuredMesh::count_sites( + const ParticleBank* bank, int64_t length, bool* outside) const { // Determine shape of array for counts std::size_t m = this->n_bins(); @@ -989,10 +987,8 @@ void RegularMesh::to_hdf5(hid_t group) const close_group(mesh_group); } -xt::xtensor -RegularMesh::count_sites(const Particle::Bank* bank, - int64_t length, - bool* outside) const +xt::xtensor RegularMesh::count_sites( + const ParticleBank* bank, int64_t length, bool* outside) const { // Determine shape of array for counts std::size_t m = this->n_bins(); diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 31294e3c9..3f2454a73 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -230,7 +230,7 @@ void MgxsInterface::read_header(const std::string& path_cross_sections) void put_mgxs_header_data_to_globals() { // Get the minimum and maximum energies - int neutron = static_cast(Particle::Type::neutron); + int neutron = static_cast(ParticleType::neutron); data::energy_min[neutron] = data::mg.energy_bins_.back(); data::energy_max[neutron] = data::mg.energy_bins_.front(); diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 832677fff..1f4eb9c81 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -327,7 +327,7 @@ void Nuclide::create_derived(const Function1D* prompt_photons, const Function1D* auto xs = xt::adapt(rx->xs_[t].value); for (const auto& p : rx->products_) { - if (p.particle_ == Particle::Type::photon) { + if (p.particle_ == ParticleType::photon) { auto pprod = xt::view(xs_[t], xt::range(j, j+n), XS_PHOTON_PROD); for (int k = 0; k < n; ++k) { double E = grid_[t].energy[k+j]; @@ -445,7 +445,7 @@ void Nuclide::create_derived(const Function1D* prompt_photons, const Function1D* void Nuclide::init_grid() { - int neutron = static_cast(Particle::Type::neutron); + int neutron = static_cast(ParticleType::neutron); double E_min = data::energy_min[neutron]; double E_max = data::energy_max[neutron]; int M = settings::n_log_bins; @@ -494,7 +494,8 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const for (int i = 1; i < rx->products_.size(); ++i) { // Skip any non-neutron products const auto& product = rx->products_[i]; - if (product.particle_ != Particle::Type::neutron) continue; + if (product.particle_ != ParticleType::neutron) + continue; // Evaluate yield if (product.emission_mode_ == EmissionMode::delayed) { diff --git a/src/output.cpp b/src/output.cpp index 458026766..ecc720421 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -147,16 +147,16 @@ void print_particle(Particle& p) { // Display particle type and ID. switch (p.type()) { - case Particle::Type::neutron: + case ParticleType::neutron: fmt::print("Neutron "); break; - case Particle::Type::photon: + case ParticleType::photon: fmt::print("Photon "); break; - case Particle::Type::electron: + case ParticleType::electron: fmt::print("Electron "); break; - case Particle::Type::positron: + case ParticleType::positron: fmt::print("Positron "); break; default: diff --git a/src/particle.cpp b/src/particle.cpp index 2e2413174..bb0a843ea 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -38,9 +38,9 @@ namespace openmc { void LocalCoord::rotate(const std::vector& rotation) { - this->r = this->r.rotate(rotation); - this->u = this->u.rotate(rotation); - this->rotated = true; + r = r.rotate(rotation); + u = u.rotate(rotation); + rotated = true; } void @@ -59,141 +59,118 @@ LocalCoord::reset() // Particle implementation //============================================================================== -Particle::Particle() +void Particle::create_secondary( + double wgt, Direction u, double E, ParticleType type) { - // Create and clear coordinate levels - coord_.resize(model::n_coord_levels); - cell_last_.resize(model::n_coord_levels); - clear(); + secondary_bank().emplace_back(); - for (int& n : n_delayed_bank_) { - n = 0; - } - - // Create microscopic cross section caches - neutron_xs_.resize(data::nuclides.size()); - photon_xs_.resize(data::elements.size()); -} - -void -Particle::clear() -{ - // Reset any coordinate levels - for (auto& level : coord_) level.reset(); - n_coord_ = 1; -} - -void -Particle::create_secondary(double wgt, Direction u, double E, Type type) -{ - secondary_bank_.emplace_back(); - - auto& bank {secondary_bank_.back()}; + auto& bank {secondary_bank().back()}; bank.particle = type; bank.wgt = wgt; - bank.r = this->r(); + bank.r = r(); bank.u = u; - bank.E = settings::run_CE ? E : g_; + bank.E = settings::run_CE ? E : g(); - n_bank_second_ += 1; + n_bank_second() += 1; } -void -Particle::from_source(const Bank* src) +void Particle::from_source(const ParticleBank* src) { // Reset some attributes - this->clear(); - alive_ = true; - surface_ = 0; - cell_born_ = C_NONE; - material_ = C_NONE; - n_collision_ = 0; - fission_ = false; - std::fill(flux_derivs_.begin(), flux_derivs_.end(), 0.0); + clear(); + alive() = true; + surface() = 0; + cell_born() = C_NONE; + material() = C_NONE; + n_collision() = 0; + fission() = false; + zero_flux_derivs(); // Copy attributes from source bank site - type_ = src->particle; - wgt_ = src->wgt; - wgt_last_ = src->wgt; - this->r() = src->r; - this->u() = src->u; - r_last_current_ = src->r; - r_last_ = src->r; - u_last_ = src->u; + type() = src->particle; + wgt() = src->wgt; + wgt_last() = src->wgt; + r() = src->r; + u() = src->u; + r_last_current() = src->r; + r_last() = src->r; + u_last() = src->u; if (settings::run_CE) { - E_ = src->E; - g_ = 0; + E() = src->E; + g() = 0; } else { - g_ = static_cast(src->E); - g_last_ = static_cast(src->E); - E_ = data::mg.energy_bin_avg_[g_]; + g() = static_cast(src->E); + g_last() = static_cast(src->E); + E() = data::mg.energy_bin_avg_[g()]; } - E_last_ = E_; + E_last() = E(); } void Particle::event_calculate_xs() { // Set the random number stream - if (type_ == Particle::Type::neutron) { - stream_ = STREAM_TRACKING; + if (type() == ParticleType::neutron) { + stream() = STREAM_TRACKING; } else { - stream_ = STREAM_PHOTON; + stream() = STREAM_PHOTON; } // Store pre-collision particle properties - wgt_last_ = wgt_; - E_last_ = E_; - u_last_ = this->u(); - r_last_ = this->r(); + wgt_last() = wgt(); + E_last() = E(); + u_last() = u(); + r_last() = r(); // Reset event variables - event_ = TallyEvent::KILL; - event_nuclide_ = NUCLIDE_NONE; - event_mt_ = REACTION_NONE; + event() = TallyEvent::KILL; + event_nuclide() = NUCLIDE_NONE; + event_mt() = REACTION_NONE; // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles - if (coord_[n_coord_ - 1].cell == C_NONE) { + if (coord(n_coord() - 1).cell == C_NONE) { if (!exhaustive_find_cell(*this)) { - this->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 - if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell; + if (cell_born() == C_NONE) + cell_born() = coord(n_coord() - 1).cell; } // Write particle track. - if (write_track_) write_particle_track(*this); + if (write_track()) + write_particle_track(*this); if (settings::check_overlaps) check_cell_overlap(*this); // Calculate microscopic and macroscopic cross sections - if (material_ != MATERIAL_VOID) { + if (material() != MATERIAL_VOID) { if (settings::run_CE) { - if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) { + if (material() != material_last() || sqrtkT() != sqrtkT_last()) { // If the material is the same as the last material and the // temperature hasn't changed, we don't need to lookup cross // sections again. - model::materials[material_]->calculate_xs(*this); + model::materials[material()]->calculate_xs(*this); } } else { // Get the MG data; unlike the CE case above, we have to re-calculate // cross sections for every collision since the cross sections may // be angle-dependent - data::mg.macro_xs_[material_].calculate_xs(*this); + data::mg.macro_xs_[material()].calculate_xs(*this); // Update the particle's group while we know we are multi-group - g_last_ = g_; + g_last() = g(); } } else { - macro_xs_.total = 0.0; - macro_xs_.absorption = 0.0; - macro_xs_.fission = 0.0; - macro_xs_.nu_fission = 0.0; + macro_xs().total = 0.0; + macro_xs().absorption = 0.0; + macro_xs().fission = 0.0; + macro_xs().nu_fission = 0.0; } } @@ -201,24 +178,23 @@ void Particle::event_advance() { // Find the distance to the nearest boundary - boundary_ = distance_to_boundary(*this); + boundary() = distance_to_boundary(*this); // Sample a distance to collision - if (type_ == Particle::Type::electron || - type_ == Particle::Type::positron) { - collision_distance_ = 0.0; - } else if (macro_xs_.total == 0.0) { - collision_distance_ = INFINITY; + if (type() == ParticleType::electron || type() == ParticleType::positron) { + collision_distance() = 0.0; + } else if (macro_xs().total == 0.0) { + collision_distance() = INFINITY; } else { - collision_distance_ = -std::log(prn(this->current_seed())) / macro_xs_.total; + collision_distance() = -std::log(prn(current_seed())) / macro_xs().total; } // Select smaller of the two distances - double distance = std::min(boundary_.distance, collision_distance_); + double distance = std::min(boundary().distance, collision_distance()); // Advance particle - for (int j = 0; j < n_coord_; ++j) { - coord_[j].r += distance * coord_[j].u; + for (int j = 0; j < n_coord(); ++j) { + coord(j).r += distance * coord(j).u; } // Score track-length tallies @@ -228,8 +204,8 @@ Particle::event_advance() // Score track-length estimate of k-eff if (settings::run_mode == RunMode::EIGENVALUE && - type_ == Particle::Type::neutron) { - keff_tally_tracklength_ += wgt_ * distance * macro_xs_.nu_fission; + type() == ParticleType::neutron) { + keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission; } // Score flux derivative accumulators for differential tallies. @@ -242,25 +218,25 @@ void Particle::event_cross_surface() { // Set surface that particle is on and adjust coordinate levels - surface_ = boundary_.surface_index; - n_coord_ = boundary_.coord_level; + surface() = boundary().surface_index; + n_coord() = boundary().coord_level; // Saving previous cell data - for (int j = 0; j < n_coord_; ++j) { - cell_last_[j] = coord_[j].cell; + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell; } - n_coord_last_ = n_coord_; + n_coord_last() = n_coord(); - if (boundary_.lattice_translation[0] != 0 || - boundary_.lattice_translation[1] != 0 || - boundary_.lattice_translation[2] != 0) { + if (boundary().lattice_translation[0] != 0 || + boundary().lattice_translation[1] != 0 || + boundary().lattice_translation[2] != 0) { // Particle crosses lattice boundary - cross_lattice(*this, boundary_); - event_ = TallyEvent::LATTICE; + cross_lattice(*this, boundary()); + event() = TallyEvent::LATTICE; } else { // Particle crosses surface - this->cross_surface(); - event_ = TallyEvent::SURFACE; + cross_surface(); + event() = TallyEvent::SURFACE; } // Score cell to cell partial currents if (!model::active_surface_tallies.empty()) { @@ -273,9 +249,8 @@ Particle::event_collide() { // Score collision estimate of keff if (settings::run_mode == RunMode::EIGENVALUE && - type_ == Particle::Type::neutron) { - keff_tally_collision_ += wgt_ * macro_xs_.nu_fission - / macro_xs_.total; + type() == ParticleType::neutron) { + keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total; } // Score surface current tallies -- this has to be done before the collision @@ -286,7 +261,7 @@ Particle::event_collide() score_surface_tally(*this, model::active_meshsurf_tallies); // Clear surface component - surface_ = 0; + surface() = 0; if (settings::run_CE) { collision(*this); @@ -307,32 +282,32 @@ Particle::event_collide() } // Reset banked weight during collision - n_bank_ = 0; - n_bank_second_ = 0; - wgt_bank_ = 0.0; - for (int& v : n_delayed_bank_) v = 0; + n_bank() = 0; + n_bank_second() = 0; + wgt_bank() = 0.0; + zero_delayed_bank(); // Reset fission logical - fission_ = false; + fission() = false; // Save coordinates for tallying purposes - r_last_current_ = this->r(); + r_last_current() = r(); // Set last material to none since cross sections will need to be // re-evaluated - material_last_ = C_NONE; + material_last() = C_NONE; // Set all directions to base level -- right now, after a collision, only // the base level directions are changed - for (int j = 0; j < n_coord_ - 1; ++j) { - if (coord_[j + 1].rotated) { + for (int j = 0; j < n_coord() - 1; ++j) { + if (coord(j + 1).rotated) { // If next level is rotated, apply rotation matrix - const auto& m {model::cells[coord_[j].cell]->rotation_}; - const auto& u {coord_[j].u}; - coord_[j + 1].u = u.rotate(m); + const auto& m {model::cells[coord(j).cell]->rotation_}; + const auto& u {coord(j).u}; + coord(j + 1).u = u.rotate(m); } else { // Otherwise, copy this level's direction - coord_[j+1].u = coord_[j].u; + coord(j + 1).u = coord(j).u; } } @@ -344,24 +319,26 @@ void Particle::event_revive_from_secondary() { // If particle has too many events, display warning and kill it - ++n_event_; - if (n_event_ == MAX_EVENTS) { - warning("Particle " + std::to_string(id_) + - " underwent maximum number of events."); - alive_ = false; + ++n_event(); + if (n_event() == MAX_EVENTS) { + warning("Particle " + std::to_string(id()) + + " underwent maximum number of events."); + alive() = false; } // Check for secondary particles if this particle is dead - if (!alive_) { + if (!alive()) { // If no secondary particles, break out of event loop - if (secondary_bank_.empty()) return; + if (secondary_bank().empty()) + return; - this->from_source(&secondary_bank_.back()); - secondary_bank_.pop_back(); - n_event_ = 0; + from_source(&secondary_bank().back()); + secondary_bank().pop_back(); + n_event() = 0; // Enter new particle in particle track file - if (write_track_) add_particle_track(*this); + if (write_track()) + add_particle_track(*this); } } @@ -373,32 +350,32 @@ Particle::event_death() #endif // Finish particle track output. - if (write_track_) { + if (write_track()) { write_particle_track(*this); finalize_particle_track(*this); } // Contribute tally reduction variables to global accumulator #pragma omp atomic - global_tally_absorption += keff_tally_absorption_; - #pragma omp atomic - global_tally_collision += keff_tally_collision_; - #pragma omp atomic - global_tally_tracklength += keff_tally_tracklength_; - #pragma omp atomic - global_tally_leakage += keff_tally_leakage_; + global_tally_absorption += keff_tally_absorption(); +#pragma omp atomic + global_tally_collision += keff_tally_collision(); +#pragma omp atomic + global_tally_tracklength += keff_tally_tracklength(); +#pragma omp atomic + global_tally_leakage += keff_tally_leakage(); // Reset particle tallies once accumulated - keff_tally_absorption_ = 0.0; - keff_tally_collision_ = 0.0; - keff_tally_tracklength_ = 0.0; - keff_tally_leakage_ = 0.0; + keff_tally_absorption() = 0.0; + keff_tally_collision() = 0.0; + keff_tally_tracklength() = 0.0; + keff_tally_leakage() = 0.0; // Record the number of progeny created by this particle. // This data will be used to efficiently sort the fission bank. if (settings::run_mode == RunMode::EIGENVALUE) { - int64_t offset = id_ - 1 - simulation::work_index[mpi::rank]; - simulation::progeny_per_particle[offset] = n_progeny_; + int64_t offset = id() - 1 - simulation::work_index[mpi::rank]; + simulation::progeny_per_particle[offset] = n_progeny(); } } @@ -406,24 +383,24 @@ Particle::event_death() void Particle::cross_surface() { - int i_surface = std::abs(surface_); + int i_surface = std::abs(surface()); // TODO: off-by-one const auto& surf {model::surfaces[i_surface - 1].get()}; - if (settings::verbosity >= 10 || trace_) { + if (settings::verbosity >= 10 || trace()) { write_message(1, " Crossing surface {}", surf->id_); } if (surf->surf_source_ && simulation::current_batch == settings::n_batches) { - Particle::Bank site; - site.r = this->r(); - site.u = this->u(); - site.E = this->E_; - site.wgt = this->wgt_; - site.delayed_group = this->delayed_group_; + ParticleBank site; + site.r = r(); + site.u = u(); + site.E = E(); + site.wgt = wgt(); + site.delayed_group = delayed_group(); site.surf_id = surf->id_; - site.particle = this->type_; - site.parent_id = this->id_; - site.progeny_id = this->n_progeny_; + site.particle = type(); + site.parent_id = id(); + site.progeny_id = n_progeny(); int64_t idx = simulation::surf_source_bank.thread_safe_append(site); } @@ -438,18 +415,19 @@ Particle::cross_surface() #ifdef DAGMC if (settings::dagmc) { - auto cellp = dynamic_cast(model::cells[cell_last_[0]].get()); + auto cellp = dynamic_cast(model::cells[cell_last(0)].get()); // TODO: off-by-one - auto surfp = dynamic_cast(model::surfaces[std::abs(surface_) - 1].get()); + auto surfp = + dynamic_cast(model::surfaces[std::abs(surface()) - 1].get()); int32_t i_cell = next_cell(cellp, surfp) - 1; // save material and temp - material_last_ = material_; - sqrtkT_last_ = sqrtkT_; + material_last() = material(); + sqrtkT_last() = sqrtkT(); // set new cell value - coord_[0].cell = i_cell; - cell_instance_ = 0; - material_ = model::cells[i_cell]->material_[0]; - sqrtkT_ = model::cells[i_cell]->sqrtkT_[0]; + coord(0).cell = i_cell; + cell_instance() = 0; + material() = model::cells[i_cell]->material_[0]; + sqrtkT() = model::cells[i_cell]->sqrtkT_[0]; return; } #endif @@ -461,8 +439,8 @@ Particle::cross_surface() // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS // Remove lower coordinate levels and assignment of surface - surface_ = 0; - n_coord_ = 1; + surface() = 0; + n_coord() = 1; bool found = exhaustive_find_cell(*this); if (settings::run_mode != RunMode::PLOTTING && (!found)) { @@ -471,16 +449,16 @@ Particle::cross_surface() // the particle is really traveling tangent to a surface, if we move it // forward a tiny bit it should fix the problem. - n_coord_ = 1; - this->r() += TINY_BIT * this->u(); + n_coord() = 1; + r() += TINY_BIT * u(); // Couldn't find next cell anywhere! This probably means there is an actual // undefined region in the geometry. if (!exhaustive_find_cell(*this)) { - this->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; } } @@ -490,7 +468,7 @@ void Particle::cross_vacuum_bc(const Surface& surf) { // Kill the particle - alive_ = false; + alive() = false; // Score any surface current tallies -- note that the particle is moved // forward slightly so that if the mesh boundary is on the surface, it is @@ -500,15 +478,15 @@ Particle::cross_vacuum_bc(const Surface& surf) // TODO: Find a better solution to score surface currents than // physically moving the particle forward slightly - this->r() += TINY_BIT * this->u(); + r() += TINY_BIT * u(); score_surface_tally(*this, model::active_meshsurf_tallies); } // Score to global leakage tally - keff_tally_leakage_ += wgt_; + keff_tally_leakage() += wgt(); // Display message - if (settings::verbosity >= 10 || trace_) { + if (settings::verbosity >= 10 || trace()) { write_message(1, " Leaked out of surface {}", surf.id_); } } @@ -517,9 +495,9 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u) { // Do not handle reflective boundary conditions on lower universes - if (n_coord_ != 1) { - this->mark_as_lost("Cannot reflect particle " + std::to_string(id_) + - " off surface in a lower universe."); + if (n_coord() != 1) { + mark_as_lost("Cannot reflect particle " + std::to_string(id()) + + " off surface in a lower universe."); return; } @@ -536,36 +514,36 @@ Particle::cross_reflective_bc(const Surface& surf, Direction new_u) if (!model::active_meshsurf_tallies.empty()) { Position r {this->r()}; - this->r() -= TINY_BIT * this->u(); + this->r() -= TINY_BIT * u(); score_surface_tally(*this, model::active_meshsurf_tallies); this->r() = r; } // Set the new particle direction - this->u() = new_u; + u() = new_u; // Reassign particle's cell and surface - coord_[0].cell = cell_last_[n_coord_last_ - 1]; - surface_ = -surface_; + coord(0).cell = cell_last(n_coord_last() - 1); + surface() = -surface(); // If a reflective surface is coincident with a lattice or universe // boundary, it is necessary to redetermine the particle's coordinates in // the lower universes. // (unless we're using a dagmc model, which has exactly one universe) if (!settings::dagmc) { - n_coord_ = 1; + n_coord() = 1; if (!neighbor_list_find_cell(*this)) { - this->mark_as_lost("Couldn't find particle after reflecting from surface " - + std::to_string(surf.id_) + "."); + mark_as_lost("Couldn't find particle after reflecting from surface " + + std::to_string(surf.id_) + "."); return; } } // Set previous coordinate going slightly past surface crossing - r_last_current_ = this->r() + TINY_BIT*this->u(); + r_last_current() = r() + TINY_BIT * u(); // Diagnostic message - if (settings::verbosity >= 10 || trace_) { + if (settings::verbosity >= 10 || trace()) { write_message(1, " Reflected from surface {}", surf.id_); } } @@ -575,8 +553,9 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r, Direction new_u, int new_surface) { // Do not handle periodic boundary conditions on lower universes - if (n_coord_ != 1) { - this->mark_as_lost("Cannot transfer particle " + std::to_string(id_) + + if (n_coord() != 1) { + mark_as_lost( + "Cannot transfer particle " + std::to_string(id()) + " across surface in a lower universe. Boundary conditions must be " "applied to root universe."); return; @@ -587,7 +566,7 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r, // case the surface crossing is coincident with a mesh boundary if (!model::active_meshsurf_tallies.empty()) { Position r {this->r()}; - this->r() -= TINY_BIT * this->u(); + this->r() -= TINY_BIT * u(); score_surface_tally(*this, model::active_meshsurf_tallies); this->r() = r; } @@ -597,23 +576,25 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r, u() = new_u; // Reassign particle's surface - surface_ = new_surface; + surface() = new_surface; // Figure out what cell particle is in now - n_coord_ = 1; + n_coord() = 1; if (!neighbor_list_find_cell(*this)) { - 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."); + 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; } // Set previous coordinate going slightly past surface crossing - r_last_current_ = this->r() + TINY_BIT*this->u(); + r_last_current() = r() + TINY_BIT * u(); // Diagnostic message - if (settings::verbosity >= 10 || trace_) { + if (settings::verbosity >= 10 || trace()) { write_message(1, " Hit periodic boundary on surface {}", surf.id_); } } @@ -626,8 +607,8 @@ Particle::mark_as_lost(const char* message) write_restart(); // Increment number of lost particles - alive_ = false; - #pragma omp atomic + alive() = false; +#pragma omp atomic simulation::n_lost_particles += 1; // Count the total number of simulated particles (on this processor) @@ -650,9 +631,9 @@ Particle::write_restart() const // Set up file name auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output, - simulation::current_batch, id_); + simulation::current_batch, id()); - #pragma omp critical (WriteParticleRestart) +#pragma omp critical (WriteParticleRestart) { // Create file hid_t file_id = file_open(filename, 'w'); @@ -683,10 +664,10 @@ Particle::write_restart() const default: break; } - write_dataset(file_id, "id", id_); - write_dataset(file_id, "type", static_cast(type_)); + write_dataset(file_id, "id", id()); + write_dataset(file_id, "type", static_cast(type())); - int64_t i = current_work_; + int64_t i = current_work(); if (settings::run_mode == RunMode::EIGENVALUE) { //take source data from primary bank for eigenvalue simulation write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt); @@ -711,31 +692,31 @@ Particle::write_restart() const } // #pragma omp critical } -std::string particle_type_to_str(Particle::Type type) +std::string particle_type_to_str(ParticleType type) { switch (type) { - case Particle::Type::neutron: - return "neutron"; - case Particle::Type::photon: - return "photon"; - case Particle::Type::electron: - return "electron"; - case Particle::Type::positron: - return "positron"; + case ParticleType::neutron: + return "neutron"; + case ParticleType::photon: + return "photon"; + case ParticleType::electron: + return "electron"; + case ParticleType::positron: + return "positron"; } UNREACHABLE(); } -Particle::Type str_to_particle_type(std::string str) +ParticleType str_to_particle_type(std::string str) { if (str == "neutron") { - return Particle::Type::neutron; + return ParticleType::neutron; } else if (str == "photon") { - return Particle::Type::photon; + return ParticleType::photon; } else if (str == "electron") { - return Particle::Type::electron; + return ParticleType::electron; } else if (str == "positron") { - return Particle::Type::positron; + return ParticleType::positron; } else { throw std::invalid_argument{fmt::format("Invalid particle name: {}", str)}; } diff --git a/src/particle_data.cpp b/src/particle_data.cpp new file mode 100644 index 000000000..8f4db5a8a --- /dev/null +++ b/src/particle_data.cpp @@ -0,0 +1,35 @@ +#include "openmc/particle_data.h" + +#include "openmc/geometry.h" +#include "openmc/nuclide.h" +#include "openmc/photon.h" +#include "openmc/tallies/derivative.h" +#include "openmc/tallies/filter.h" +#include "openmc/tallies/tally.h" + +namespace openmc { + +ParticleData::ParticleData() +{ + // Create and clear coordinate levels + coord_.resize(model::n_coord_levels); + cell_last_.resize(model::n_coord_levels); + clear(); + + zero_delayed_bank(); + + // Every particle starts with no accumulated flux derivative. + if (!model::active_tallies.empty()) { + flux_derivs_.resize(model::tally_derivs.size()); + zero_flux_derivs(); + } + + // Allocate space for tally filter matches + filter_matches_.resize(model::tally_filters.size()); + + // Create microscopic cross section caches + neutron_xs_.resize(data::nuclides.size()); + photon_xs_.resize(data::elements.size()); +} + +} // namespace openmc diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 9d9be30ef..32b64154a 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -44,7 +44,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) read_dataset(file_id, "id", p.id()); int type; read_dataset(file_id, "type", type); - p.type() = static_cast(type); + p.type() = static_cast(type); read_dataset(file_id, "weight", p.wgt()); read_dataset(file_id, "energy", p.E()); read_dataset(file_id, "xyz", p.r()); @@ -116,16 +116,6 @@ void run_particle_restart() if (p.write_track()) add_particle_track(p); - // Every particle starts with no accumulated flux derivative. - if (!model::active_tallies.empty()) { - p.flux_derivs().resize(model::tally_derivs.size(), 0.0); - std::fill(p.flux_derivs().begin(), p.flux_derivs().end(), 0.0); - } - - // Allocate space for tally filter matches (TODO shouldn't this be in the - // particle constructor, instead?) - p.filter_matches().resize(model::tally_filters.size()); - // Transport neutron transport_history_based_single_particle(p); diff --git a/src/photon.cpp b/src/photon.cpp index 4722a2fac..e1d262302 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -231,7 +231,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) close_group(rgroup); // Truncate the bremsstrahlung data at the cutoff energy - int photon = static_cast(Particle::Type::photon); + int photon = static_cast(ParticleType::photon); const auto& E {electron_energy}; double cutoff = settings::energy_cutoff[photon]; if (cutoff > E(0)) { @@ -669,7 +669,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); double E = shell.binding_energy; - p.create_secondary(p.wgt(), u, E, Particle::Type::photon); + p.create_secondary(p.wgt(), u, E, ParticleType::photon); return; } @@ -701,7 +701,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Non-radiative transition -- Auger/Coster-Kronig effect // Create auger electron - p.create_secondary(p.wgt(), u, E, Particle::Type::electron); + p.create_secondary(p.wgt(), u, E, ParticleType::electron); // Fill hole left by emitted auger electron int i_hole = shell_map_.at(secondary); @@ -711,7 +711,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Radiative transition -- get X-ray energy // Create fluorescent photon - p.create_secondary(p.wgt(), u, E, Particle::Type::photon); + p.create_secondary(p.wgt(), u, E, ParticleType::photon); } // Fill hole created by electron transitioning to the photoelectron hole diff --git a/src/physics.cpp b/src/physics.cpp index 5b58e0aa0..c268c3dc5 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -40,16 +40,16 @@ void collision(Particle& p) // Sample reaction for the material the particle is in switch (p.type()) { - case Particle::Type::neutron: + case ParticleType::neutron: sample_neutron_reaction(p); break; - case Particle::Type::photon: + case ParticleType::photon: sample_photon_reaction(p); break; - case Particle::Type::electron: + case ParticleType::electron: sample_electron_reaction(p); break; - case Particle::Type::positron: + case ParticleType::positron: sample_positron_reaction(p); break; } @@ -66,11 +66,11 @@ void collision(Particle& p) std::string msg; if (p.event() == TallyEvent::KILL) { msg = fmt::format(" Killed. Energy = {} eV.", p.E()); - } else if (p.type() == Particle::Type::neutron) { + } else if (p.type() == ParticleType::neutron) { msg = fmt::format(" {} with {}. Energy = {} eV.", reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_, p.E()); - } else if (p.type() == Particle::Type::photon) { + } else if (p.type() == ParticleType::photon) { msg = fmt::format(" {} with {}. Energy = {} eV.", reaction_name(p.event_mt()), to_element(data::nuclides[p.event_nuclide()]->name_), p.E()); @@ -187,9 +187,9 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) for (int i = 0; i < nu; ++i) { // Initialize fission site object with particle data - Particle::Bank site; + ParticleBank site; site.r = p.r(); - site.particle = Particle::Type::neutron; + site.particle = ParticleType::neutron; site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; @@ -221,7 +221,7 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Write fission particles to nuBank p.nu_bank().emplace_back(); - Particle::NuBank* nu_bank_entry = &p.nu_bank().back(); + NuBank* nu_bank_entry = &p.nu_bank().back(); nu_bank_entry->wgt = site.wgt; nu_bank_entry->E = site.E; nu_bank_entry->delayed_group = site.delayed_group; @@ -251,7 +251,7 @@ void sample_photon_reaction(Particle& p) // Kill photon if below energy cutoff -- an extra check is made here because // photons with energy below the cutoff may have been produced by neutrons // reactions or atomic relaxation - int photon = static_cast(Particle::Type::photon); + int photon = static_cast(ParticleType::photon); if (p.E() < settings::energy_cutoff[photon]) { p.E() = 0.0; p.alive() = false; @@ -300,12 +300,12 @@ void sample_photon_reaction(Particle& p) // Create Compton electron double phi = 2.0*PI*prn(p.current_seed()); double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b; - int electron = static_cast(Particle::Type::electron); + int electron = static_cast(ParticleType::electron); if (E_electron >= settings::energy_cutoff[electron]) { double mu_electron = (alpha - alpha_out*mu) / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu); Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed()); - p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); } // TODO: Compton subshell data does not match atomic relaxation data @@ -366,7 +366,7 @@ void sample_photon_reaction(Particle& p) u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); // Create secondary electron - p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons @@ -391,11 +391,11 @@ void sample_photon_reaction(Particle& p) // Create secondary electron Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed()); - p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); // Create secondary positron u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed()); - p.create_secondary(p.wgt(), u, E_positron, Particle::Type::positron); + p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron); p.event() = TallyEvent::ABSORB; p.event_mt() = PAIR_PROD; @@ -436,8 +436,8 @@ void sample_positron_reaction(Particle& p) u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); // Create annihilation photon pair traveling in opposite directions - p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, Particle::Type::photon); - p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, Particle::Type::photon); + p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon); + p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon); p.E() = 0.0; p.alive() = false; @@ -569,7 +569,7 @@ void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product + f*(rx->xs_[i_temp].value[i_grid - threshold + 1])); for (int j = 0; j < rx->products_.size(); ++j) { - if (rx->products_[j].particle_ == Particle::Type::photon) { + if (rx->products_[j].particle_ == ParticleType::photon) { // For fission, artificially increase the photon yield to account // for delayed photons double f = 1.0; @@ -1014,7 +1014,8 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_ return vt * rotate_angle(u, mu, nullptr, seed); } -void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Particle::Bank* site, uint64_t* seed) +void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, + ParticleBank* site, uint64_t* seed) { // Sample cosine of angle -- fission neutrons are always emitted // isotropically. Sometimes in ACE data, fission reactions actually have @@ -1066,7 +1067,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part rx.products_[group].sample(E_in, site->E, mu, seed); // resample if energy is greater than maximum neutron energy - constexpr int neutron = static_cast(Particle::Type::neutron); + constexpr int neutron = static_cast(ParticleType::neutron); if (site->E < data::energy_max[neutron]) break; // check for large number of resamples @@ -1091,7 +1092,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part rx.products_[0].sample(E_in, site->E, mu, seed); // resample if energy is greater than maximum neutron energy - constexpr int neutron = static_cast(Particle::Type::neutron); + constexpr int neutron = static_cast(ParticleType::neutron); if (site->E < data::energy_max[neutron]) break; // check for large number of resamples @@ -1146,7 +1147,7 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) if (std::floor(yield) == yield) { // If yield is integral, create exactly that many secondary particles for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { - p.create_secondary(p.wgt(), p.u(), p.E(), Particle::Type::neutron); + p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron); } } else { // Otherwise, change weight of particle based on yield @@ -1191,8 +1192,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide) } // Create the secondary photon - p.create_secondary(wgt, u, E, Particle::Type::photon); - + p.create_secondary(wgt, u, E, ParticleType::photon); } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index bd4966dd8..286440838 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -13,6 +13,7 @@ #include "openmc/math_functions.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" +#include "openmc/particle.h" #include "openmc/physics_common.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" @@ -126,9 +127,9 @@ create_fission_sites(Particle& p) for (int i = 0; i < nu; ++i) { // Initialize fission site object with particle data - Particle::Bank site; + ParticleBank site; site.r = p.r(); - site.particle = Particle::Type::neutron; + site.particle = ParticleType::neutron; site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; @@ -179,7 +180,7 @@ create_fission_sites(Particle& p) // Write fission particles to nuBank p.nu_bank().emplace_back(); - Particle::NuBank* nu_bank_entry = &p.nu_bank().back(); + NuBank* nu_bank_entry = &p.nu_bank().back(); nu_bank_entry->wgt = site.wgt; nu_bank_entry->E = site.E; nu_bank_entry->delayed_group = site.delayed_group; diff --git a/src/reaction.cpp b/src/reaction.cpp index 3e4d3b981..0d5a80004 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -73,7 +73,7 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) // mark fission reactions so that we avoid the angle sampling. if (is_fission(mt_)) { for (auto& p : products_) { - if (p.particle_ == Particle::Type::neutron) { + if (p.particle_ == ParticleType::neutron) { for (auto& d : p.distribution_) { auto d_ = dynamic_cast(d.get()); if (d_) d_->fission() = true; diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 6a3b231cc..5d80fb7ba 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -42,7 +42,7 @@ ReactionProduct::ReactionProduct(hid_t group) if (emission_mode_ == EmissionMode::delayed) { if (attribute_exists(group, "decay_rate")) { read_attribute(group, "decay_rate", decay_rate_); - } else if (particle_ == Particle::Type::neutron) { + } else if (particle_ == ParticleType::neutron) { warning(fmt::format("Decay rate doesn't exist for delayed neutron " "emission ({}).", object_name(group))); } diff --git a/src/simulation.cpp b/src/simulation.cpp index d5fdfccf4..5b9fc1c4d 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -511,11 +511,6 @@ void initialize_history(Particle& p, int64_t index_source) #pragma omp atomic simulation::total_weight += p.wgt(); - initialize_history_partial(p); -} - -void initialize_history_partial(Particle& p) -{ // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { p.invalidate_neutron_xs(); @@ -524,16 +519,6 @@ void initialize_history_partial(Particle& p) // Prepare to write out particle track. if (p.write_track()) add_particle_track(p); - - // Every particle starts with no accumulated flux derivative. - if (!model::active_tallies.empty()) - { - p.flux_derivs().resize(model::tally_derivs.size(), 0.0); - std::fill(p.flux_derivs().begin(), p.flux_derivs().end(), 0.0); - } - - // Allocate space for tally filter matches - p.filter_matches().resize(model::tally_filters.size()); } int overall_generation() @@ -573,7 +558,7 @@ void initialize_data() data::energy_min = {0.0, 0.0}; for (const auto& nuc : data::nuclides) { if (nuc->grid_.size() >= 1) { - int neutron = static_cast(Particle::Type::neutron); + int neutron = static_cast(ParticleType::neutron); data::energy_min[neutron] = std::max(data::energy_min[neutron], nuc->grid_[0].energy.front()); data::energy_max[neutron] = std::min(data::energy_max[neutron], @@ -584,7 +569,7 @@ void initialize_data() if (settings::photon_transport) { for (const auto& elem : data::elements) { if (elem->energy_.size() >= 1) { - int photon = static_cast(Particle::Type::photon); + int photon = static_cast(ParticleType::photon); int n = elem->energy_.size(); data::energy_min[photon] = std::max(data::energy_min[photon], std::exp(elem->energy_(1))); @@ -597,7 +582,7 @@ void initialize_data() // Determine if minimum/maximum energy for bremsstrahlung is greater/less // than the current minimum/maximum if (data::ttb_e_grid.size() >= 1) { - int photon = static_cast(Particle::Type::photon); + int photon = static_cast(ParticleType::photon); int n_e = data::ttb_e_grid.size(); data::energy_min[photon] = std::max(data::energy_min[photon], std::exp(data::ttb_e_grid(1))); @@ -613,7 +598,7 @@ void initialize_data() // grid has not been allocated if (nuc->grid_.size() > 0) { double max_E = nuc->grid_[0].energy.back(); - int neutron = static_cast(Particle::Type::neutron); + int neutron = static_cast(ParticleType::neutron); if (max_E == data::energy_max[neutron]) { write_message(7, "Maximum neutron transport energy: {} eV for {}", data::energy_max[neutron], nuc->name_); @@ -630,7 +615,7 @@ void initialize_data() for (auto& nuc : data::nuclides) { nuc->init_grid(); } - int neutron = static_cast(Particle::Type::neutron); + int neutron = static_cast(ParticleType::neutron); simulation::log_spacing = std::log(data::energy_max[neutron] / data::energy_min[neutron]) / settings::n_log_bins; } diff --git a/src/source.cpp b/src/source.cpp index 2bbe7c321..89131f359 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -56,9 +56,9 @@ IndependentSource::IndependentSource(pugi::xml_node node) if (check_for_node(node, "particle")) { auto temp_str = get_node_value(node, "particle", true, true); if (temp_str == "neutron") { - particle_ = Particle::Type::neutron; + particle_ = ParticleType::neutron; } else if (temp_str == "photon") { - particle_ = Particle::Type::photon; + particle_ = ParticleType::photon; settings::photon_transport = true; } else { fatal_error(std::string("Unknown source particle type: ") + temp_str); @@ -141,9 +141,9 @@ IndependentSource::IndependentSource(pugi::xml_node node) } } -Particle::Bank IndependentSource::sample(uint64_t* seed) const +ParticleBank IndependentSource::sample(uint64_t* seed) const { - Particle::Bank site; + ParticleBank site; // Set weight to one by default site.wgt = 1.0; @@ -263,7 +263,7 @@ FileSource::FileSource(std::string path) file_close(file_id); } -Particle::Bank FileSource::sample(uint64_t* seed) const +ParticleBank FileSource::sample(uint64_t* seed) const { size_t i_site = sites_.size()*prn(seed); return sites_[i_site]; @@ -349,7 +349,7 @@ void initialize_source() } } -Particle::Bank sample_external_source(uint64_t* seed) +ParticleBank sample_external_source(uint64_t* seed) { // Determine total source strength double total_strength = 0.0; @@ -368,7 +368,7 @@ Particle::Bank sample_external_source(uint64_t* seed) } // Sample source site from i-th source distribution - Particle::Bank site {model::external_sources[i]->sample(seed)}; + ParticleBank site {model::external_sources[i]->sample(seed)}; // If running in MG, convert site.E to group if (!settings::run_CE) { diff --git a/src/state_point.cpp b/src/state_point.cpp index d1663c0b9..1db9fa38b 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -512,14 +512,17 @@ hid_t h5banktype() { // - openmc/statepoint.py // - docs/source/io_formats/statepoint.rst // - docs/source/io_formats/source.rst - hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct Particle::Bank)); - H5Tinsert(banktype, "r", HOFFSET(Particle::Bank, r), postype); - H5Tinsert(banktype, "u", HOFFSET(Particle::Bank, u), postype); - H5Tinsert(banktype, "E", HOFFSET(Particle::Bank, E), H5T_NATIVE_DOUBLE); - H5Tinsert(banktype, "wgt", HOFFSET(Particle::Bank, wgt), H5T_NATIVE_DOUBLE); - H5Tinsert(banktype, "delayed_group", HOFFSET(Particle::Bank, delayed_group), H5T_NATIVE_INT); - H5Tinsert(banktype, "surf_id", HOFFSET(Particle::Bank, surf_id), H5T_NATIVE_INT); - H5Tinsert(banktype, "particle", HOFFSET(Particle::Bank, particle), H5T_NATIVE_INT); + hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct ParticleBank)); + H5Tinsert(banktype, "r", HOFFSET(ParticleBank, r), postype); + H5Tinsert(banktype, "u", HOFFSET(ParticleBank, u), postype); + H5Tinsert(banktype, "E", HOFFSET(ParticleBank, E), H5T_NATIVE_DOUBLE); + H5Tinsert(banktype, "wgt", HOFFSET(ParticleBank, wgt), H5T_NATIVE_DOUBLE); + H5Tinsert(banktype, "delayed_group", HOFFSET(ParticleBank, delayed_group), + H5T_NATIVE_INT); + H5Tinsert( + banktype, "surf_id", HOFFSET(ParticleBank, surf_id), H5T_NATIVE_INT); + H5Tinsert( + banktype, "particle", HOFFSET(ParticleBank, particle), H5T_NATIVE_INT); H5Tclose(postype); return banktype; @@ -595,9 +598,9 @@ write_source_bank(hid_t group_id, bool surf_source_bank) // Set vectors for source bank and starting bank index of each process std::vector* bank_index = &simulation::work_index; - std::vector* source_bank = &simulation::source_bank; + std::vector* source_bank = &simulation::source_bank; std::vector surf_source_index_vector; - std::vector surf_source_bank_vector; + std::vector surf_source_bank_vector; // Reset dataspace sizes and vectors for surface source bank if (surf_source_bank) { @@ -653,7 +656,8 @@ write_source_bank(hid_t group_id, bool surf_source_bank) // Save source bank sites since the array is overwritten below #ifdef OPENMC_MPI - std::vector temp_source {source_bank->begin(), source_bank->end()}; + std::vector temp_source { + source_bank->begin(), source_bank->end()}; #endif for (int i = 0; i < mpi::n_procs; ++i) { @@ -710,7 +714,8 @@ std::string dtype_member_names(hid_t dtype_id) return names; } -void read_source_bank(hid_t group_id, std::vector& sites, bool distribute) +void read_source_bank( + hid_t group_id, std::vector& sites, bool distribute) { hid_t banktype = h5banktype(); diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 1a950c8d6..00e2d7db1 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -109,7 +109,7 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, // perturbated variable. const auto& deriv {model::tally_derivs[tally.deriv_]}; - const auto flux_deriv = p.flux_derivs()[tally.deriv_]; + const auto flux_deriv = p.flux_derivs(tally.deriv_); // Handle special cases where we know that d_c/d_p must be zero. if (score_bin == SCORE_FLUX) { @@ -558,7 +558,7 @@ score_track_derivative(Particle& p, double distance) for (auto idx = 0; idx < model::tally_derivs.size(); idx++) { const auto& deriv = model::tally_derivs[idx]; - auto& flux_deriv = p.flux_derivs()[idx]; + auto& flux_deriv = p.flux_derivs(idx); if (deriv.diff_material != material.id_) continue; switch (deriv.variable) { @@ -606,7 +606,7 @@ void score_collision_derivative(Particle& p) for (auto idx = 0; idx < model::tally_derivs.size(); idx++) { const auto& deriv = model::tally_derivs[idx]; - auto& flux_deriv = p.flux_derivs()[idx]; + auto& flux_deriv = p.flux_derivs(idx); if (deriv.diff_material != material.id_) continue; diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index eba56421b..b5f0660f6 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -11,16 +11,15 @@ ParticleFilter::from_xml(pugi::xml_node node) { auto particles = get_node_array(node, "bins"); - // Convert to vector of Particle::Type - std::vector types; + // Convert to vector of ParticleType + std::vector types; for (auto& p : particles) { types.push_back(str_to_particle_type(p)); } this->set_particles(types); } -void -ParticleFilter::set_particles(gsl::span particles) +void ParticleFilter::set_particles(gsl::span particles) { // Clear existing particles particles_.clear(); diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index dd67e02ad..a4dc835a9 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -196,7 +196,7 @@ Tally::Tally(pugi::xml_node node) const auto& f = model::tally_filters[particle_filter_index].get(); auto pf = dynamic_cast(f); for (auto p : pf->particles()) { - if (p != Particle::Type::neutron) { + if (p != ParticleType::neutron) { warning(fmt::format("Particle filter other than NEUTRON used with " "photon transport turned off. All tallies for particle type {}" " will have no scores", static_cast(p))); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index c643d014d..c2f51bb9d 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -547,7 +547,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // Get the pre-collision energy of the particle. auto E = p.E_last(); - using Type = Particle::Type; + using Type = ParticleType; for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; @@ -2181,7 +2181,7 @@ void score_analog_tally_ce(Particle& p) // Note that the heating score does NOT use the flux and will be non-zero for // electrons/positrons. double flux = - (p.type() == Particle::Type::neutron || p.type() == Particle::Type::photon) + (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) ? 1.0 : 0.0; @@ -2356,8 +2356,7 @@ void score_collision_tally(Particle& p) { // Determine the collision estimate of the flux double flux = 0.0; - if (p.type() == Particle::Type::neutron || - p.type() == Particle::Type::photon) { + if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) { if (!settings::survival_biasing) { flux = p.wgt_last() / p.macro_xs().total; } else { diff --git a/tests/regression_tests/source_dlopen/source_sampling.cpp b/tests/regression_tests/source_dlopen/source_sampling.cpp index 8eaf6265c..d212f9675 100644 --- a/tests/regression_tests/source_dlopen/source_sampling.cpp +++ b/tests/regression_tests/source_dlopen/source_sampling.cpp @@ -1,17 +1,17 @@ #include #include +#include "openmc/particle_data.h" #include "openmc/random_lcg.h" #include "openmc/source.h" -#include "openmc/particle.h" class CustomSource : public openmc::Source { - openmc::Particle::Bank sample(uint64_t *seed) const + openmc::ParticleBank sample(uint64_t* seed) const { - openmc::Particle::Bank particle; + openmc::ParticleBank particle; // wgt - particle.particle = openmc::Particle::Type::neutron; + particle.particle = openmc::ParticleType::neutron; particle.wgt = 1.0; // position diff --git a/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp b/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp index da7d70e37..ba5958dd5 100644 --- a/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp +++ b/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp @@ -1,16 +1,16 @@ +#include "openmc/particle_data.h" #include "openmc/source.h" -#include "openmc/particle.h" class CustomSource : public openmc::Source { public: CustomSource(double energy) : energy_(energy) { } // Samples from an instance of this class. - openmc::Particle::Bank sample(uint64_t* seed) const + openmc::ParticleBank sample(uint64_t* seed) const { - openmc::Particle::Bank particle; + openmc::ParticleBank particle; // wgt - particle.particle = openmc::Particle::Type::neutron; + particle.particle = openmc::ParticleType::neutron; particle.wgt = 1.0; // position particle.r.x = 0.0;