crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1202 from paulromano/particle-improvements

Browse source 
Improvements to Particle class
This commit is contained in:
Sterling Harper 2019-03-23 20:29:51 -04:00 • committed by GitHub
commit 99e5d0ba12
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
43 changed files with 749 additions and 897 deletions
Download patch file Download diff file Expand all files Collapse all files

4
CMakeLists.txt
View file

@ -20,9 +20,6 @@ option(optimize "Turn on all compiler optimization flags" OFF)
option(coverage "Compile with coverage analysis flags" OFF)
option(dagmc "Enable support for DAGMC (CAD) geometry" OFF)
# Maximum number of nested coordinates levels
set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels")
#===============================================================================
# MPI for distributed-memory parallelism
#===============================================================================
@ -260,7 +257,6 @@ target_include_directories(libopenmc
# Set compile flags
target_compile_options(libopenmc PRIVATE ${cxxflags})
target_compile_definitions(libopenmc PUBLIC -DMAX_COORD=${maxcoord})
if (HDF5_IS_PARALLEL)
target_compile_definitions(libopenmc PRIVATE -DPHDF5)
endif()

25
docs/source/devguide/styleguide.rst
View file

@ -12,6 +12,11 @@ adding new code in OpenMC.
C++
---
Indentation
-----------
Use two spaces per indentation level.
Miscellaneous
-------------
@ -126,6 +131,15 @@ single declaration to avoid confusion:
Curly braces
------------
For a class declaration, the opening brace should be on the same line that
lists the name of the class.
.. code-block:: C++
class Matrix {
...
};
For a function definition, the opening and closing braces should each be on
their own lines. This helps distinguish function code from the argument list.
If the entire function fits on one or two lines, then the braces can be on the
@ -210,11 +224,18 @@ Use of third-party Python packages should be limited to numpy_, scipy_,
matplotlib_, pandas_, and h5py_. Use of other third-party packages must be
implemented as optional dependencies rather than required dependencies.
Prefer pathlib_ when working with filesystem paths over functions in the os_
module or other standard-library modules. Functions that accept arguments that
represent a filesystem path should work with both strings and Path_ objects.
.. _C++ Core Guidelines: http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines
.. _PEP8: https://www.python.org/dev/peps/pep-0008/
.. _numpydoc: https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt
.. _numpydoc: https://numpydoc.readthedocs.io/en/latest/format.html
.. _numpy: http://www.numpy.org/
.. _scipy: https://www.scipy.org/
.. _matplotlib: https://matplotlib.org/
.. _pandas: https://pandas.pydata.org/
.. _h5py: http://www.h5py.org/
.. _h5py: https://www.h5py.org/
.. _pathlib: https://docs.python.org/3/library/pathlib.html
.. _os: https://docs.python.org/3/library/os.html
.. _Path: https://docs.python.org/3/library/pathlib.html#pathlib.Path

9
docs/source/io_formats/settings.rst
View file

@ -730,15 +730,6 @@ sections.
*Default*: 10 K
---------------------
``<threads>`` Element
---------------------
The ``<threads>`` element indicates the number of OpenMP threads to be used for
a simulation. It has no attributes and accepts a positive integer value.
*Default*: None (Determined by environment variable :envvar:`OMP_NUM_THREADS`)
.. _trace:
-------------------

2
docs/source/methods/photon_physics.rst
View file

@ -5,7 +5,7 @@ Photon Physics
==============
Photons, being neutral particles, behave much in the same manner as neutrons,
traveling in straight lines and experiencing occasional collisions which change
traveling in straight lines and experiencing occasional collisions that change
their energy and direction. Photons undergo four basic interactions as they pass
through matter: coherent (Rayleigh) scattering, incoherent (Compton) scattering,
photoelectric effect, and pair/triplet production. Photons with energy in the

3
docs/source/usersguide/install.rst
View file

@ -243,9 +243,6 @@ coverage
Compile and link code instrumented for coverage analysis. This is typically
used in conjunction with gcov_.
maxcoord
Maximum number of nested coordinate levels in geometry. Defaults to 10.
To set any of these options (e.g. turning on debug mode), the following form
should be used:

5
include/openmc/bank.h
View file

@ -20,15 +20,14 @@ namespace openmc {
namespace simulation {
extern "C" int64_t n_bank;
extern std::vector<Particle::Bank> source_bank;
extern std::vector<Particle::Bank> fission_bank;
extern std::vector<Particle::Bank> secondary_bank;
#ifdef _OPENMP
extern std::vector<Particle::Bank> master_fission_bank;
#endif
#pragma omp threadprivate(fission_bank, n_bank)
#pragma omp threadprivate(fission_bank, secondary_bank)
} // namespace simulation

2
include/openmc/constants.h
View file

@ -403,8 +403,6 @@ constexpr int LEAKAGE {3};
// Miscellaneous
constexpr int C_NONE {-1};
constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
constexpr int ERROR_INT {-2147483647}; // TODO: use <numeric_limits> when F90
// interop is gone
// Interpolation rules
enum class Interpolation {

29
include/openmc/geometry.h
View file

@ -1,6 +1,8 @@
#ifndef OPENMC_GEOMETRY_H
#define OPENMC_GEOMETRY_H
#include <array>
#include <cmath>
#include <cstdint>
#include <vector>
@ -15,18 +17,29 @@ namespace openmc {
namespace model {
extern "C" int root_universe;
extern int root_universe; //!< Index of root universe
extern int n_coord_levels; //!< Number of CSG coordinate levels
extern std::vector<int64_t> overlap_check_count;
} // namespace model
//==============================================================================
// Information about nearest boundary crossing
//==============================================================================
struct BoundaryInfo {
double distance {INFINITY}; //!< distance to nearest boundary
int surface_index {0}; //!< if boundary is surface, index in surfaces vector
int coord_level; //!< coordinate level after crossing boundary
std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
};
//==============================================================================
//! Check for overlapping cells at a particle's position.
//==============================================================================
extern "C" bool
check_cell_overlap(Particle* p);
bool check_cell_overlap(Particle* p);
//==============================================================================
//! Locate a particle in the geometry tree and set its geometry data fields.
@ -40,23 +53,19 @@ check_cell_overlap(Particle* p);
//! valid geometry coordinate stack.
//==============================================================================
extern "C" bool
find_cell(Particle* p, bool use_neighbor_lists);
bool find_cell(Particle* p, bool use_neighbor_lists);
//==============================================================================
//! Move a particle into a new lattice tile.
//==============================================================================
extern "C" void
cross_lattice(Particle* p, int lattice_translation[3]);
void cross_lattice(Particle* p, const BoundaryInfo& boundary);
//==============================================================================
//! Find the next boundary a particle will intersect.
//==============================================================================
extern "C" void
distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
int lattice_translation[3], int* next_level);
BoundaryInfo distance_to_boundary(Particle* p);
} // namespace openmc

9
include/openmc/geometry_aux.h
View file

@ -57,7 +57,7 @@ void finalize_geometry(std::vector<std::vector<double>>& nuc_temps,
//! \return The index of the root universe.
//==============================================================================
extern "C" int32_t find_root_universe();
int32_t find_root_universe();
//==============================================================================
//! Populate all data structures needed for distribcells.
@ -74,7 +74,7 @@ void prepare_distribcell();
//! the root universe).
//==============================================================================
extern "C" void count_cell_instances(int32_t univ_indx);
void count_cell_instances(int32_t univ_indx);
//==============================================================================
//! Recursively search through universes and count universe instances.
@ -84,8 +84,7 @@ extern "C" void count_cell_instances(int32_t univ_indx);
//! search_univ.
//==============================================================================
extern "C" int
count_universe_instances(int32_t search_univ, int32_t target_univ_id);
int count_universe_instances(int32_t search_univ, int32_t target_univ_id);
//==============================================================================
//! Build a character array representing the path to a distribcell instance.
@ -107,7 +106,7 @@ distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset);
//! \return The number of coordinate levels.
//==============================================================================
extern "C" int maximum_levels(int32_t univ);
int maximum_levels(int32_t univ);
//==============================================================================
//! Deallocates global vectors and maps for cells, universes, and lattices.

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

7
include/openmc/mesh.h
View file

@ -96,14 +96,11 @@ public:
//! Count number of bank sites in each mesh bin / energy bin
//
//! \param[in] n Number of bank sites
//! \param[in] bank Array of bank sites
//! \param[in] n_energy Number of energies
//! \param[in] energies Array of energies
//! \param[out] Whether any bank sites are outside the mesh
//! \return Array indicating number of sites in each mesh/energy bin
xt::xarray<double> count_sites(int64_t n, const Particle::Bank* bank,
int n_energy, const double* energies, bool* outside) const;
xt::xarray<double> count_sites(const std::vector<Particle::Bank>& bank,
bool* outside) const;
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions

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
//==============================================================================

214
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 {
@ -24,15 +26,14 @@ namespace openmc {
// use to store the bins for delayed group tallies.
constexpr int MAX_DELAYED_GROUPS {8};
// Maximum number of secondary particles created
constexpr int MAX_SECONDARY {1000};
// Maximum number of lost particles
constexpr int MAX_LOST_PARTICLES {10};
// 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 +53,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
};
//==============================================================================
// MACROXS contains cached macroscopic cross sections for the material a
// particle is traveling through
//==============================================================================
struct MacroXS {
double total; //!< macroscopic total xs
double absorption; //!< macroscopic absorption xs
double fission; //!< macroscopic fission xs
double nu_fission; //!< macroscopic production xs
double photon_prod; //!< macroscopic photon production xs
// Photon cross sections
double coherent; //!< macroscopic coherent xs
double incoherent; //!< macroscopic incoherent xs
double photoelectric; //!< macroscopic photoelectric xs
double pair_production; //!< macroscopic pair production xs
};
//============================================================================
//! 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,19 +150,87 @@ 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) const;
//! 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
// Cross section caches
std::vector<NuclideMicroXS> neutron_xs_; //!< Microscopic neutron cross sections
std::vector<ElementMicroXS> photon_xs_; //!< Microscopic photon cross sections
MacroXS macro_xs_; //!< Macroscopic cross sections
int64_t id_; //!< Unique ID
Type type_ {Type::neutron}; //!< Particle type (n, p, e, etc.)
int n_coord_ {1}; //!< number of current coordinate levels
int cell_instance_; //!< offset for distributed properties
LocalCoord coord_[MAX_COORD]; //!< coordinates for all levels
std::vector<LocalCoord> coord_; //!< coordinates for all levels
// Particle coordinates before crossing a surface
int n_coord_last_ {1}; //!< number of current coordinates
int cell_last_[MAX_COORD]; //!< coordinates for all levels
std::vector<int> cell_last_; //!< coordinates for all levels
// Energy data
double E_; //!< post-collision energy in eV
@ -135,65 +280,6 @@ public:
// Track output
bool write_track_ {false};
// 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

15
include/openmc/physics.h
View file

@ -7,6 +7,8 @@
#include "openmc/position.h"
#include "openmc/reaction.h"
#include <vector>
namespace openmc {
//==============================================================================
@ -47,24 +49,23 @@ int sample_nuclide(const Particle* p);
//! Determine the average total, prompt, and delayed neutrons produced from
//! fission and creates appropriate bank sites.
void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
Particle::Bank* bank_array, int64_t* bank_size, int64_t bank_capacity);
std::vector<Particle::Bank>& bank);
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

9
include/openmc/physics_mg.h
View file

@ -8,6 +8,8 @@
#include "openmc/particle.h"
#include "openmc/nuclide.h"
#include <vector>
namespace openmc {
//! \brief samples particle behavior after a collision event.
@ -31,12 +33,9 @@ scatter(Particle* p);
//! \brief Determines the average total, prompt and delayed neutrons produced
//! from fission and creates the appropriate bank sites.
//! \param p Particle to operate on
//! \param bank_array The particle bank to populate
//! \param size_bank Number of particles currently in the bank
//! \param bank_array_size Allocated size of the bank
//! \param bank The particle bank to populate
void
create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size);
create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank);
//! \brief Handles an absorption event
//! \param p Particle to operate on

8
include/openmc/simulation.h
View file

@ -37,19 +37,15 @@ extern "C" int restart_batch; //!< batch at which a restart job resumed
extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
extern "C" int total_gen; //!< total number of generations simulated
extern double total_weight; //!< Total source weight in a batch
extern "C" int64_t work; //!< number of particles per process
extern int64_t work_per_rank; //!< number of particles per MPI rank
extern std::vector<double> k_generation;
extern std::vector<int64_t> work_index;
// Threadprivate variables
extern "C" bool trace; //!< flag to show debug information
#ifdef _OPENMP
extern "C" int n_threads; //!< number of OpenMP threads
extern "C" int thread_id; //!< ID of a given thread
#endif
#pragma omp threadprivate(current_work, thread_id, trace)
#pragma omp threadprivate(current_work, trace)
} // namespace simulation

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 {

19
openmc/settings.py
View file

@ -130,8 +130,6 @@ class Settings(object):
range. 'multipole' is a boolean indicating whether or not the windowed
multipole method should be used to evaluate resolved resonance cross
sections.
threads : int
Number of OpenMP threads
trace : tuple or list
Show detailed information about a single particle, indicated by three
integers: the batch number, generation number, and particle number
@ -197,7 +195,6 @@ class Settings(object):
self._statepoint = {}
self._sourcepoint = {}
self._threads = None
self._no_reduce = None
self._verbosity = None
@ -312,10 +309,6 @@ class Settings(object):
def statepoint(self):
return self._statepoint
@property
def threads(self):
return self._threads
@property
def no_reduce(self):
return self._no_reduce
@ -623,12 +616,6 @@ class Settings(object):
self._temperature = temperature
@threads.setter
def threads(self, threads):
cv.check_type('number of threads', threads, Integral)
cv.check_greater_than('number of threads', threads, 0)
self._threads = threads
@trace.setter
def trace(self, trace):
cv.check_type('trace', trace, Iterable, Integral)
@ -885,11 +872,6 @@ class Settings(object):
else:
element.text = str(value)
def _create_threads_subelement(self, root):
if self._threads is not None:
element = ET.SubElement(root, "threads")
element.text = str(self._threads)
def _create_trace_subelement(self, root):
if self._trace is not None:
element = ET.SubElement(root, "trace")
@ -980,7 +962,6 @@ class Settings(object):
self._create_entropy_mesh_subelement(root_element)
self._create_trigger_subelement(root_element)
self._create_no_reduce_subelement(root_element)
self._create_threads_subelement(root_element)
self._create_verbosity_subelement(root_element)
self._create_tabular_legendre_subelements(root_element)
self._create_temperature_subelements(root_element)

5
src/bank.cpp
View file

@ -16,10 +16,9 @@ namespace openmc {
namespace simulation {
int64_t n_bank;
std::vector<Particle::Bank> source_bank;
std::vector<Particle::Bank> fission_bank;
std::vector<Particle::Bank> secondary_bank;
#ifdef _OPENMP
std::vector<Particle::Bank> master_fission_bank;
#endif
@ -33,7 +32,7 @@ std::vector<Particle::Bank> master_fission_bank;
void free_memory_bank()
{
simulation::source_bank.clear();
#pragma omp parallel
#pragma omp parallel
{
simulation::fission_bank.clear();
}

68
src/eigenvalue.cpp
View file

@ -20,9 +20,14 @@
#include "openmc/timer.h"
#include "openmc/tallies/tally.h"
#ifdef _OPENMP
#include <omp.h>
#endif
#include <algorithm> // for min
#include <array>
#include <cmath> // for sqrt, abs, pow
#include <iterator> // for back_inserter
#include <string>
namespace openmc {
@ -81,21 +86,22 @@ void synchronize_bank()
#ifdef OPENMC_MPI
int64_t start = 0;
MPI_Exscan(&simulation::n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
int64_t n_bank = simulation::fission_bank.size();
MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
// While we would expect the value of start on rank 0 to be 0, the MPI
// standard says that the receive buffer on rank 0 is undefined and not
// significant
if (mpi::rank == 0) start = 0;
int64_t finish = start + simulation::n_bank;
int64_t finish = start + simulation::fission_bank.size();
int64_t total = finish;
MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
#else
int64_t start = 0;
int64_t finish = simulation::n_bank;
int64_t total = simulation::n_bank;
int64_t finish = simulation::fission_bank.size();
int64_t total = simulation::fission_bank.size();
#endif
// If there are not that many particles per generation, it's possible that no
@ -103,7 +109,7 @@ void synchronize_bank()
// extra logic to treat this circumstance, we really want to ensure the user
// runs enough particles to avoid this in the first place.
if (simulation::n_bank == 0) {
if (simulation::fission_bank.empty()) {
fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
}
@ -132,23 +138,23 @@ void synchronize_bank()
// Allocate temporary source bank
int64_t index_temp = 0;
std::vector<Particle::Bank> temp_sites(3*simulation::work);
std::vector<Particle::Bank> temp_sites(3*simulation::work_per_rank);
for (int64_t i = 0; i < simulation::n_bank; ++i) {
for (const auto& site : simulation::fission_bank) {
// If there are less than n_particles particles banked, automatically add
// int(n_particles/total) sites to temp_sites. For example, if you need
// 1000 and 300 were banked, this would add 3 source sites per banked site
// and the remaining 100 would be randomly sampled.
if (total < settings::n_particles) {
for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
temp_sites[index_temp] = simulation::fission_bank[i];
temp_sites[index_temp] = site;
++index_temp;
}
}
// Randomly sample sites needed
if (prn() < p_sample) {
temp_sites[index_temp] = simulation::fission_bank[i];
temp_sites[index_temp] = site;
++index_temp;
}
}
@ -189,7 +195,7 @@ void synchronize_bank()
// fission bank
sites_needed = settings::n_particles - finish;
for (int i = 0; i < sites_needed; ++i) {
int i_bank = simulation::n_bank - sites_needed + i;
int i_bank = simulation::fission_bank.size() - sites_needed + i;
temp_sites[index_temp] = simulation::fission_bank[i_bank];
++index_temp;
}
@ -346,38 +352,32 @@ void calculate_average_keff()
#ifdef _OPENMP
void join_bank_from_threads()
{
// Initialize the total number of fission bank sites
int64_t total = 0;
int n_threads = omp_get_max_threads();
#pragma omp parallel
#pragma omp parallel
{
// Copy thread fission bank sites to one shared copy
#pragma omp for ordered schedule(static)
for (int i = 0; i < simulation::n_threads; ++i) {
#pragma omp ordered
#pragma omp for ordered schedule(static)
for (int i = 0; i < n_threads; ++i) {
#pragma omp ordered
{
std::copy(
&simulation::fission_bank[0],
&simulation::fission_bank[0] + simulation::n_bank,
&simulation::master_fission_bank[total]
simulation::fission_bank.cbegin(),
simulation::fission_bank.cend(),
std::back_inserter(simulation::master_fission_bank)
);
total += simulation::n_bank;
}
}
// Make sure all threads have made it to this point
#pragma omp barrier
#pragma omp barrier
// Now copy the shared fission bank sites back to the master thread's copy.
if (simulation::thread_id == 0) {
simulation::n_bank = total;
std::copy(
&simulation::master_fission_bank[0],
&simulation::master_fission_bank[0] + simulation::n_bank,
&simulation::fission_bank[0]
);
if (omp_get_thread_num() == 0) {
simulation::fission_bank = simulation::master_fission_bank;
simulation::master_fission_bank.clear();
} else {
simulation::n_bank = 0;
simulation::fission_bank.clear();
}
}
}
@ -537,8 +537,8 @@ void shannon_entropy()
// Get source weight in each mesh bin
bool sites_outside;
xt::xtensor<double, 1> p = m->count_sites(simulation::n_bank,
simulation::fission_bank.data(), 0, nullptr, &sites_outside);
xt::xtensor<double, 1> p = m->count_sites(simulation::fission_bank,
&sites_outside);
// display warning message if there were sites outside entropy box
if (sites_outside) {
@ -577,8 +577,8 @@ void ufs_count_sites()
} else {
// count number of source sites in each ufs mesh cell
bool sites_outside;
simulation::source_frac = m->count_sites(simulation::work,
simulation::source_bank.data(), 0, nullptr, &sites_outside);
simulation::source_frac = m->count_sites(simulation::source_bank,
&sites_outside);
// Check for sites outside of the mesh
if (mpi::master && sites_outside) {
@ -597,7 +597,7 @@ void ufs_count_sites()
// Since the total starting weight is not equal to n_particles, we need to
// renormalize the weight of the source sites
for (int i = 0; i < simulation::work; ++i) {
for (int i = 0; i < simulation::work_per_rank; ++i) {
simulation::source_bank[i].wgt *= settings::n_particles / total;
}
}

59
src/geometry.cpp
View file

@ -22,6 +22,7 @@ namespace openmc {
namespace model {
int root_universe {-1};
int n_coord_levels;
std::vector<int64_t> overlap_check_count;
@ -31,8 +32,7 @@ std::vector<int64_t> overlap_check_count;
// Non-member functions
//==============================================================================
extern "C" bool
check_cell_overlap(Particle* p)
bool check_cell_overlap(Particle* p)
{
int n_coord = p->n_coord_;
@ -245,7 +245,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
//==============================================================================
extern "C" bool
bool
find_cell(Particle* p, bool use_neighbor_lists)
{
// Determine universe (if not yet set, use root universe).
@ -257,7 +257,7 @@ find_cell(Particle* p, bool use_neighbor_lists)
}
// Reset all the deeper coordinate levels.
for (int i = p->n_coord_; i < MAX_COORD; i++) {
for (int i = p->n_coord_; i < p->coord_.size(); i++) {
p->coord_[i].reset();
}
@ -287,8 +287,8 @@ find_cell(Particle* p, bool use_neighbor_lists)
//==============================================================================
extern "C" void
cross_lattice(Particle* p, int lattice_translation[3])
void
cross_lattice(Particle* p, const BoundaryInfo& boundary)
{
auto& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice]};
@ -302,9 +302,9 @@ cross_lattice(Particle* p, int lattice_translation[3])
}
// Set the lattice indices.
p->coord_[p->n_coord_-1].lattice_x += lattice_translation[0];
p->coord_[p->n_coord_-1].lattice_y += lattice_translation[1];
p->coord_[p->n_coord_-1].lattice_z += lattice_translation[2];
p->coord_[p->n_coord_-1].lattice_x += boundary.lattice_translation[0];
p->coord_[p->n_coord_-1].lattice_y += boundary.lattice_translation[1];
p->coord_[p->n_coord_-1].lattice_z += boundary.lattice_translation[2];
std::array<int, 3> i_xyz {p->coord_[p->n_coord_-1].lattice_x,
p->coord_[p->n_coord_-1].lattice_y,
p->coord_[p->n_coord_-1].lattice_z};
@ -345,16 +345,11 @@ cross_lattice(Particle* p, int lattice_translation[3])
//==============================================================================
extern "C" void
distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
int lattice_translation[3], int* next_level)
BoundaryInfo distance_to_boundary(Particle* p)
{
*dist = INFINITY;
BoundaryInfo info;
double d_lat = INFINITY;
double d_surf = INFINITY;
lattice_translation[0] = 0;
lattice_translation[1] = 0;
lattice_translation[2] = 0;
int32_t level_surf_cross;
std::array<int, 3> level_lat_trans;
@ -401,43 +396,43 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
// If the boundary on this coordinate level is coincident with a boundary on
// a higher level then we need to make sure that the higher level boundary
// is selected. This logic must consider floating point precision.
double& d = info.distance;
if (d_surf < d_lat) {
if (*dist == INFINITY || ((*dist) - d_surf)/(*dist) >= FP_REL_PRECISION) {
*dist = d_surf;
if (d == INFINITY || (d - d_surf)/d >= FP_REL_PRECISION) {
d = d_surf;
// If the cell is not simple, it is possible that both the negative and
// positive half-space were given in the region specification. Thus, we
// have to explicitly check which half-space the particle would be
// traveling into if the surface is crossed
if (c.simple_) {
*surface_crossed = level_surf_cross;
info.surface_index = level_surf_cross;
} else {
Position r_hit = r + d_surf * u;
Surface& surf {*model::surfaces[std::abs(level_surf_cross)-1]};
Direction norm = surf.normal(r_hit);
if (u.dot(norm) > 0) {
*surface_crossed = std::abs(level_surf_cross);
info.surface_index = std::abs(level_surf_cross);
} else {
*surface_crossed = -std::abs(level_surf_cross);
info.surface_index = -std::abs(level_surf_cross);
}
}
lattice_translation[0] = 0;
lattice_translation[1] = 0;
lattice_translation[2] = 0;
*next_level = i + 1;
info.lattice_translation[0] = 0;
info.lattice_translation[1] = 0;
info.lattice_translation[2] = 0;
info.coord_level = i + 1;
}
} else {
if (*dist == INFINITY || ((*dist) - d_lat)/(*dist) >= FP_REL_PRECISION) {
*dist = d_lat;
*surface_crossed = F90_NONE;
lattice_translation[0] = level_lat_trans[0];
lattice_translation[1] = level_lat_trans[1];
lattice_translation[2] = level_lat_trans[2];
*next_level = i + 1;
if (d == INFINITY || (d - d_lat)/d >= FP_REL_PRECISION) {
d = d_lat;
info.surface_index = 0;
info.lattice_translation = level_lat_trans;
info.coord_level = i + 1;
}
}
}
return info;
}
//==============================================================================

10
src/geometry_aux.cpp
View file

@ -217,14 +217,8 @@ void finalize_geometry(std::vector<std::vector<double>>& nuc_temps,
// Determine desired temperatures for each nuclide and S(a,b) table
get_temperatures(nuc_temps, thermal_temps);
// Check to make sure there are not too many nested coordinate levels in the
// geometry since the coordinate list is statically allocated for performance
// reasons
if (maximum_levels(model::root_universe) > MAX_COORD) {
fatal_error("Too many nested coordinate levels in the geometry. "
"Try increasing the maximum number of coordinate levels by "
"providing the CMake -Dmaxcoord= option.");
}
// Determine number of nested coordinate levels in the geometry
model::n_coord_levels = maximum_levels(model::root_universe);
}
//==============================================================================

6
src/initialize.cpp
View file

@ -196,13 +196,13 @@ parse_command_line(int argc, char* argv[])
#ifdef _OPENMP
// Read and set number of OpenMP threads
simulation::n_threads = std::stoi(argv[i]);
if (simulation::n_threads < 1) {
int n_threads = std::stoi(argv[i]);
if (n_threads < 1) {
std::string msg {"Number of threads must be positive."};
strcpy(openmc_err_msg, msg.c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
omp_set_num_threads(simulation::n_threads);
omp_set_num_threads(n_threads);
#else
if (mpi::master)
warning("Ignoring number of threads specified on command line.");

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.macro_xs_.total = 0.0;
p.macro_xs_.absorption = 0.0;
p.macro_xs_.fission = 0.0;
p.macro_xs_.nu_fission = 0.0;
if (p.type_ == Particle::Type::neutron) {
this->calculate_neutron_xs(p);
@ -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.neutron_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.macro_xs_.total += atom_density * micro.total;
p.macro_xs_.absorption += atom_density * micro.absorption;
p.macro_xs_.fission += atom_density * micro.fission;
p.macro_xs_.nu_fission += atom_density * micro.nu_fission;
}
}
void Material::calculate_photon_xs(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.macro_xs_.coherent = 0.0;
p.macro_xs_.incoherent = 0.0;
p.macro_xs_.photoelectric = 0.0;
p.macro_xs_.pair_production = 0.0;
// Add contribution from each nuclide in material
for (int i = 0; i < nuclide_.size(); ++i) {
@ -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.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.macro_xs_.total += atom_density * micro.total;
p.macro_xs_.coherent += atom_density * micro.coherent;
p.macro_xs_.incoherent += atom_density * micro.incoherent;
p.macro_xs_.photoelectric += atom_density * micro.photoelectric;
p.macro_xs_.pair_production += atom_density * micro.pair_production;
}
}

31
src/mesh.cpp
View file

@ -732,25 +732,21 @@ void RegularMesh::to_hdf5(hid_t group) const
close_group(mesh_group);
}
xt::xarray<double> RegularMesh::count_sites(int64_t n, const Particle::Bank* bank,
int n_energy, const double* energies, bool* outside) const
xt::xarray<double>
RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
bool* outside) const
{
// Determine shape of array for counts
std::size_t m = xt::prod(shape_)();
std::vector<std::size_t> shape;
if (n_energy > 0) {
shape = {m, static_cast<std::size_t>(n_energy - 1)};
} else {
shape = {m};
}
std::vector<std::size_t> shape = {m};
// Create array of zeros
xt::xarray<double> cnt {shape, 0.0};
bool outside_ = false;
for (int64_t i = 0; i < n; ++i) {
for (const auto& site : bank) {
// determine scoring bin for entropy mesh
int mesh_bin = get_bin(bank[i].r);
int mesh_bin = get_bin(site.r);
// if outside mesh, skip particle
if (mesh_bin < 0) {
@ -758,19 +754,8 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Particle::Bank* ban
continue;
}
if (n_energy > 0) {
double E = bank[i].E;
if (E >= energies[0] && E <= energies[n_energy - 1]) {
// determine energy bin
int e_bin = lower_bound_index(energies, energies + n_energy, E);
// Add to appropriate bin
cnt(mesh_bin, e_bin) += bank[i].wgt;
}
} else {
// Add to appropriate bin
cnt(mesh_bin) += bank[i].wgt;
}
// Add to appropriate bin
cnt(mesh_bin) += site.wgt;
}
// Create copy of count data

127
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
//==============================================================================
@ -459,12 +454,15 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const
return (*fission_rx_[0]->products_[0].yield_)(E);
}
}
#ifdef __GNUC__
__builtin_unreachable();
#endif
}
void Nuclide::calculate_elastic_xs() 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.neutron_xs_[i_nuclide_]};
int i_temp = micro.index_temp;
int i_grid = micro.index_grid;
double f = micro.interp_factor;
@ -498,42 +496,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.neutron_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 +544,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 +587,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 +598,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 +655,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 +677,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.neutron_xs_[i_nuclide_]};
// Set flag that S(a,b) treatment should be used for scattering
micro.index_sab = i_sab;
@ -717,14 +714,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 +732,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.neutron_xs_[i_nuclide_];
micro.use_ptable = true;
// Create a shorthand for the URR data
@ -745,7 +742,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 +770,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 +781,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 +835,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 +855,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;
}
}

4
src/output.cpp
View file

@ -11,7 +11,9 @@
#include <unordered_map>
#include <utility> // for pair
#ifdef _OPENMP
#include <omp.h>
#endif
#include "xtensor/xview.hpp"
#include "openmc/capi.h"
@ -189,7 +191,7 @@ extern "C" void print_particle(Particle* p)
}
// Display miscellaneous info.
if (p->surface_ != ERROR_INT) {
if (p->surface_ != 0) {
const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]};
std::cout << " Surface = " << std::copysign(surf.id_, p->surface_) << "\n";
}

94
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"
@ -49,12 +50,18 @@ LocalCoord::reset()
Particle::Particle()
{
// Clear coordinate lists
// Create and clear coordinate levels
coord_.resize(model::n_coord_levels);
cell_last_.resize(model::n_coord_levels);
clear();
for (int& n : n_delayed_bank_) {
n = 0;
}
// Create microscopic cross section caches
neutron_xs_.resize(data::nuclides.size());
photon_xs_.resize(data::elements.size());
}
void
@ -66,19 +73,16 @@ Particle::clear()
}
void
Particle::create_secondary(Direction u, double E, Type type)
Particle::create_secondary(Direction u, double E, Type type) const
{
if (n_secondary_ == MAX_SECONDARY) {
fatal_error("Too many secondary particles created.");
}
simulation::secondary_bank.emplace_back();
int64_t n = n_secondary_;
secondary_bank_[n].particle = type;
secondary_bank_[n].wgt = wgt_;
secondary_bank_[n].r = this->r();
secondary_bank_[n].u = u;
secondary_bank_[n].E = settings::run_CE ? E : g_;
++n_secondary_;
auto& bank {simulation::secondary_bank.back()};
bank.particle = type;
bank.wgt = wgt_;
bank.r = this->r();
bank.u = u;
bank.E = settings::run_CE ? E : g_;
}
void
@ -130,9 +134,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 : neutron_xs_) micro.last_E = 0.0;
}
// Prepare to write out particle track.
@ -186,41 +188,35 @@ 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);
macro_xs_.total, macro_xs_.absorption, macro_xs_.nu_fission);
// Finally, update the particle group while we have already checked
// for if multi-group
g_last_ = g_;
}
} else {
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;
macro_xs_.total = 0.0;
macro_xs_.absorption = 0.0;
macro_xs_.fission = 0.0;
macro_xs_.nu_fission = 0.0;
}
// Find the distance to the nearest boundary
double d_boundary;
int surface_crossed;
int lattice_translation[3];
int next_level;
distance_to_boundary(this, &d_boundary, &surface_crossed,
lattice_translation, &next_level);
auto boundary = distance_to_boundary(this);
// Sample a distance to collision
double d_collision;
if (type_ == Particle::Type::electron ||
type_ == Particle::Type::positron) {
d_collision = 0.0;
} else if (simulation::material_xs.total == 0.0) {
} else if (macro_xs_.total == 0.0) {
d_collision = INFINITY;
} else {
d_collision = -std::log(prn()) / simulation::material_xs.total;
d_collision = -std::log(prn()) / macro_xs_.total;
}
// Select smaller of the two distances
double distance = std::min(d_boundary, d_collision);
double distance = std::min(boundary.distance, d_collision);
// Advance particle
for (int j = 0; j < n_coord_; ++j) {
@ -235,7 +231,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 * macro_xs_.nu_fission;
}
// Score flux derivative accumulators for differential tallies.
@ -243,11 +239,13 @@ Particle::transport()
score_track_derivative(this, distance);
}
if (d_collision > d_boundary) {
if (d_collision > boundary.distance) {
// ====================================================================
// PARTICLE CROSSES SURFACE
if (next_level > 0) n_coord_ = next_level;
// Set surface that particle is on and adjust coordinate levels
surface_ = boundary.surface_index;
n_coord_ = boundary.coord_level;
// Saving previous cell data
for (int j = 0; j < n_coord_; ++j) {
@ -255,15 +253,14 @@ Particle::transport()
}
n_coord_last_ = n_coord_;
if (lattice_translation[0] != 0 || lattice_translation[1] != 0 ||
lattice_translation[2] != 0) {
if (boundary.lattice_translation[0] != 0 ||
boundary.lattice_translation[1] != 0 ||
boundary.lattice_translation[2] != 0) {
// Particle crosses lattice boundary
surface_ = ERROR_INT;
cross_lattice(this, lattice_translation);
cross_lattice(this, boundary);
event_ = EVENT_LATTICE;
} else {
// Particle crosses surface
surface_ = surface_crossed;
this->cross_surface();
event_ = EVENT_SURFACE;
}
@ -278,8 +275,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_ * macro_xs_.nu_fission
/ macro_xs_.total;
}
// Score surface current tallies -- this has to be done before the collision
@ -290,7 +287,7 @@ Particle::transport()
score_surface_tally(this, model::active_meshsurf_tallies);
// Clear surface component
surface_ = ERROR_INT;
surface_ = 0;
if (settings::run_CE) {
collision(this);
@ -356,10 +353,10 @@ Particle::transport()
// Check for secondary particles if this particle is dead
if (!alive_) {
// If no secondary particles, break out of event loop
if (n_secondary_ == 0) break;
if (simulation::secondary_bank.empty()) break;
this->from_source(&secondary_bank_[n_secondary_ - 1]);
--n_secondary_;
this->from_source(&simulation::secondary_bank.back());
simulation::secondary_bank.pop_back();
n_event = 0;
// Enter new particle in particle track file
@ -555,7 +552,7 @@ Particle::cross_surface()
// COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
// Remove lower coordinate levels and assignment of surface
surface_ = ERROR_INT;
surface_ = 0;
n_coord_ = 1;
bool found = find_cell(this, false);
@ -589,11 +586,12 @@ Particle::mark_as_lost(const char* message)
// Increment number of lost particles
alive_ = false;
#pragma omp atomic
#pragma omp atomic
simulation::n_lost_particles += 1;
// Count the total number of simulated particles (on this processor)
auto n = simulation::current_batch * settings::gen_per_batch * simulation::work;
auto n = simulation::current_batch * settings::gen_per_batch *
simulation::work_per_rank;
// Abort the simulation if the maximum number of lost particles has been
// reached
@ -614,7 +612,7 @@ Particle::write_restart() const
filename << settings::path_output << "particle_" << simulation::current_batch
<< '_' << id_ << ".h5";
#pragma omp critical (WriteParticleRestart)
#pragma omp critical (WriteParticleRestart)
{
// Create file
hid_t file_id = file_open(filename.str(), 'w');

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

176
src/physics.cpp
View file

@ -89,14 +89,20 @@ 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());
create_fission_sites(p, i_nuclide, rx, simulation::fission_bank);
} else if (settings::run_mode == RUN_MODE_FIXEDSOURCE &&
settings::create_fission_neutrons) {
create_fission_sites(p, i_nuclide, rx, p->secondary_bank_,
&p->n_secondary_, MAX_SECONDARY);
create_fission_sites(p, i_nuclide, rx, simulation::secondary_bank);
// Make sure particle population doesn't grow out of control for
// subcritical multiplication problems.
if (simulation::secondary_bank.size() >= 10000) {
fatal_error("The secondary particle bank appears to be growing without "
"bound. You are likely running a subcritical multiplication problem "
"with k-effective close to or greater than one.");
}
}
}
@ -110,7 +116,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->neutron_xs_[i_nuclide].absorption > 0.0) {
absorption(p, i_nuclide);
} else {
p->wgt_absorb_ = 0.0;
@ -137,63 +143,44 @@ void sample_neutron_reaction(Particle* p)
void
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
Particle::Bank* bank_array, int64_t* size_bank, int64_t bank_capacity)
std::vector<Particle::Bank>& bank)
{
// TODO: Heat generation from fission
// If uniform fission source weighting is turned on, we increase or decrease
// the expected number of fission sites produced
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
// 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->neutron_xs_[
i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].total;
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
if (prn() <= (nu_t - nu)) ++nu;
// Check for the bank size getting hit. For fixed source calculations, this
// is a fatal error; for eigenvalue calculations, it just means that k-eff
// was too high for a single batch.
if (*size_bank + nu > bank_capacity) {
if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
throw std::runtime_error{"Secondary particle bank size limit reached."
" If you are running a subcritical multiplication problem,"
" k-effective may be too close to one."};
} else {
if (mpi::master) {
warning("Maximum number of sites in fission bank reached. This can"
" result in irreproducible results using different numbers of"
" processes/threads.");
}
}
}
// Begin banking the source neutrons
// First, if our bank is full then don't continue
if (nu == 0 || *size_bank == bank_capacity) return;
if (nu == 0) return;
// Initialize the counter of delayed neutrons encountered for each delayed
// group.
double nu_d[MAX_DELAYED_GROUPS] = {0.};
p->fission_ = true;
for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) {
for (int i = 0; i < nu; ++i) {
// Create new bank site and get reference to last element
bank.emplace_back();
auto& site {bank.back()};
// Bank source neutrons by copying the particle data
bank_array[i].r = p->r();
// Set that the bank particle is a neutron
bank_array[i].particle = Particle::Type::neutron;
// Set the weight of the fission bank site
bank_array[i].wgt = 1. / weight;
site.r = p->r();
site.particle = Particle::Type::neutron;
site.wgt = 1. / weight;
// Sample delayed group and angle/energy for fission reaction
sample_fission_neutron(i_nuclide, rx, p->E_, &bank_array[i]);
sample_fission_neutron(i_nuclide, rx, p->E_, &site);
// Set the delayed group on the particle as well
p->delayed_group_ = bank_array[i].delayed_group;
p->delayed_group_ = site.delayed_group;
// Increment the number of neutrons born delayed
if (p->delayed_group_ > 0) {
@ -201,9 +188,6 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
}
}
// Increment number of bank sites
*size_bank = std::min(*size_bank + nu, bank_capacity);
// Store the total weight banked for analog fission tallies
p->n_bank_ = nu;
p->wgt_bank_ = nu / weight;
@ -227,7 +211,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->photon_xs_[i_element]};
const auto& element {data::elements[i_element]};
// Calculate photon energy over electron rest mass equivalent
@ -408,7 +392,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->macro_xs_.total;
// Get pointers to nuclide/density arrays
const auto& mat {model::materials[p->material_]};
@ -421,7 +405,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->neutron_xs_[i_nuclide].total;
if (prob >= cutoff) return i_nuclide;
}
@ -433,7 +417,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->macro_xs_.total;
// Get pointers to elements, densities
const auto& mat {model::materials[p->material_]};
@ -454,7 +438,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->photon_xs_[i_element].total;
// Increment probability to compare to cutoff
prob += sigma;
@ -464,7 +448,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 +456,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->neutron_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->neutron_xs_[i_nuclide].index_temp;
int i_grid = p->neutron_xs_[i_nuclide].index_grid;
double f = p->neutron_xs_[i_nuclide].interp_factor;
double cutoff = prn() * p->neutron_xs_[i_nuclide].fission;
double prob = 0.0;
// Loop through each partial fission reaction type
@ -509,13 +493,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->neutron_xs_[i_nuclide].index_temp;
int i_grid = p->neutron_xs_[i_nuclide].index_grid;
double f = p->neutron_xs_[i_nuclide].interp_factor;
double cutoff = prn() * p->neutron_xs_[i_nuclide].photon_prod;
double prob = 0.0;
// Loop through each reaction type
@ -534,7 +518,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 +532,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->neutron_xs_[i_nuclide].absorption /
p->neutron_xs_[i_nuclide].total;
// Adjust weight of particle by probability of absorption
p->wgt_ -= p->wgt_absorb_;
@ -557,17 +541,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->neutron_xs_[
i_nuclide].nu_fission / p->neutron_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->neutron_xs_[i_nuclide].absorption >
prn() * p->neutron_xs_[i_nuclide].total) {
// Score absorption estimate of keff
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
global_tally_absorption += p->wgt_ * simulation::micro_xs[
i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption;
global_tally_absorption += p->wgt_ * p->neutron_xs_[
i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
}
p->alive_ = false;
@ -584,7 +568,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->neutron_xs_[i_nuclide]};
int i_temp = micro.index_temp;
int i_grid = micro.index_grid;
double f = micro.interp_factor;
@ -596,7 +580,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 +592,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 +603,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 +661,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->neutron_xs_[i_nuclide].use_ptable) {
v_t = sample_target_velocity(nuc.get(), p->E_, p->u(), v_n,
p->neutron_xs_[i_nuclide].elastic, kT);
}
// Velocity of center-of-mass
@ -712,7 +695,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 +711,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->neutron_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,
@ -907,8 +890,9 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
} // case RVS, DBRC
} // switch (sampling_method)
throw std::runtime_error{"Unable to sample target velocity for "
"elastic scattering."};
#ifdef __GNUC__
__builtin_unreachable();
#endif
}
Direction
@ -1105,8 +1089,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->neutron_xs_[i_nuclide].photon_prod /
p->neutron_xs_[i_nuclide].total;
int y = static_cast<int>(y_t);
if (prn() <= y_t - y) ++y;
@ -1115,7 +1099,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];

