From 8a8b7b2c0ac9d6c0896cfae560a65476e8cd4df9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Mar 2019 09:12:37 -0500 Subject: [PATCH 01/14] Move cross section caches into Particle class --- include/openmc/material.h | 6 +- include/openmc/nuclide.h | 82 +--------- include/openmc/particle.h | 202 ++++++++++++++++++------- include/openmc/photon.h | 23 +-- include/openmc/physics.h | 11 +- include/openmc/tallies/tally_scoring.h | 6 +- include/openmc/thermal.h | 2 +- src/material.cpp | 49 +++--- src/nuclide.cpp | 124 ++++++++------- src/particle.cpp | 30 ++-- src/particle_restart.cpp | 14 -- src/photon.cpp | 12 +- src/physics.cpp | 111 +++++++------- src/physics_mg.cpp | 18 +-- src/simulation.cpp | 14 -- src/tallies/derivative.cpp | 142 ++++++++--------- src/tallies/tally_scoring.cpp | 184 +++++++++++----------- 17 files changed, 495 insertions(+), 535 deletions(-) diff --git a/include/openmc/material.h b/include/openmc/material.h index 932945ffce..3109f6d787 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -47,7 +47,7 @@ public: explicit Material(pugi::xml_node material_node); // Methods - void calculate_xs(const Particle& p) const; + void calculate_xs(Particle& p) const; //! Assign thermal scattering tables to specific nuclides within the material //! so the code knows when to apply bound thermal scattering data @@ -104,8 +104,8 @@ private: //! Normalize density void normalize_density(); - void calculate_neutron_xs(const Particle& p) const; - void calculate_photon_xs(const Particle& p) const; + void calculate_neutron_xs(Particle& p) const; + void calculate_photon_xs(Particle& p) const; }; //============================================================================== diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index c8b34d08e2..4be898a64e 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -13,6 +13,7 @@ #include "openmc/constants.h" #include "openmc/endf.h" +#include "openmc/particle.h" #include "openmc/reaction.h" #include "openmc/reaction_product.h" #include "openmc/urr.h" @@ -20,69 +21,6 @@ namespace openmc { -//============================================================================== -// Constants -//============================================================================== - -constexpr double CACHE_INVALID {-1.0}; - -//============================================================================== -//! 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)) -}; - -//============================================================================== -// MATERIALMACROXS contains cached macroscopic cross sections for the material a -// particle is traveling through -//============================================================================== - -struct MaterialMacroXS { - 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 -}; - //============================================================================== // Data for a nuclide //============================================================================== @@ -102,14 +40,13 @@ public: //! Initialize logarithmic grid for energy searches void init_grid(); - void calculate_xs(int i_sab, double E, int i_log_union, - double sqrtkT, double sab_frac); + void calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle& p); - void calculate_sab_xs(int i_sab, double E, double sqrtkT, double sab_frac); + void calculate_sab_xs(int i_sab, double sab_frac, Particle& p); // Methods double nu(double E, EmissionMode mode, int group=0) const; - void calculate_elastic_xs() const; + void calculate_elastic_xs(Particle& p) const; //! Determines the microscopic 0K elastic cross section at a trial relative //! energy used in resonance scattering @@ -117,7 +54,7 @@ public: //! \brief Determines cross sections in the unresolved resonance range //! from probability tables. - void calculate_urr_xs(int i_temp, double E) const; + void calculate_urr_xs(int i_temp, Particle& p) const; // Data members std::string name_; //!< Name of nuclide, e.g. "U235" @@ -194,15 +131,6 @@ extern std::unordered_map nuclide_map; } // namespace data -namespace simulation { - -// Cross section caches -extern NuclideMicroXS* micro_xs; -extern MaterialMacroXS material_xs; -#pragma omp threadprivate(micro_xs, material_xs) - -} // namespace simulation - //============================================================================== // Non-member functions //============================================================================== diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 2cb30a44fb..e971c30da3 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -6,9 +6,11 @@ #include #include +#include // for unique_ptr #include #include +#include "openmc/constants.h" #include "openmc/position.h" namespace openmc { @@ -33,6 +35,8 @@ constexpr int MAX_LOST_PARTICLES {10}; // Maximum number of lost particles, relative to the total number of particles constexpr double REL_MAX_LOST_PARTICLES {1.0e-6}; +constexpr double CACHE_INVALID {-1.0}; + //============================================================================== // Class declarations //============================================================================== @@ -52,12 +56,88 @@ struct LocalCoord { void reset(); }; +//============================================================================== +//! 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 +}; + +//============================================================================== +// MATERIALMACROXS contains cached macroscopic cross sections for the material a +// particle is traveling through +//============================================================================== + +struct MaterialMacroXS { + 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 +}; + //============================================================================ //! State of a particle being transported through geometry //============================================================================ class Particle { public: + //========================================================================== + // Aliases and type definitions + //! Particle types enum class Type { neutron, photon, electron, positron @@ -73,9 +153,76 @@ public: Type particle; }; + //========================================================================== // Constructors + Particle(); + //========================================================================== + // Methods and accessors + + // 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; } + + //! resets all coordinate levels for the particle + void clear(); + + //! create a secondary particle + // + //! stores the current phase space attributes of the particle in the + //! secondary bank and increments the number of sites in the secondary bank. + //! \param u Direction of the secondary particle + //! \param E Energy of the secondary particle in [eV] + //! \param type Particle type + void create_secondary(Direction u, double E, Type type); + + //! initialize from a source site + // + //! initializes a particle from data stored in a source site. The source + //! 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); + + //! Transport a particle from birth to death + void transport(); + + //! Cross a surface and handle boundary conditions + void cross_surface(); + + //! mark a particle as lost and create a particle restart file + //! \param message A warning message to display + void mark_as_lost(const char* message); + + void mark_as_lost(const std::string& message) + {mark_as_lost(message.c_str());} + + void mark_as_lost(const std::stringstream& message) + {mark_as_lost(message.str());} + + //! create a particle restart HDF5 file + void write_restart() const; + + //========================================================================== + // Data members + + std::vector micro_xs_; + std::vector micro_photon_xs_; + MaterialMacroXS material_xs_; + int64_t id_; //!< Unique ID Type type_ {Type::neutron}; //!< Particle type (n, p, e, etc.) @@ -139,61 +286,6 @@ public: // Secondary particles created int64_t n_secondary_ {}; Bank secondary_bank_[MAX_SECONDARY]; - - // 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; } - - //! resets all coordinate levels for the particle - void clear(); - - //! create a secondary particle - // - //! stores the current phase space attributes of the particle in the - //! secondary bank and increments the number of sites in the secondary bank. - //! \param u Direction of the secondary particle - //! \param E Energy of the secondary particle in [eV] - //! \param type Particle type - void create_secondary(Direction u, double E, Type type); - - //! initialize from a source site - // - //! initializes a particle from data stored in a source site. The source - //! 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); - - //! Transport a particle from birth to death - void transport(); - - //! Cross a surface and handle boundary conditions - void cross_surface(); - - //! mark a particle as lost and create a particle restart file - //! \param message A warning message to display - void mark_as_lost(const char* message); - - void mark_as_lost(const std::string& message) - {mark_as_lost(message.c_str());} - - void mark_as_lost(const std::stringstream& message) - {mark_as_lost(message.str());} - - //! create a particle restart HDF5 file - void write_restart() const; }; } // namespace openmc diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 7c20a716d6..902d242321 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -42,7 +42,7 @@ public: PhotonInteraction(hid_t group, int i_element); // Methods - void calculate_xs(double E) const; + void calculate_xs(Particle& p) const; void compton_scatter(double alpha, bool doppler, double* alpha_out, double* mu, int* i_shell) const; @@ -98,22 +98,6 @@ private: void compton_doppler(double alpha, double mu, double* E_out, int* i_shell) const; }; -//============================================================================== -//! 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 -}; - //============================================================================== // Non-member functions //============================================================================== @@ -136,11 +120,6 @@ extern std::unordered_map element_map; } // namespace data -namespace simulation { -extern ElementMicroXS* micro_photon_xs; -#pragma omp threadprivate(micro_photon_xs) -} // namespace simulation - } // namespace openmc #endif // OPENMC_PHOTON_H diff --git a/include/openmc/physics.h b/include/openmc/physics.h index a7dd20aad3..5c8d64f5f9 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -51,20 +51,19 @@ void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, int sample_element(Particle* p); -Reaction* sample_fission(int i_nuclide, double E); +Reaction* sample_fission(int i_nuclide, const Particle* p); -void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product); +void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_product); void absorption(Particle* p, int i_nuclide); void scatter(Particle*, int i_nuclide); //! Treats the elastic scattering of a neutron with a target. -void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E, - Direction& u, double& mu_lab); +void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, + Particle* p); -void sab_scatter(int i_nuclide, int i_sab, double& E, - Direction& u, double& mu); +void sab_scatter(int i_nuclide, int i_sab, Particle* p); //! samples the target velocity. The constant cross section free gas model is //! the default method. Methods for correctly accounting for the energy diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index af03b1aa23..2c5a29c176 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -59,14 +59,14 @@ private: //! since collisions do not occur in voids. // //! \param p The particle being tracked -void score_collision_tally(const Particle* p); +void score_collision_tally(Particle* p); //! Score tallies based on a simple count of events (for continuous energy). // //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked -void score_analog_tally_ce(const Particle* p); +void score_analog_tally_ce(Particle* p); //! Score tallies based on a simple count of events (for multigroup). // @@ -83,7 +83,7 @@ void score_analog_tally_mg(const Particle* p); // //! \param p The particle being tracked //! \param distance The distance in [cm] traveled by the particle -void score_tracklength_tally(const Particle* p, double distance); +void score_tracklength_tally(Particle* p, double distance); //! Score surface or mesh-surface tallies for particle currents. // diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index 624cf31bd6..b737f005b3 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -10,7 +10,7 @@ #include "xtensor/xtensor.hpp" #include "openmc/hdf5_interface.h" -#include "openmc/nuclide.h" +#include "openmc/particle.h" namespace openmc { diff --git a/src/material.cpp b/src/material.cpp index da68a536dd..0d2fe4bbc3 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -729,13 +729,13 @@ void Material::init_nuclide_index() } } -void Material::calculate_xs(const Particle& p) const +void Material::calculate_xs(Particle& p) const { // Set all material macroscopic cross sections to zero - simulation::material_xs.total = 0.0; - simulation::material_xs.absorption = 0.0; - simulation::material_xs.fission = 0.0; - simulation::material_xs.nu_fission = 0.0; + p.material_xs_.total = 0.0; + p.material_xs_.absorption = 0.0; + p.material_xs_.fission = 0.0; + p.material_xs_.nu_fission = 0.0; if (p.type_ == Particle::Type::neutron) { this->calculate_neutron_xs(p); @@ -744,7 +744,7 @@ void Material::calculate_xs(const Particle& p) const } } -void Material::calculate_neutron_xs(const Particle& p) const +void Material::calculate_neutron_xs(Particle& p) const { // Find energy index on energy grid int neutron = static_cast(Particle::Type::neutron); @@ -792,13 +792,12 @@ void Material::calculate_neutron_xs(const Particle& p) const int i_nuclide = nuclide_[i]; // Calculate microscopic cross section for this nuclide - const auto& micro {simulation::micro_xs[i_nuclide]}; + const auto& micro {p.micro_xs_[i_nuclide]}; if (p.E_ != micro.last_E || p.sqrtkT_ != micro.last_sqrtkT || i_sab != micro.index_sab || sab_frac != micro.sab_frac) { - data::nuclides[i_nuclide]->calculate_xs(i_sab, p.E_, i_grid, - p.sqrtkT_, sab_frac); + data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p); } // ====================================================================== @@ -808,19 +807,19 @@ void Material::calculate_neutron_xs(const Particle& p) const double atom_density = atom_density_(i); // Add contributions to cross sections - simulation::material_xs.total += atom_density * micro.total; - simulation::material_xs.absorption += atom_density * micro.absorption; - simulation::material_xs.fission += atom_density * micro.fission; - simulation::material_xs.nu_fission += atom_density * micro.nu_fission; + p.material_xs_.total += atom_density * micro.total; + p.material_xs_.absorption += atom_density * micro.absorption; + p.material_xs_.fission += atom_density * micro.fission; + p.material_xs_.nu_fission += atom_density * micro.nu_fission; } } -void Material::calculate_photon_xs(const Particle& p) const +void Material::calculate_photon_xs(Particle& p) const { - simulation::material_xs.coherent = 0.0; - simulation::material_xs.incoherent = 0.0; - simulation::material_xs.photoelectric = 0.0; - simulation::material_xs.pair_production = 0.0; + p.material_xs_.coherent = 0.0; + p.material_xs_.incoherent = 0.0; + p.material_xs_.photoelectric = 0.0; + p.material_xs_.pair_production = 0.0; // Add contribution from each nuclide in material for (int i = 0; i < nuclide_.size(); ++i) { @@ -831,9 +830,9 @@ void Material::calculate_photon_xs(const Particle& p) const int i_element = element_[i]; // Calculate microscopic cross section for this nuclide - const auto& micro {simulation::micro_photon_xs[i_element]}; + const auto& micro {p.micro_photon_xs_[i_element]}; if (p.E_ != micro.last_E) { - data::elements[i_element].calculate_xs(p.E_); + data::elements[i_element].calculate_xs(p); } // ======================================================================== @@ -843,11 +842,11 @@ void Material::calculate_photon_xs(const Particle& p) const double atom_density = atom_density_(i); // Add contributions to material macroscopic cross sections - simulation::material_xs.total += atom_density * micro.total; - simulation::material_xs.coherent += atom_density * micro.coherent; - simulation::material_xs.incoherent += atom_density * micro.incoherent; - simulation::material_xs.photoelectric += atom_density * micro.photoelectric; - simulation::material_xs.pair_production += atom_density * micro.pair_production; + p.material_xs_.total += atom_density * micro.total; + p.material_xs_.coherent += atom_density * micro.coherent; + p.material_xs_.incoherent += atom_density * micro.incoherent; + p.material_xs_.photoelectric += atom_density * micro.photoelectric; + p.material_xs_.pair_production += atom_density * micro.pair_production; } } diff --git a/src/nuclide.cpp b/src/nuclide.cpp index cd2461186a..87965dfc3e 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -34,11 +34,6 @@ std::vector> nuclides; std::unordered_map nuclide_map; } // namespace data -namespace simulation { -NuclideMicroXS* micro_xs; -MaterialMacroXS material_xs; -} // namespace simulation - //============================================================================== // Nuclide implementation //============================================================================== @@ -461,10 +456,10 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const } } -void Nuclide::calculate_elastic_xs() const +void Nuclide::calculate_elastic_xs(Particle& p) const { // Get temperature index, grid index, and interpolation factor - auto& micro = simulation::micro_xs[i_nuclide_]; + auto& micro {p.micro_xs_[i_nuclide_]}; int i_temp = micro.index_temp; int i_grid = micro.index_grid; double f = micro.interp_factor; @@ -498,42 +493,41 @@ double Nuclide::elastic_xs_0K(double E) const return (1.0 - f)*elastic_0K_[i_grid] + f*elastic_0K_[i_grid + 1]; } -void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, - double sqrtkT, double sab_frac) +void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle& p) { - auto& micro_xs = simulation::micro_xs[i_nuclide_]; + auto& micro {p.micro_xs_[i_nuclide_]}; // Initialize cached cross sections to zero - micro_xs.elastic = CACHE_INVALID; - micro_xs.thermal = 0.0; - micro_xs.thermal_elastic = 0.0; + micro.elastic = CACHE_INVALID; + micro.thermal = 0.0; + micro.thermal_elastic = 0.0; // Check to see if there is multipole data present at this energy bool use_mp = false; if (multipole_) { - use_mp = (E >= multipole_->E_min_ && E <= multipole_->E_max_); + use_mp = (p.E_ >= multipole_->E_min_ && p.E_ <= multipole_->E_max_); } // Evaluate multipole or interpolate if (use_mp) { // Call multipole kernel double sig_s, sig_a, sig_f; - std::tie(sig_s, sig_a, sig_f) = multipole_->evaluate(E, sqrtkT); + std::tie(sig_s, sig_a, sig_f) = multipole_->evaluate(p.E_, p.sqrtkT_); - micro_xs.total = sig_s + sig_a; - micro_xs.elastic = sig_s; - micro_xs.absorption = sig_a; - micro_xs.fission = sig_f; - micro_xs.nu_fission = fissionable_ ? - sig_f * this->nu(E, EmissionMode::total) : 0.0; + micro.total = sig_s + sig_a; + micro.elastic = sig_s; + micro.absorption = sig_a; + micro.fission = sig_f; + micro.nu_fission = fissionable_ ? + sig_f * this->nu(p.E_, EmissionMode::total) : 0.0; if (simulation::need_depletion_rx) { // Only non-zero reaction is (n,gamma) - micro_xs.reaction[0] = sig_a - sig_f; + micro.reaction[0] = sig_a - sig_f; // Set all other reaction cross sections to zero for (int i = 1; i < DEPLETION_RX.size(); ++i) { - micro_xs.reaction[i] = 0.0; + micro.reaction[i] = 0.0; } } @@ -547,13 +541,13 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, // resonance range, so the value here does not matter. index_temp is // set to -1 to force a segfault in case a developer messes up and tries // to use it with multipole. - micro_xs.index_temp = -1; - micro_xs.index_grid = -1; - micro_xs.interp_factor = 0.0; + micro.index_temp = -1; + micro.index_grid = -1; + micro.interp_factor = 0.0; } else { // Find the appropriate temperature index. - double kT = sqrtkT*sqrtkT; + double kT = p.sqrtkT_*p.sqrtkT_; double f; int i_temp = -1; switch (settings::temperature_method) { @@ -590,9 +584,9 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, const auto& xs {xs_[i_temp]}; int i_grid; - if (E < grid.energy.front()) { + if (p.E_ < grid.energy.front()) { i_grid = 0; - } else if (E > grid.energy.back()) { + } else if (p.E_ > grid.energy.back()) { i_grid = grid.energy.size() - 2; } else { // Determine bounding indices based on which equal log-spaced @@ -601,49 +595,49 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, int i_high = grid.grid_index[i_log_union + 1] + 1; // Perform binary search over reduced range - i_grid = i_low + lower_bound_index(&grid.energy[i_low], &grid.energy[i_high], E); + i_grid = i_low + lower_bound_index(&grid.energy[i_low], &grid.energy[i_high], p.E_); } // check for rare case where two energy points are the same if (grid.energy[i_grid] == grid.energy[i_grid + 1]) ++i_grid; // calculate interpolation factor - f = (E - grid.energy[i_grid]) / + f = (p.E_ - grid.energy[i_grid]) / (grid.energy[i_grid + 1]- grid.energy[i_grid]); - micro_xs.index_temp = i_temp; - micro_xs.index_grid = i_grid; - micro_xs.interp_factor = f; + micro.index_temp = i_temp; + micro.index_grid = i_grid; + micro.interp_factor = f; // Calculate microscopic nuclide total cross section - micro_xs.total = (1.0 - f)*xs(i_grid, XS_TOTAL) + micro.total = (1.0 - f)*xs(i_grid, XS_TOTAL) + f*xs(i_grid + 1, XS_TOTAL); // Calculate microscopic nuclide absorption cross section - micro_xs.absorption = (1.0 - f)*xs(i_grid, XS_ABSORPTION) + micro.absorption = (1.0 - f)*xs(i_grid, XS_ABSORPTION) + f*xs(i_grid + 1, XS_ABSORPTION); if (fissionable_) { // Calculate microscopic nuclide total cross section - micro_xs.fission = (1.0 - f)*xs(i_grid, XS_FISSION) + micro.fission = (1.0 - f)*xs(i_grid, XS_FISSION) + f*xs(i_grid + 1, XS_FISSION); // Calculate microscopic nuclide nu-fission cross section - micro_xs.nu_fission = (1.0 - f)*xs(i_grid, XS_NU_FISSION) + micro.nu_fission = (1.0 - f)*xs(i_grid, XS_NU_FISSION) + f*xs(i_grid + 1, XS_NU_FISSION); } else { - micro_xs.fission = 0.0; - micro_xs.nu_fission = 0.0; + micro.fission = 0.0; + micro.nu_fission = 0.0; } // Calculate microscopic nuclide photon production cross section - micro_xs.photon_prod = (1.0 - f)*xs(i_grid, XS_PHOTON_PROD) + micro.photon_prod = (1.0 - f)*xs(i_grid, XS_PHOTON_PROD) + f*xs(i_grid + 1, XS_PHOTON_PROD); // Depletion-related reactions if (simulation::need_depletion_rx) { // Initialize all reaction cross sections to zero - for (double& xs_i : micro_xs.reaction) { + for (double& xs_i : micro.reaction) { xs_i = 0.0; } @@ -658,14 +652,14 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, // Physics says that (n,gamma) is not a threshold reaction, so we don't // need to specifically check its threshold index if (j == 0) { - micro_xs.reaction[0] = (1.0 - f)*rx_xs[i_grid] + micro.reaction[0] = (1.0 - f)*rx_xs[i_grid] + f*rx_xs[i_grid + 1]; continue; } int threshold = rx->xs_[i_temp].threshold; if (i_grid >= threshold) { - micro_xs.reaction[j] = (1.0 - f)*rx_xs[i_grid - threshold] + + micro.reaction[j] = (1.0 - f)*rx_xs[i_grid - threshold] + f*rx_xs[i_grid - threshold + 1]; } else if (j >= 3) { // One can show that the the threshold for (n,(x+1)n) is always @@ -680,35 +674,35 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, } // Initialize sab treatment to false - micro_xs.index_sab = C_NONE; - micro_xs.sab_frac = 0.0; + micro.index_sab = C_NONE; + micro.sab_frac = 0.0; // Initialize URR probability table treatment to false - micro_xs.use_ptable = false; + micro.use_ptable = false; // If there is S(a,b) data for this nuclide, we need to set the sab_scatter // and sab_elastic cross sections and correct the total and elastic cross // sections. - if (i_sab >= 0) this->calculate_sab_xs(i_sab, E, sqrtkT, sab_frac); + if (i_sab >= 0) this->calculate_sab_xs(i_sab, sab_frac, p); // If the particle is in the unresolved resonance range and there are // probability tables, we need to determine cross sections from the table if (settings::urr_ptables_on && urr_present_ && !use_mp) { - int n = urr_data_[micro_xs.index_temp].n_energy_; - if ((E > urr_data_[micro_xs.index_temp].energy_(0)) && - (E < urr_data_[micro_xs.index_temp].energy_(n-1))) { - this->calculate_urr_xs(micro_xs.index_temp, E); + int n = urr_data_[micro.index_temp].n_energy_; + if ((p.E_ > urr_data_[micro.index_temp].energy_(0)) && + (p.E_ < urr_data_[micro.index_temp].energy_(n-1))) { + this->calculate_urr_xs(micro.index_temp, p); } } - micro_xs.last_E = E; - micro_xs.last_sqrtkT = sqrtkT; + micro.last_E = p.E_; + micro.last_sqrtkT = p.sqrtkT_; } -void Nuclide::calculate_sab_xs(int i_sab, double E, double sqrtkT, double sab_frac) +void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p) { - auto& micro {simulation::micro_xs[i_nuclide_]}; + auto& micro {p.micro_xs_[i_nuclide_]}; // Set flag that S(a,b) treatment should be used for scattering micro.index_sab = i_sab; @@ -717,14 +711,14 @@ void Nuclide::calculate_sab_xs(int i_sab, double E, double sqrtkT, double sab_fr int i_temp; double elastic; double inelastic; - data::thermal_scatt[i_sab]->calculate_xs(E, sqrtkT, &i_temp, &elastic, &inelastic); + data::thermal_scatt[i_sab]->calculate_xs(p.E_, p.sqrtkT_, &i_temp, &elastic, &inelastic); // Store the S(a,b) cross sections. micro.thermal = sab_frac * (elastic + inelastic); micro.thermal_elastic = sab_frac * elastic; // Calculate free atom elastic cross section - this->calculate_elastic_xs(); + this->calculate_elastic_xs(p); // Correct total and elastic cross sections micro.total = micro.total + micro.thermal - sab_frac*micro.elastic; @@ -735,9 +729,9 @@ void Nuclide::calculate_sab_xs(int i_sab, double E, double sqrtkT, double sab_fr micro.sab_frac = sab_frac; } -void Nuclide::calculate_urr_xs(int i_temp, double E) const +void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const { - auto& micro = simulation::micro_xs[i_nuclide_]; + auto& micro = p.micro_xs_[i_nuclide_]; micro.use_ptable = true; // Create a shorthand for the URR data @@ -745,7 +739,7 @@ void Nuclide::calculate_urr_xs(int i_temp, double E) const // Determine the energy table int i_energy = 0; - while (E >= urr.energy_(i_energy + 1)) {++i_energy;}; + while (p.E_ >= urr.energy_(i_energy + 1)) {++i_energy;}; // Sample the probability table using the cumulative distribution @@ -773,7 +767,7 @@ void Nuclide::calculate_urr_xs(int i_temp, double E) const double f; if (urr.interp_ == Interpolation::lin_lin) { // Determine the interpolation factor on the table - f = (E - urr.energy_(i_energy)) / + f = (p.E_ - urr.energy_(i_energy)) / (urr.energy_(i_energy + 1) - urr.energy_(i_energy)); elastic = (1. - f) * urr.prob_(i_energy, URR_ELASTIC, i_low) + @@ -784,7 +778,7 @@ void Nuclide::calculate_urr_xs(int i_temp, double E) const f * urr.prob_(i_energy + 1, URR_N_GAMMA, i_up); } else if (urr.interp_ == Interpolation::log_log) { // Determine interpolation factor on the table - f = std::log(E / urr.energy_(i_energy)) / + f = std::log(p.E_ / urr.energy_(i_energy)) / std::log(urr.energy_(i_energy + 1) / urr.energy_(i_energy)); // Calculate the elastic cross section/factor @@ -838,7 +832,7 @@ void Nuclide::calculate_urr_xs(int i_temp, double E) const // Multiply by smooth cross-section if needed if (urr.multiply_smooth_) { - calculate_elastic_xs(); + calculate_elastic_xs(p); elastic *= micro.elastic; capture *= (micro.absorption - micro.fission); fission *= micro.fission; @@ -858,7 +852,7 @@ void Nuclide::calculate_urr_xs(int i_temp, double E) const // Determine nu-fission cross-section if (fissionable_) { - micro.nu_fission = nu(E, EmissionMode::total) * micro.fission; + micro.nu_fission = nu(p.E_, EmissionMode::total) * micro.fission; } } diff --git a/src/particle.cpp b/src/particle.cpp index ca7fccfe78..f5084cd1c2 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -15,6 +15,7 @@ #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" +#include "openmc/photon.h" #include "openmc/physics.h" #include "openmc/physics_mg.h" #include "openmc/random_lcg.h" @@ -55,6 +56,10 @@ Particle::Particle() for (int& n : n_delayed_bank_) { n = 0; } + + // Create microscopic cross section caches + micro_xs_.resize(data::nuclides.size()); + micro_photon_xs_.resize(data::elements.size()); } void @@ -130,9 +135,7 @@ Particle::transport() // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { - for (int i = 0; i < data::nuclides.size(); ++i) { - simulation::micro_xs[i].last_E = 0.0; - } + for (auto& micro : micro_xs_) micro.last_E = 0.0; } // Prepare to write out particle track. @@ -186,18 +189,17 @@ Particle::transport() } else { // Get the MG data calculate_xs_c(material_, g_, sqrtkT_, this->u_local(), - simulation::material_xs.total, simulation::material_xs.absorption, - simulation::material_xs.nu_fission); + material_xs_.total, material_xs_.absorption, material_xs_.nu_fission); // Finally, update the particle group while we have already checked // for if multi-group g_last_ = g_; } } else { - simulation::material_xs.total = 0.0; - simulation::material_xs.absorption = 0.0; - simulation::material_xs.fission = 0.0; - simulation::material_xs.nu_fission = 0.0; + material_xs_.total = 0.0; + material_xs_.absorption = 0.0; + material_xs_.fission = 0.0; + material_xs_.nu_fission = 0.0; } // Find the distance to the nearest boundary @@ -213,10 +215,10 @@ Particle::transport() if (type_ == Particle::Type::electron || type_ == Particle::Type::positron) { d_collision = 0.0; - } else if (simulation::material_xs.total == 0.0) { + } else if (material_xs_.total == 0.0) { d_collision = INFINITY; } else { - d_collision = -std::log(prn()) / simulation::material_xs.total; + d_collision = -std::log(prn()) / material_xs_.total; } // Select smaller of the two distances @@ -235,7 +237,7 @@ Particle::transport() // Score track-length estimate of k-eff if (settings::run_mode == RUN_MODE_EIGENVALUE && type_ == Particle::Type::neutron) { - global_tally_tracklength += wgt_ * distance * simulation::material_xs.nu_fission; + global_tally_tracklength += wgt_ * distance * material_xs_.nu_fission; } // Score flux derivative accumulators for differential tallies. @@ -278,8 +280,8 @@ Particle::transport() // Score collision estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE && type_ == Particle::Type::neutron) { - global_tally_collision += wgt_ * simulation::material_xs.nu_fission - / simulation::material_xs.total; + global_tally_collision += wgt_ * material_xs_.nu_fission + / material_xs_.total; } // Score surface current tallies -- this has to be done before the collision diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 0cface319a..0bc677cfd8 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -71,13 +71,6 @@ void run_particle_restart() // Set verbosity high settings::verbosity = 10; - // Create cross section caches - #pragma omp parallel - { - simulation::micro_xs = new NuclideMicroXS[data::nuclides.size()]; - simulation::micro_photon_xs = new ElementMicroXS[data::elements.size()]; - } - // Initialize the particle to be tracked Particle p; @@ -105,13 +98,6 @@ void run_particle_restart() // Write output if particle made it print_particle(&p); - - // Clear cross section caches - #pragma omp parallel - { - delete[] simulation::micro_xs; - delete[] simulation::micro_photon_xs; - } } } // namespace openmc diff --git a/src/photon.cpp b/src/photon.cpp index d977d62ac3..9841e9e3bf 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -32,10 +32,6 @@ std::unordered_map element_map; } // namespace data -namespace simulation { -ElementMicroXS* micro_photon_xs; -} // namespace simulation - //============================================================================== // PhotonInteraction implementation //============================================================================== @@ -433,12 +429,12 @@ void PhotonInteraction::compton_doppler(double alpha, double mu, *i_shell = shell; } -void PhotonInteraction::calculate_xs(double E) const +void PhotonInteraction::calculate_xs(Particle& p) const { // Perform binary search on the element energy grid in order to determine // which points to interpolate between int n_grid = energy_.size(); - double log_E = std::log(E); + double log_E = std::log(p.E_); int i_grid; if (log_E <= energy_[0]) { i_grid = 0; @@ -456,7 +452,7 @@ void PhotonInteraction::calculate_xs(double E) const // calculate interpolation factor double f = (log_E - energy_(i_grid)) / (energy_(i_grid+1) - energy_(i_grid)); - auto& xs {simulation::micro_photon_xs[i_element_]}; + auto& xs {p.micro_photon_xs_[i_element_]}; xs.index_grid = i_grid; xs.interp_factor = f; @@ -490,7 +486,7 @@ void PhotonInteraction::calculate_xs(double E) const // Calculate microscopic total cross section xs.total = xs.coherent + xs.incoherent + xs.photoelectric + xs.pair_production; - xs.last_E = E; + xs.last_E = p.E_; } double PhotonInteraction::rayleigh_scatter(double alpha) const diff --git a/src/physics.cpp b/src/physics.cpp index 858aae8f4b..029ab18c86 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -89,7 +89,7 @@ void sample_neutron_reaction(Particle* p) const auto& nuc {data::nuclides[i_nuclide]}; if (nuc->fissionable_) { - Reaction* rx = sample_fission(i_nuclide, p->E_); + Reaction* rx = sample_fission(i_nuclide, p); if (settings::run_mode == RUN_MODE_EIGENVALUE) { create_fission_sites(p, i_nuclide, rx, simulation::fission_bank.data(), &simulation::n_bank, simulation::fission_bank.size()); @@ -110,7 +110,7 @@ void sample_neutron_reaction(Particle* p) // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs - if (simulation::micro_xs[i_nuclide].absorption > 0.0) { + if (p->micro_xs_[i_nuclide].absorption > 0.0) { absorption(p, i_nuclide); } else { p->wgt_absorb_ = 0.0; @@ -146,8 +146,8 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p->wgt_ / simulation::keff * weight * simulation::micro_xs[ - i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].total; + double nu_t = p->wgt_ / simulation::keff * weight * p->micro_xs_[ + i_nuclide].nu_fission / p->micro_xs_[i_nuclide].total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -227,7 +227,7 @@ void sample_photon_reaction(Particle* p) // Sample element within material int i_element = sample_element(p); p->event_nuclide_ = i_element; - const auto& micro {simulation::micro_photon_xs[i_element]}; + const auto& micro {p->micro_photon_xs_[i_element]}; const auto& element {data::elements[i_element]}; // Calculate photon energy over electron rest mass equivalent @@ -408,7 +408,7 @@ void sample_positron_reaction(Particle* p) int sample_nuclide(const Particle* p) { // Sample cumulative distribution function - double cutoff = prn() * simulation::material_xs.total; + double cutoff = prn() * p->material_xs_.total; // Get pointers to nuclide/density arrays const auto& mat {model::materials[p->material_]}; @@ -421,7 +421,7 @@ int sample_nuclide(const Particle* p) double atom_density = mat->atom_density_[i]; // Increment probability to compare to cutoff - prob += atom_density * simulation::micro_xs[i_nuclide].total; + prob += atom_density * p->micro_xs_[i_nuclide].total; if (prob >= cutoff) return i_nuclide; } @@ -433,7 +433,7 @@ int sample_nuclide(const Particle* p) int sample_element(Particle* p) { // Sample cumulative distribution function - double cutoff = prn() * simulation::material_xs.total; + double cutoff = prn() * p->material_xs_.total; // Get pointers to elements, densities const auto& mat {model::materials[p->material_]}; @@ -454,7 +454,7 @@ int sample_element(Particle* p) double atom_density = mat->atom_density_[i]; // Determine microscopic cross section - double sigma = atom_density * simulation::micro_photon_xs[i_element].total; + double sigma = atom_density * p->micro_photon_xs_[i_element].total; // Increment probability to compare to cutoff prob += sigma; @@ -464,7 +464,7 @@ int sample_element(Particle* p) return i_element; } -Reaction* sample_fission(int i_nuclide, double E) +Reaction* sample_fission(int i_nuclide, const Particle* p) { // Get pointer to nuclide const auto& nuc {data::nuclides[i_nuclide]}; @@ -472,23 +472,23 @@ Reaction* sample_fission(int i_nuclide, double E) // If we're in the URR, by default use the first fission reaction. We also // default to the first reaction if we know that there are no partial fission // reactions - if (simulation::micro_xs[i_nuclide].use_ptable || !nuc->has_partial_fission_) { + if (p->micro_xs_[i_nuclide].use_ptable || !nuc->has_partial_fission_) { return nuc->fission_rx_[0]; } // Check to see if we are in a windowed multipole range. WMP only supports // the first fission reaction. if (nuc->multipole_) { - if (E >= nuc->multipole_->E_min_ && E <= nuc->multipole_->E_max_) { + if (p->E_ >= nuc->multipole_->E_min_ && p->E_ <= nuc->multipole_->E_max_) { return nuc->fission_rx_[0]; } } // Get grid index and interpolatoin factor and sample fission cdf - int i_temp = simulation::micro_xs[i_nuclide].index_temp; - int i_grid = simulation::micro_xs[i_nuclide].index_grid; - double f = simulation::micro_xs[i_nuclide].interp_factor; - double cutoff = prn() * simulation::micro_xs[i_nuclide].fission; + int i_temp = p->micro_xs_[i_nuclide].index_temp; + int i_grid = p->micro_xs_[i_nuclide].index_grid; + double f = p->micro_xs_[i_nuclide].interp_factor; + double cutoff = prn() * p->micro_xs_[i_nuclide].fission; double prob = 0.0; // Loop through each partial fission reaction type @@ -509,13 +509,13 @@ Reaction* sample_fission(int i_nuclide, double E) throw std::runtime_error{"No fission reaction was sampled for " + nuc->name_}; } -void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product) +void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_product) { // Get grid index and interpolation factor and sample photon production cdf - int i_temp = simulation::micro_xs[i_nuclide].index_temp; - int i_grid = simulation::micro_xs[i_nuclide].index_grid; - double f = simulation::micro_xs[i_nuclide].interp_factor; - double cutoff = prn() * simulation::micro_xs[i_nuclide].photon_prod; + int i_temp = p->micro_xs_[i_nuclide].index_temp; + int i_grid = p->micro_xs_[i_nuclide].index_grid; + double f = p->micro_xs_[i_nuclide].interp_factor; + double cutoff = prn() * p->micro_xs_[i_nuclide].photon_prod; double prob = 0.0; // Loop through each reaction type @@ -534,7 +534,7 @@ void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product) for (int j = 0; j < rx->products_.size(); ++j) { if (rx->products_[j].particle_ == Particle::Type::photon) { // add to cumulative probability - prob += (*rx->products_[j].yield_)(E) * xs; + prob += (*rx->products_[j].yield_)(p->E_) * xs; *i_rx = i; *i_product = j; @@ -548,8 +548,8 @@ void absorption(Particle* p, int i_nuclide) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->wgt_absorb_ = p->wgt_ * simulation::micro_xs[i_nuclide].absorption / - simulation::micro_xs[i_nuclide].total; + p->wgt_absorb_ = p->wgt_ * p->micro_xs_[i_nuclide].absorption / + p->micro_xs_[i_nuclide].total; // Adjust weight of particle by probability of absorption p->wgt_ -= p->wgt_absorb_; @@ -557,17 +557,17 @@ void absorption(Particle* p, int i_nuclide) // Score implicit absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { - global_tally_absorption += p->wgt_absorb_ * simulation::micro_xs[ - i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption; + global_tally_absorption += p->wgt_absorb_ * p->micro_xs_[ + i_nuclide].nu_fission / p->micro_xs_[i_nuclide].absorption; } } else { // See if disappearance reaction happens - if (simulation::micro_xs[i_nuclide].absorption > - prn() * simulation::micro_xs[i_nuclide].total) { + if (p->micro_xs_[i_nuclide].absorption > + prn() * p->micro_xs_[i_nuclide].total) { // Score absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { - global_tally_absorption += p->wgt_ * simulation::micro_xs[ - i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption; + global_tally_absorption += p->wgt_ * p->micro_xs_[ + i_nuclide].nu_fission / p->micro_xs_[i_nuclide].absorption; } p->alive_ = false; @@ -584,7 +584,7 @@ void scatter(Particle* p, int i_nuclide) // Get pointer to nuclide and grid index/interpolation factor const auto& nuc {data::nuclides[i_nuclide]}; - const auto& micro {simulation::micro_xs[i_nuclide]}; + const auto& micro {p->micro_xs_[i_nuclide]}; int i_temp = micro.index_temp; int i_grid = micro.index_grid; double f = micro.interp_factor; @@ -596,7 +596,7 @@ void scatter(Particle* p, int i_nuclide) // Calculate elastic cross section if it wasn't precalculated if (micro.elastic == CACHE_INVALID) { - nuc->calculate_elastic_xs(); + nuc->calculate_elastic_xs(*p); } double prob = micro.elastic - micro.thermal; @@ -608,8 +608,7 @@ void scatter(Particle* p, int i_nuclide) double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp]; // Perform collision physics for elastic scattering - elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, - p->E_, p->u(), p->mu_); + elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p); p->event_mt_ = ELASTIC; sampled = true; @@ -620,7 +619,7 @@ void scatter(Particle* p, int i_nuclide) // ======================================================================= // S(A,B) SCATTERING - sab_scatter(i_nuclide, micro.index_sab, p->E_, p->u(), p->mu_); + sab_scatter(i_nuclide, micro.index_sab, p); p->event_mt_ = ELASTIC; sampled = true; @@ -678,23 +677,23 @@ void scatter(Particle* p, int i_nuclide) } } -void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E, - Direction& u, double& mu_lab) +void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, + Particle* p) { // get pointer to nuclide const auto& nuc {data::nuclides[i_nuclide]}; - double vel = std::sqrt(E); + double vel = std::sqrt(p->E_); double awr = nuc->awr_; // Neutron velocity in LAB - Direction v_n = vel*u; + Direction v_n = vel*p->u(); // Sample velocity of target nucleus Direction v_t {}; - if (!simulation::micro_xs[i_nuclide].use_ptable) { - v_t = sample_target_velocity(nuc.get(), E, u, v_n, - simulation::micro_xs[i_nuclide].elastic, kT); + if (!p->micro_xs_[i_nuclide].use_ptable) { + v_t = sample_target_velocity(nuc.get(), p->E_, p->u(), v_n, + p->micro_xs_[i_nuclide].elastic, kT); } // Velocity of center-of-mass @@ -712,7 +711,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E, auto& d = rx.products_[0].distribution_[0]; auto d_ = dynamic_cast(d.get()); if (d_) { - mu_cm = d_->angle().sample(E); + mu_cm = d_->angle().sample(p->E_); } else { mu_cm = 2.0*prn() - 1.0; } @@ -728,35 +727,35 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E, // Transform back to LAB frame v_n += v_cm; - E = v_n.dot(v_n); - vel = std::sqrt(E); + p->E_ = v_n.dot(v_n); + vel = std::sqrt(p->E_); // compute cosine of scattering angle in LAB frame by taking dot product of // neutron's pre- and post-collision angle - mu_lab = u.dot(v_n) / vel; + p->mu_ = p->u().dot(v_n) / vel; // Set energy and direction of particle in LAB frame - u = v_n / vel; + p->u() = v_n / vel; // Because of floating-point roundoff, it may be possible for mu_lab to be // outside of the range [-1,1). In these cases, we just set mu_lab to exactly // -1 or 1 - if (std::abs(mu_lab) > 1.0) mu_lab = std::copysign(1.0, mu_lab); + if (std::abs(p->mu_) > 1.0) p->mu_ = std::copysign(1.0, p->mu_); } -void sab_scatter(int i_nuclide, int i_sab, double& E, Direction& u, double& mu) +void sab_scatter(int i_nuclide, int i_sab, Particle* p) { // Determine temperature index - const auto& micro {simulation::micro_xs[i_nuclide]}; + const auto& micro {p->micro_xs_[i_nuclide]}; int i_temp = micro.index_temp_sab; // Sample energy and angle double E_out; - data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, E, &E_out, &mu); + data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p->E_, &E_out, &p->mu_); // Set energy to outgoing, change direction of particle - E = E_out; - u = rotate_angle(u, mu, nullptr); + p->E_ = E_out; + p->u() = rotate_angle(p->u(), p->mu_, nullptr); } Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, @@ -1105,8 +1104,8 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) void sample_secondary_photons(Particle* p, int i_nuclide) { // Sample the number of photons produced - double y_t = p->wgt_ * simulation::micro_xs[i_nuclide].photon_prod / - simulation::micro_xs[i_nuclide].total; + double y_t = p->wgt_ * p->micro_xs_[i_nuclide].photon_prod / + p->micro_xs_[i_nuclide].total; int y = static_cast(y_t); if (prn() <= y_t - y) ++y; @@ -1115,7 +1114,7 @@ void sample_secondary_photons(Particle* p, int i_nuclide) // Sample the reaction and product int i_rx; int i_product; - sample_photon_product(i_nuclide, p->E_, &i_rx, &i_product); + sample_photon_product(i_nuclide, p, &i_rx, &i_product); // Sample the outgoing energy and angle auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index fa6365a3a5..9565a627c6 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -60,7 +60,7 @@ sample_reaction(Particle* p) // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs. - if (simulation::material_xs.absorption > 0.) { + if (p->material_xs_.absorption > 0.) { absorption(p); } else { p->wgt_absorb_ = 0.; @@ -113,7 +113,7 @@ create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank // Determine the expected number of neutrons produced double nu_t = p->wgt_ / simulation::keff * weight * - simulation::material_xs.nu_fission / simulation::material_xs.total; + p->material_xs_.nu_fission / p->material_xs_.total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -204,7 +204,7 @@ absorption(Particle* p) if (settings::survival_biasing) { // Determine weight absorbed in survival biasing p->wgt_absorb_ = p->wgt_ * - simulation::material_xs.absorption / simulation::material_xs.total; + p->material_xs_.absorption / p->material_xs_.total; // Adjust weight of particle by the probability of absorption p->wgt_ -= p->wgt_absorb_; @@ -213,14 +213,14 @@ absorption(Particle* p) // Score implicit absorpion estimate of keff #pragma omp atomic global_tally_absorption += p->wgt_absorb_ * - simulation::material_xs.nu_fission / - simulation::material_xs.absorption; + p->material_xs_.nu_fission / + p->material_xs_.absorption; } else { - if (simulation::material_xs.absorption > - prn() * simulation::material_xs.total) { + if (p->material_xs_.absorption > + prn() * p->material_xs_.total) { #pragma omp atomic - global_tally_absorption += p->wgt_ * simulation::material_xs.nu_fission / - simulation::material_xs.absorption; + global_tally_absorption += p->wgt_ * p->material_xs_.nu_fission / + p->material_xs_.absorption; p->alive_ = false; p->event_ = EVENT_ABSORB; } diff --git a/src/simulation.cpp b/src/simulation.cpp index 7d30ccea9c..b4db5e0e6d 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -78,13 +78,6 @@ int openmc_simulation_init() mat->init_nuclide_index(); } - // Create cross section caches - #pragma omp parallel - { - simulation::micro_xs = new NuclideMicroXS[data::nuclides.size()]; - simulation::micro_photon_xs = new ElementMicroXS[data::elements.size()]; - } - // Reset global variables -- this is done before loading state point (as that // will potentially populate k_generation and entropy) simulation::current_batch = 0; @@ -133,13 +126,6 @@ int openmc_simulation_finalize() mat->mat_nuclide_index_.clear(); } - // Clear cross section caches - #pragma omp parallel - { - delete[] simulation::micro_xs; - delete[] simulation::micro_photon_xs; - } - // Increment total number of generations simulation::total_gen += simulation::current_batch*settings::gen_per_batch; diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 315cc081e2..411995fd78 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -221,12 +221,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (i_nuclide == -1 && simulation::material_xs.total > 0.0) { + if (i_nuclide == -1 && p->material_xs_.total > 0.0) { score *= flux_deriv - + simulation::micro_xs[deriv.diff_nuclide].total - / simulation::material_xs.total; + + p->micro_xs_[deriv.diff_nuclide].total + / p->material_xs_.total; } else if (i_nuclide == deriv.diff_nuclide - && simulation::micro_xs[i_nuclide].total) { + && p->micro_xs_[i_nuclide].total) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -234,13 +234,13 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_SCATTER: - if (i_nuclide == -1 && (simulation::material_xs.total - - simulation::material_xs.absorption) > 0.0) { + if (i_nuclide == -1 && (p->material_xs_.total + - p->material_xs_.absorption) > 0.0) { score *= flux_deriv - + (simulation::micro_xs[deriv.diff_nuclide].total - - simulation::micro_xs[deriv.diff_nuclide].absorption) - / (simulation::material_xs.total - - simulation::material_xs.absorption); + + (p->micro_xs_[deriv.diff_nuclide].total + - p->micro_xs_[deriv.diff_nuclide].absorption) + / (p->material_xs_.total + - p->material_xs_.absorption); } else if (i_nuclide == deriv.diff_nuclide) { score *= flux_deriv + 1. / atom_density; } else { @@ -249,12 +249,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_ABSORPTION: - if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) { + if (i_nuclide == -1 && p->material_xs_.absorption > 0.0) { score *= flux_deriv - + simulation::micro_xs[deriv.diff_nuclide].absorption - / simulation::material_xs.absorption; + + p->micro_xs_[deriv.diff_nuclide].absorption + / p->material_xs_.absorption; } else if (i_nuclide == deriv.diff_nuclide - && simulation::micro_xs[i_nuclide].absorption) { + && p->micro_xs_[i_nuclide].absorption) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -262,12 +262,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_FISSION: - if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) { + if (i_nuclide == -1 && p->material_xs_.fission > 0.0) { score *= flux_deriv - + simulation::micro_xs[deriv.diff_nuclide].fission - / simulation::material_xs.fission; + + p->micro_xs_[deriv.diff_nuclide].fission + / p->material_xs_.fission; } else if (i_nuclide == deriv.diff_nuclide - && simulation::micro_xs[i_nuclide].fission) { + && p->micro_xs_[i_nuclide].fission) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -275,12 +275,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_NU_FISSION: - if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) { + if (i_nuclide == -1 && p->material_xs_.nu_fission > 0.0) { score *= flux_deriv - + simulation::micro_xs[deriv.diff_nuclide].nu_fission - / simulation::material_xs.nu_fission; + + p->micro_xs_[deriv.diff_nuclide].nu_fission + / p->material_xs_.nu_fission; } else if (i_nuclide == deriv.diff_nuclide - && simulation::micro_xs[i_nuclide].nu_fission) { + && p->micro_xs_[i_nuclide].nu_fission) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -331,64 +331,64 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (simulation::micro_xs[p->event_nuclide_].total) { + if (p->micro_xs_[p->event_nuclide_].total) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i) - / simulation::material_xs.total; + / p->material_xs_.total; } else { score *= flux_deriv; } break; case SCORE_SCATTER: - if (simulation::micro_xs[p->event_nuclide_].total - - simulation::micro_xs[p->event_nuclide_].absorption) { + if (p->micro_xs_[p->event_nuclide_].total + - p->micro_xs_[p->event_nuclide_].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_s * material.atom_density_(i) - / (simulation::material_xs.total - - simulation::material_xs.absorption); + / (p->material_xs_.total + - p->material_xs_.absorption); } else { score *= flux_deriv; } break; case SCORE_ABSORPTION: - if (simulation::micro_xs[p->event_nuclide_].absorption) { + if (p->micro_xs_[p->event_nuclide_].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_a * material.atom_density_(i) - / simulation::material_xs.absorption; + / p->material_xs_.absorption; } else { score *= flux_deriv; } break; case SCORE_FISSION: - if (simulation::micro_xs[p->event_nuclide_].fission) { + if (p->micro_xs_[p->event_nuclide_].fission) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_f * material.atom_density_(i) - / simulation::material_xs.fission; + / p->material_xs_.fission; } else { score *= flux_deriv; } break; case SCORE_NU_FISSION: - if (simulation::micro_xs[p->event_nuclide_].fission) { - double nu = simulation::micro_xs[p->event_nuclide_].nu_fission - / simulation::micro_xs[p->event_nuclide_].fission; + if (p->micro_xs_[p->event_nuclide_].fission) { + double nu = p->micro_xs_[p->event_nuclide_].nu_fission + / p->micro_xs_[p->event_nuclide_].fission; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + nu * dsig_f * material.atom_density_(i) - / simulation::material_xs.nu_fission; + / p->material_xs_.nu_fission; } else { score *= flux_deriv; } @@ -413,141 +413,141 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (i_nuclide == -1 && simulation::material_xs.total > 0.0) { + if (i_nuclide == -1 && p->material_xs_.total > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && simulation::micro_xs[i_nuc].total) { + && p->micro_xs_[i_nuc].total) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / simulation::material_xs.total; - } else if (simulation::micro_xs[i_nuclide].total) { + score *= flux_deriv + cum_dsig / p->material_xs_.total; + } else if (p->micro_xs_[i_nuclide].total) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv - + (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total; + + (dsig_s + dsig_a) / p->micro_xs_[i_nuclide].total; } else { score *= flux_deriv; } break; case SCORE_SCATTER: - if (i_nuclide == -1 && (simulation::material_xs.total - - simulation::material_xs.absorption)) { + if (i_nuclide == -1 && (p->material_xs_.total + - p->material_xs_.absorption)) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && (simulation::micro_xs[i_nuc].total - - simulation::micro_xs[i_nuc].absorption)) { + && (p->micro_xs_[i_nuc].total + - p->micro_xs_[i_nuc].absorption)) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_s * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / (simulation::material_xs.total - - simulation::material_xs.absorption); - } else if (simulation::micro_xs[i_nuclide].total - - simulation::micro_xs[i_nuclide].absorption) { + score *= flux_deriv + cum_dsig / (p->material_xs_.total + - p->material_xs_.absorption); + } else if (p->micro_xs_[i_nuclide].total + - p->micro_xs_[i_nuclide].absorption) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); - score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total - - simulation::micro_xs[i_nuclide].absorption); + score *= flux_deriv + dsig_s / (p->micro_xs_[i_nuclide].total + - p->micro_xs_[i_nuclide].absorption); } else { score *= flux_deriv; } break; case SCORE_ABSORPTION: - if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) { + if (i_nuclide == -1 && p->material_xs_.absorption > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && simulation::micro_xs[i_nuc].absorption) { + && p->micro_xs_[i_nuc].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_a * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / simulation::material_xs.absorption; - } else if (simulation::micro_xs[i_nuclide].absorption) { + score *= flux_deriv + cum_dsig / p->material_xs_.absorption; + } else if (p->micro_xs_[i_nuclide].absorption) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv - + dsig_a / simulation::micro_xs[i_nuclide].absorption; + + dsig_a / p->micro_xs_[i_nuclide].absorption; } else { score *= flux_deriv; } break; case SCORE_FISSION: - if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) { + if (i_nuclide == -1 && p->material_xs_.fission > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && simulation::micro_xs[i_nuc].fission) { + && p->micro_xs_[i_nuc].fission) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_f * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / simulation::material_xs.fission; - } else if (simulation::micro_xs[i_nuclide].fission) { + score *= flux_deriv + cum_dsig / p->material_xs_.fission; + } else if (p->micro_xs_[i_nuclide].fission) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv - + dsig_f / simulation::micro_xs[i_nuclide].fission; + + dsig_f / p->micro_xs_[i_nuclide].fission; } else { score *= flux_deriv; } break; case SCORE_NU_FISSION: - if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) { + if (i_nuclide == -1 && p->material_xs_.nu_fission > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && simulation::micro_xs[i_nuc].fission) { - double nu = simulation::micro_xs[i_nuc].nu_fission - / simulation::micro_xs[i_nuc].fission; + && p->micro_xs_[i_nuc].fission) { + double nu = p->micro_xs_[i_nuc].nu_fission + / p->micro_xs_[i_nuc].fission; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += nu * dsig_f * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / simulation::material_xs.nu_fission; - } else if (simulation::micro_xs[i_nuclide].fission) { + score *= flux_deriv + cum_dsig / p->material_xs_.nu_fission; + } else if (p->micro_xs_[i_nuclide].fission) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv - + dsig_f / simulation::micro_xs[i_nuclide].fission; + + dsig_f / p->micro_xs_[i_nuclide].fission; } else { score *= flux_deriv; } @@ -582,7 +582,7 @@ score_track_derivative(const Particle* p, double distance) // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist // (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist - deriv.flux_deriv -= distance * simulation::material_xs.total + deriv.flux_deriv -= distance * p->material_xs_.total / material.density_gpcc_; break; @@ -591,7 +591,7 @@ score_track_derivative(const Particle* p, double distance) // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist deriv.flux_deriv -= distance - * simulation::micro_xs[deriv.diff_nuclide].total; + * p->micro_xs_[deriv.diff_nuclide].total; break; case DIFF_TEMPERATURE: @@ -660,7 +660,7 @@ void score_collision_derivative(const Particle* p) // phi is proportional to Sigma_s // (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s // (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s - const auto& micro_xs {simulation::micro_xs[i_nuc]}; + const auto& micro_xs {p->micro_xs_[i_nuc]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index ed4d092545..57e810dfd2 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -317,7 +317,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) //! is not used for analog tallies. void -score_general_ce(const Particle* p, int i_tally, int start_index, +score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, int i_nuclide, double atom_density, double flux) { auto& tally {*model::tallies[i_tally]}; @@ -349,7 +349,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (p->type_ == Particle::Type::neutron || p->type_ == Particle::Type::photon) { - score *= flux / simulation::material_xs.total; + score *= flux / p->material_xs_.total; } else { score = 0.; } @@ -373,9 +373,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { if (i_nuclide >= 0) { - score = simulation::micro_xs[i_nuclide].total * atom_density * flux; + score = p->micro_xs_[i_nuclide].total * atom_density * flux; } else { - score = simulation::material_xs.total * flux; + score = p->material_xs_.total * flux; } } break; @@ -393,7 +393,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { score = p->wgt_last_; } - score *= flux / simulation::material_xs.total; + score *= flux / p->material_xs_.total; } else { score = flux; } @@ -411,11 +411,11 @@ score_general_ce(const Particle* p, int i_tally, int start_index, score = p->wgt_last_ * flux; } else { if (i_nuclide >= 0) { - score = (simulation::micro_xs[i_nuclide].total - - simulation::micro_xs[i_nuclide].absorption) * atom_density * flux; + score = (p->micro_xs_[i_nuclide].total + - p->micro_xs_[i_nuclide].absorption) * atom_density * flux; } else { - score = (simulation::material_xs.total - - simulation::material_xs.absorption) * flux; + score = (p->material_xs_.total + - p->material_xs_.absorption) * flux; } } break; @@ -458,26 +458,26 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } } else { if (i_nuclide >= 0) { - score = simulation::micro_xs[i_nuclide].absorption * atom_density + score = p->micro_xs_[i_nuclide].absorption * atom_density * flux; } else { - score = simulation::material_xs.absorption * flux; + score = p->material_xs_.absorption * flux; } } break; case SCORE_FISSION: - if (simulation::material_xs.absorption == 0) continue; + if (p->material_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + if (p->micro_xs_[p->event_nuclide_].absorption > 0) { score = p->wgt_absorb_ - * simulation::micro_xs[p->event_nuclide_].fission - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].fission + / p->micro_xs_[p->event_nuclide_].absorption * flux; } else { score = 0.; } @@ -488,21 +488,21 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // weight entering the collision as the estimate for the fission // reaction rate score = p->wgt_last_ - * simulation::micro_xs[p->event_nuclide_].fission - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].fission + / p->micro_xs_[p->event_nuclide_].absorption * flux; } } else { if (i_nuclide >= 0) { - score = simulation::micro_xs[i_nuclide].fission * atom_density * flux; + score = p->micro_xs_[i_nuclide].fission * atom_density * flux; } else { - score = simulation::material_xs.fission * flux; + score = p->material_xs_.fission * flux; } } break; case SCORE_NU_FISSION: - if (simulation::material_xs.absorption == 0) continue; + if (p->material_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { @@ -516,10 +516,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + if (p->micro_xs_[p->event_nuclide_].absorption > 0) { score = p->wgt_absorb_ - * simulation::micro_xs[p->event_nuclide_].nu_fission - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].nu_fission + / p->micro_xs_[p->event_nuclide_].absorption * flux; } else { score = 0.; } @@ -535,17 +535,17 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } } else { if (i_nuclide >= 0) { - score = simulation::micro_xs[i_nuclide].nu_fission * atom_density + score = p->micro_xs_[i_nuclide].nu_fission * atom_density * flux; } else { - score = simulation::material_xs.nu_fission * flux; + score = p->material_xs_.nu_fission * flux; } } break; case SCORE_PROMPT_NU_FISSION: - if (simulation::material_xs.absorption == 0) continue; + if (p->material_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { @@ -559,12 +559,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + if (p->micro_xs_[p->event_nuclide_].absorption > 0) { score = p->wgt_absorb_ - * simulation::micro_xs[p->event_nuclide_].fission + * p->micro_xs_[p->event_nuclide_].fission * data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::prompt) - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + / p->micro_xs_[p->event_nuclide_].absorption * flux; } else { score = 0.; } @@ -583,7 +583,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } } else { if (i_nuclide >= 0) { - score = simulation::micro_xs[i_nuclide].fission + score = p->micro_xs_[i_nuclide].fission * data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::prompt) * atom_density * flux; @@ -595,7 +595,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += simulation::micro_xs[j_nuclide].fission + score += p->micro_xs_[j_nuclide].fission * data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::prompt) * atom_density * flux; @@ -607,7 +607,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_DELAYED_NU_FISSION: - if (simulation::material_xs.absorption == 0) continue; + if (p->material_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { @@ -621,7 +621,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - if (simulation::micro_xs[p->event_nuclide_].absorption > 0 + if (p->micro_xs_[p->event_nuclide_].absorption > 0 && data::nuclides[p->event_nuclide_]->fissionable_) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; @@ -634,8 +634,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = p->wgt_absorb_ * yield - * simulation::micro_xs[p->event_nuclide_].fission - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].fission + / p->micro_xs_[p->event_nuclide_].absorption * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } @@ -645,10 +645,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs score = p->wgt_absorb_ - * simulation::micro_xs[p->event_nuclide_].fission + * p->micro_xs_[p->event_nuclide_].fission * data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed) - / simulation::micro_xs[p->event_nuclide_].absorption *flux; + / p->micro_xs_[p->event_nuclide_].absorption *flux; } } } else { @@ -694,7 +694,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto d = filt.groups_[d_bin]; auto yield = data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = simulation::micro_xs[i_nuclide].fission * yield + score = p->micro_xs_[i_nuclide].fission * yield * atom_density * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } @@ -702,7 +702,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { // If the delayed group filter is not present, compute the score // by multiplying the delayed-nu-fission macro xs by the flux - score = simulation::micro_xs[i_nuclide].fission + score = p->micro_xs_[i_nuclide].fission * data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed) * atom_density * flux; @@ -723,7 +723,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto d = filt.groups_[d_bin]; auto yield = data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = simulation::micro_xs[j_nuclide].fission * yield + score = p->micro_xs_[j_nuclide].fission * yield * atom_density * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -738,7 +738,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += simulation::micro_xs[j_nuclide].fission + score += p->micro_xs_[j_nuclide].fission * data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed) * atom_density * flux; @@ -751,14 +751,14 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_DECAY_RATE: - if (simulation::material_xs.absorption == 0) continue; + if (p->material_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (simulation::micro_xs[p->event_nuclide_].absorption > 0 + if (p->micro_xs_[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; if (tally.delayedgroup_filter_ != C_NONE) { @@ -773,8 +773,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; score = p->wgt_absorb_ * yield - * simulation::micro_xs[p->event_nuclide_].fission - / simulation::micro_xs[p->event_nuclide_].absorption + * p->micro_xs_[p->event_nuclide_].fission + / p->micro_xs_[p->event_nuclide_].absorption * rate * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -796,8 +796,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; score += rate * p->wgt_absorb_ - * simulation::micro_xs[p->event_nuclide_].fission * yield - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].fission * yield + / p->micro_xs_[p->event_nuclide_].absorption * flux; } } } @@ -854,7 +854,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = simulation::micro_xs[i_nuclide].fission * yield * flux + score = p->micro_xs_[i_nuclide].fission * yield * flux * atom_density * rate; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } @@ -870,7 +870,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += simulation::micro_xs[i_nuclide].fission * flux + score += p->micro_xs_[i_nuclide].fission * flux * yield * atom_density * rate; } } @@ -894,7 +894,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = simulation::micro_xs[j_nuclide].fission * yield + score = p->micro_xs_[j_nuclide].fission * yield * flux * atom_density * rate; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -922,7 +922,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += simulation::micro_xs[j_nuclide].fission + score += p->micro_xs_[j_nuclide].fission * yield * atom_density * flux * rate; } } @@ -935,7 +935,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_KAPPA_FISSION: - if (simulation::material_xs.absorption == 0.) continue; + if (p->material_xs_.absorption == 0.) continue; score = 0.; // Kappa-fission values are determined from the Q-value listed for the // fission cross section. @@ -945,12 +945,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (simulation::micro_xs[p->event_nuclide_].absorption > 0 + if (p->micro_xs_[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; score = p->wgt_absorb_ * rxn.q_value_ - * simulation::micro_xs[p->event_nuclide_].fission - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].fission + / p->micro_xs_[p->event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events @@ -959,12 +959,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // weight entering the collision as the estimate for the fission // reaction rate const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (simulation::micro_xs[p->event_nuclide_].absorption > 0 + if (p->micro_xs_[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; score = p->wgt_last_ * rxn.q_value_ - * simulation::micro_xs[p->event_nuclide_].fission - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].fission + / p->micro_xs_[p->event_nuclide_].absorption * flux; } } } else { @@ -972,7 +972,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, const auto& nuc {*data::nuclides[i_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = rxn.q_value_ * simulation::micro_xs[i_nuclide].fission + score = rxn.q_value_ * p->micro_xs_[i_nuclide].fission * atom_density * flux; } } else { @@ -984,7 +984,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, const auto& nuc {*data::nuclides[j_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score += rxn.q_value_ * simulation::micro_xs[j_nuclide].fission + score += rxn.q_value_ * p->micro_xs_[j_nuclide].fission * atom_density * flux; } } @@ -1007,9 +1007,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, score = p->wgt_last_ * flux; } else { if (i_nuclide >= 0) { - if (simulation::micro_xs[i_nuclide].elastic == CACHE_INVALID) - data::nuclides[i_nuclide]->calculate_elastic_xs(); - score = simulation::micro_xs[i_nuclide].elastic * atom_density * flux; + if (p->micro_xs_[i_nuclide].elastic == CACHE_INVALID) + data::nuclides[i_nuclide]->calculate_elastic_xs(*p); + score = p->micro_xs_[i_nuclide].elastic * atom_density * flux; } else { score = 0.; if (p->material_ != MATERIAL_VOID) { @@ -1017,9 +1017,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - if (simulation::micro_xs[j_nuclide].elastic == CACHE_INVALID) - data::nuclides[j_nuclide]->calculate_elastic_xs(); - score += simulation::micro_xs[j_nuclide].elastic * atom_density + if (p->micro_xs_[j_nuclide].elastic == CACHE_INVALID) + data::nuclides[j_nuclide]->calculate_elastic_xs(*p); + score += p->micro_xs_[j_nuclide].elastic * atom_density * flux; } } @@ -1030,7 +1030,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_FISS_Q_PROMPT: case SCORE_FISS_Q_RECOV: //continue; - if (simulation::material_xs.absorption == 0.) continue; + if (p->material_xs_.absorption == 0.) continue; score = 0.; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing) { @@ -1038,7 +1038,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + if (p->micro_xs_[p->event_nuclide_].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) @@ -1048,8 +1048,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, q_value = (*nuc.fission_q_recov_)(p->E_last_); } score = p->wgt_absorb_ * q_value - * simulation::micro_xs[p->event_nuclide_].fission - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].fission + / p->micro_xs_[p->event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events @@ -1058,7 +1058,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // weight entering the collision as the estimate for the fission // reaction rate const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + if (p->micro_xs_[p->event_nuclide_].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) @@ -1068,8 +1068,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, q_value = (*nuc.fission_q_recov_)(p->E_last_); } score = p->wgt_last_ * q_value - * simulation::micro_xs[p->event_nuclide_].fission - / simulation::micro_xs[p->event_nuclide_].absorption * flux; + * p->micro_xs_[p->event_nuclide_].fission + / p->micro_xs_[p->event_nuclide_].absorption * flux; } } } else { @@ -1083,7 +1083,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (nuc.fission_q_recov_) q_value = (*nuc.fission_q_recov_)(p->E_last_); } - score = q_value * simulation::micro_xs[i_nuclide].fission + score = q_value * p->micro_xs_[i_nuclide].fission * atom_density * flux; } else { if (p->material_ != MATERIAL_VOID) { @@ -1100,7 +1100,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (nuc.fission_q_recov_) q_value = (*nuc.fission_q_recov_)(p->E_last_); } - score += q_value * simulation::micro_xs[j_nuclide].fission + score += q_value * p->micro_xs_[j_nuclide].fission * atom_density * flux; } } @@ -1130,7 +1130,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case N_4N: m = 5; break; } if (i_nuclide >= 0) { - score = simulation::micro_xs[i_nuclide].reaction[m] * atom_density + score = p->micro_xs_[i_nuclide].reaction[m] * atom_density * flux; } else { score = 0.; @@ -1139,7 +1139,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += simulation::micro_xs[j_nuclide].reaction[m] + score += p->micro_xs_[j_nuclide].reaction[m] * atom_density * flux; } } @@ -1164,10 +1164,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = simulation::micro_xs[i_nuclide].index_temp; + auto i_temp = p->micro_xs_[i_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = simulation::micro_xs[i_nuclide].index_grid; - auto f = simulation::micro_xs[i_nuclide].interp_factor; + auto i_grid = p->micro_xs_[i_nuclide].index_grid; + auto f = p->micro_xs_[i_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { score = ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1184,10 +1184,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = simulation::micro_xs[j_nuclide].index_temp; + auto i_temp = p->micro_xs_[j_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = simulation::micro_xs[j_nuclide].index_grid; - auto f = simulation::micro_xs[j_nuclide].interp_factor; + auto i_grid = p->micro_xs_[j_nuclide].index_grid; + auto f = p->micro_xs_[j_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { score += ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1291,7 +1291,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } else { score = p->wgt_last_; } - score *= flux / simulation::material_xs.total; + score *= flux / p->material_xs_.total; } else { score = flux; } @@ -1320,7 +1320,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = get_nuclide_xs(i_nuclide, MG_GET_XS_TOTAL, p_g) * atom_density * flux; } else { - score = simulation::material_xs.total * flux; + score = p->material_xs_.total * flux; } } break; @@ -1438,7 +1438,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = atom_density * flux * get_nuclide_xs(i_nuclide, MG_GET_XS_ABSORPTION, p_g); } else { - score = simulation::material_xs.absorption * flux; + score = p->material_xs_.absorption * flux; } } break; @@ -1920,7 +1920,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, //! Tally rates for when the user requests a tally on all nuclides. void -score_all_nuclides(const Particle* p, int i_tally, double flux, +score_all_nuclides(Particle* p, int i_tally, double flux, int filter_index) { const Tally& tally {*model::tallies[i_tally]}; @@ -1955,7 +1955,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux, } } -void score_analog_tally_ce(const Particle* p) +void score_analog_tally_ce(Particle* p) { for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2059,7 +2059,7 @@ void score_analog_tally_mg(const Particle* p) } void -score_tracklength_tally(const Particle* p, double distance) +score_tracklength_tally(Particle* p, double distance) { // Determine the tracklength estimate of the flux double flux = p->wgt_ * distance; @@ -2122,14 +2122,14 @@ score_tracklength_tally(const Particle* p, double distance) match.bins_present_ = false; } -void score_collision_tally(const Particle* p) +void score_collision_tally(Particle* p) { // Determine the collision estimate of the flux double flux; if (!settings::survival_biasing) { - flux = p->wgt_last_ / simulation::material_xs.total; + flux = p->wgt_last_ / p->material_xs_.total; } else { - flux = (p->wgt_last_ + p->wgt_absorb_) / simulation::material_xs.total; + flux = (p->wgt_last_ + p->wgt_absorb_) / p->material_xs_.total; } for (auto i_tally : model::active_collision_tallies) { From 52fc93e6e038df2455b80adf268fa5f3af60b173 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Mar 2019 09:18:54 -0500 Subject: [PATCH 02/14] Rename cross section caches --- include/openmc/particle.h | 11 ++- src/material.cpp | 38 ++++---- src/nuclide.cpp | 8 +- src/particle.cpp | 26 +++--- src/photon.cpp | 2 +- src/physics.cpp | 62 ++++++------- src/physics_mg.cpp | 18 ++-- src/tallies/derivative.cpp | 142 ++++++++++++++-------------- src/tallies/tally_scoring.cpp | 170 +++++++++++++++++----------------- 9 files changed, 239 insertions(+), 238 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index e971c30da3..2130b6454a 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -111,11 +111,11 @@ struct ElementMicroXS { }; //============================================================================== -// MATERIALMACROXS contains cached macroscopic cross sections for the material a +// MACROXS contains cached macroscopic cross sections for the material a // particle is traveling through //============================================================================== -struct MaterialMacroXS { +struct MacroXS { double total; //!< macroscopic total xs double absorption; //!< macroscopic absorption xs double fission; //!< macroscopic fission xs @@ -219,9 +219,10 @@ public: //========================================================================== // Data members - std::vector micro_xs_; - std::vector micro_photon_xs_; - MaterialMacroXS material_xs_; + // 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.) diff --git a/src/material.cpp b/src/material.cpp index 0d2fe4bbc3..94a6214836 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -732,10 +732,10 @@ void Material::init_nuclide_index() void Material::calculate_xs(Particle& p) const { // Set all material macroscopic cross sections to zero - p.material_xs_.total = 0.0; - p.material_xs_.absorption = 0.0; - p.material_xs_.fission = 0.0; - p.material_xs_.nu_fission = 0.0; + p.macro_xs_.total = 0.0; + p.macro_xs_.absorption = 0.0; + p.macro_xs_.fission = 0.0; + p.macro_xs_.nu_fission = 0.0; if (p.type_ == Particle::Type::neutron) { this->calculate_neutron_xs(p); @@ -792,7 +792,7 @@ void Material::calculate_neutron_xs(Particle& p) const int i_nuclide = nuclide_[i]; // Calculate microscopic cross section for this nuclide - const auto& micro {p.micro_xs_[i_nuclide]}; + const auto& micro {p.neutron_xs_[i_nuclide]}; if (p.E_ != micro.last_E || p.sqrtkT_ != micro.last_sqrtkT || i_sab != micro.index_sab @@ -807,19 +807,19 @@ void Material::calculate_neutron_xs(Particle& p) const double atom_density = atom_density_(i); // Add contributions to cross sections - p.material_xs_.total += atom_density * micro.total; - p.material_xs_.absorption += atom_density * micro.absorption; - p.material_xs_.fission += atom_density * micro.fission; - p.material_xs_.nu_fission += atom_density * micro.nu_fission; + p.macro_xs_.total += atom_density * micro.total; + p.macro_xs_.absorption += atom_density * micro.absorption; + p.macro_xs_.fission += atom_density * micro.fission; + p.macro_xs_.nu_fission += atom_density * micro.nu_fission; } } void Material::calculate_photon_xs(Particle& p) const { - p.material_xs_.coherent = 0.0; - p.material_xs_.incoherent = 0.0; - p.material_xs_.photoelectric = 0.0; - p.material_xs_.pair_production = 0.0; + p.macro_xs_.coherent = 0.0; + p.macro_xs_.incoherent = 0.0; + p.macro_xs_.photoelectric = 0.0; + p.macro_xs_.pair_production = 0.0; // Add contribution from each nuclide in material for (int i = 0; i < nuclide_.size(); ++i) { @@ -830,7 +830,7 @@ void Material::calculate_photon_xs(Particle& p) const int i_element = element_[i]; // Calculate microscopic cross section for this nuclide - const auto& micro {p.micro_photon_xs_[i_element]}; + const auto& micro {p.photon_xs_[i_element]}; if (p.E_ != micro.last_E) { data::elements[i_element].calculate_xs(p); } @@ -842,11 +842,11 @@ void Material::calculate_photon_xs(Particle& p) const double atom_density = atom_density_(i); // Add contributions to material macroscopic cross sections - p.material_xs_.total += atom_density * micro.total; - p.material_xs_.coherent += atom_density * micro.coherent; - p.material_xs_.incoherent += atom_density * micro.incoherent; - p.material_xs_.photoelectric += atom_density * micro.photoelectric; - p.material_xs_.pair_production += atom_density * micro.pair_production; + p.macro_xs_.total += atom_density * micro.total; + p.macro_xs_.coherent += atom_density * micro.coherent; + p.macro_xs_.incoherent += atom_density * micro.incoherent; + p.macro_xs_.photoelectric += atom_density * micro.photoelectric; + p.macro_xs_.pair_production += atom_density * micro.pair_production; } } diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 87965dfc3e..cfb5a669f2 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -459,7 +459,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const void Nuclide::calculate_elastic_xs(Particle& p) const { // Get temperature index, grid index, and interpolation factor - auto& micro {p.micro_xs_[i_nuclide_]}; + auto& micro {p.neutron_xs_[i_nuclide_]}; int i_temp = micro.index_temp; int i_grid = micro.index_grid; double f = micro.interp_factor; @@ -495,7 +495,7 @@ double Nuclide::elastic_xs_0K(double E) const void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle& p) { - auto& micro {p.micro_xs_[i_nuclide_]}; + auto& micro {p.neutron_xs_[i_nuclide_]}; // Initialize cached cross sections to zero micro.elastic = CACHE_INVALID; @@ -702,7 +702,7 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p) { - auto& micro {p.micro_xs_[i_nuclide_]}; + auto& micro {p.neutron_xs_[i_nuclide_]}; // Set flag that S(a,b) treatment should be used for scattering micro.index_sab = i_sab; @@ -731,7 +731,7 @@ void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p) void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const { - auto& micro = p.micro_xs_[i_nuclide_]; + auto& micro = p.neutron_xs_[i_nuclide_]; micro.use_ptable = true; // Create a shorthand for the URR data diff --git a/src/particle.cpp b/src/particle.cpp index f5084cd1c2..88ff251dd5 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -58,8 +58,8 @@ Particle::Particle() } // Create microscopic cross section caches - micro_xs_.resize(data::nuclides.size()); - micro_photon_xs_.resize(data::elements.size()); + neutron_xs_.resize(data::nuclides.size()); + photon_xs_.resize(data::elements.size()); } void @@ -135,7 +135,7 @@ Particle::transport() // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { - for (auto& micro : micro_xs_) micro.last_E = 0.0; + for (auto& micro : neutron_xs_) micro.last_E = 0.0; } // Prepare to write out particle track. @@ -189,17 +189,17 @@ Particle::transport() } else { // Get the MG data calculate_xs_c(material_, g_, sqrtkT_, this->u_local(), - material_xs_.total, material_xs_.absorption, material_xs_.nu_fission); + macro_xs_.total, macro_xs_.absorption, macro_xs_.nu_fission); // Finally, update the particle group while we have already checked // for if multi-group g_last_ = g_; } } else { - material_xs_.total = 0.0; - material_xs_.absorption = 0.0; - material_xs_.fission = 0.0; - material_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; } // Find the distance to the nearest boundary @@ -215,10 +215,10 @@ Particle::transport() if (type_ == Particle::Type::electron || type_ == Particle::Type::positron) { d_collision = 0.0; - } else if (material_xs_.total == 0.0) { + } else if (macro_xs_.total == 0.0) { d_collision = INFINITY; } else { - d_collision = -std::log(prn()) / material_xs_.total; + d_collision = -std::log(prn()) / macro_xs_.total; } // Select smaller of the two distances @@ -237,7 +237,7 @@ Particle::transport() // Score track-length estimate of k-eff if (settings::run_mode == RUN_MODE_EIGENVALUE && type_ == Particle::Type::neutron) { - global_tally_tracklength += wgt_ * distance * material_xs_.nu_fission; + global_tally_tracklength += wgt_ * distance * macro_xs_.nu_fission; } // Score flux derivative accumulators for differential tallies. @@ -280,8 +280,8 @@ Particle::transport() // Score collision estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE && type_ == Particle::Type::neutron) { - global_tally_collision += wgt_ * material_xs_.nu_fission - / material_xs_.total; + global_tally_collision += wgt_ * macro_xs_.nu_fission + / macro_xs_.total; } // Score surface current tallies -- this has to be done before the collision diff --git a/src/photon.cpp b/src/photon.cpp index 9841e9e3bf..79345b90e0 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -452,7 +452,7 @@ void PhotonInteraction::calculate_xs(Particle& p) const // calculate interpolation factor double f = (log_E - energy_(i_grid)) / (energy_(i_grid+1) - energy_(i_grid)); - auto& xs {p.micro_photon_xs_[i_element_]}; + auto& xs {p.photon_xs_[i_element_]}; xs.index_grid = i_grid; xs.interp_factor = f; diff --git a/src/physics.cpp b/src/physics.cpp index 029ab18c86..4cc9b12a03 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -110,7 +110,7 @@ void sample_neutron_reaction(Particle* p) // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs - if (p->micro_xs_[i_nuclide].absorption > 0.0) { + if (p->neutron_xs_[i_nuclide].absorption > 0.0) { absorption(p, i_nuclide); } else { p->wgt_absorb_ = 0.0; @@ -146,8 +146,8 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p->wgt_ / simulation::keff * weight * p->micro_xs_[ - i_nuclide].nu_fission / p->micro_xs_[i_nuclide].total; + double nu_t = p->wgt_ / simulation::keff * weight * p->neutron_xs_[ + i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -227,7 +227,7 @@ void sample_photon_reaction(Particle* p) // Sample element within material int i_element = sample_element(p); p->event_nuclide_ = i_element; - const auto& micro {p->micro_photon_xs_[i_element]}; + const auto& micro {p->photon_xs_[i_element]}; const auto& element {data::elements[i_element]}; // Calculate photon energy over electron rest mass equivalent @@ -408,7 +408,7 @@ void sample_positron_reaction(Particle* p) int sample_nuclide(const Particle* p) { // Sample cumulative distribution function - double cutoff = prn() * p->material_xs_.total; + double cutoff = prn() * p->macro_xs_.total; // Get pointers to nuclide/density arrays const auto& mat {model::materials[p->material_]}; @@ -421,7 +421,7 @@ int sample_nuclide(const Particle* p) double atom_density = mat->atom_density_[i]; // Increment probability to compare to cutoff - prob += atom_density * p->micro_xs_[i_nuclide].total; + prob += atom_density * p->neutron_xs_[i_nuclide].total; if (prob >= cutoff) return i_nuclide; } @@ -433,7 +433,7 @@ int sample_nuclide(const Particle* p) int sample_element(Particle* p) { // Sample cumulative distribution function - double cutoff = prn() * p->material_xs_.total; + double cutoff = prn() * p->macro_xs_.total; // Get pointers to elements, densities const auto& mat {model::materials[p->material_]}; @@ -454,7 +454,7 @@ int sample_element(Particle* p) double atom_density = mat->atom_density_[i]; // Determine microscopic cross section - double sigma = atom_density * p->micro_photon_xs_[i_element].total; + double sigma = atom_density * p->photon_xs_[i_element].total; // Increment probability to compare to cutoff prob += sigma; @@ -472,7 +472,7 @@ Reaction* sample_fission(int i_nuclide, const Particle* p) // If we're in the URR, by default use the first fission reaction. We also // default to the first reaction if we know that there are no partial fission // reactions - if (p->micro_xs_[i_nuclide].use_ptable || !nuc->has_partial_fission_) { + if (p->neutron_xs_[i_nuclide].use_ptable || !nuc->has_partial_fission_) { return nuc->fission_rx_[0]; } @@ -485,10 +485,10 @@ Reaction* sample_fission(int i_nuclide, const Particle* p) } // Get grid index and interpolatoin factor and sample fission cdf - int i_temp = p->micro_xs_[i_nuclide].index_temp; - int i_grid = p->micro_xs_[i_nuclide].index_grid; - double f = p->micro_xs_[i_nuclide].interp_factor; - double cutoff = prn() * p->micro_xs_[i_nuclide].fission; + int i_temp = p->neutron_xs_[i_nuclide].index_temp; + int i_grid = p->neutron_xs_[i_nuclide].index_grid; + double f = p->neutron_xs_[i_nuclide].interp_factor; + double cutoff = prn() * p->neutron_xs_[i_nuclide].fission; double prob = 0.0; // Loop through each partial fission reaction type @@ -512,10 +512,10 @@ Reaction* sample_fission(int i_nuclide, const Particle* p) void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_product) { // Get grid index and interpolation factor and sample photon production cdf - int i_temp = p->micro_xs_[i_nuclide].index_temp; - int i_grid = p->micro_xs_[i_nuclide].index_grid; - double f = p->micro_xs_[i_nuclide].interp_factor; - double cutoff = prn() * p->micro_xs_[i_nuclide].photon_prod; + int i_temp = p->neutron_xs_[i_nuclide].index_temp; + int i_grid = p->neutron_xs_[i_nuclide].index_grid; + double f = p->neutron_xs_[i_nuclide].interp_factor; + double cutoff = prn() * p->neutron_xs_[i_nuclide].photon_prod; double prob = 0.0; // Loop through each reaction type @@ -548,8 +548,8 @@ void absorption(Particle* p, int i_nuclide) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->wgt_absorb_ = p->wgt_ * p->micro_xs_[i_nuclide].absorption / - p->micro_xs_[i_nuclide].total; + p->wgt_absorb_ = p->wgt_ * p->neutron_xs_[i_nuclide].absorption / + p->neutron_xs_[i_nuclide].total; // Adjust weight of particle by probability of absorption p->wgt_ -= p->wgt_absorb_; @@ -557,17 +557,17 @@ void absorption(Particle* p, int i_nuclide) // Score implicit absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { - global_tally_absorption += p->wgt_absorb_ * p->micro_xs_[ - i_nuclide].nu_fission / p->micro_xs_[i_nuclide].absorption; + global_tally_absorption += p->wgt_absorb_ * p->neutron_xs_[ + i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption; } } else { // See if disappearance reaction happens - if (p->micro_xs_[i_nuclide].absorption > - prn() * p->micro_xs_[i_nuclide].total) { + if (p->neutron_xs_[i_nuclide].absorption > + prn() * p->neutron_xs_[i_nuclide].total) { // Score absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { - global_tally_absorption += p->wgt_ * p->micro_xs_[ - i_nuclide].nu_fission / p->micro_xs_[i_nuclide].absorption; + global_tally_absorption += p->wgt_ * p->neutron_xs_[ + i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption; } p->alive_ = false; @@ -584,7 +584,7 @@ void scatter(Particle* p, int i_nuclide) // Get pointer to nuclide and grid index/interpolation factor const auto& nuc {data::nuclides[i_nuclide]}; - const auto& micro {p->micro_xs_[i_nuclide]}; + const auto& micro {p->neutron_xs_[i_nuclide]}; int i_temp = micro.index_temp; int i_grid = micro.index_grid; double f = micro.interp_factor; @@ -691,9 +691,9 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, // Sample velocity of target nucleus Direction v_t {}; - if (!p->micro_xs_[i_nuclide].use_ptable) { + if (!p->neutron_xs_[i_nuclide].use_ptable) { v_t = sample_target_velocity(nuc.get(), p->E_, p->u(), v_n, - p->micro_xs_[i_nuclide].elastic, kT); + p->neutron_xs_[i_nuclide].elastic, kT); } // Velocity of center-of-mass @@ -746,7 +746,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, void sab_scatter(int i_nuclide, int i_sab, Particle* p) { // Determine temperature index - const auto& micro {p->micro_xs_[i_nuclide]}; + const auto& micro {p->neutron_xs_[i_nuclide]}; int i_temp = micro.index_temp_sab; // Sample energy and angle @@ -1104,8 +1104,8 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) void sample_secondary_photons(Particle* p, int i_nuclide) { // Sample the number of photons produced - double y_t = p->wgt_ * p->micro_xs_[i_nuclide].photon_prod / - p->micro_xs_[i_nuclide].total; + double y_t = p->wgt_ * p->neutron_xs_[i_nuclide].photon_prod / + p->neutron_xs_[i_nuclide].total; int y = static_cast(y_t); if (prn() <= y_t - y) ++y; diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 9565a627c6..affe0ec197 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -60,7 +60,7 @@ sample_reaction(Particle* p) // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs. - if (p->material_xs_.absorption > 0.) { + if (p->macro_xs_.absorption > 0.) { absorption(p); } else { p->wgt_absorb_ = 0.; @@ -113,7 +113,7 @@ create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank // Determine the expected number of neutrons produced double nu_t = p->wgt_ / simulation::keff * weight * - p->material_xs_.nu_fission / p->material_xs_.total; + p->macro_xs_.nu_fission / p->macro_xs_.total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -204,7 +204,7 @@ absorption(Particle* p) if (settings::survival_biasing) { // Determine weight absorbed in survival biasing p->wgt_absorb_ = p->wgt_ * - p->material_xs_.absorption / p->material_xs_.total; + p->macro_xs_.absorption / p->macro_xs_.total; // Adjust weight of particle by the probability of absorption p->wgt_ -= p->wgt_absorb_; @@ -213,14 +213,14 @@ absorption(Particle* p) // Score implicit absorpion estimate of keff #pragma omp atomic global_tally_absorption += p->wgt_absorb_ * - p->material_xs_.nu_fission / - p->material_xs_.absorption; + p->macro_xs_.nu_fission / + p->macro_xs_.absorption; } else { - if (p->material_xs_.absorption > - prn() * p->material_xs_.total) { + if (p->macro_xs_.absorption > + prn() * p->macro_xs_.total) { #pragma omp atomic - global_tally_absorption += p->wgt_ * p->material_xs_.nu_fission / - p->material_xs_.absorption; + global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission / + p->macro_xs_.absorption; p->alive_ = false; p->event_ = EVENT_ABSORB; } diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 411995fd78..290d093843 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -221,12 +221,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (i_nuclide == -1 && p->material_xs_.total > 0.0) { + if (i_nuclide == -1 && p->macro_xs_.total > 0.0) { score *= flux_deriv - + p->micro_xs_[deriv.diff_nuclide].total - / p->material_xs_.total; + + p->neutron_xs_[deriv.diff_nuclide].total + / p->macro_xs_.total; } else if (i_nuclide == deriv.diff_nuclide - && p->micro_xs_[i_nuclide].total) { + && p->neutron_xs_[i_nuclide].total) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -234,13 +234,13 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_SCATTER: - if (i_nuclide == -1 && (p->material_xs_.total - - p->material_xs_.absorption) > 0.0) { + if (i_nuclide == -1 && (p->macro_xs_.total + - p->macro_xs_.absorption) > 0.0) { score *= flux_deriv - + (p->micro_xs_[deriv.diff_nuclide].total - - p->micro_xs_[deriv.diff_nuclide].absorption) - / (p->material_xs_.total - - p->material_xs_.absorption); + + (p->neutron_xs_[deriv.diff_nuclide].total + - p->neutron_xs_[deriv.diff_nuclide].absorption) + / (p->macro_xs_.total + - p->macro_xs_.absorption); } else if (i_nuclide == deriv.diff_nuclide) { score *= flux_deriv + 1. / atom_density; } else { @@ -249,12 +249,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_ABSORPTION: - if (i_nuclide == -1 && p->material_xs_.absorption > 0.0) { + if (i_nuclide == -1 && p->macro_xs_.absorption > 0.0) { score *= flux_deriv - + p->micro_xs_[deriv.diff_nuclide].absorption - / p->material_xs_.absorption; + + p->neutron_xs_[deriv.diff_nuclide].absorption + / p->macro_xs_.absorption; } else if (i_nuclide == deriv.diff_nuclide - && p->micro_xs_[i_nuclide].absorption) { + && p->neutron_xs_[i_nuclide].absorption) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -262,12 +262,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_FISSION: - if (i_nuclide == -1 && p->material_xs_.fission > 0.0) { + if (i_nuclide == -1 && p->macro_xs_.fission > 0.0) { score *= flux_deriv - + p->micro_xs_[deriv.diff_nuclide].fission - / p->material_xs_.fission; + + p->neutron_xs_[deriv.diff_nuclide].fission + / p->macro_xs_.fission; } else if (i_nuclide == deriv.diff_nuclide - && p->micro_xs_[i_nuclide].fission) { + && p->neutron_xs_[i_nuclide].fission) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -275,12 +275,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_NU_FISSION: - if (i_nuclide == -1 && p->material_xs_.nu_fission > 0.0) { + if (i_nuclide == -1 && p->macro_xs_.nu_fission > 0.0) { score *= flux_deriv - + p->micro_xs_[deriv.diff_nuclide].nu_fission - / p->material_xs_.nu_fission; + + p->neutron_xs_[deriv.diff_nuclide].nu_fission + / p->macro_xs_.nu_fission; } else if (i_nuclide == deriv.diff_nuclide - && p->micro_xs_[i_nuclide].nu_fission) { + && p->neutron_xs_[i_nuclide].nu_fission) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -331,64 +331,64 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (p->micro_xs_[p->event_nuclide_].total) { + if (p->neutron_xs_[p->event_nuclide_].total) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i) - / p->material_xs_.total; + / p->macro_xs_.total; } else { score *= flux_deriv; } break; case SCORE_SCATTER: - if (p->micro_xs_[p->event_nuclide_].total - - p->micro_xs_[p->event_nuclide_].absorption) { + if (p->neutron_xs_[p->event_nuclide_].total + - p->neutron_xs_[p->event_nuclide_].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_s * material.atom_density_(i) - / (p->material_xs_.total - - p->material_xs_.absorption); + / (p->macro_xs_.total + - p->macro_xs_.absorption); } else { score *= flux_deriv; } break; case SCORE_ABSORPTION: - if (p->micro_xs_[p->event_nuclide_].absorption) { + if (p->neutron_xs_[p->event_nuclide_].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_a * material.atom_density_(i) - / p->material_xs_.absorption; + / p->macro_xs_.absorption; } else { score *= flux_deriv; } break; case SCORE_FISSION: - if (p->micro_xs_[p->event_nuclide_].fission) { + if (p->neutron_xs_[p->event_nuclide_].fission) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_f * material.atom_density_(i) - / p->material_xs_.fission; + / p->macro_xs_.fission; } else { score *= flux_deriv; } break; case SCORE_NU_FISSION: - if (p->micro_xs_[p->event_nuclide_].fission) { - double nu = p->micro_xs_[p->event_nuclide_].nu_fission - / p->micro_xs_[p->event_nuclide_].fission; + if (p->neutron_xs_[p->event_nuclide_].fission) { + double nu = p->neutron_xs_[p->event_nuclide_].nu_fission + / p->neutron_xs_[p->event_nuclide_].fission; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + nu * dsig_f * material.atom_density_(i) - / p->material_xs_.nu_fission; + / p->macro_xs_.nu_fission; } else { score *= flux_deriv; } @@ -413,141 +413,141 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (i_nuclide == -1 && p->material_xs_.total > 0.0) { + if (i_nuclide == -1 && p->macro_xs_.total > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && p->micro_xs_[i_nuc].total) { + && p->neutron_xs_[i_nuc].total) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / p->material_xs_.total; - } else if (p->micro_xs_[i_nuclide].total) { + score *= flux_deriv + cum_dsig / p->macro_xs_.total; + } else if (p->neutron_xs_[i_nuclide].total) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv - + (dsig_s + dsig_a) / p->micro_xs_[i_nuclide].total; + + (dsig_s + dsig_a) / p->neutron_xs_[i_nuclide].total; } else { score *= flux_deriv; } break; case SCORE_SCATTER: - if (i_nuclide == -1 && (p->material_xs_.total - - p->material_xs_.absorption)) { + if (i_nuclide == -1 && (p->macro_xs_.total + - p->macro_xs_.absorption)) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && (p->micro_xs_[i_nuc].total - - p->micro_xs_[i_nuc].absorption)) { + && (p->neutron_xs_[i_nuc].total + - p->neutron_xs_[i_nuc].absorption)) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_s * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / (p->material_xs_.total - - p->material_xs_.absorption); - } else if (p->micro_xs_[i_nuclide].total - - p->micro_xs_[i_nuclide].absorption) { + score *= flux_deriv + cum_dsig / (p->macro_xs_.total + - p->macro_xs_.absorption); + } else if (p->neutron_xs_[i_nuclide].total + - p->neutron_xs_[i_nuclide].absorption) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); - score *= flux_deriv + dsig_s / (p->micro_xs_[i_nuclide].total - - p->micro_xs_[i_nuclide].absorption); + score *= flux_deriv + dsig_s / (p->neutron_xs_[i_nuclide].total + - p->neutron_xs_[i_nuclide].absorption); } else { score *= flux_deriv; } break; case SCORE_ABSORPTION: - if (i_nuclide == -1 && p->material_xs_.absorption > 0.0) { + if (i_nuclide == -1 && p->macro_xs_.absorption > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && p->micro_xs_[i_nuc].absorption) { + && p->neutron_xs_[i_nuc].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_a * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / p->material_xs_.absorption; - } else if (p->micro_xs_[i_nuclide].absorption) { + score *= flux_deriv + cum_dsig / p->macro_xs_.absorption; + } else if (p->neutron_xs_[i_nuclide].absorption) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv - + dsig_a / p->micro_xs_[i_nuclide].absorption; + + dsig_a / p->neutron_xs_[i_nuclide].absorption; } else { score *= flux_deriv; } break; case SCORE_FISSION: - if (i_nuclide == -1 && p->material_xs_.fission > 0.0) { + if (i_nuclide == -1 && p->macro_xs_.fission > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && p->micro_xs_[i_nuc].fission) { + && p->neutron_xs_[i_nuc].fission) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_f * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / p->material_xs_.fission; - } else if (p->micro_xs_[i_nuclide].fission) { + score *= flux_deriv + cum_dsig / p->macro_xs_.fission; + } else if (p->neutron_xs_[i_nuclide].fission) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv - + dsig_f / p->micro_xs_[i_nuclide].fission; + + dsig_f / p->neutron_xs_[i_nuclide].fission; } else { score *= flux_deriv; } break; case SCORE_NU_FISSION: - if (i_nuclide == -1 && p->material_xs_.nu_fission > 0.0) { + if (i_nuclide == -1 && p->macro_xs_.nu_fission > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; if (multipole_in_range(&nuc, p->E_last_) - && p->micro_xs_[i_nuc].fission) { - double nu = p->micro_xs_[i_nuc].nu_fission - / p->micro_xs_[i_nuc].fission; + && p->neutron_xs_[i_nuc].fission) { + double nu = p->neutron_xs_[i_nuc].nu_fission + / p->neutron_xs_[i_nuc].fission; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += nu * dsig_f * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / p->material_xs_.nu_fission; - } else if (p->micro_xs_[i_nuclide].fission) { + score *= flux_deriv + cum_dsig / p->macro_xs_.nu_fission; + } else if (p->neutron_xs_[i_nuclide].fission) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv - + dsig_f / p->micro_xs_[i_nuclide].fission; + + dsig_f / p->neutron_xs_[i_nuclide].fission; } else { score *= flux_deriv; } @@ -582,7 +582,7 @@ score_track_derivative(const Particle* p, double distance) // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist // (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist - deriv.flux_deriv -= distance * p->material_xs_.total + deriv.flux_deriv -= distance * p->macro_xs_.total / material.density_gpcc_; break; @@ -591,7 +591,7 @@ score_track_derivative(const Particle* p, double distance) // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist deriv.flux_deriv -= distance - * p->micro_xs_[deriv.diff_nuclide].total; + * p->neutron_xs_[deriv.diff_nuclide].total; break; case DIFF_TEMPERATURE: @@ -660,7 +660,7 @@ void score_collision_derivative(const Particle* p) // phi is proportional to Sigma_s // (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s // (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s - const auto& micro_xs {p->micro_xs_[i_nuc]}; + const auto& micro_xs {p->neutron_xs_[i_nuc]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 57e810dfd2..de851ea3cb 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -349,7 +349,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, if (p->type_ == Particle::Type::neutron || p->type_ == Particle::Type::photon) { - score *= flux / p->material_xs_.total; + score *= flux / p->macro_xs_.total; } else { score = 0.; } @@ -373,9 +373,9 @@ score_general_ce(Particle* p, int i_tally, int start_index, } else { if (i_nuclide >= 0) { - score = p->micro_xs_[i_nuclide].total * atom_density * flux; + score = p->neutron_xs_[i_nuclide].total * atom_density * flux; } else { - score = p->material_xs_.total * flux; + score = p->macro_xs_.total * flux; } } break; @@ -393,7 +393,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, } else { score = p->wgt_last_; } - score *= flux / p->material_xs_.total; + score *= flux / p->macro_xs_.total; } else { score = flux; } @@ -411,11 +411,11 @@ score_general_ce(Particle* p, int i_tally, int start_index, score = p->wgt_last_ * flux; } else { if (i_nuclide >= 0) { - score = (p->micro_xs_[i_nuclide].total - - p->micro_xs_[i_nuclide].absorption) * atom_density * flux; + score = (p->neutron_xs_[i_nuclide].total + - p->neutron_xs_[i_nuclide].absorption) * atom_density * flux; } else { - score = (p->material_xs_.total - - p->material_xs_.absorption) * flux; + score = (p->macro_xs_.total + - p->macro_xs_.absorption) * flux; } } break; @@ -458,26 +458,26 @@ score_general_ce(Particle* p, int i_tally, int start_index, } } else { if (i_nuclide >= 0) { - score = p->micro_xs_[i_nuclide].absorption * atom_density + score = p->neutron_xs_[i_nuclide].absorption * atom_density * flux; } else { - score = p->material_xs_.absorption * flux; + score = p->macro_xs_.absorption * flux; } } break; case SCORE_FISSION: - if (p->material_xs_.absorption == 0) continue; + if (p->macro_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - if (p->micro_xs_[p->event_nuclide_].absorption > 0) { + if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { score = p->wgt_absorb_ - * p->micro_xs_[p->event_nuclide_].fission - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].fission + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } else { score = 0.; } @@ -488,21 +488,21 @@ score_general_ce(Particle* p, int i_tally, int start_index, // weight entering the collision as the estimate for the fission // reaction rate score = p->wgt_last_ - * p->micro_xs_[p->event_nuclide_].fission - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].fission + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } } else { if (i_nuclide >= 0) { - score = p->micro_xs_[i_nuclide].fission * atom_density * flux; + score = p->neutron_xs_[i_nuclide].fission * atom_density * flux; } else { - score = p->material_xs_.fission * flux; + score = p->macro_xs_.fission * flux; } } break; case SCORE_NU_FISSION: - if (p->material_xs_.absorption == 0) continue; + if (p->macro_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { @@ -516,10 +516,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - if (p->micro_xs_[p->event_nuclide_].absorption > 0) { + if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { score = p->wgt_absorb_ - * p->micro_xs_[p->event_nuclide_].nu_fission - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].nu_fission + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } else { score = 0.; } @@ -535,17 +535,17 @@ score_general_ce(Particle* p, int i_tally, int start_index, } } else { if (i_nuclide >= 0) { - score = p->micro_xs_[i_nuclide].nu_fission * atom_density + score = p->neutron_xs_[i_nuclide].nu_fission * atom_density * flux; } else { - score = p->material_xs_.nu_fission * flux; + score = p->macro_xs_.nu_fission * flux; } } break; case SCORE_PROMPT_NU_FISSION: - if (p->material_xs_.absorption == 0) continue; + if (p->macro_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { @@ -559,12 +559,12 @@ score_general_ce(Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - if (p->micro_xs_[p->event_nuclide_].absorption > 0) { + if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { score = p->wgt_absorb_ - * p->micro_xs_[p->event_nuclide_].fission + * p->neutron_xs_[p->event_nuclide_].fission * data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::prompt) - / p->micro_xs_[p->event_nuclide_].absorption * flux; + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } else { score = 0.; } @@ -583,7 +583,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, } } else { if (i_nuclide >= 0) { - score = p->micro_xs_[i_nuclide].fission + score = p->neutron_xs_[i_nuclide].fission * data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::prompt) * atom_density * flux; @@ -595,7 +595,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p->micro_xs_[j_nuclide].fission + score += p->neutron_xs_[j_nuclide].fission * data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::prompt) * atom_density * flux; @@ -607,7 +607,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_DELAYED_NU_FISSION: - if (p->material_xs_.absorption == 0) continue; + if (p->macro_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { @@ -621,7 +621,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - if (p->micro_xs_[p->event_nuclide_].absorption > 0 + if (p->neutron_xs_[p->event_nuclide_].absorption > 0 && data::nuclides[p->event_nuclide_]->fissionable_) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; @@ -634,8 +634,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto yield = data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = p->wgt_absorb_ * yield - * p->micro_xs_[p->event_nuclide_].fission - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].fission + / p->neutron_xs_[p->event_nuclide_].absorption * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } @@ -645,10 +645,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs score = p->wgt_absorb_ - * p->micro_xs_[p->event_nuclide_].fission + * p->neutron_xs_[p->event_nuclide_].fission * data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed) - / p->micro_xs_[p->event_nuclide_].absorption *flux; + / p->neutron_xs_[p->event_nuclide_].absorption *flux; } } } else { @@ -694,7 +694,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto d = filt.groups_[d_bin]; auto yield = data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p->micro_xs_[i_nuclide].fission * yield + score = p->neutron_xs_[i_nuclide].fission * yield * atom_density * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } @@ -702,7 +702,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, } else { // If the delayed group filter is not present, compute the score // by multiplying the delayed-nu-fission macro xs by the flux - score = p->micro_xs_[i_nuclide].fission + score = p->neutron_xs_[i_nuclide].fission * data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed) * atom_density * flux; @@ -723,7 +723,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto d = filt.groups_[d_bin]; auto yield = data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p->micro_xs_[j_nuclide].fission * yield + score = p->neutron_xs_[j_nuclide].fission * yield * atom_density * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -738,7 +738,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p->micro_xs_[j_nuclide].fission + score += p->neutron_xs_[j_nuclide].fission * data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed) * atom_density * flux; @@ -751,14 +751,14 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_DECAY_RATE: - if (p->material_xs_.absorption == 0) continue; + if (p->macro_xs_.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->micro_xs_[p->event_nuclide_].absorption > 0 + if (p->neutron_xs_[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; if (tally.delayedgroup_filter_ != C_NONE) { @@ -773,8 +773,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; score = p->wgt_absorb_ * yield - * p->micro_xs_[p->event_nuclide_].fission - / p->micro_xs_[p->event_nuclide_].absorption + * p->neutron_xs_[p->event_nuclide_].fission + / p->neutron_xs_[p->event_nuclide_].absorption * rate * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -796,8 +796,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; score += rate * p->wgt_absorb_ - * p->micro_xs_[p->event_nuclide_].fission * yield - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].fission * yield + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } } } @@ -854,7 +854,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p->micro_xs_[i_nuclide].fission * yield * flux + score = p->neutron_xs_[i_nuclide].fission * yield * flux * atom_density * rate; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } @@ -870,7 +870,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += p->micro_xs_[i_nuclide].fission * flux + score += p->neutron_xs_[i_nuclide].fission * flux * yield * atom_density * rate; } } @@ -894,7 +894,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p->micro_xs_[j_nuclide].fission * yield + score = p->neutron_xs_[j_nuclide].fission * yield * flux * atom_density * rate; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -922,7 +922,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += p->micro_xs_[j_nuclide].fission + score += p->neutron_xs_[j_nuclide].fission * yield * atom_density * flux * rate; } } @@ -935,7 +935,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_KAPPA_FISSION: - if (p->material_xs_.absorption == 0.) continue; + if (p->macro_xs_.absorption == 0.) continue; score = 0.; // Kappa-fission values are determined from the Q-value listed for the // fission cross section. @@ -945,12 +945,12 @@ score_general_ce(Particle* p, int i_tally, int start_index, // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->micro_xs_[p->event_nuclide_].absorption > 0 + if (p->neutron_xs_[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; score = p->wgt_absorb_ * rxn.q_value_ - * p->micro_xs_[p->event_nuclide_].fission - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].fission + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events @@ -959,12 +959,12 @@ score_general_ce(Particle* p, int i_tally, int start_index, // weight entering the collision as the estimate for the fission // reaction rate const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->micro_xs_[p->event_nuclide_].absorption > 0 + if (p->neutron_xs_[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; score = p->wgt_last_ * rxn.q_value_ - * p->micro_xs_[p->event_nuclide_].fission - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].fission + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } } } else { @@ -972,7 +972,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, const auto& nuc {*data::nuclides[i_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = rxn.q_value_ * p->micro_xs_[i_nuclide].fission + score = rxn.q_value_ * p->neutron_xs_[i_nuclide].fission * atom_density * flux; } } else { @@ -984,7 +984,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, const auto& nuc {*data::nuclides[j_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score += rxn.q_value_ * p->micro_xs_[j_nuclide].fission + score += rxn.q_value_ * p->neutron_xs_[j_nuclide].fission * atom_density * flux; } } @@ -1007,9 +1007,9 @@ score_general_ce(Particle* p, int i_tally, int start_index, score = p->wgt_last_ * flux; } else { if (i_nuclide >= 0) { - if (p->micro_xs_[i_nuclide].elastic == CACHE_INVALID) + if (p->neutron_xs_[i_nuclide].elastic == CACHE_INVALID) data::nuclides[i_nuclide]->calculate_elastic_xs(*p); - score = p->micro_xs_[i_nuclide].elastic * atom_density * flux; + score = p->neutron_xs_[i_nuclide].elastic * atom_density * flux; } else { score = 0.; if (p->material_ != MATERIAL_VOID) { @@ -1017,9 +1017,9 @@ score_general_ce(Particle* p, int i_tally, int start_index, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - if (p->micro_xs_[j_nuclide].elastic == CACHE_INVALID) + if (p->neutron_xs_[j_nuclide].elastic == CACHE_INVALID) data::nuclides[j_nuclide]->calculate_elastic_xs(*p); - score += p->micro_xs_[j_nuclide].elastic * atom_density + score += p->neutron_xs_[j_nuclide].elastic * atom_density * flux; } } @@ -1030,7 +1030,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_FISS_Q_PROMPT: case SCORE_FISS_Q_RECOV: //continue; - if (p->material_xs_.absorption == 0.) continue; + if (p->macro_xs_.absorption == 0.) continue; score = 0.; if (tally.estimator_ == ESTIMATOR_ANALOG) { if (settings::survival_biasing) { @@ -1038,7 +1038,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->micro_xs_[p->event_nuclide_].absorption > 0) { + if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) @@ -1048,8 +1048,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, q_value = (*nuc.fission_q_recov_)(p->E_last_); } score = p->wgt_absorb_ * q_value - * p->micro_xs_[p->event_nuclide_].fission - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].fission + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events @@ -1058,7 +1058,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, // weight entering the collision as the estimate for the fission // reaction rate const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->micro_xs_[p->event_nuclide_].absorption > 0) { + if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) @@ -1068,8 +1068,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, q_value = (*nuc.fission_q_recov_)(p->E_last_); } score = p->wgt_last_ * q_value - * p->micro_xs_[p->event_nuclide_].fission - / p->micro_xs_[p->event_nuclide_].absorption * flux; + * p->neutron_xs_[p->event_nuclide_].fission + / p->neutron_xs_[p->event_nuclide_].absorption * flux; } } } else { @@ -1083,7 +1083,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, if (nuc.fission_q_recov_) q_value = (*nuc.fission_q_recov_)(p->E_last_); } - score = q_value * p->micro_xs_[i_nuclide].fission + score = q_value * p->neutron_xs_[i_nuclide].fission * atom_density * flux; } else { if (p->material_ != MATERIAL_VOID) { @@ -1100,7 +1100,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, if (nuc.fission_q_recov_) q_value = (*nuc.fission_q_recov_)(p->E_last_); } - score += q_value * p->micro_xs_[j_nuclide].fission + score += q_value * p->neutron_xs_[j_nuclide].fission * atom_density * flux; } } @@ -1130,7 +1130,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, case N_4N: m = 5; break; } if (i_nuclide >= 0) { - score = p->micro_xs_[i_nuclide].reaction[m] * atom_density + score = p->neutron_xs_[i_nuclide].reaction[m] * atom_density * flux; } else { score = 0.; @@ -1139,7 +1139,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p->micro_xs_[j_nuclide].reaction[m] + score += p->neutron_xs_[j_nuclide].reaction[m] * atom_density * flux; } } @@ -1164,10 +1164,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->micro_xs_[i_nuclide].index_temp; + auto i_temp = p->neutron_xs_[i_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->micro_xs_[i_nuclide].index_grid; - auto f = p->micro_xs_[i_nuclide].interp_factor; + auto i_grid = p->neutron_xs_[i_nuclide].index_grid; + auto f = p->neutron_xs_[i_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { score = ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1184,10 +1184,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->micro_xs_[j_nuclide].index_temp; + auto i_temp = p->neutron_xs_[j_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->micro_xs_[j_nuclide].index_grid; - auto f = p->micro_xs_[j_nuclide].interp_factor; + auto i_grid = p->neutron_xs_[j_nuclide].index_grid; + auto f = p->neutron_xs_[j_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { score += ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1291,7 +1291,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } else { score = p->wgt_last_; } - score *= flux / p->material_xs_.total; + score *= flux / p->macro_xs_.total; } else { score = flux; } @@ -1320,7 +1320,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = get_nuclide_xs(i_nuclide, MG_GET_XS_TOTAL, p_g) * atom_density * flux; } else { - score = p->material_xs_.total * flux; + score = p->macro_xs_.total * flux; } } break; @@ -1438,7 +1438,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = atom_density * flux * get_nuclide_xs(i_nuclide, MG_GET_XS_ABSORPTION, p_g); } else { - score = p->material_xs_.absorption * flux; + score = p->macro_xs_.absorption * flux; } } break; @@ -2127,9 +2127,9 @@ void score_collision_tally(Particle* p) // Determine the collision estimate of the flux double flux; if (!settings::survival_biasing) { - flux = p->wgt_last_ / p->material_xs_.total; + flux = p->wgt_last_ / p->macro_xs_.total; } else { - flux = (p->wgt_last_ + p->wgt_absorb_) / p->material_xs_.total; + flux = (p->wgt_last_ + p->wgt_absorb_) / p->macro_xs_.total; } for (auto i_tally : model::active_collision_tallies) { From 429f186223e146ee3eb54be9df2917f7b92961b8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 28 Oct 2018 16:11:33 -0600 Subject: [PATCH 03/14] Fix grammatical error in photon docs --- docs/source/methods/photon_physics.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index 23809f6745..2d9a62ef37 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -5,7 +5,7 @@ Photon Physics ============== Photons, being neutral particles, behave much in the same manner as neutrons, -traveling in straight lines and experiencing occasional collisions which change +traveling in straight lines and experiencing occasional collisions that change their energy and direction. Photons undergo four basic interactions as they pass through matter: coherent (Rayleigh) scattering, incoherent (Compton) scattering, photoelectric effect, and pair/triplet production. Photons with energy in the From d70cfda07c29958a58380ca7d4c015bc5172c137 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Mar 2019 11:34:33 -0500 Subject: [PATCH 04/14] A few style guide additions --- docs/source/devguide/styleguide.rst | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/docs/source/devguide/styleguide.rst b/docs/source/devguide/styleguide.rst index 0f7669a244..cbf665cc8e 100644 --- a/docs/source/devguide/styleguide.rst +++ b/docs/source/devguide/styleguide.rst @@ -12,6 +12,11 @@ adding new code in OpenMC. C++ --- +Indentation +----------- + +Use two spaces per indentation level. + Miscellaneous ------------- @@ -126,6 +131,15 @@ single declaration to avoid confusion: Curly braces ------------ +For a class declaration, the opening brace should be on the same line that +lists the name of the class. + +.. code-block:: C++ + + class Matrix { + ... + }; + For a function definition, the opening and closing braces should each be on their own lines. This helps distinguish function code from the argument list. If the entire function fits on one or two lines, then the braces can be on the @@ -210,11 +224,18 @@ Use of third-party Python packages should be limited to numpy_, scipy_, matplotlib_, pandas_, and h5py_. Use of other third-party packages must be implemented as optional dependencies rather than required dependencies. +Prefer pathlib_ when working with filesystem paths over functions in the os_ +module or other standard-library modules. Functions that accept arguments that +represent a filesystem path should work with both strings and Path_ objects. + .. _C++ Core Guidelines: http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines .. _PEP8: https://www.python.org/dev/peps/pep-0008/ -.. _numpydoc: https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt +.. _numpydoc: https://numpydoc.readthedocs.io/en/latest/format.html .. _numpy: http://www.numpy.org/ .. _scipy: https://www.scipy.org/ .. _matplotlib: https://matplotlib.org/ .. _pandas: https://pandas.pydata.org/ -.. _h5py: http://www.h5py.org/ +.. _h5py: https://www.h5py.org/ +.. _pathlib: https://docs.python.org/3/library/pathlib.html +.. _os: https://docs.python.org/3/library/os.html +.. _Path: https://docs.python.org/3/library/pathlib.html#pathlib.Path From 86d187a0a71cfdeee7fd5de4ace49764ec1f1725 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Mar 2019 12:26:53 -0500 Subject: [PATCH 05/14] Turn Particle.coord_ into a std::vector (remove MAX_COORD) --- CMakeLists.txt | 4 ---- docs/source/usersguide/install.rst | 3 --- include/openmc/geometry.h | 3 ++- include/openmc/geometry_aux.h | 9 ++++----- include/openmc/particle.h | 4 ++-- src/geometry.cpp | 3 ++- src/geometry_aux.cpp | 10 ++-------- src/particle.cpp | 4 +++- 8 files changed, 15 insertions(+), 25 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f943d09697..930b1e0bb8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -20,9 +20,6 @@ option(optimize "Turn on all compiler optimization flags" OFF) option(coverage "Compile with coverage analysis flags" OFF) option(dagmc "Enable support for DAGMC (CAD) geometry" OFF) -# Maximum number of nested coordinates levels -set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels") - #=============================================================================== # MPI for distributed-memory parallelism #=============================================================================== @@ -260,7 +257,6 @@ target_include_directories(libopenmc # Set compile flags target_compile_options(libopenmc PRIVATE ${cxxflags}) -target_compile_definitions(libopenmc PUBLIC -DMAX_COORD=${maxcoord}) if (HDF5_IS_PARALLEL) target_compile_definitions(libopenmc PRIVATE -DPHDF5) endif() diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index d44ee95940..96b4f5024c 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -243,9 +243,6 @@ coverage Compile and link code instrumented for coverage analysis. This is typically used in conjunction with gcov_. -maxcoord - Maximum number of nested coordinate levels in geometry. Defaults to 10. - To set any of these options (e.g. turning on debug mode), the following form should be used: diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index a2da2c5f22..9bb1a9ae7f 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -15,7 +15,8 @@ namespace openmc { namespace model { -extern "C" int root_universe; +extern int root_universe; //!< Index of root universe +extern int n_coord_levels; //!< Number of CSG coordinate levels extern std::vector overlap_check_count; diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index 5898c61f31..ffcdcf885f 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -57,7 +57,7 @@ void finalize_geometry(std::vector>& nuc_temps, //! \return The index of the root universe. //============================================================================== -extern "C" int32_t find_root_universe(); +int32_t find_root_universe(); //============================================================================== //! Populate all data structures needed for distribcells. @@ -74,7 +74,7 @@ void prepare_distribcell(); //! the root universe). //============================================================================== -extern "C" void count_cell_instances(int32_t univ_indx); +void count_cell_instances(int32_t univ_indx); //============================================================================== //! Recursively search through universes and count universe instances. @@ -84,8 +84,7 @@ extern "C" void count_cell_instances(int32_t univ_indx); //! search_univ. //============================================================================== -extern "C" int -count_universe_instances(int32_t search_univ, int32_t target_univ_id); +int count_universe_instances(int32_t search_univ, int32_t target_univ_id); //============================================================================== //! Build a character array representing the path to a distribcell instance. @@ -107,7 +106,7 @@ distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset); //! \return The number of coordinate levels. //============================================================================== -extern "C" int maximum_levels(int32_t univ); +int maximum_levels(int32_t univ); //============================================================================== //! Deallocates global vectors and maps for cells, universes, and lattices. diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 2130b6454a..ef1e076bb6 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -229,11 +229,11 @@ public: int n_coord_ {1}; //!< number of current coordinate levels int cell_instance_; //!< offset for distributed properties - LocalCoord coord_[MAX_COORD]; //!< coordinates for all levels + std::vector coord_; //!< coordinates for all levels // Particle coordinates before crossing a surface int n_coord_last_ {1}; //!< number of current coordinates - int cell_last_[MAX_COORD]; //!< coordinates for all levels + std::vector cell_last_; //!< coordinates for all levels // Energy data double E_; //!< post-collision energy in eV diff --git a/src/geometry.cpp b/src/geometry.cpp index 3c1c11b0e8..cde2cdc1b6 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -22,6 +22,7 @@ namespace openmc { namespace model { int root_universe {-1}; +int n_coord_levels; std::vector overlap_check_count; @@ -257,7 +258,7 @@ find_cell(Particle* p, bool use_neighbor_lists) } // Reset all the deeper coordinate levels. - for (int i = p->n_coord_; i < MAX_COORD; i++) { + for (int i = p->n_coord_; i < p->coord_.size(); i++) { p->coord_[i].reset(); } diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index c560743f86..e19b5f946c 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -217,14 +217,8 @@ void finalize_geometry(std::vector>& nuc_temps, // Determine desired temperatures for each nuclide and S(a,b) table get_temperatures(nuc_temps, thermal_temps); - // Check to make sure there are not too many nested coordinate levels in the - // geometry since the coordinate list is statically allocated for performance - // reasons - if (maximum_levels(model::root_universe) > MAX_COORD) { - fatal_error("Too many nested coordinate levels in the geometry. " - "Try increasing the maximum number of coordinate levels by " - "providing the CMake -Dmaxcoord= option."); - } + // Determine number of nested coordinate levels in the geometry + model::n_coord_levels = maximum_levels(model::root_universe); } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index 88ff251dd5..c6f73f4f12 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -50,7 +50,9 @@ LocalCoord::reset() Particle::Particle() { - // Clear coordinate lists + // Create and clear coordinate levels + coord_.resize(model::n_coord_levels); + cell_last_.resize(model::n_coord_levels); clear(); for (int& n : n_delayed_bank_) { From d980b4a5ea3e6050031825f6318b55c52a613c18 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Mar 2019 07:37:33 -0500 Subject: [PATCH 06/14] Make secondary_bank_ a vector and improve interfaces --- include/openmc/bank.h | 4 +-- include/openmc/mesh.h | 7 ++-- include/openmc/particle.h | 3 +- include/openmc/physics.h | 4 ++- include/openmc/physics_mg.h | 9 +++-- src/bank.cpp | 2 -- src/eigenvalue.cpp | 56 ++++++++++++++---------------- src/mesh.cpp | 31 +++++------------ src/particle.cpp | 23 ++++++------- src/physics.cpp | 52 ++++++++-------------------- src/physics_mg.cpp | 65 ++++++++++------------------------- src/simulation.cpp | 12 +++---- src/tallies/tally_scoring.cpp | 6 ++-- 13 files changed, 95 insertions(+), 179 deletions(-) diff --git a/include/openmc/bank.h b/include/openmc/bank.h index b7fab9ac66..4597369c78 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -20,15 +20,13 @@ namespace openmc { namespace simulation { -extern "C" int64_t n_bank; - extern std::vector source_bank; extern std::vector fission_bank; #ifdef _OPENMP extern std::vector master_fission_bank; #endif -#pragma omp threadprivate(fission_bank, n_bank) +#pragma omp threadprivate(fission_bank) } // namespace simulation diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3fc7db7816..f53bba224d 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -96,14 +96,11 @@ public: //! Count number of bank sites in each mesh bin / energy bin // - //! \param[in] n Number of bank sites //! \param[in] bank Array of bank sites - //! \param[in] n_energy Number of energies - //! \param[in] energies Array of energies //! \param[out] Whether any bank sites are outside the mesh //! \return Array indicating number of sites in each mesh/energy bin - xt::xarray count_sites(int64_t n, const Particle::Bank* bank, - int n_energy, const double* energies, bool* outside) const; + xt::xarray count_sites(const std::vector& bank, + bool* outside) const; int id_ {-1}; //!< User-specified ID int n_dimension_; //!< Number of dimensions diff --git a/include/openmc/particle.h b/include/openmc/particle.h index ef1e076bb6..d52c67ea55 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -285,8 +285,7 @@ public: bool write_track_ {false}; // Secondary particles created - int64_t n_secondary_ {}; - Bank secondary_bank_[MAX_SECONDARY]; + std::vector secondary_bank_; }; } // namespace openmc diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 5c8d64f5f9..107380eb84 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -7,6 +7,8 @@ #include "openmc/position.h" #include "openmc/reaction.h" +#include + namespace openmc { //============================================================================== @@ -47,7 +49,7 @@ int sample_nuclide(const Particle* p); //! Determine the average total, prompt, and delayed neutrons produced from //! fission and creates appropriate bank sites. void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, - Particle::Bank* bank_array, int64_t* bank_size, int64_t bank_capacity); + std::vector& bank); int sample_element(Particle* p); diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index cdf1f7089c..f246cb6ecb 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -8,6 +8,8 @@ #include "openmc/particle.h" #include "openmc/nuclide.h" +#include + namespace openmc { //! \brief samples particle behavior after a collision event. @@ -31,12 +33,9 @@ scatter(Particle* p); //! \brief Determines the average total, prompt and delayed neutrons produced //! from fission and creates the appropriate bank sites. //! \param p Particle to operate on -//! \param bank_array The particle bank to populate -//! \param size_bank Number of particles currently in the bank -//! \param bank_array_size Allocated size of the bank +//! \param bank The particle bank to populate void -create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank, - int64_t bank_array_size); +create_fission_sites(Particle* p, std::vector& bank); //! \brief Handles an absorption event //! \param p Particle to operate on diff --git a/src/bank.cpp b/src/bank.cpp index 8b7b436889..fc48a77615 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -16,8 +16,6 @@ namespace openmc { namespace simulation { -int64_t n_bank; - std::vector source_bank; std::vector fission_bank; #ifdef _OPENMP diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index c995a02d7d..1701955295 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -23,6 +23,7 @@ #include // for min #include #include // for sqrt, abs, pow +#include // for back_inserter #include namespace openmc { @@ -81,21 +82,22 @@ void synchronize_bank() #ifdef OPENMC_MPI int64_t start = 0; - MPI_Exscan(&simulation::n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm); + int64_t n_bank = simulation::fission_bank.size(); + MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm); // While we would expect the value of start on rank 0 to be 0, the MPI // standard says that the receive buffer on rank 0 is undefined and not // significant if (mpi::rank == 0) start = 0; - int64_t finish = start + simulation::n_bank; + int64_t finish = start + simulation::fission_bank.size(); int64_t total = finish; MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm); #else int64_t start = 0; - int64_t finish = simulation::n_bank; - int64_t total = simulation::n_bank; + int64_t finish = simulation::fission_bank.size(); + int64_t total = simulation::fission_bank.size(); #endif // If there are not that many particles per generation, it's possible that no @@ -103,7 +105,7 @@ void synchronize_bank() // extra logic to treat this circumstance, we really want to ensure the user // runs enough particles to avoid this in the first place. - if (simulation::n_bank == 0) { + if (simulation::fission_bank.empty()) { fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank)); } @@ -134,21 +136,21 @@ void synchronize_bank() int64_t index_temp = 0; std::vector temp_sites(3*simulation::work); - for (int64_t i = 0; i < simulation::n_bank; ++i) { + for (const auto& site : simulation::fission_bank) { // If there are less than n_particles particles banked, automatically add // int(n_particles/total) sites to temp_sites. For example, if you need // 1000 and 300 were banked, this would add 3 source sites per banked site // and the remaining 100 would be randomly sampled. if (total < settings::n_particles) { for (int64_t j = 1; j <= settings::n_particles / total; ++j) { - temp_sites[index_temp] = simulation::fission_bank[i]; + temp_sites[index_temp] = site; ++index_temp; } } // Randomly sample sites needed if (prn() < p_sample) { - temp_sites[index_temp] = simulation::fission_bank[i]; + temp_sites[index_temp] = site; ++index_temp; } } @@ -189,7 +191,7 @@ void synchronize_bank() // fission bank sites_needed = settings::n_particles - finish; for (int i = 0; i < sites_needed; ++i) { - int i_bank = simulation::n_bank - sites_needed + i; + int i_bank = simulation::fission_bank.size() - sites_needed + i; temp_sites[index_temp] = simulation::fission_bank[i_bank]; ++index_temp; } @@ -346,38 +348,30 @@ void calculate_average_keff() #ifdef _OPENMP void join_bank_from_threads() { - // Initialize the total number of fission bank sites - int64_t total = 0; - -#pragma omp parallel + #pragma omp parallel { // Copy thread fission bank sites to one shared copy -#pragma omp for ordered schedule(static) + #pragma omp for ordered schedule(static) for (int i = 0; i < simulation::n_threads; ++i) { -#pragma omp ordered + #pragma omp ordered { std::copy( - &simulation::fission_bank[0], - &simulation::fission_bank[0] + simulation::n_bank, - &simulation::master_fission_bank[total] + simulation::fission_bank.cbegin(), + simulation::fission_bank.cend(), + std::back_inserter(simulation::master_fission_bank) ); - total += simulation::n_bank; } } // Make sure all threads have made it to this point -#pragma omp barrier + #pragma omp barrier // Now copy the shared fission bank sites back to the master thread's copy. if (simulation::thread_id == 0) { - simulation::n_bank = total; - std::copy( - &simulation::master_fission_bank[0], - &simulation::master_fission_bank[0] + simulation::n_bank, - &simulation::fission_bank[0] - ); + simulation::fission_bank = simulation::master_fission_bank; + simulation::master_fission_bank.clear(); } else { - simulation::n_bank = 0; + simulation::fission_bank.clear(); } } } @@ -537,8 +531,8 @@ void shannon_entropy() // Get source weight in each mesh bin bool sites_outside; - xt::xtensor p = m->count_sites(simulation::n_bank, - simulation::fission_bank.data(), 0, nullptr, &sites_outside); + xt::xtensor p = m->count_sites(simulation::fission_bank, + &sites_outside); // display warning message if there were sites outside entropy box if (sites_outside) { @@ -577,8 +571,8 @@ void ufs_count_sites() } else { // count number of source sites in each ufs mesh cell bool sites_outside; - simulation::source_frac = m->count_sites(simulation::work, - simulation::source_bank.data(), 0, nullptr, &sites_outside); + simulation::source_frac = m->count_sites(simulation::source_bank, + &sites_outside); // Check for sites outside of the mesh if (mpi::master && sites_outside) { diff --git a/src/mesh.cpp b/src/mesh.cpp index b75a0aad1f..ddfa272a6d 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -732,25 +732,21 @@ void RegularMesh::to_hdf5(hid_t group) const close_group(mesh_group); } -xt::xarray RegularMesh::count_sites(int64_t n, const Particle::Bank* bank, - int n_energy, const double* energies, bool* outside) const +xt::xarray +RegularMesh::count_sites(const std::vector& bank, + bool* outside) const { // Determine shape of array for counts std::size_t m = xt::prod(shape_)(); - std::vector shape; - if (n_energy > 0) { - shape = {m, static_cast(n_energy - 1)}; - } else { - shape = {m}; - } + std::vector shape = {m}; // Create array of zeros xt::xarray cnt {shape, 0.0}; bool outside_ = false; - for (int64_t i = 0; i < n; ++i) { + for (const auto& site : bank) { // determine scoring bin for entropy mesh - int mesh_bin = get_bin(bank[i].r); + int mesh_bin = get_bin(site.r); // if outside mesh, skip particle if (mesh_bin < 0) { @@ -758,19 +754,8 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Particle::Bank* ban continue; } - if (n_energy > 0) { - double E = bank[i].E; - if (E >= energies[0] && E <= energies[n_energy - 1]) { - // determine energy bin - int e_bin = lower_bound_index(energies, energies + n_energy, E); - - // Add to appropriate bin - cnt(mesh_bin, e_bin) += bank[i].wgt; - } - } else { - // Add to appropriate bin - cnt(mesh_bin) += bank[i].wgt; - } + // Add to appropriate bin + cnt(mesh_bin) += site.wgt; } // Create copy of count data diff --git a/src/particle.cpp b/src/particle.cpp index c6f73f4f12..92045e2cdb 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -75,17 +75,14 @@ Particle::clear() void Particle::create_secondary(Direction u, double E, Type type) { - if (n_secondary_ == MAX_SECONDARY) { - fatal_error("Too many secondary particles created."); - } + secondary_bank_.emplace_back(); - int64_t n = n_secondary_; - secondary_bank_[n].particle = type; - secondary_bank_[n].wgt = wgt_; - secondary_bank_[n].r = this->r(); - secondary_bank_[n].u = u; - secondary_bank_[n].E = settings::run_CE ? E : g_; - ++n_secondary_; + auto& bank {secondary_bank_.back()}; + bank.particle = type; + bank.wgt = wgt_; + bank.r = this->r(); + bank.u = u; + bank.E = settings::run_CE ? E : g_; } void @@ -360,10 +357,10 @@ Particle::transport() // Check for secondary particles if this particle is dead if (!alive_) { // If no secondary particles, break out of event loop - if (n_secondary_ == 0) break; + if (secondary_bank_.empty()) break; - this->from_source(&secondary_bank_[n_secondary_ - 1]); - --n_secondary_; + this->from_source(&secondary_bank_.back()); + secondary_bank_.pop_back(); n_event = 0; // Enter new particle in particle track file diff --git a/src/physics.cpp b/src/physics.cpp index 4cc9b12a03..7d44aa17c6 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -91,12 +91,10 @@ void sample_neutron_reaction(Particle* p) if (nuc->fissionable_) { Reaction* rx = sample_fission(i_nuclide, p); if (settings::run_mode == RUN_MODE_EIGENVALUE) { - create_fission_sites(p, i_nuclide, rx, simulation::fission_bank.data(), - &simulation::n_bank, simulation::fission_bank.size()); + create_fission_sites(p, i_nuclide, rx, simulation::fission_bank); } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE && settings::create_fission_neutrons) { - create_fission_sites(p, i_nuclide, rx, p->secondary_bank_, - &p->n_secondary_, MAX_SECONDARY); + create_fission_sites(p, i_nuclide, rx, p->secondary_bank_); } } @@ -137,10 +135,8 @@ void sample_neutron_reaction(Particle* p) void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, - Particle::Bank* bank_array, int64_t* size_bank, int64_t bank_capacity) + std::vector& bank) { - // TODO: Heat generation from fission - // If uniform fission source weighting is turned on, we increase or decrease // the expected number of fission sites produced double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; @@ -153,47 +149,30 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, int nu = static_cast(nu_t); if (prn() <= (nu_t - nu)) ++nu; - // Check for the bank size getting hit. For fixed source calculations, this - // is a fatal error; for eigenvalue calculations, it just means that k-eff - // was too high for a single batch. - if (*size_bank + nu > bank_capacity) { - if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { - throw std::runtime_error{"Secondary particle bank size limit reached." - " If you are running a subcritical multiplication problem," - " k-effective may be too close to one."}; - } else { - if (mpi::master) { - warning("Maximum number of sites in fission bank reached. This can" - " result in irreproducible results using different numbers of" - " processes/threads."); - } - } - } - // Begin banking the source neutrons // First, if our bank is full then don't continue - if (nu == 0 || *size_bank == bank_capacity) return; + if (nu == 0) return; // Initialize the counter of delayed neutrons encountered for each delayed // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; p->fission_ = true; - for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) { + for (int i = 0; i < nu; ++i) { + // Create new bank site and get reference to last element + bank.emplace_back(); + auto& site {bank.back()}; + // Bank source neutrons by copying the particle data - bank_array[i].r = p->r(); - - // Set that the bank particle is a neutron - bank_array[i].particle = Particle::Type::neutron; - - // Set the weight of the fission bank site - bank_array[i].wgt = 1. / weight; + site.r = p->r(); + site.particle = Particle::Type::neutron; + site.wgt = 1. / weight; // Sample delayed group and angle/energy for fission reaction - sample_fission_neutron(i_nuclide, rx, p->E_, &bank_array[i]); + sample_fission_neutron(i_nuclide, rx, p->E_, &site); // Set the delayed group on the particle as well - p->delayed_group_ = bank_array[i].delayed_group; + p->delayed_group_ = site.delayed_group; // Increment the number of neutrons born delayed if (p->delayed_group_ > 0) { @@ -201,9 +180,6 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, } } - // Increment number of bank sites - *size_bank = std::min(*size_bank + nu, bank_capacity); - // Store the total weight banked for analog fission tallies p->n_bank_ = nu; p->wgt_bank_ = nu / weight; diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index affe0ec197..2a6b700f03 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -48,13 +48,10 @@ sample_reaction(Particle* p) if (model::materials[p->material_]->fissionable_) { if (settings::run_mode == RUN_MODE_EIGENVALUE) { - create_fission_sites( - p, simulation::fission_bank.data(), &simulation::n_bank, - simulation::fission_bank.size()); + create_fission_sites(p, simulation::fission_bank); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { - create_fission_sites(p, p->secondary_bank_, &(p->n_secondary_), - MAX_SECONDARY); + create_fission_sites(p, p->secondary_bank_); } } @@ -102,11 +99,8 @@ scatter(Particle* p) } void -create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank, - int64_t bank_array_size) +create_fission_sites(Particle* p, std::vector& bank) { - // TODO: Heat generation from fission - // If uniform fission source weighting is turned on, we increase or decrease // the expected number of fission sites produced double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; @@ -121,62 +115,42 @@ create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank nu++; } - // Check for the bank size getting hit. For fixed source calculations, this - // is a fatal error; for eigenvalue calculations, it just means that k-eff - // was too high for a single batch. - if (*size_bank + nu > bank_array_size) { - if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { - throw std::runtime_error{"Secondary particle bank size limit reached." - " If you are running a subcritical multiplication problem," - " k-effective may be too close to one."}; - } else { - if (mpi::master) { - std::stringstream msg; - msg << "Maximum number of sites in fission bank reached. This can" - " result in irreproducible results using different numbers of" - " processes/threads."; - warning(msg); - } - } - } - // Begin banking the source neutrons // First, if our bank is full then don't continue - if ((nu == 0) || (*size_bank == bank_array_size)) return; + if (nu == 0) return; // Initialize the counter of delayed neutrons encountered for each delayed // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; p->fission_ = true; - for (size_t i = static_cast(*size_bank); - i < static_cast(std::min(*size_bank + nu, bank_array_size)); i++) { + for (int i = 0; i < nu; ++i) { + // Create new bank site and get reference to last element + bank.emplace_back(); + auto& site {bank.back()}; + // Bank source neutrons by copying the particle data - bank_array[i].r = p->r(); - - // Set that the bank particle is a neutron - bank_array[i].particle = Particle::Type::neutron; - - // Set the weight of the fission bank site - bank_array[i].wgt = 1. / weight; + site.r = p->r(); + site.particle = Particle::Type::neutron; + site.wgt = 1. / weight; // Sample the cosine of the angle, assuming fission neutrons are emitted // isotropically - double mu = 2. * prn() - 1.; + double mu = 2.*prn() - 1.; // Sample the azimuthal angle uniformly in [0, 2.pi) double phi = 2. * PI * prn(); - bank_array[i].u.x = mu; - bank_array[i].u.y = std::sqrt(1. - mu * mu) * std::cos(phi); - bank_array[i].u.z = std::sqrt(1. - mu * mu) * std::sin(phi); + site.u.x = mu; + site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi); + site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi); // Sample secondary energy distribution for the fission reaction and set // the energy in the fission bank int dg; int gout; data::macro_xs[p->material_].sample_fission_energy(p->g_ - 1, dg, gout); - bank_array[i].E = gout + 1; - bank_array[i].delayed_group = dg + 1; + site.E = gout + 1; + site.delayed_group = dg + 1; // Set the delayed group on the particle as well p->delayed_group_ = dg + 1; @@ -187,9 +161,6 @@ create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank } } - // Increment number of bank sites - *size_bank = std::min(*size_bank + nu, bank_array_size); - // Store the total weight banked for analog fission tallies p->n_bank_ = nu; p->wgt_bank_ = nu / weight; diff --git a/src/simulation.cpp b/src/simulation.cpp index b4db5e0e6d..fca0c85282 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -279,14 +279,14 @@ void allocate_banks() { simulation::thread_id = omp_get_thread_num(); if (simulation::thread_id == 0) { - simulation::fission_bank.resize(3*simulation::work); + simulation::fission_bank.reserve(3*simulation::work); } else { - simulation::fission_bank.resize(3*simulation::work / simulation::n_threads); + simulation::fission_bank.reserve(3*simulation::work / simulation::n_threads); } } - simulation::master_fission_bank.resize(3*simulation::work); + simulation::master_fission_bank.reserve(3*simulation::work); #else - simulation::fission_bank.resize(3*simulation::work); + simulation::fission_bank.reserve(3*simulation::work); #endif } } @@ -386,8 +386,8 @@ void finalize_batch() void initialize_generation() { if (settings::run_mode == RUN_MODE_EIGENVALUE) { - // Reset number of fission bank sites - simulation::n_bank = 0; + // Clear out the fission bank + simulation::fission_bank.clear(); // Count source sites if using uniform fission source weighting if (settings::ufs_on) ufs_count_sites(); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index de851ea3cb..7e235fa458 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -192,7 +192,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) // loop over number of particles banked for (auto i = 0; i < p->n_bank_; ++i) { - auto i_bank = simulation::n_bank - p->n_bank_ + i; + auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; // get the delayed group @@ -813,7 +813,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, // contribution to the fission bank to the score. score = 0.; for (auto i = 0; i < p->n_bank_; ++i) { - auto i_bank = simulation::n_bank - p->n_bank_ + i; + auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; auto g = bank.delayed_group; if (g != 0) { @@ -1785,7 +1785,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // contribution to the fission bank to the score. score = 0.; for (auto i = 0; i < p->n_bank_; ++i) { - auto i_bank = simulation::n_bank - p->n_bank_ + i; + auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; auto g = bank.delayed_group; if (g != 0) { From 49bf4523c37e7a449eb133cc00e020c4cc784a0a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Mar 2019 10:11:39 -0500 Subject: [PATCH 07/14] Remove thread_id and n_threads global variables --- docs/source/io_formats/settings.rst | 9 --------- include/openmc/simulation.h | 6 +----- openmc/settings.py | 19 ------------------- src/bank.cpp | 2 +- src/eigenvalue.cpp | 10 ++++++++-- src/initialize.cpp | 6 +++--- src/output.cpp | 2 ++ src/particle.cpp | 4 ++-- src/physics_mg.cpp | 15 ++++++--------- src/random_lcg.cpp | 2 +- src/settings.cpp | 16 ++-------------- src/simulation.cpp | 26 +++++++++++--------------- tests/unit_tests/test_settings.py | 1 - 13 files changed, 37 insertions(+), 81 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index e5c0ca1600..9c68a9730b 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -730,15 +730,6 @@ sections. *Default*: 10 K ---------------------- -```` Element ---------------------- - -The ```` element indicates the number of OpenMP threads to be used for -a simulation. It has no attributes and accepts a positive integer value. - - *Default*: None (Determined by environment variable :envvar:`OMP_NUM_THREADS`) - .. _trace: ------------------- diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 3eb5b1f054..8d43108f14 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -44,12 +44,8 @@ extern std::vector work_index; // Threadprivate variables extern "C" bool trace; //!< flag to show debug information -#ifdef _OPENMP -extern "C" int n_threads; //!< number of OpenMP threads -extern "C" int thread_id; //!< ID of a given thread -#endif -#pragma omp threadprivate(current_work, thread_id, trace) +#pragma omp threadprivate(current_work, trace) } // namespace simulation diff --git a/openmc/settings.py b/openmc/settings.py index 4be2eaaa31..f348bc6cb0 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -130,8 +130,6 @@ class Settings(object): range. 'multipole' is a boolean indicating whether or not the windowed multipole method should be used to evaluate resolved resonance cross sections. - threads : int - Number of OpenMP threads trace : tuple or list Show detailed information about a single particle, indicated by three integers: the batch number, generation number, and particle number @@ -197,7 +195,6 @@ class Settings(object): self._statepoint = {} self._sourcepoint = {} - self._threads = None self._no_reduce = None self._verbosity = None @@ -312,10 +309,6 @@ class Settings(object): def statepoint(self): return self._statepoint - @property - def threads(self): - return self._threads - @property def no_reduce(self): return self._no_reduce @@ -623,12 +616,6 @@ class Settings(object): self._temperature = temperature - @threads.setter - def threads(self, threads): - cv.check_type('number of threads', threads, Integral) - cv.check_greater_than('number of threads', threads, 0) - self._threads = threads - @trace.setter def trace(self, trace): cv.check_type('trace', trace, Iterable, Integral) @@ -885,11 +872,6 @@ class Settings(object): else: element.text = str(value) - def _create_threads_subelement(self, root): - if self._threads is not None: - element = ET.SubElement(root, "threads") - element.text = str(self._threads) - def _create_trace_subelement(self, root): if self._trace is not None: element = ET.SubElement(root, "trace") @@ -980,7 +962,6 @@ class Settings(object): self._create_entropy_mesh_subelement(root_element) self._create_trigger_subelement(root_element) self._create_no_reduce_subelement(root_element) - self._create_threads_subelement(root_element) self._create_verbosity_subelement(root_element) self._create_tabular_legendre_subelements(root_element) self._create_temperature_subelements(root_element) diff --git a/src/bank.cpp b/src/bank.cpp index fc48a77615..a80cfbf375 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -31,7 +31,7 @@ std::vector master_fission_bank; void free_memory_bank() { simulation::source_bank.clear(); -#pragma omp parallel + #pragma omp parallel { simulation::fission_bank.clear(); } diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 1701955295..f5ed493fb8 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -20,6 +20,10 @@ #include "openmc/timer.h" #include "openmc/tallies/tally.h" +#ifdef _OPENMP +#include +#endif + #include // for min #include #include // for sqrt, abs, pow @@ -348,11 +352,13 @@ void calculate_average_keff() #ifdef _OPENMP void join_bank_from_threads() { + int n_threads = omp_get_max_threads(); + #pragma omp parallel { // Copy thread fission bank sites to one shared copy #pragma omp for ordered schedule(static) - for (int i = 0; i < simulation::n_threads; ++i) { + for (int i = 0; i < n_threads; ++i) { #pragma omp ordered { std::copy( @@ -367,7 +373,7 @@ void join_bank_from_threads() #pragma omp barrier // Now copy the shared fission bank sites back to the master thread's copy. - if (simulation::thread_id == 0) { + if (omp_get_thread_num() == 0) { simulation::fission_bank = simulation::master_fission_bank; simulation::master_fission_bank.clear(); } else { diff --git a/src/initialize.cpp b/src/initialize.cpp index 3aa31184f3..61390651f5 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -196,13 +196,13 @@ parse_command_line(int argc, char* argv[]) #ifdef _OPENMP // Read and set number of OpenMP threads - simulation::n_threads = std::stoi(argv[i]); - if (simulation::n_threads < 1) { + int n_threads = std::stoi(argv[i]); + if (n_threads < 1) { std::string msg {"Number of threads must be positive."}; strcpy(openmc_err_msg, msg.c_str()); return OPENMC_E_INVALID_ARGUMENT; } - omp_set_num_threads(simulation::n_threads); + omp_set_num_threads(n_threads); #else if (mpi::master) warning("Ignoring number of threads specified on command line."); diff --git a/src/output.cpp b/src/output.cpp index 92ad14e9d1..82aaca9891 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -11,7 +11,9 @@ #include #include // for pair +#ifdef _OPENMP #include +#endif #include "xtensor/xview.hpp" #include "openmc/capi.h" diff --git a/src/particle.cpp b/src/particle.cpp index 92045e2cdb..ab0f6effba 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -590,7 +590,7 @@ Particle::mark_as_lost(const char* message) // Increment number of lost particles alive_ = false; -#pragma omp atomic + #pragma omp atomic simulation::n_lost_particles += 1; // Count the total number of simulated particles (on this processor) @@ -615,7 +615,7 @@ Particle::write_restart() const filename << settings::path_output << "particle_" << simulation::current_batch << '_' << id_ << ".h5"; -#pragma omp critical (WriteParticleRestart) + #pragma omp critical (WriteParticleRestart) { // Create file hid_t file_id = file_open(filename.str(), 'w'); diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 2a6b700f03..84dfa9cdea 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -174,22 +174,19 @@ absorption(Particle* p) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->wgt_absorb_ = p->wgt_ * - p->macro_xs_.absorption / p->macro_xs_.total; + p->wgt_absorb_ = p->wgt_ * p->macro_xs_.absorption / p->macro_xs_.total; // Adjust weight of particle by the probability of absorption p->wgt_ -= p->wgt_absorb_; p->wgt_last_ = p->wgt_; // Score implicit absorpion estimate of keff -#pragma omp atomic - global_tally_absorption += p->wgt_absorb_ * - p->macro_xs_.nu_fission / - p->macro_xs_.absorption; + #pragma omp atomic + global_tally_absorption += p->wgt_absorb_ * p->macro_xs_.nu_fission / + p->macro_xs_.absorption; } else { - if (p->macro_xs_.absorption > - prn() * p->macro_xs_.total) { -#pragma omp atomic + if (p->macro_xs_.absorption > prn() * p->macro_xs_.total) { + #pragma omp atomic global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission / p->macro_xs_.absorption; p->alive_ = false; diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index cf5d7e34bd..5002b42794 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -141,7 +141,7 @@ extern "C" void openmc_set_seed(int64_t new_seed) { seed = new_seed; -#pragma omp parallel + #pragma omp parallel { for (int i = 0; i < N_STREAMS; i++) { prn_seed[i] = seed + i; diff --git a/src/settings.cpp b/src/settings.cpp index dfa5ebb7e5..4c041a080d 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -383,20 +383,8 @@ void read_settings_xml() // Number of OpenMP threads if (check_for_node(root, "threads")) { -#ifdef _OPENMP - if (simulation::n_threads == 0) { - simulation::n_threads = std::stoi(get_node_value(root, "threads")); - if (simulation::n_threads < 1) { - std::stringstream msg; - msg << "Invalid number of threads: " << simulation::n_threads; - fatal_error(msg); - } - omp_set_num_threads(simulation::n_threads); - } -#else - if (mpi::master) warning("OpenMC was not compiled with OpenMP support; " - "ignoring number of threads."); -#endif + if (mpi::master) warning("The element has been deprecated. Use " + "the OMP_NUM_THREADS environment variable to set the number of threads."); } // ========================================================================== diff --git a/src/simulation.cpp b/src/simulation.cpp index fca0c85282..ea313e94c0 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -20,7 +20,9 @@ #include "openmc/tallies/tally.h" #include "openmc/tallies/trigger.h" +#ifdef _OPENMP #include +#endif #include "xtensor/xview.hpp" #include @@ -59,7 +61,7 @@ int openmc_simulation_init() calculate_work(); // Allocate array for matching filter bins -#pragma omp parallel + #pragma omp parallel { simulation::filter_matches.resize(model::tally_filters.size()); } @@ -136,7 +138,7 @@ int openmc_simulation_finalize() // Write tally results to tallies.out if (settings::output_tallies && mpi::master) write_tallies(); -#pragma omp parallel + #pragma omp parallel { simulation::filter_matches.clear(); } @@ -186,7 +188,7 @@ int openmc_next_batch(int* status) // ==================================================================== // LOOP OVER PARTICLES -#pragma omp parallel for schedule(runtime) + #pragma omp parallel for schedule(runtime) for (int64_t i_work = 1; i_work <= simulation::work; ++i_work) { simulation::current_work = i_work; @@ -250,10 +252,6 @@ std::vector work_index; // Threadprivate variables bool trace; //!< flag to show debug information -#ifdef _OPENMP -int n_threads {-1}; //!< number of OpenMP threads -int thread_id; //!< ID of a given thread -#endif } // namespace simulation @@ -273,15 +271,13 @@ void allocate_banks() // a generation, there is also a 'master_fission_bank' that is used to // collect the sites from each thread. - simulation::n_threads = omp_get_max_threads(); - -#pragma omp parallel + #pragma omp parallel { - simulation::thread_id = omp_get_thread_num(); - if (simulation::thread_id == 0) { + if (omp_get_thread_num() == 0) { simulation::fission_bank.reserve(3*simulation::work); } else { - simulation::fission_bank.reserve(3*simulation::work / simulation::n_threads); + int n_threads = omp_get_num_threads(); + simulation::fission_bank.reserve(3*simulation::work / n_threads); } } simulation::master_fission_bank.reserve(3*simulation::work); @@ -403,9 +399,9 @@ void finalize_generation() auto& gt = simulation::global_tallies; // Update global tallies with the omp private accumulation variables -#pragma omp parallel + #pragma omp parallel { -#pragma omp critical(increment_global_tallies) + #pragma omp critical(increment_global_tallies) { if (settings::run_mode == RUN_MODE_EIGENVALUE) { gt(K_COLLISION, RESULT_VALUE) += global_tally_collision; diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index 697eb095fa..e42f6240f5 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -36,7 +36,6 @@ def test_export_to_xml(run_in_tmpdir): s.tabular_legendre = {'enable': True, 'num_points': 50} s.temperature = {'default': 293.6, 'method': 'interpolation', 'multipole': True, 'range': (200., 1000.)} - s.threads = 8 s.trace = (10, 1, 20) s.track = [1, 1, 1, 2, 1, 1] s.ufs_mesh = mesh From 279104cfed48518373efb59186758e96b4acb84e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Mar 2019 10:55:05 -0500 Subject: [PATCH 08/14] Move secondary bank outside of Particle class --- include/openmc/bank.h | 3 ++- include/openmc/particle.h | 5 +---- src/bank.cpp | 1 + src/particle.cpp | 12 ++++++------ src/physics.cpp | 2 +- src/physics_mg.cpp | 2 +- 6 files changed, 12 insertions(+), 13 deletions(-) diff --git a/include/openmc/bank.h b/include/openmc/bank.h index 4597369c78..93e2be669f 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -22,11 +22,12 @@ namespace simulation { extern std::vector source_bank; extern std::vector fission_bank; +extern std::vector secondary_bank; #ifdef _OPENMP extern std::vector master_fission_bank; #endif -#pragma omp threadprivate(fission_bank) +#pragma omp threadprivate(fission_bank, secondary_bank) } // namespace simulation diff --git a/include/openmc/particle.h b/include/openmc/particle.h index d52c67ea55..c6156e481d 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -187,7 +187,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(Direction u, double E, Type type); + void create_secondary(Direction u, double E, Type type) const; //! initialize from a source site // @@ -283,9 +283,6 @@ public: // Track output bool write_track_ {false}; - - // Secondary particles created - std::vector secondary_bank_; }; } // namespace openmc diff --git a/src/bank.cpp b/src/bank.cpp index a80cfbf375..172ecf84f3 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -18,6 +18,7 @@ namespace simulation { std::vector source_bank; std::vector fission_bank; +std::vector secondary_bank; #ifdef _OPENMP std::vector master_fission_bank; #endif diff --git a/src/particle.cpp b/src/particle.cpp index ab0f6effba..8ca07ca82a 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -73,11 +73,11 @@ Particle::clear() } void -Particle::create_secondary(Direction u, double E, Type type) +Particle::create_secondary(Direction u, double E, Type type) const { - secondary_bank_.emplace_back(); + simulation::secondary_bank.emplace_back(); - auto& bank {secondary_bank_.back()}; + auto& bank {simulation::secondary_bank.back()}; bank.particle = type; bank.wgt = wgt_; bank.r = this->r(); @@ -357,10 +357,10 @@ Particle::transport() // Check for secondary particles if this particle is dead if (!alive_) { // If no secondary particles, break out of event loop - if (secondary_bank_.empty()) break; + if (simulation::secondary_bank.empty()) break; - this->from_source(&secondary_bank_.back()); - secondary_bank_.pop_back(); + this->from_source(&simulation::secondary_bank.back()); + simulation::secondary_bank.pop_back(); n_event = 0; // Enter new particle in particle track file diff --git a/src/physics.cpp b/src/physics.cpp index 7d44aa17c6..2d48855c77 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -94,7 +94,7 @@ void sample_neutron_reaction(Particle* p) create_fission_sites(p, i_nuclide, rx, simulation::fission_bank); } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE && settings::create_fission_neutrons) { - create_fission_sites(p, i_nuclide, rx, p->secondary_bank_); + create_fission_sites(p, i_nuclide, rx, simulation::secondary_bank); } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 84dfa9cdea..f99d15aac7 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -51,7 +51,7 @@ sample_reaction(Particle* p) create_fission_sites(p, simulation::fission_bank); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { - create_fission_sites(p, p->secondary_bank_); + create_fission_sites(p, simulation::secondary_bank); } } From c049557c8575bf76dbf0346633e7707fad76e94b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Mar 2019 14:38:24 -0500 Subject: [PATCH 09/14] Get rid of MAX_SECONDARY constant --- include/openmc/particle.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index c6156e481d..8ec53de447 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -26,9 +26,6 @@ namespace openmc { // use to store the bins for delayed group tallies. constexpr int MAX_DELAYED_GROUPS {8}; -// Maximum number of secondary particles created -constexpr int MAX_SECONDARY {1000}; - // Maximum number of lost particles constexpr int MAX_LOST_PARTICLES {10}; From 492d984525a665f97baf02326a8b407d884883fc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Mar 2019 13:26:04 -0500 Subject: [PATCH 10/14] Change name of work -> work_per_rank --- include/openmc/simulation.h | 2 +- src/eigenvalue.cpp | 4 ++-- src/particle.cpp | 3 ++- src/simulation.cpp | 16 ++++++++-------- src/source.cpp | 4 ++-- src/state_point.cpp | 8 ++++---- 6 files changed, 19 insertions(+), 18 deletions(-) diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 8d43108f14..fef6ebbcd1 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -37,7 +37,7 @@ extern "C" int restart_batch; //!< batch at which a restart job resumed extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? extern "C" int total_gen; //!< total number of generations simulated extern double total_weight; //!< Total source weight in a batch -extern "C" int64_t work; //!< number of particles per process +extern int64_t work_per_rank; //!< number of particles per MPI rank extern std::vector k_generation; extern std::vector work_index; diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index f5ed493fb8..08a4d3f91c 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -138,7 +138,7 @@ void synchronize_bank() // Allocate temporary source bank int64_t index_temp = 0; - std::vector temp_sites(3*simulation::work); + std::vector temp_sites(3*simulation::work_per_rank); for (const auto& site : simulation::fission_bank) { // If there are less than n_particles particles banked, automatically add @@ -597,7 +597,7 @@ void ufs_count_sites() // Since the total starting weight is not equal to n_particles, we need to // renormalize the weight of the source sites - for (int i = 0; i < simulation::work; ++i) { + for (int i = 0; i < simulation::work_per_rank; ++i) { simulation::source_bank[i].wgt *= settings::n_particles / total; } } diff --git a/src/particle.cpp b/src/particle.cpp index 8ca07ca82a..d2b845bf41 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -594,7 +594,8 @@ Particle::mark_as_lost(const char* message) simulation::n_lost_particles += 1; // Count the total number of simulated particles (on this processor) - auto n = simulation::current_batch * settings::gen_per_batch * simulation::work; + auto n = simulation::current_batch * settings::gen_per_batch * + simulation::work_per_rank; // Abort the simulation if the maximum number of lost particles has been // reached diff --git a/src/simulation.cpp b/src/simulation.cpp index ea313e94c0..fc7a788517 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -189,7 +189,7 @@ int openmc_next_batch(int* status) // LOOP OVER PARTICLES #pragma omp parallel for schedule(runtime) - for (int64_t i_work = 1; i_work <= simulation::work; ++i_work) { + for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) { simulation::current_work = i_work; // grab source particle from bank @@ -245,7 +245,7 @@ int restart_batch; bool satisfy_triggers {false}; int total_gen {0}; double total_weight; -int64_t work; +int64_t work_per_rank; std::vector k_generation; std::vector work_index; @@ -262,7 +262,7 @@ bool trace; //!< flag to show debug information void allocate_banks() { // Allocate source bank - simulation::source_bank.resize(simulation::work); + simulation::source_bank.resize(simulation::work_per_rank); if (settings::run_mode == RUN_MODE_EIGENVALUE) { #ifdef _OPENMP @@ -274,15 +274,15 @@ void allocate_banks() #pragma omp parallel { if (omp_get_thread_num() == 0) { - simulation::fission_bank.reserve(3*simulation::work); + simulation::fission_bank.reserve(3*simulation::work_per_rank); } else { int n_threads = omp_get_num_threads(); - simulation::fission_bank.reserve(3*simulation::work / n_threads); + simulation::fission_bank.reserve(3*simulation::work_per_rank / n_threads); } } - simulation::master_fission_bank.reserve(3*simulation::work); + simulation::master_fission_bank.reserve(3*simulation::work_per_rank); #else - simulation::fission_bank.reserve(3*simulation::work); + simulation::fission_bank.reserve(3*simulation::work_per_rank); #endif } } @@ -507,7 +507,7 @@ void calculate_work() int64_t work_i = i < remainder ? min_work + 1 : min_work; // Set number of particles - if (mpi::rank == i) simulation::work = work_i; + if (mpi::rank == i) simulation::work_per_rank = work_i; // Set index into source bank for rank i i_bank += work_i; diff --git a/src/source.cpp b/src/source.cpp index 36d66394ca..0baa21b170 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -264,7 +264,7 @@ void initialize_source() } else { // Generation source sites from specified distribution in user input - for (int64_t i = 0; i < simulation::work; ++i) { + for (int64_t i = 0; i < simulation::work_per_rank; ++i) { // initialize random number seed int64_t id = simulation::total_gen*settings::n_particles + simulation::work_index[mpi::rank] + i + 1; @@ -330,7 +330,7 @@ void free_memory_source() void fill_source_bank_fixedsource() { if (settings::path_source.empty()) { - for (int64_t i = 0; i < simulation::work; ++i) { + for (int64_t i = 0; i < simulation::work_per_rank; ++i) { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation()) * settings::n_particles + simulation::work_index[mpi::rank] + i + 1; diff --git a/src/state_point.cpp b/src/state_point.cpp index d657316f97..54f1cb38c5 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -537,7 +537,7 @@ write_source_bank(hid_t group_id) H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); // Create another data space but for each proc individually - hsize_t count[] {static_cast(simulation::work)}; + hsize_t count[] {static_cast(simulation::work_per_rank)}; hid_t memspace = H5Screate_simple(1, count, nullptr); // Select hyperslab for this dataspace @@ -569,7 +569,7 @@ write_source_bank(hid_t group_id) // Save source bank sites since the souce_bank array is overwritten below #ifdef OPENMC_MPI std::vector temp_source {simulation::source_bank.begin(), - simulation::source_bank.begin() + simulation::work}; + simulation::source_bank.begin() + simulation::work_per_rank}; #endif for (int i = 0; i < mpi::n_procs; ++i) { @@ -607,7 +607,7 @@ write_source_bank(hid_t group_id) #endif } else { #ifdef OPENMC_MPI - MPI_Send(simulation::source_bank.data(), simulation::work, mpi::bank, + MPI_Send(simulation::source_bank.data(), simulation::work_per_rank, mpi::bank, 0, mpi::rank, mpi::intracomm); #endif } @@ -625,7 +625,7 @@ void read_source_bank(hid_t group_id) hid_t dset = H5Dopen(group_id, "source_bank", H5P_DEFAULT); // Create another data space but for each proc individually - hsize_t dims[] {static_cast(simulation::work)}; + hsize_t dims[] {static_cast(simulation::work_per_rank)}; hid_t memspace = H5Screate_simple(1, dims, nullptr); // Make sure source bank is big enough From 0fbe024e11c373e985fd5220353e57bdbca8c2dc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Mar 2019 14:47:42 -0500 Subject: [PATCH 11/14] Get rid of ERROR_INT and use 0 for "no surface" --- include/openmc/constants.h | 2 -- src/output.cpp | 2 +- src/particle.cpp | 6 +++--- 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 87388c91a5..d9307abe20 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -403,8 +403,6 @@ constexpr int LEAKAGE {3}; // Miscellaneous constexpr int C_NONE {-1}; constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE -constexpr int ERROR_INT {-2147483647}; // TODO: use when F90 - // interop is gone // Interpolation rules enum class Interpolation { diff --git a/src/output.cpp b/src/output.cpp index 82aaca9891..df0976d94a 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -191,7 +191,7 @@ extern "C" void print_particle(Particle* p) } // Display miscellaneous info. - if (p->surface_ != ERROR_INT) { + if (p->surface_ != 0) { const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]}; std::cout << " Surface = " << std::copysign(surf.id_, p->surface_) << "\n"; } diff --git a/src/particle.cpp b/src/particle.cpp index d2b845bf41..b307aec8d1 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -259,7 +259,7 @@ Particle::transport() if (lattice_translation[0] != 0 || lattice_translation[1] != 0 || lattice_translation[2] != 0) { // Particle crosses lattice boundary - surface_ = ERROR_INT; + surface_ = 0; cross_lattice(this, lattice_translation); event_ = EVENT_LATTICE; } else { @@ -291,7 +291,7 @@ Particle::transport() score_surface_tally(this, model::active_meshsurf_tallies); // Clear surface component - surface_ = ERROR_INT; + surface_ = 0; if (settings::run_CE) { collision(this); @@ -556,7 +556,7 @@ Particle::cross_surface() // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS // Remove lower coordinate levels and assignment of surface - surface_ = ERROR_INT; + surface_ = 0; n_coord_ = 1; bool found = find_cell(this, false); From 4971779f39ab84ae5826e0835ce9dea56a53a0ee Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Mar 2019 16:07:00 -0500 Subject: [PATCH 12/14] Use __builtin_reachable to tell gcc about unreachable code --- src/nuclide.cpp | 3 +++ src/physics.cpp | 5 +++-- src/tallies/filter_particle.cpp | 3 +++ 3 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index cfb5a669f2..8a6f1aa57b 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -454,6 +454,9 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const return (*fission_rx_[0]->products_[0].yield_)(E); } } +#ifdef __GNUC__ + __builtin_unreachable(); +#endif } void Nuclide::calculate_elastic_xs(Particle& p) const diff --git a/src/physics.cpp b/src/physics.cpp index 2d48855c77..3edaf4545e 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -882,8 +882,9 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, } // case RVS, DBRC } // switch (sampling_method) - throw std::runtime_error{"Unable to sample target velocity for " - "elastic scattering."}; +#ifdef __GNUC__ + __builtin_unreachable(); +#endif } Direction diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index a2902785ae..1d4b849375 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -50,6 +50,9 @@ ParticleFilter::text_label(int bin) const case Particle::Type::positron: return "Particle: positron"; } +#ifdef __GNUC__ + __builtin_unreachable(); +#endif } } // namespace openmc From 2e94189b1bcfa60c70fd6f06912c0757c72f6cdd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Mar 2019 22:33:54 -0500 Subject: [PATCH 13/14] Make sure subcrit multiplication problems abort gracefully when k > 1 --- src/physics.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/physics.cpp b/src/physics.cpp index 3edaf4545e..1cfd316580 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -95,6 +95,14 @@ void sample_neutron_reaction(Particle* p) } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE && settings::create_fission_neutrons) { create_fission_sites(p, i_nuclide, rx, simulation::secondary_bank); + + // Make sure particle population doesn't grow out of control for + // subcritical multiplication problems. + if (simulation::secondary_bank.size() >= 10000) { + fatal_error("The secondary particle bank appears to be growing without " + "bound. You are likely running a subcritical multiplication problem " + "with k-effective close to or greater than one."); + } } } From c60ecd7e8f688927904e818535c2d30d420b76e0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 21 Mar 2019 07:43:42 -0500 Subject: [PATCH 14/14] Have distance_to_boundary return a struct --- include/openmc/geometry.h | 26 +++++++++++------- src/geometry.cpp | 56 +++++++++++++++++---------------------- src/particle.cpp | 24 +++++++---------- 3 files changed, 52 insertions(+), 54 deletions(-) diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 9bb1a9ae7f..5feca2992a 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -1,6 +1,8 @@ #ifndef OPENMC_GEOMETRY_H #define OPENMC_GEOMETRY_H +#include +#include #include #include @@ -22,12 +24,22 @@ extern std::vector overlap_check_count; } // namespace model +//============================================================================== +// 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 +}; + //============================================================================== //! Check for overlapping cells at a particle's position. //============================================================================== -extern "C" bool -check_cell_overlap(Particle* p); +bool check_cell_overlap(Particle* p); //============================================================================== //! Locate a particle in the geometry tree and set its geometry data fields. @@ -41,23 +53,19 @@ check_cell_overlap(Particle* p); //! valid geometry coordinate stack. //============================================================================== -extern "C" bool -find_cell(Particle* p, bool use_neighbor_lists); +bool find_cell(Particle* p, bool use_neighbor_lists); //============================================================================== //! Move a particle into a new lattice tile. //============================================================================== -extern "C" void -cross_lattice(Particle* p, int lattice_translation[3]); +void cross_lattice(Particle* p, const BoundaryInfo& boundary); //============================================================================== //! Find the next boundary a particle will intersect. //============================================================================== -extern "C" void -distance_to_boundary(Particle* p, double* dist, int* surface_crossed, - int lattice_translation[3], int* next_level); +BoundaryInfo distance_to_boundary(Particle* p); } // namespace openmc diff --git a/src/geometry.cpp b/src/geometry.cpp index cde2cdc1b6..6eaa61fe24 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -32,8 +32,7 @@ std::vector overlap_check_count; // Non-member functions //============================================================================== -extern "C" bool -check_cell_overlap(Particle* p) +bool check_cell_overlap(Particle* p) { int n_coord = p->n_coord_; @@ -246,7 +245,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) //============================================================================== -extern "C" bool +bool find_cell(Particle* p, bool use_neighbor_lists) { // Determine universe (if not yet set, use root universe). @@ -288,8 +287,8 @@ find_cell(Particle* p, bool use_neighbor_lists) //============================================================================== -extern "C" void -cross_lattice(Particle* p, int lattice_translation[3]) +void +cross_lattice(Particle* p, const BoundaryInfo& boundary) { auto& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice]}; @@ -303,9 +302,9 @@ cross_lattice(Particle* p, int lattice_translation[3]) } // Set the lattice indices. - p->coord_[p->n_coord_-1].lattice_x += lattice_translation[0]; - p->coord_[p->n_coord_-1].lattice_y += lattice_translation[1]; - p->coord_[p->n_coord_-1].lattice_z += lattice_translation[2]; + p->coord_[p->n_coord_-1].lattice_x += boundary.lattice_translation[0]; + p->coord_[p->n_coord_-1].lattice_y += boundary.lattice_translation[1]; + p->coord_[p->n_coord_-1].lattice_z += boundary.lattice_translation[2]; std::array i_xyz {p->coord_[p->n_coord_-1].lattice_x, p->coord_[p->n_coord_-1].lattice_y, p->coord_[p->n_coord_-1].lattice_z}; @@ -346,16 +345,11 @@ cross_lattice(Particle* p, int lattice_translation[3]) //============================================================================== -extern "C" void -distance_to_boundary(Particle* p, double* dist, int* surface_crossed, - int lattice_translation[3], int* next_level) +BoundaryInfo distance_to_boundary(Particle* p) { - *dist = INFINITY; + BoundaryInfo info; double d_lat = INFINITY; double d_surf = INFINITY; - lattice_translation[0] = 0; - lattice_translation[1] = 0; - lattice_translation[2] = 0; int32_t level_surf_cross; std::array level_lat_trans; @@ -402,43 +396,43 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, // If the boundary on this coordinate level is coincident with a boundary on // a higher level then we need to make sure that the higher level boundary // is selected. This logic must consider floating point precision. + double& d = info.distance; if (d_surf < d_lat) { - if (*dist == INFINITY || ((*dist) - d_surf)/(*dist) >= FP_REL_PRECISION) { - *dist = d_surf; + if (d == INFINITY || (d - d_surf)/d >= FP_REL_PRECISION) { + d = d_surf; // If the cell is not simple, it is possible that both the negative and // positive half-space were given in the region specification. Thus, we // have to explicitly check which half-space the particle would be // traveling into if the surface is crossed if (c.simple_) { - *surface_crossed = level_surf_cross; + info.surface_index = level_surf_cross; } else { Position r_hit = r + d_surf * u; Surface& surf {*model::surfaces[std::abs(level_surf_cross)-1]}; Direction norm = surf.normal(r_hit); if (u.dot(norm) > 0) { - *surface_crossed = std::abs(level_surf_cross); + info.surface_index = std::abs(level_surf_cross); } else { - *surface_crossed = -std::abs(level_surf_cross); + info.surface_index = -std::abs(level_surf_cross); } } - lattice_translation[0] = 0; - lattice_translation[1] = 0; - lattice_translation[2] = 0; - *next_level = i + 1; + info.lattice_translation[0] = 0; + info.lattice_translation[1] = 0; + info.lattice_translation[2] = 0; + info.coord_level = i + 1; } } else { - if (*dist == INFINITY || ((*dist) - d_lat)/(*dist) >= FP_REL_PRECISION) { - *dist = d_lat; - *surface_crossed = F90_NONE; - lattice_translation[0] = level_lat_trans[0]; - lattice_translation[1] = level_lat_trans[1]; - lattice_translation[2] = level_lat_trans[2]; - *next_level = i + 1; + if (d == INFINITY || (d - d_lat)/d >= FP_REL_PRECISION) { + d = d_lat; + info.surface_index = 0; + info.lattice_translation = level_lat_trans; + info.coord_level = i + 1; } } } + return info; } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index b307aec8d1..4e0066b20e 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -202,12 +202,7 @@ Particle::transport() } // Find the distance to the nearest boundary - double d_boundary; - int surface_crossed; - int lattice_translation[3]; - int next_level; - distance_to_boundary(this, &d_boundary, &surface_crossed, - lattice_translation, &next_level); + auto boundary = distance_to_boundary(this); // Sample a distance to collision double d_collision; @@ -221,7 +216,7 @@ Particle::transport() } // Select smaller of the two distances - double distance = std::min(d_boundary, d_collision); + double distance = std::min(boundary.distance, d_collision); // Advance particle for (int j = 0; j < n_coord_; ++j) { @@ -244,11 +239,13 @@ Particle::transport() score_track_derivative(this, distance); } - if (d_collision > d_boundary) { + if (d_collision > boundary.distance) { // ==================================================================== // PARTICLE CROSSES SURFACE - if (next_level > 0) n_coord_ = next_level; + // Set surface that particle is on and adjust coordinate levels + surface_ = boundary.surface_index; + n_coord_ = boundary.coord_level; // Saving previous cell data for (int j = 0; j < n_coord_; ++j) { @@ -256,15 +253,14 @@ Particle::transport() } n_coord_last_ = n_coord_; - if (lattice_translation[0] != 0 || lattice_translation[1] != 0 || - lattice_translation[2] != 0) { + if (boundary.lattice_translation[0] != 0 || + boundary.lattice_translation[1] != 0 || + boundary.lattice_translation[2] != 0) { // Particle crosses lattice boundary - surface_ = 0; - cross_lattice(this, lattice_translation); + cross_lattice(this, boundary); event_ = EVENT_LATTICE; } else { // Particle crosses surface - surface_ = surface_crossed; this->cross_surface(); event_ = EVENT_SURFACE; }