crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Move cross section caches into Particle class

Browse source
This commit is contained in:
Paul Romano 2019-03-18 09:12:37 -05:00
parent 522b6be8eb
commit 8a8b7b2c0a
17 changed files with 495 additions and 535 deletions
Download patch file Download diff file Expand all files Collapse all files

6
include/openmc/material.h
View file

@ -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;
};
//==============================================================================

82
include/openmc/nuclide.h
View file

@ -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<std::string, int> 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
//==============================================================================

202
include/openmc/particle.h
View file

@ -6,9 +6,11 @@
#include <array>
#include <cstdint>
#include <memory> // for unique_ptr
#include <sstream>
#include <string>
#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<NuclideMicroXS> micro_xs_;
std::vector<ElementMicroXS> 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

23
include/openmc/photon.h
View file

@ -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<std::string, int> 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

11
include/openmc/physics.h
View file

@ -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

6
include/openmc/tallies/tally_scoring.h
View file

@ -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.
//

2
include/openmc/thermal.h
View file

@ -10,7 +10,7 @@
#include "xtensor/xtensor.hpp"
#include "openmc/hdf5_interface.h"
#include "openmc/nuclide.h"
#include "openmc/particle.h"
namespace openmc {

49
src/material.cpp
View file

@ -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<int>(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;
}
}

124
src/nuclide.cpp
View file

@ -34,11 +34,6 @@ std::vector<std::unique_ptr<Nuclide>> nuclides;
std::unordered_map<std::string, int> 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;
}
}

30
src/particle.cpp
View file

@ -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

14
src/particle_restart.cpp
View file

@ -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

12
src/photon.cpp
View file

@ -32,10 +32,6 @@ std::unordered_map<std::string, int> 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

111
src/physics.cpp
View file

@ -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<int>(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<UncorrelatedAngleEnergy*>(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<int>(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];

18
src/physics_mg.cpp
View file

@ -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<int>(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;
}

14
src/simulation.cpp
View file

@ -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;

142
src/tallies/derivative.cpp
View file

@ -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_);

184
src/tallies/tally_scoring.cpp
View file

@ -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) {