88
src/physics_mg.cpp
View file

@ -48,19 +48,16 @@ sample_reaction(Particle* p)
if (model::materials[p->material_]->fissionable_) {
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
create_fission_sites(
p, simulation::fission_bank.data(), &simulation::n_bank,
simulation::fission_bank.size());
create_fission_sites(p, simulation::fission_bank);
} else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
(settings::create_fission_neutrons)) {
create_fission_sites(p, p->secondary_bank_, &(p->n_secondary_),
MAX_SECONDARY);
create_fission_sites(p, simulation::secondary_bank);
}
}
// 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->macro_xs_.absorption > 0.) {
absorption(p);
} else {
p->wgt_absorb_ = 0.;
@ -102,18 +99,15 @@ scatter(Particle* p)
}
void
create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size)
create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
{
// TODO: Heat generation from fission
// If uniform fission source weighting is turned on, we increase or decrease
// the expected number of fission sites produced
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
// 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->macro_xs_.nu_fission / p->macro_xs_.total;
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
@ -121,62 +115,42 @@ create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank
nu++;
}
// Check for the bank size getting hit. For fixed source calculations, this
// is a fatal error; for eigenvalue calculations, it just means that k-eff
// was too high for a single batch.
if (*size_bank + nu > bank_array_size) {
if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
throw std::runtime_error{"Secondary particle bank size limit reached."
" If you are running a subcritical multiplication problem,"
" k-effective may be too close to one."};
} else {
if (mpi::master) {
std::stringstream msg;
msg << "Maximum number of sites in fission bank reached. This can"
" result in irreproducible results using different numbers of"
" processes/threads.";
warning(msg);
}
}
}
// Begin banking the source neutrons
// First, if our bank is full then don't continue
if ((nu == 0) || (*size_bank == bank_array_size)) return;
if (nu == 0) return;
// Initialize the counter of delayed neutrons encountered for each delayed
// group.
double nu_d[MAX_DELAYED_GROUPS] = {0.};
p->fission_ = true;
for (size_t i = static_cast<size_t>(*size_bank);
i < static_cast<size_t>(std::min(*size_bank + nu, bank_array_size)); i++) {
for (int i = 0; i < nu; ++i) {
// Create new bank site and get reference to last element
bank.emplace_back();
auto& site {bank.back()};
// Bank source neutrons by copying the particle data
bank_array[i].r = p->r();
// Set that the bank particle is a neutron
bank_array[i].particle = Particle::Type::neutron;
// Set the weight of the fission bank site
bank_array[i].wgt = 1. / weight;
site.r = p->r();
site.particle = Particle::Type::neutron;
site.wgt = 1. / weight;
// Sample the cosine of the angle, assuming fission neutrons are emitted
// isotropically
double mu = 2. * prn() - 1.;
double mu = 2.*prn() - 1.;
// Sample the azimuthal angle uniformly in [0, 2.pi)
double phi = 2. * PI * prn();
bank_array[i].u.x = mu;
bank_array[i].u.y = std::sqrt(1. - mu * mu) * std::cos(phi);
bank_array[i].u.z = std::sqrt(1. - mu * mu) * std::sin(phi);
site.u.x = mu;
site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi);
site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi);
// Sample secondary energy distribution for the fission reaction and set
// the energy in the fission bank
int dg;
int gout;
data::macro_xs[p->material_].sample_fission_energy(p->g_ - 1, dg, gout);
bank_array[i].E = gout + 1;
bank_array[i].delayed_group = dg + 1;
site.E = gout + 1;
site.delayed_group = dg + 1;
// Set the delayed group on the particle as well
p->delayed_group_ = dg + 1;
@ -187,9 +161,6 @@ create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank
}
}
// Increment number of bank sites
*size_bank = std::min(*size_bank + nu, bank_array_size);
// Store the total weight banked for analog fission tallies
p->n_bank_ = nu;
p->wgt_bank_ = nu / weight;
@ -203,24 +174,21 @@ 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->wgt_absorb_ = p->wgt_ * p->macro_xs_.absorption / p->macro_xs_.total;
// Adjust weight of particle by the probability of absorption
p->wgt_ -= p->wgt_absorb_;
p->wgt_last_ = p->wgt_;
// Score implicit absorpion estimate of keff
#pragma omp atomic
global_tally_absorption += p->wgt_absorb_ *
simulation::material_xs.nu_fission /
simulation::material_xs.absorption;
#pragma omp atomic
global_tally_absorption += p->wgt_absorb_ * p->macro_xs_.nu_fission /
p->macro_xs_.absorption;
} else {
if (simulation::material_xs.absorption >
prn() * simulation::material_xs.total) {
#pragma omp atomic
global_tally_absorption += p->wgt_ * simulation::material_xs.nu_fission /
simulation::material_xs.absorption;
if (p->macro_xs_.absorption > prn() * p->macro_xs_.total) {
#pragma omp atomic
global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission /
p->macro_xs_.absorption;
p->alive_ = false;
p->event_ = EVENT_ABSORB;
}

2
src/random_lcg.cpp
View file

@ -141,7 +141,7 @@ extern "C" void
openmc_set_seed(int64_t new_seed)
{
seed = new_seed;
#pragma omp parallel
#pragma omp parallel
{
for (int i = 0; i < N_STREAMS; i++) {
prn_seed[i] = seed + i;

16
src/settings.cpp
View file

@ -383,20 +383,8 @@ void read_settings_xml()
// Number of OpenMP threads
if (check_for_node(root, "threads")) {
#ifdef _OPENMP
if (simulation::n_threads == 0) {
simulation::n_threads = std::stoi(get_node_value(root, "threads"));
if (simulation::n_threads < 1) {
std::stringstream msg;
msg << "Invalid number of threads: " << simulation::n_threads;
fatal_error(msg);
}
omp_set_num_threads(simulation::n_threads);
}
#else
if (mpi::master) warning("OpenMC was not compiled with OpenMP support; "
"ignoring number of threads.");
#endif
if (mpi::master) warning("The <threads> element has been deprecated. Use "
"the OMP_NUM_THREADS environment variable to set the number of threads.");
}
// ==========================================================================

58
src/simulation.cpp
View file

@ -20,7 +20,9 @@
#include "openmc/tallies/tally.h"
#include "openmc/tallies/trigger.h"
#ifdef _OPENMP
#include <omp.h>
#endif
#include "xtensor/xview.hpp"
#include <algorithm>
@ -59,7 +61,7 @@ int openmc_simulation_init()
calculate_work();
// Allocate array for matching filter bins
#pragma omp parallel
#pragma omp parallel
{
simulation::filter_matches.resize(model::tally_filters.size());
}
@ -78,13 +80,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 +128,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;
@ -150,7 +138,7 @@ int openmc_simulation_finalize()
// Write tally results to tallies.out
if (settings::output_tallies && mpi::master) write_tallies();
#pragma omp parallel
#pragma omp parallel
{
simulation::filter_matches.clear();
}
@ -200,8 +188,8 @@ int openmc_next_batch(int* status)
// ====================================================================
// LOOP OVER PARTICLES
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work; ++i_work) {
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
simulation::current_work = i_work;
// grab source particle from bank
@ -257,17 +245,13 @@ int restart_batch;
bool satisfy_triggers {false};
int total_gen {0};
double total_weight;
int64_t work;
int64_t work_per_rank;
std::vector<double> k_generation;
std::vector<int64_t> work_index;
// Threadprivate variables
bool trace; //!< flag to show debug information
#ifdef _OPENMP
int n_threads {-1}; //!< number of OpenMP threads
int thread_id; //!< ID of a given thread
#endif
} // namespace simulation
@ -278,7 +262,7 @@ int thread_id; //!< ID of a given thread
void allocate_banks()
{
// Allocate source bank
simulation::source_bank.resize(simulation::work);
simulation::source_bank.resize(simulation::work_per_rank);
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
#ifdef _OPENMP
@ -287,20 +271,18 @@ void allocate_banks()
// a generation, there is also a 'master_fission_bank' that is used to
// collect the sites from each thread.
simulation::n_threads = omp_get_max_threads();
#pragma omp parallel
#pragma omp parallel
{
simulation::thread_id = omp_get_thread_num();
if (simulation::thread_id == 0) {
simulation::fission_bank.resize(3*simulation::work);
if (omp_get_thread_num() == 0) {
simulation::fission_bank.reserve(3*simulation::work_per_rank);
} else {
simulation::fission_bank.resize(3*simulation::work / simulation::n_threads);
int n_threads = omp_get_num_threads();
simulation::fission_bank.reserve(3*simulation::work_per_rank / n_threads);
}
}
simulation::master_fission_bank.resize(3*simulation::work);
simulation::master_fission_bank.reserve(3*simulation::work_per_rank);
#else
simulation::fission_bank.resize(3*simulation::work);
simulation::fission_bank.reserve(3*simulation::work_per_rank);
#endif
}
}
@ -400,8 +382,8 @@ void finalize_batch()
void initialize_generation()
{
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
// Reset number of fission bank sites
simulation::n_bank = 0;
// Clear out the fission bank
simulation::fission_bank.clear();
// Count source sites if using uniform fission source weighting
if (settings::ufs_on) ufs_count_sites();
@ -417,9 +399,9 @@ void finalize_generation()
auto& gt = simulation::global_tallies;
// Update global tallies with the omp private accumulation variables
#pragma omp parallel
#pragma omp parallel
{
#pragma omp critical(increment_global_tallies)
#pragma omp critical(increment_global_tallies)
{
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
gt(K_COLLISION, RESULT_VALUE) += global_tally_collision;
@ -525,7 +507,7 @@ void calculate_work()
int64_t work_i = i < remainder ? min_work + 1 : min_work;
// Set number of particles
if (mpi::rank == i) simulation::work = work_i;
if (mpi::rank == i) simulation::work_per_rank = work_i;
// Set index into source bank for rank i
i_bank += work_i;

4
src/source.cpp
View file

@ -264,7 +264,7 @@ void initialize_source()
} else {
// Generation source sites from specified distribution in user input
for (int64_t i = 0; i < simulation::work; ++i) {
for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
// initialize random number seed
int64_t id = simulation::total_gen*settings::n_particles +
simulation::work_index[mpi::rank] + i + 1;
@ -330,7 +330,7 @@ void free_memory_source()
void fill_source_bank_fixedsource()
{
if (settings::path_source.empty()) {
for (int64_t i = 0; i < simulation::work; ++i) {
for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
// initialize random number seed
int64_t id = (simulation::total_gen + overall_generation()) *
settings::n_particles + simulation::work_index[mpi::rank] + i + 1;

8
src/state_point.cpp
View file

@ -537,7 +537,7 @@ write_source_bank(hid_t group_id)
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
// Create another data space but for each proc individually
hsize_t count[] {static_cast<hsize_t>(simulation::work)};
hsize_t count[] {static_cast<hsize_t>(simulation::work_per_rank)};
hid_t memspace = H5Screate_simple(1, count, nullptr);
// Select hyperslab for this dataspace
@ -569,7 +569,7 @@ write_source_bank(hid_t group_id)
// Save source bank sites since the souce_bank array is overwritten below
#ifdef OPENMC_MPI
std::vector<Particle::Bank> temp_source {simulation::source_bank.begin(),
simulation::source_bank.begin() + simulation::work};
simulation::source_bank.begin() + simulation::work_per_rank};
#endif
for (int i = 0; i < mpi::n_procs; ++i) {
@ -607,7 +607,7 @@ write_source_bank(hid_t group_id)
#endif
} else {
#ifdef OPENMC_MPI
MPI_Send(simulation::source_bank.data(), simulation::work, mpi::bank,
MPI_Send(simulation::source_bank.data(), simulation::work_per_rank, mpi::bank,
0, mpi::rank, mpi::intracomm);
#endif
}
@ -625,7 +625,7 @@ void read_source_bank(hid_t group_id)
hid_t dset = H5Dopen(group_id, "source_bank", H5P_DEFAULT);
// Create another data space but for each proc individually
hsize_t dims[] {static_cast<hsize_t>(simulation::work)};
hsize_t dims[] {static_cast<hsize_t>(simulation::work_per_rank)};
hid_t memspace = H5Screate_simple(1, dims, nullptr);
// Make sure source bank is big enough

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->macro_xs_.total > 0.0) {
score *= flux_deriv
+ simulation::micro_xs[deriv.diff_nuclide].total
/ simulation::material_xs.total;
+ p->neutron_xs_[deriv.diff_nuclide].total
/ p->macro_xs_.total;
} else if (i_nuclide == deriv.diff_nuclide
&& simulation::micro_xs[i_nuclide].total) {
&& p->neutron_xs_[i_nuclide].total) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
@ -234,13 +234,13 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
break;
case SCORE_SCATTER:
if (i_nuclide == -1 && (simulation::material_xs.total
- simulation::material_xs.absorption) > 0.0) {
if (i_nuclide == -1 && (p->macro_xs_.total
- p->macro_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->neutron_xs_[deriv.diff_nuclide].total
- p->neutron_xs_[deriv.diff_nuclide].absorption)
/ (p->macro_xs_.total
- p->macro_xs_.absorption);
} else if (i_nuclide == deriv.diff_nuclide) {
score *= flux_deriv + 1. / atom_density;
} else {
@ -249,12 +249,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
break;
case SCORE_ABSORPTION:
if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
if (i_nuclide == -1 && p->macro_xs_.absorption > 0.0) {
score *= flux_deriv
+ simulation::micro_xs[deriv.diff_nuclide].absorption
/ simulation::material_xs.absorption;
+ p->neutron_xs_[deriv.diff_nuclide].absorption
/ p->macro_xs_.absorption;
} else if (i_nuclide == deriv.diff_nuclide
&& simulation::micro_xs[i_nuclide].absorption) {
&& p->neutron_xs_[i_nuclide].absorption) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
@ -262,12 +262,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
break;
case SCORE_FISSION:
if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
if (i_nuclide == -1 && p->macro_xs_.fission > 0.0) {
score *= flux_deriv
+ simulation::micro_xs[deriv.diff_nuclide].fission
/ simulation::material_xs.fission;
+ p->neutron_xs_[deriv.diff_nuclide].fission
/ p->macro_xs_.fission;
} else if (i_nuclide == deriv.diff_nuclide
&& simulation::micro_xs[i_nuclide].fission) {
&& p->neutron_xs_[i_nuclide].fission) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
@ -275,12 +275,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
break;
case SCORE_NU_FISSION:
if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
if (i_nuclide == -1 && p->macro_xs_.nu_fission > 0.0) {
score *= flux_deriv
+ simulation::micro_xs[deriv.diff_nuclide].nu_fission
/ simulation::material_xs.nu_fission;
+ p->neutron_xs_[deriv.diff_nuclide].nu_fission
/ p->macro_xs_.nu_fission;
} else if (i_nuclide == deriv.diff_nuclide
&& simulation::micro_xs[i_nuclide].nu_fission) {
&& p->neutron_xs_[i_nuclide].nu_fission) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
@ -331,64 +331,64 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
switch (score_bin) {
case SCORE_TOTAL:
if (simulation::micro_xs[p->event_nuclide_].total) {
if (p->neutron_xs_[p->event_nuclide_].total) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i)
/ simulation::material_xs.total;
/ p->macro_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->neutron_xs_[p->event_nuclide_].total
- p->neutron_xs_[p->event_nuclide_].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_s * material.atom_density_(i)
/ (simulation::material_xs.total
- simulation::material_xs.absorption);
/ (p->macro_xs_.total
- p->macro_xs_.absorption);
} else {
score *= flux_deriv;
}
break;
case SCORE_ABSORPTION:
if (simulation::micro_xs[p->event_nuclide_].absorption) {
if (p->neutron_xs_[p->event_nuclide_].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_a * material.atom_density_(i)
/ simulation::material_xs.absorption;
/ p->macro_xs_.absorption;
} else {
score *= flux_deriv;
}
break;
case SCORE_FISSION:
if (simulation::micro_xs[p->event_nuclide_].fission) {
if (p->neutron_xs_[p->event_nuclide_].fission) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_f * material.atom_density_(i)
/ simulation::material_xs.fission;
/ p->macro_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->neutron_xs_[p->event_nuclide_].fission) {
double nu = p->neutron_xs_[p->event_nuclide_].nu_fission
/ p->neutron_xs_[p->event_nuclide_].fission;
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + nu * dsig_f * material.atom_density_(i)
/ simulation::material_xs.nu_fission;
/ p->macro_xs_.nu_fission;
} else {
score *= flux_deriv;
}
@ -413,141 +413,141 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
switch (score_bin) {
case SCORE_TOTAL:
if (i_nuclide == -1 && simulation::material_xs.total > 0.0) {
if (i_nuclide == -1 && p->macro_xs_.total > 0.0) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->E_last_)
&& simulation::micro_xs[i_nuc].total) {
&& p->neutron_xs_[i_nuc].total) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / simulation::material_xs.total;
} else if (simulation::micro_xs[i_nuclide].total) {
score *= flux_deriv + cum_dsig / p->macro_xs_.total;
} else if (p->neutron_xs_[i_nuclide].total) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv
+ (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total;
+ (dsig_s + dsig_a) / p->neutron_xs_[i_nuclide].total;
} else {
score *= flux_deriv;
}
break;
case SCORE_SCATTER:
if (i_nuclide == -1 && (simulation::material_xs.total
- simulation::material_xs.absorption)) {
if (i_nuclide == -1 && (p->macro_xs_.total
- p->macro_xs_.absorption)) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->E_last_)
&& (simulation::micro_xs[i_nuc].total
- simulation::micro_xs[i_nuc].absorption)) {
&& (p->neutron_xs_[i_nuc].total
- p->neutron_xs_[i_nuc].absorption)) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += dsig_s * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / (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->macro_xs_.total
- p->macro_xs_.absorption);
} else if (p->neutron_xs_[i_nuclide].total
- p->neutron_xs_[i_nuclide].absorption) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total
- simulation::micro_xs[i_nuclide].absorption);
score *= flux_deriv + dsig_s / (p->neutron_xs_[i_nuclide].total
- p->neutron_xs_[i_nuclide].absorption);
} else {
score *= flux_deriv;
}
break;
case SCORE_ABSORPTION:
if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
if (i_nuclide == -1 && p->macro_xs_.absorption > 0.0) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->E_last_)
&& simulation::micro_xs[i_nuc].absorption) {
&& p->neutron_xs_[i_nuc].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += dsig_a * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / simulation::material_xs.absorption;
} else if (simulation::micro_xs[i_nuclide].absorption) {
score *= flux_deriv + cum_dsig / p->macro_xs_.absorption;
} else if (p->neutron_xs_[i_nuclide].absorption) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv
+ dsig_a / simulation::micro_xs[i_nuclide].absorption;
+ dsig_a / p->neutron_xs_[i_nuclide].absorption;
} else {
score *= flux_deriv;
}
break;
case SCORE_FISSION:
if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
if (i_nuclide == -1 && p->macro_xs_.fission > 0.0) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->E_last_)
&& simulation::micro_xs[i_nuc].fission) {
&& p->neutron_xs_[i_nuc].fission) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += dsig_f * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / simulation::material_xs.fission;
} else if (simulation::micro_xs[i_nuclide].fission) {
score *= flux_deriv + cum_dsig / p->macro_xs_.fission;
} else if (p->neutron_xs_[i_nuclide].fission) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv
+ dsig_f / simulation::micro_xs[i_nuclide].fission;
+ dsig_f / p->neutron_xs_[i_nuclide].fission;
} else {
score *= flux_deriv;
}
break;
case SCORE_NU_FISSION:
if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
if (i_nuclide == -1 && p->macro_xs_.nu_fission > 0.0) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->E_last_)
&& simulation::micro_xs[i_nuc].fission) {
double nu = simulation::micro_xs[i_nuc].nu_fission
/ simulation::micro_xs[i_nuc].fission;
&& p->neutron_xs_[i_nuc].fission) {
double nu = p->neutron_xs_[i_nuc].nu_fission
/ p->neutron_xs_[i_nuc].fission;
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += nu * dsig_f * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / simulation::material_xs.nu_fission;
} else if (simulation::micro_xs[i_nuclide].fission) {
score *= flux_deriv + cum_dsig / p->macro_xs_.nu_fission;
} else if (p->neutron_xs_[i_nuclide].fission) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv
+ dsig_f / simulation::micro_xs[i_nuclide].fission;
+ dsig_f / p->neutron_xs_[i_nuclide].fission;
} else {
score *= flux_deriv;
}
@ -582,7 +582,7 @@ score_track_derivative(const Particle* p, double distance)
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
// (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
deriv.flux_deriv -= distance * simulation::material_xs.total
deriv.flux_deriv -= distance * p->macro_xs_.total
/ material.density_gpcc_;
break;
@ -591,7 +591,7 @@ score_track_derivative(const Particle* p, double distance)
// (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
// (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
deriv.flux_deriv -= distance
* simulation::micro_xs[deriv.diff_nuclide].total;
* p->neutron_xs_[deriv.diff_nuclide].total;
break;
case DIFF_TEMPERATURE:
@ -660,7 +660,7 @@ void score_collision_derivative(const Particle* p)
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
const auto& micro_xs {simulation::micro_xs[i_nuc]};
const auto& micro_xs {p->neutron_xs_[i_nuc]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);

