From 8a8b7b2c0ac9d6c0896cfae560a65476e8cd4df9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Mar 2019 09:12:37 -0500 Subject: [PATCH] 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 932945ffc..3109f6d78 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 c8b34d08e..4be898a64 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 2cb30a44f..e971c30da 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 7c20a716d..902d24232 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 a7dd20aad..5c8d64f5f 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 af03b1aa2..2c5a29c17 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 624cf31bd..b737f005b 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 da68a536d..0d2fe4bbc 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 cd2461186..87965dfc3 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 ca7fccfe7..f5084cd1c 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 0cface319..0bc677cfd 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 d977d62ac..9841e9e3b 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 858aae8f4..029ab18c8 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 fa6365a3a..9565a627c 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 7d30ccea9..b4db5e0e6 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 315cc081e..411995fd7 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 ed4d09254..57e810dfd 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) {