3
src/tallies/filter_particle.cpp
View file

@ -50,6 +50,9 @@ ParticleFilter::text_label(int bin) const
case Particle::Type::positron:
return "Particle: positron";
}
#ifdef __GNUC__
__builtin_unreachable();
#endif
}
} // namespace openmc

190
src/tallies/tally_scoring.cpp
View file

@ -192,7 +192,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
// loop over number of particles banked
for (auto i = 0; i < p->n_bank_; ++i) {
auto i_bank = simulation::n_bank - p->n_bank_ + i;
auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i;
const auto& bank = simulation::fission_bank[i_bank];
// get the delayed group
@ -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->macro_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->neutron_xs_[i_nuclide].total * atom_density * flux;
} else {
score = simulation::material_xs.total * flux;
score = p->macro_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->macro_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->neutron_xs_[i_nuclide].total
- p->neutron_xs_[i_nuclide].absorption) * atom_density * flux;
} else {
score = (simulation::material_xs.total
- simulation::material_xs.absorption) * flux;
score = (p->macro_xs_.total
- p->macro_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->neutron_xs_[i_nuclide].absorption * atom_density
* flux;
} else {
score = simulation::material_xs.absorption * flux;
score = p->macro_xs_.absorption * flux;
}
}
break;
case SCORE_FISSION:
if (simulation::material_xs.absorption == 0) continue;
if (p->macro_xs_.absorption == 0) continue;
if (tally.estimator_ == ESTIMATOR_ANALOG) {
if (settings::survival_biasing) {
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// fission
if (simulation::micro_xs[p->event_nuclide_].absorption > 0) {
if (p->neutron_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->neutron_xs_[p->event_nuclide_].fission
/ p->neutron_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->neutron_xs_[p->event_nuclide_].fission
/ p->neutron_xs_[p->event_nuclide_].absorption * flux;
}
} else {
if (i_nuclide >= 0) {
score = simulation::micro_xs[i_nuclide].fission * atom_density * flux;
score = p->neutron_xs_[i_nuclide].fission * atom_density * flux;
} else {
score = simulation::material_xs.fission * flux;
score = p->macro_xs_.fission * flux;
}
}
break;
case SCORE_NU_FISSION:
if (simulation::material_xs.absorption == 0) continue;
if (p->macro_xs_.absorption == 0) continue;
if (tally.estimator_ == ESTIMATOR_ANALOG) {
if (settings::survival_biasing || p->fission_) {
if (tally.energyout_filter_ != C_NONE) {
@ -516,10 +516,10 @@ score_general_ce(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->neutron_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->neutron_xs_[p->event_nuclide_].nu_fission
/ p->neutron_xs_[p->event_nuclide_].absorption * flux;
} else {
score = 0.;
}
@ -535,17 +535,17 @@ score_general_ce(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->neutron_xs_[i_nuclide].nu_fission * atom_density
* flux;
} else {
score = simulation::material_xs.nu_fission * flux;
score = p->macro_xs_.nu_fission * flux;
}
}
break;
case SCORE_PROMPT_NU_FISSION:
if (simulation::material_xs.absorption == 0) continue;
if (p->macro_xs_.absorption == 0) continue;
if (tally.estimator_ == ESTIMATOR_ANALOG) {
if (settings::survival_biasing || p->fission_) {
if (tally.energyout_filter_ != C_NONE) {
@ -559,12 +559,12 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].absorption > 0) {
score = p->wgt_absorb_
* simulation::micro_xs[p->event_nuclide_].fission
* p->neutron_xs_[p->event_nuclide_].fission
* data::nuclides[p->event_nuclide_]
->nu(E, ReactionProduct::EmissionMode::prompt)
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
/ p->neutron_xs_[p->event_nuclide_].absorption * flux;
} else {
score = 0.;
}
@ -583,7 +583,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
}
} else {
if (i_nuclide >= 0) {
score = simulation::micro_xs[i_nuclide].fission
score = p->neutron_xs_[i_nuclide].fission
* data::nuclides[i_nuclide]
->nu(E, ReactionProduct::EmissionMode::prompt)
* atom_density * flux;
@ -595,7 +595,7 @@ score_general_ce(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->neutron_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->macro_xs_.absorption == 0) continue;
if (tally.estimator_ == ESTIMATOR_ANALOG) {
if (settings::survival_biasing || p->fission_) {
if (tally.energyout_filter_ != C_NONE) {
@ -621,7 +621,7 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].absorption > 0
&& data::nuclides[p->event_nuclide_]->fissionable_) {
if (tally.delayedgroup_filter_ != C_NONE) {
auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@ -634,8 +634,8 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].fission
/ p->neutron_xs_[p->event_nuclide_].absorption * flux;
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
@ -645,10 +645,10 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].fission
* data::nuclides[p->event_nuclide_]
->nu(E, ReactionProduct::EmissionMode::delayed)
/ simulation::micro_xs[p->event_nuclide_].absorption *flux;
/ p->neutron_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->neutron_xs_[i_nuclide].fission * yield
* atom_density * flux;
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
@ -702,7 +702,7 @@ score_general_ce(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->neutron_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->neutron_xs_[j_nuclide].fission * yield
* atom_density * flux;
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
@ -738,7 +738,7 @@ score_general_ce(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->neutron_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->macro_xs_.absorption == 0) continue;
if (tally.estimator_ == ESTIMATOR_ANALOG) {
if (settings::survival_biasing) {
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// delayed-nu-fission
const auto& nuc {*data::nuclides[p->event_nuclide_]};
if (simulation::micro_xs[p->event_nuclide_].absorption > 0
if (p->neutron_xs_[p->event_nuclide_].absorption > 0
&& nuc.fissionable_) {
const auto& rxn {*nuc.fission_rx_[0]};
if (tally.delayedgroup_filter_ != C_NONE) {
@ -773,8 +773,8 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].fission
/ p->neutron_xs_[p->event_nuclide_].absorption
* rate * flux;
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
@ -796,8 +796,8 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].fission * yield
/ p->neutron_xs_[p->event_nuclide_].absorption * flux;
}
}
}
@ -813,7 +813,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// contribution to the fission bank to the score.
score = 0.;
for (auto i = 0; i < p->n_bank_; ++i) {
auto i_bank = simulation::n_bank - p->n_bank_ + i;
auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i;
const auto& bank = simulation::fission_bank[i_bank];
auto g = bank.delayed_group;
if (g != 0) {
@ -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->neutron_xs_[i_nuclide].fission * yield * flux
* atom_density * rate;
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
@ -870,7 +870,7 @@ score_general_ce(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->neutron_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->neutron_xs_[j_nuclide].fission * yield
* flux * atom_density * rate;
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
@ -922,7 +922,7 @@ score_general_ce(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->neutron_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->macro_xs_.absorption == 0.) continue;
score = 0.;
// Kappa-fission values are determined from the Q-value listed for the
// fission cross section.
@ -945,12 +945,12 @@ score_general_ce(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->neutron_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->neutron_xs_[p->event_nuclide_].fission
/ p->neutron_xs_[p->event_nuclide_].absorption * flux;
}
} else {
// Skip any non-absorption events
@ -959,12 +959,12 @@ score_general_ce(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->neutron_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->neutron_xs_[p->event_nuclide_].fission
/ p->neutron_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->neutron_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->neutron_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->neutron_xs_[i_nuclide].elastic == CACHE_INVALID)
data::nuclides[i_nuclide]->calculate_elastic_xs(*p);
score = p->neutron_xs_[i_nuclide].elastic * atom_density * flux;
} else {
score = 0.;
if (p->material_ != MATERIAL_VOID) {
@ -1017,9 +1017,9 @@ score_general_ce(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->neutron_xs_[j_nuclide].elastic == CACHE_INVALID)
data::nuclides[j_nuclide]->calculate_elastic_xs(*p);
score += p->neutron_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->macro_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->neutron_xs_[p->event_nuclide_].absorption > 0) {
double q_value = 0.;
if (score_bin == SCORE_FISS_Q_PROMPT) {
if (nuc.fission_q_prompt_)
@ -1048,8 +1048,8 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].fission
/ p->neutron_xs_[p->event_nuclide_].absorption * flux;
}
} else {
// Skip any non-absorption events
@ -1058,7 +1058,7 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].absorption > 0) {
double q_value = 0.;
if (score_bin == SCORE_FISS_Q_PROMPT) {
if (nuc.fission_q_prompt_)
@ -1068,8 +1068,8 @@ score_general_ce(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->neutron_xs_[p->event_nuclide_].fission
/ p->neutron_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->neutron_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->neutron_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->neutron_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->neutron_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->neutron_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->neutron_xs_[i_nuclide].index_grid;
auto f = p->neutron_xs_[i_nuclide].interp_factor;
const auto& xs {rxn.xs_[i_temp]};
if (i_grid >= xs.threshold) {
score = ((1.0 - f) * xs.value[i_grid-xs.threshold]
@ -1184,10 +1184,10 @@ score_general_ce(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->neutron_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->neutron_xs_[j_nuclide].index_grid;
auto f = p->neutron_xs_[j_nuclide].interp_factor;
const auto& xs {rxn.xs_[i_temp]};
if (i_grid >= xs.threshold) {
score += ((1.0 - f) * xs.value[i_grid-xs.threshold]
@ -1291,7 +1291,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
} else {
score = p->wgt_last_;
}
score *= flux / simulation::material_xs.total;
score *= flux / p->macro_xs_.total;
} else {
score = flux;
}
@ -1320,7 +1320,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
score = get_nuclide_xs(i_nuclide, MG_GET_XS_TOTAL, p_g)
* atom_density * flux;
} else {
score = simulation::material_xs.total * flux;
score = p->macro_xs_.total * flux;
}
}
break;
@ -1438,7 +1438,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
score = atom_density * flux
* get_nuclide_xs(i_nuclide, MG_GET_XS_ABSORPTION, p_g);
} else {
score = simulation::material_xs.absorption * flux;
score = p->macro_xs_.absorption * flux;
}
}
break;
@ -1785,7 +1785,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// contribution to the fission bank to the score.
score = 0.;
for (auto i = 0; i < p->n_bank_; ++i) {
auto i_bank = simulation::n_bank - p->n_bank_ + i;
auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i;
const auto& bank = simulation::fission_bank[i_bank];
auto g = bank.delayed_group;
if (g != 0) {
@ -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->macro_xs_.total;
} else {
flux = (p->wgt_last_ + p->wgt_absorb_) / simulation::material_xs.total;
flux = (p->wgt_last_ + p->wgt_absorb_) / p->macro_xs_.total;
}
for (auto i_tally : model::active_collision_tallies) {

1
tests/unit_tests/test_settings.py
View file

@ -36,7 +36,6 @@ def test_export_to_xml(run_in_tmpdir):
s.tabular_legendre = {'enable': True, 'num_points': 50}
s.temperature = {'default': 293.6, 'method': 'interpolation',
'multipole': True, 'range': (200., 1000.)}
s.threads = 8
s.trace = (10, 1, 20)
s.track = [1, 1, 1, 2, 1, 1]
s.ufs_mesh = mesh