crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-25 12:35:29 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'develop' into add_enums2

Browse source
This commit is contained in:
Gavin Ridley 2020-01-21 13:13:50 -05:00 • committed by GitHub
commit 7df7fc280b
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
46 changed files with 1014 additions and 694 deletions
Download patch file Download diff file Expand all files Collapse all files

19
include/openmc/bank.h
View file

@ -7,11 +7,6 @@
#include "openmc/particle.h"
#include "openmc/position.h"
// Without an explicit instantiation of vector<Bank>, the Intel compiler
// will complain about the threadprivate directive on filter_matches. Note that
// this has to happen *outside* of the openmc namespace
extern template class std::vector<openmc::Particle::Bank>;
namespace openmc {
//==============================================================================
@ -21,13 +16,11 @@ namespace openmc {
namespace simulation {
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, secondary_bank)
extern std::unique_ptr<Particle::Bank[]> fission_bank;
extern int64_t fission_bank_length;
extern int64_t fission_bank_max;
extern std::vector<int64_t> progeny_per_particle;
} // namespace simulation
@ -35,8 +28,12 @@ extern std::vector<Particle::Bank> master_fission_bank;
// Non-member functions
//==============================================================================
void sort_fission_bank();
void free_memory_bank();
void init_fission_bank(int64_t max);
} // namespace openmc
#endif // OPENMC_BANK_H

6
include/openmc/cell.h
View file

@ -113,7 +113,7 @@ public:
//! Find the oncoming boundary of this cell.
virtual std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface) const = 0;
distance(Position r, Direction u, int32_t on_surface, Particle* p) const = 0;
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! \param group_id An HDF5 group id.
@ -204,7 +204,7 @@ public:
contains(Position r, Direction u, int32_t on_surface) const;
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface) const;
distance(Position r, Direction u, int32_t on_surface, Particle* p) const;
void to_hdf5(hid_t group_id) const;
@ -248,7 +248,7 @@ public:
bool contains(Position r, Direction u, int32_t on_surface) const;
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface) const;
distance(Position r, Direction u, int32_t on_surface, Particle* p) const;
BoundingBox bounding_box() const;

8
include/openmc/dagmc.h
View file

@ -14,14 +14,6 @@ extern "C" const bool dagmc_enabled;
namespace openmc {
namespace simulation {
extern moab::DagMC::RayHistory history; //!< facet history for DagMC particles
extern Direction last_dir; //!< last direction passed to DagMC's ray_fire
#pragma omp threadprivate(history, last_dir)
}
namespace model {
extern moab::DagMC* DAG;
}

8
include/openmc/eigenvalue.h
View file

@ -42,14 +42,6 @@ void calculate_generation_keff();
//! generations. It also broadcasts the value from the master process.
void calculate_average_keff();
#ifdef _OPENMP
//! Join threadprivate fission banks into a single fission bank
//!
//! Note that this operation is necessarily sequential to preserve the order of
//! the bank when using varying numbers of threads.
void join_bank_from_threads();
#endif
//! Calculates a minimum variance estimate of k-effective
//!
//! The minimum variance estimate is based on a linear combination of the

11
include/openmc/geometry.h
View file

@ -24,17 +24,6 @@ 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 two distances by coincidence tolerance
//==============================================================================

2
include/openmc/mesh.h
View file

@ -161,7 +161,7 @@ public:
//! \param[in] bank Array of bank sites
//! \param[out] Whether any bank sites are outside the mesh
//! \return Array indicating number of sites in each mesh/energy bin
xt::xtensor<double, 1> count_sites(const std::vector<Particle::Bank>& bank,
xt::xtensor<double, 1> count_sites(const Particle::Bank* bank, int64_t length,
bool* outside) const;
// Data members

79
include/openmc/particle.h
View file

@ -14,6 +14,12 @@
#include "openmc/constants.h"
#include "openmc/position.h"
#include "openmc/random_lcg.h"
#include "openmc/tallies/filter_match.h"
#ifdef DAGMC
#include "DagMC.hpp"
#endif
namespace openmc {
@ -131,6 +137,17 @@ struct MacroXS {
double pair_production; //!< macroscopic pair production xs
};
//==============================================================================
// Information about nearest boundary crossing
//==============================================================================
struct BoundaryInfo {
double distance {INFINITY}; //!< distance to nearest boundary
int surface_index {0}; //!< if boundary is surface, index in surfaces vector
int coord_level; //!< coordinate level after crossing boundary
std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
};
//============================================================================
//! State of a particle being transported through geometry
//============================================================================
@ -146,6 +163,12 @@ public:
};
//! Saved ("banked") state of a particle
//! NOTE: This structure's MPI type is built in initialize_mpi() of
//! initialize.cpp. Any changes made to the struct here must also be
//! made when building the Bank MPI type in initialize_mpi().
//! NOTE: This structure is also used on the python side, and is defined
//! in lib/core.py. Changes made to the type here must also be made to the
//! python defintion.
struct Bank {
Position r;
Direction u;
@ -153,6 +176,15 @@ public:
double wgt;
int delayed_group;
Type particle;
int64_t parent_id;
int64_t progeny_id;
};
//! Saved ("banked") state of a particle, for nu-fission tallying
struct NuBank {
double E; //!< particle energy
double wgt; //!< particle weight
int delayed_group; //!< particle delayed group
};
//==========================================================================
@ -199,8 +231,13 @@ public:
//! \param src Source site data
void from_source(const Bank* src);
//! Transport a particle from birth to death
void transport();
// Coarse-grained particle events
void event_calculate_xs();
void event_advance();
void event_cross_surface();
void event_collide();
void event_revive_from_secondary();
void event_death();
//! Cross a surface and handle boundary conditions
void cross_surface();
@ -276,10 +313,13 @@ public:
//!< sites banked
// Indices for various arrays
int surface_ {0}; //!< index for surface particle is on
int surface_ {0}; //!< index for surface particle is on
int cell_born_ {-1}; //!< index for cell particle was born in
int material_ {-1}; //!< index for current material
int material_last_ {-1}; //!< index for last material
// Boundary information
BoundaryInfo boundary_;
// Temperature of current cell
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
@ -294,6 +334,39 @@ public:
// Current PRNG state
uint64_t seeds_[N_STREAMS]; // current seeds
int stream_; // current RNG stream
// Secondary particle bank
std::vector<Particle::Bank> secondary_bank_;
int64_t current_work_; // current work index
std::vector<double> flux_derivs_; // for derivatives for this particle
std::vector<FilterMatch> filter_matches_; // tally filter matches
std::vector<std::vector<Position>> tracks_; // tracks for outputting to file
std::vector<NuBank> nu_bank_; // bank of most recently fissioned particles
// Global tally accumulators
double keff_tally_absorption_ {0.0};
double keff_tally_collision_ {0.0};
double keff_tally_tracklength_ {0.0};
double keff_tally_leakage_ {0.0};
bool trace_ {false}; //!< flag to show debug information
double collision_distance_; // distance to particle's next closest collision
int n_event_ {0}; // number of events executed in this particle's history
// DagMC state variables
#ifdef DAGMC
moab::DagMC::RayHistory history_;
Direction last_dir_;
#endif
int64_t n_progeny_ {0}; // Number of progeny produced by this particle
};
} // namespace openmc

3
include/openmc/physics.h
View file

@ -48,8 +48,7 @@ int sample_nuclide(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,
std::vector<Particle::Bank>& bank);
void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx);
int sample_element(Particle* p);

3
include/openmc/physics_mg.h
View file

@ -33,9 +33,8 @@ 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 The particle bank to populate
void
create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank);
create_fission_sites(Particle* p);
//! \brief Handles an absorption event
//! \param p Particle to operate on

17
include/openmc/simulation.h
View file

@ -24,7 +24,6 @@ namespace simulation {
extern "C" int current_batch; //!< current batch
extern "C" int current_gen; //!< current fission generation
extern "C" int64_t current_work; //!< index in source back of current particle
extern "C" bool initialized; //!< has simulation been initialized?
extern "C" double keff; //!< average k over batches
extern "C" double keff_std; //!< standard deviation of average k
@ -46,11 +45,6 @@ extern const RegularMesh* ufs_mesh;
extern std::vector<double> k_generation;
extern std::vector<int64_t> work_index;
// Threadprivate variables
extern "C" bool trace; //!< flag to show debug information
#pragma omp threadprivate(current_work, trace)
} // namespace simulation
//==============================================================================
@ -69,8 +63,12 @@ void initialize_batch();
//! Initialize a fission generation
void initialize_generation();
//! Full initialization of a particle history
void initialize_history(Particle* p, int64_t index_source);
//! Helper function for initialize_history() that is called independently elsewhere
void initialize_history_partial(Particle* p);
//! Finalize a batch
//!
//! Handles synchronization and accumulation of tallies, calculation of Shannon
@ -90,6 +88,13 @@ void broadcast_results();
void free_memory_simulation();
//! Simulate a single particle history (and all generated secondary particles,
//! if enabled), from birth to death
void transport_history_based_single_particle(Particle& p);
//! Simulate all particle histories using history-based parallelism
void transport_history_based();
} // namespace openmc
#endif // OPENMC_SIMULATION_H

6
include/openmc/surface.h
View file

@ -11,6 +11,7 @@
#include "pugixml.hpp"
#include "openmc/constants.h"
#include "openmc/particle.h"
#include "openmc/position.h"
#include "dagmc.h"
@ -112,8 +113,9 @@ public:
//! Determine the direction of a ray reflected from the surface.
//! \param[in] r The point at which the ray is incident.
//! \param[in] u Incident direction of the ray
//! \param[inout] p Pointer to the particle
//! \return Outgoing direction of the ray
virtual Direction reflect(Position r, Direction u) const;
virtual Direction reflect(Position r, Direction u, Particle* p) const;
virtual Direction diffuse_reflect(Position r, Direction u,
uint64_t* seed) const;
@ -170,7 +172,7 @@ public:
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
Direction reflect(Position r, Direction u) const;
Direction reflect(Position r, Direction u, Particle* p) const;
void to_hdf5(hid_t group_id) const;

13
include/openmc/tallies/derivative.h
View file

@ -27,7 +27,6 @@ struct TallyDerivative {
int id; //!< User-defined identifier
int diff_material; //!< Material this derivative is applied to
int diff_nuclide; //!< Nuclide this material is applied to
double flux_deriv; //!< Derivative of the current particle's weight
TallyDerivative() {}
explicit TallyDerivative(pugi::xml_node node);
@ -54,16 +53,13 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
//! further tallies are scored.
//
//! \param p The particle being tracked
void score_collision_derivative(const Particle* p);
void score_collision_derivative(Particle* p);
//! Adjust diff tally flux derivatives for a particle tracking event.
//
//! \param p The particle being tracked
//! \param distance The distance in [cm] traveled by the particle
void score_track_derivative(const Particle* p, double distance);
//! Set the flux derivatives on differential tallies to zero.
void zero_flux_derivs();
void score_track_derivative(Particle* p, double distance);
} // namespace openmc
@ -71,15 +67,10 @@ void zero_flux_derivs();
// Global variables
//==============================================================================
// Explicit vector template specialization declaration of threadprivate variable
// outside of the openmc namespace for the picky Intel compiler.
extern template class std::vector<openmc::TallyDerivative>;
namespace openmc {
namespace model {
extern std::vector<TallyDerivative> tally_derivs;
#pragma omp threadprivate(tally_derivs)
extern std::unordered_map<int, int> tally_deriv_map;
} // namespace model

30
include/openmc/tallies/filter.h
View file

@ -12,31 +12,10 @@
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/particle.h"
#include "openmc/tallies/filter_match.h"
#include "pugixml.hpp"
namespace openmc {
//==============================================================================
//! Stores bins and weights for filtered tally events.
//==============================================================================
class FilterMatch
{
public:
std::vector<int> bins_;
std::vector<double> weights_;
int i_bin_;
bool bins_present_ {false};
};
} // namespace openmc
// Without an explicit instantiation of vector<FilterMatch>, the Intel compiler
// will complain about the threadprivate directive on filter_matches. Note that
// this has to happen *outside* of the openmc namespace
extern template class std::vector<openmc::FilterMatch>;
namespace openmc {
//==============================================================================
@ -127,13 +106,6 @@ private:
// Global variables
//==============================================================================
namespace simulation {
extern std::vector<FilterMatch> filter_matches;
#pragma omp threadprivate(filter_matches)
} // namespace simulation
namespace model {
extern "C" int32_t n_filters;
extern std::vector<std::unique_ptr<Filter>> tally_filters;

21
include/openmc/tallies/filter_match.h Normal file
View file

@ -0,0 +1,21 @@
#ifndef OPENMC_TALLIES_FILTERMATCH_H
#define OPENMC_TALLIES_FILTERMATCH_H
namespace openmc {
//==============================================================================
//! Stores bins and weights for filtered tally events.
//==============================================================================
class FilterMatch
{
public:
std::vector<int> bins_;
std::vector<double> weights_;
int i_bin_;
bool bins_present_ {false};
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTERMATCH_H

7
include/openmc/tallies/tally.h
View file

@ -160,17 +160,10 @@ namespace simulation {
extern "C" int32_t n_realizations;
}
// It is possible to protect accumulate operations on global tallies by using an
// atomic update. However, when multiple threads accumulate to the same global
// tally, it can cause a higher cache miss rate due to invalidation. Thus, we
// use threadprivate variables to accumulate global tallies and then reduce at
// the end of a generation.
extern double global_tally_absorption;
extern double global_tally_collision;
extern double global_tally_tracklength;
extern double global_tally_leakage;
#pragma omp threadprivate(global_tally_absorption, global_tally_collision, \
global_tally_tracklength, global_tally_leakage)
//==============================================================================
// Non-member functions

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

@ -2,6 +2,7 @@
#define OPENMC_TALLIES_TALLY_SCORING_H
#include "openmc/particle.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
namespace openmc {
@ -22,12 +23,13 @@ class FilterBinIter
public:
//! Construct an iterator over bins that match a given particle's state.
FilterBinIter(const Tally& tally, const Particle* p);
FilterBinIter(const Tally& tally, Particle* p);
//! Construct an iterator over all filter bin combinations.
//
//! \param end if true, the returned iterator indicates the end of a loop.
FilterBinIter(const Tally& tally, bool end);
FilterBinIter(const Tally& tally, bool end,
std::vector<FilterMatch>* particle_filter_matches);
bool operator==(const FilterBinIter& other) const
{return index_ == other.index_;}
@ -39,6 +41,8 @@ public:
int index_ {1};
double weight_ {1.};
std::vector<FilterMatch>& filter_matches_;
private:
void compute_index_weight();
@ -73,7 +77,7 @@ void score_analog_tally_ce(Particle* p);
//! Analog tallies are triggered at every collision, not every event.
//
//! \param p The particle being tracked
void score_analog_tally_mg(const Particle* p);
void score_analog_tally_mg(Particle* p);
//! Score tallies using a tracklength estimate of the flux.
//
@ -89,7 +93,7 @@ void score_tracklength_tally(Particle* p, double distance);
//
//! \param p The particle being tracked
//! \param tallies A vector of tallies to score to
void score_surface_tally(const Particle* p, const std::vector<int>& tallies);
void score_surface_tally(Particle* p, const std::vector<int>& tallies);
} // namespace openmc

1
include/openmc/timer.h
View file

@ -21,6 +21,7 @@ extern Timer time_finalize;
extern Timer time_inactive;
extern Timer time_initialize;
extern Timer time_read_xs;
extern Timer time_sample_source;
extern Timer time_tallies;
extern Timer time_total;
extern Timer time_transport;

6
include/openmc/track_output.h
View file

@ -9,9 +9,9 @@ namespace openmc {
// Non-member functions
//==============================================================================
void add_particle_track();
void write_particle_track(const Particle& p);
void finalize_particle_track(const Particle& p);
void add_particle_track(Particle& p);
void write_particle_track(Particle& p);
void finalize_particle_track(Particle& p);
} // namespace openmc

12
openmc/data/reaction.py
View file

@ -1073,11 +1073,15 @@ class Reaction(EqualityMixin):
if i_reaction < ace.nxs[5] + 1:
# Check if angular distribution data exist
loc = int(ace.xss[ace.jxs[8] + i_reaction])
if loc <= 0:
# Angular distribution is either given as part of a product
# angle-energy distribution or is not given at all (in which
# case isotropic scattering is assumed)
if loc < 0:
# Angular distribution is given as part of a product
# angle-energy distribution
angle_dist = None
elif loc == 0:
# Angular distribution is isotropic
energy = [0.0, grid[-1]]
mu = Uniform(-1., 1.)
angle_dist = AngleDistribution(energy, [mu, mu])
else:
angle_dist = AngleDistribution.from_ace(ace, ace.jxs[9], loc)

40
openmc/geometry.py
View file

@ -77,15 +77,22 @@ class Geometry(object):
if universe.id in volume_calc.volumes:
universe.add_volume_information(volume_calc)
def export_to_xml(self, path='geometry.xml'):
def export_to_xml(self, path='geometry.xml', remove_surfs=False):
"""Export geometry to an XML file.
Parameters
----------
path : str
Path to file to write. Defaults to 'geometry.xml'.
remove_surfs : bool
Whether or not to remove redundant surfaces from the geometry when
exporting
"""
# Find and remove redundant surfaces from the geometry
if remove_surfs:
self.remove_redundant_surfaces()
# Create XML representation
root_element = ET.Element("geometry")
self.root_universe.create_xml_subelement(root_element, memo=set())
@ -382,6 +389,26 @@ class Geometry(object):
surfaces = cell.region.get_surfaces(surfaces)
return surfaces
def get_redundant_surfaces(self):
"""Return all of the topologically redundant surface ids
Returns
-------
dict
Dictionary whose keys are the ID of a redundant surface and whose
values are the topologically equivalent :class:`openmc.Surface`
that should replace it.
"""
tally = defaultdict(list)
for surf in self.get_all_surfaces().values():
coeffs = tuple(surf._coefficients[k] for k in surf._coeff_keys)
key = (surf._type,) + coeffs
tally[key].append(surf)
return {replace.id: keep
for keep, *redundant in tally.values()
for replace in redundant}
def get_materials_by_name(self, name, case_sensitive=False, matching=False):
"""Return a list of materials with matching names.
@ -579,6 +606,17 @@ class Geometry(object):
return sorted(lattices, key=lambda x: x.id)
def remove_redundant_surfaces(self):
"""Remove redundant surfaces from the geometry"""
# Get redundant surfaces
redundant_surfaces = self.get_redundant_surfaces()
# Iterate through all cells contained in the geometry
for cell in self.get_all_cells().values():
# Recursively remove redundant surfaces from regions
cell.region.remove_redundant_surfaces(redundant_surfaces)
def determine_paths(self, instances_only=False):
"""Determine paths through CSG tree for cells and materials.

4
openmc/lib/core.py
View file

@ -20,7 +20,9 @@ class _Bank(Structure):
('E', c_double),
('wgt', c_double),
('delayed_group', c_int),
('particle', c_int)]
('particle', c_int),
('parent_id', c_int64),
('progeny_id', c_int64)]
# Define input type for numpy arrays that will be passed into C++ functions

32
openmc/region.py
View file

@ -45,8 +45,7 @@ class Region(metaclass=ABCMeta):
return not self == other
def get_surfaces(self, surfaces=None):
"""
Recursively find all the surfaces referenced by a region and return them
"""Recursively find all surfaces referenced by a region and return them
Parameters
----------
@ -65,6 +64,19 @@ class Region(metaclass=ABCMeta):
surfaces = region.get_surfaces(surfaces)
return surfaces
def remove_redundant_surfaces(self, redundant_surfaces):
"""Recursively remove all redundant surfaces referenced by this region
Parameters
----------
redundant_surfaces : dict
Dictionary mapping redundant surface IDs to class:`openmc.Surface`
instances that should replace them.
"""
for region in self:
region.remove_redundant_surfaces(redundant_surfaces)
@staticmethod
def from_expression(expression, surfaces):
"""Generate a region given an infix expression.
@ -597,8 +609,7 @@ class Complement(Region):
return temp_region.bounding_box
def get_surfaces(self, surfaces=None):
"""
Recursively find and return all the surfaces referenced by the node
"""Recursively find and return all the surfaces referenced by the node
Parameters
----------
@ -617,6 +628,19 @@ class Complement(Region):
surfaces = region.get_surfaces(surfaces)
return surfaces
def remove_redundant_surfaces(self, redundant_surfaces):
"""Recursively remove all redundant surfaces referenced by this region
Parameters
----------
redundant_surfaces : dict
Dictionary mapping redundant surface IDs to class:`openmc.Surface`
instances that should replace them.
"""
for region in self.node:
region.remove_redundant_surfaces(redundant_surfaces)
def clone(self, memo=None):
"""Create a copy of this region - each of the surfaces in the
complement's node will be cloned and will have new unique IDs.

15
openmc/surface.py
View file

@ -2209,6 +2209,21 @@ class Halfspace(Region):
surfaces[self.surface.id] = self.surface
return surfaces
def remove_redundant_surfaces(self, redundant_surfaces):
"""Recursively remove all redundant surfaces referenced by this region
Parameters
----------
redundant_surfaces : dict
Dictionary mapping redundant surface IDs to surface IDs for the
:class:`openmc.Surface` instances that should replace them.
"""
surf = redundant_surfaces.get(self.surface.id)
if surf is not None:
self.surface = surf
def clone(self, memo=None):
"""Create a copy of this halfspace, with a cloned surface with a
unique ID.

100
src/bank.cpp
View file

@ -1,12 +1,11 @@
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/error.h"
#include "openmc/simulation.h"
#include "openmc/message_passing.h"
#include <cstdint>
// Explicit template instantiation definition
template class std::vector<openmc::Particle::Bank>;
namespace openmc {
@ -17,11 +16,21 @@ namespace openmc {
namespace simulation {
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
// The fission bank is allocated as a pointer, rather than a vector, as it will
// be shared by all threads in the simulation. It will be allocated to a fixed
// maximum capacity in the init_fission_bank() function. Then, Elements will be
// added to it by performing atomic incrementing captures on the fission_bank_length
// variable. A vector is not appropriate, as use of push_back() and size()
// functions would cause undefined or unintended behavior.
std::unique_ptr<Particle::Bank[]> fission_bank;
int64_t fission_bank_length {0};
int64_t fission_bank_max;
// Each entry in this vector corresponds to the number of progeny produced
// this generation for the particle located at that index. This vector is
// used to efficiently sort the fission bank after each iteration.
std::vector<int64_t> progeny_per_particle;
} // namespace simulation
@ -32,13 +41,70 @@ std::vector<Particle::Bank> master_fission_bank;
void free_memory_bank()
{
simulation::source_bank.clear();
#pragma omp parallel
{
simulation::fission_bank.clear();
simulation::fission_bank.reset();
simulation::fission_bank_length = 0;
simulation::progeny_per_particle.clear();
}
void init_fission_bank(int64_t max)
{
simulation::fission_bank_max = max;
simulation::fission_bank = std::make_unique<Particle::Bank[]>(max);
simulation::fission_bank_length = 0;
simulation::progeny_per_particle.resize(simulation::work_per_rank);
}
// Performs an O(n) sort on the fission bank, by leveraging
// the parent_id and progeny_id fields of banked particles. See the following
// paper for more details:
// "Reproducibility and Monte Carlo Eigenvalue Calculations," F.B. Brown and
// T.M. Sutton, 1992 ANS Annual Meeting, Transactions of the American Nuclear
// Society, Volume 65, Page 235.
void sort_fission_bank()
{
// Ensure we don't read off the end of the array if we ran with 0 particles
if (simulation::progeny_per_particle.size() == 0) {
return;
}
#ifdef _OPENMP
simulation::master_fission_bank.clear();
#endif
// Perform exclusive scan summation to determine starting indices in fission
// bank for each parent particle id
int64_t tmp = simulation::progeny_per_particle[0];
simulation::progeny_per_particle[0] = 0;
for (int64_t i = 1; i < simulation::progeny_per_particle.size(); i++) {
int64_t value = simulation::progeny_per_particle[i-1] + tmp;
tmp = simulation::progeny_per_particle[i];
simulation::progeny_per_particle[i] = value;
}
// TODO: C++17 introduces the exclusive_scan() function which could be
// used to replace everything above this point in this function.
// We need a scratch vector to make permutation of the fission bank into
// sorted order easy. Under normal usage conditions, the fission bank is
// over provisioned, so we can use that as scratch space.
Particle::Bank* sorted_bank;
std::vector<Particle::Bank> sorted_bank_holder;
// If there is not enough space, allocate a temporary vector and point to it
if (simulation::fission_bank_length > simulation::fission_bank_max / 2) {
sorted_bank_holder.resize(simulation::fission_bank_length);
sorted_bank = sorted_bank_holder.data();
} else { // otherwise, point sorted_bank to unused portion of the fission bank
sorted_bank = &simulation::fission_bank[simulation::fission_bank_length];
}
// Use parent and progeny indices to sort fission bank
for (int64_t i = 0; i < simulation::fission_bank_length; i++) {
const auto& site = simulation::fission_bank[i];
int64_t offset = site.parent_id - 1 - simulation::work_index[mpi::rank];
int64_t idx = simulation::progeny_per_particle[offset] + site.progeny_id;
sorted_bank[idx] = site;
}
// Copy sorted bank into the fission bank
std::copy(sorted_bank, sorted_bank + simulation::fission_bank_length,
simulation::fission_bank.get());
}
//==============================================================================
@ -59,12 +125,12 @@ extern "C" int openmc_source_bank(void** ptr, int64_t* n)
extern "C" int openmc_fission_bank(void** ptr, int64_t* n)
{
if (simulation::fission_bank.size() == 0) {
if (simulation::fission_bank_length == 0) {
set_errmsg("Fission bank has not been allocated.");
return OPENMC_E_ALLOCATE;
} else {
*ptr = simulation::fission_bank.data();
*n = simulation::fission_bank.size();
*ptr = simulation::fission_bank.get();
*n = simulation::fission_bank_length;
return 0;
}
}

12
src/cell.cpp
View file

@ -500,7 +500,7 @@ CSGCell::contains(Position r, Direction u, int32_t on_surface) const
//==============================================================================
std::pair<double, int32_t>
CSGCell::distance(Position r, Direction u, int32_t on_surface) const
CSGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
{
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
@ -786,13 +786,13 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
DAGCell::DAGCell() : Cell{} {};
std::pair<double, int32_t>
DAGCell::distance(Position r, Direction u, int32_t on_surface) const
DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
{
// if we've changed direction or we're not on a surface,
// reset the history and update last direction
if (u != simulation::last_dir || on_surface == 0) {
simulation::history.reset();
simulation::last_dir = u;
if (u != p->last_dir_ || on_surface == 0) {
p->history_.reset();
p->last_dir_ = u;
}
moab::ErrorCode rval;
@ -801,7 +801,7 @@ DAGCell::distance(Position r, Direction u, int32_t on_surface) const
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &simulation::history);
rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history_);
MB_CHK_ERR_CONT(rval);
int surf_idx;
if (hit_surf != 0) {

7
src/dagmc.cpp
View file

@ -39,13 +39,6 @@ namespace openmc {
const std::string DAGMC_FILENAME = "dagmc.h5m";
namespace simulation {
moab::DagMC::RayHistory history;
Direction last_dir;
}
namespace model {
moab::DagMC* DAG;

59
src/eigenvalue.cpp
View file

@ -20,10 +20,6 @@
#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
@ -87,7 +83,7 @@ void synchronize_bank()
#ifdef OPENMC_MPI
int64_t start = 0;
int64_t n_bank = simulation::fission_bank.size();
int64_t n_bank = simulation::fission_bank_length;
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
@ -95,14 +91,14 @@ void synchronize_bank()
// significant
if (mpi::rank == 0) start = 0;
int64_t finish = start + simulation::fission_bank.size();
int64_t finish = start + simulation::fission_bank_length;
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::fission_bank.size();
int64_t total = simulation::fission_bank.size();
int64_t finish = simulation::fission_bank_length;
int64_t total = finish;
#endif
// If there are not that many particles per generation, it's possible that no
@ -110,7 +106,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::fission_bank.empty()) {
if (simulation::fission_bank_length == 0) {
fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
}
@ -142,7 +138,9 @@ void synchronize_bank()
int64_t index_temp = 0;
std::vector<Particle::Bank> temp_sites(3*simulation::work_per_rank);
for (const auto& site : simulation::fission_bank) {
for (int64_t i = 0; i < simulation::fission_bank_length; i++ ) {
const auto& site = simulation::fission_bank[i];
// 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
@ -197,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::fission_bank.size() - sites_needed + i;
int i_bank = simulation::fission_bank_length - sites_needed + i;
temp_sites[index_temp] = simulation::fission_bank[i_bank];
++index_temp;
}
@ -351,40 +349,6 @@ void calculate_average_keff()
}
}
#ifdef _OPENMP
void join_bank_from_threads()
{
int n_threads = omp_get_max_threads();
#pragma omp parallel
{
// Copy thread fission bank sites to one shared copy
#pragma omp for ordered schedule(static)
for (int i = 0; i < n_threads; ++i) {
#pragma omp ordered
{
std::copy(
simulation::fission_bank.cbegin(),
simulation::fission_bank.cend(),
std::back_inserter(simulation::master_fission_bank)
);
}
}
// Make sure all threads have made it to this point
#pragma omp barrier
// Now copy the shared fission bank sites back to the master thread's copy.
if (omp_get_thread_num() == 0) {
simulation::fission_bank = simulation::master_fission_bank;
simulation::master_fission_bank.clear();
} else {
simulation::fission_bank.clear();
}
}
}
#endif
int openmc_get_keff(double* k_combined)
{
k_combined[0] = 0.0;
@ -542,7 +506,8 @@ void shannon_entropy()
// Get source weight in each mesh bin
bool sites_outside;
xt::xtensor<double, 1> p = simulation::entropy_mesh->count_sites(
simulation::fission_bank, &sites_outside);
simulation::fission_bank.get(), simulation::fission_bank_length,
&sites_outside);
// display warning message if there were sites outside entropy box
if (sites_outside) {
@ -581,7 +546,7 @@ void ufs_count_sites()
// count number of source sites in each ufs mesh cell
bool sites_outside;
simulation::source_frac = simulation::ufs_mesh->count_sites(
simulation::source_bank, &sites_outside);
simulation::source_bank.data(), simulation::source_bank.size(), &sites_outside);
// Check for sites outside of the mesh
if (mpi::master && sites_outside) {

6
src/geometry.cpp
View file

@ -119,7 +119,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
}
// Announce the cell that the particle is entering.
if (found && (settings::verbosity >= 10 || simulation::trace)) {
if (found && (settings::verbosity >= 10 || p->trace_)) {
std::stringstream msg;
msg << " Entering cell " << model::cells[i_cell]->id_;
write_message(msg, 1);
@ -297,7 +297,7 @@ cross_lattice(Particle* p, const BoundaryInfo& boundary)
auto& coord {p->coord_[p->n_coord_ - 1]};
auto& lat {*model::lattices[coord.lattice]};
if (settings::verbosity >= 10 || simulation::trace) {
if (settings::verbosity >= 10 || p->trace_) {
std::stringstream msg;
msg << " Crossing lattice " << lat.id_ << ". Current position ("
<< coord.lattice_x << "," << coord.lattice_y << ","
@ -370,7 +370,7 @@ BoundaryInfo distance_to_boundary(Particle* p)
Cell& c {*model::cells[coord.cell]};
// Find the oncoming surface in this cell and the distance to it.
auto surface_distance = c.distance(r, u, p->surface_);
auto surface_distance = c.distance(r, u, p->surface_, p);
d_surf = surface_distance.first;
level_surf_cross = surface_distance.second;

12
src/initialize.cpp
View file

@ -101,18 +101,20 @@ void initialize_mpi(MPI_Comm intracomm)
// Create bank datatype
Particle::Bank b;
MPI_Aint disp[6];
MPI_Aint disp[8];
MPI_Get_address(&b.r, &disp[0]);
MPI_Get_address(&b.u, &disp[1]);
MPI_Get_address(&b.E, &disp[2]);
MPI_Get_address(&b.wgt, &disp[3]);
MPI_Get_address(&b.delayed_group, &disp[4]);
MPI_Get_address(&b.particle, &disp[5]);
for (int i = 5; i >= 0; --i) disp[i] -= disp[0];
MPI_Get_address(&b.parent_id, &disp[6]);
MPI_Get_address(&b.progeny_id, &disp[7]);
for (int i = 7; i >= 0; --i) disp[i] -= disp[0];
int blocks[] {3, 3, 1, 1, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT};
MPI_Type_create_struct(6, blocks, disp, types, &mpi::bank);
int blocks[] {3, 3, 1, 1, 1, 1, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
MPI_Type_create_struct(8, blocks, disp, types, &mpi::bank);
MPI_Type_commit(&mpi::bank);
}
#endif // OPENMC_MPI

6
src/mesh.cpp
View file

@ -790,7 +790,7 @@ void RegularMesh::to_hdf5(hid_t group) const
}
xt::xtensor<double, 1>
RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
RegularMesh::count_sites(const Particle::Bank* bank, int64_t length,
bool* outside) const
{
// Determine shape of array for counts
@ -801,7 +801,9 @@ RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
xt::xarray<double> cnt {shape, 0.0};
bool outside_ = false;
for (const auto& site : bank) {
for (int64_t i = 0; i < length; i++) {
const auto& site = bank[i];
// determine scoring bin for entropy mesh
int mesh_bin = get_bin(site.r);

13
src/output.cpp
View file

@ -474,6 +474,8 @@ void print_runtime()
show_time("Time synchronizing fission bank", time_bank.elapsed(), 1);
show_time("Sampling source sites", time_bank_sample.elapsed(), 2);
show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2);
} else {
show_time("Time sampling source", time_sample_source.elapsed(), 1);
}
show_time("Time accumulating tallies", time_tallies.elapsed(), 1);
show_time("Total time for finalization", time_finalize.elapsed());
@ -678,9 +680,14 @@ write_tallies()
}
}
// Initialize Filter Matches Object
std::vector<FilterMatch> filter_matches;
// Allocate space for tally filter matches
filter_matches.resize(model::tally_filters.size());
// Loop over all filter bin combinations.
auto filter_iter = FilterBinIter(tally, false);
auto end = FilterBinIter(tally, true);
auto filter_iter = FilterBinIter(tally, false, &filter_matches);
auto end = FilterBinIter(tally, true, &filter_matches);
for (; filter_iter != end; ++filter_iter) {
auto filter_index = filter_iter.index_;
@ -691,7 +698,7 @@ write_tallies()
if (filter_index % tally.strides(i) == 0) {
auto i_filt = tally.filters(i);
const auto& filt {*model::tally_filters[i_filt]};
auto& match {simulation::filter_matches[i_filt]};
auto& match {filter_matches[i_filt]};
tallies_out << std::string(indent+1, ' ')
<< filt.text_label(match.i_bin_) << "\n";
}

507
src/particle.cpp
View file

@ -84,9 +84,9 @@ Particle::clear()
void
Particle::create_secondary(Direction u, double E, Type type)
{
simulation::secondary_bank.emplace_back();
secondary_bank_.emplace_back();
auto& bank {simulation::secondary_bank.back()};
auto& bank {secondary_bank_.back()};
bank.particle = type;
bank.wgt = wgt_;
bank.r = this->r();
@ -107,6 +107,7 @@ Particle::from_source(const Bank* src)
material_ = C_NONE;
n_collision_ = 0;
fission_ = false;
std::fill(flux_derivs_.begin(), flux_derivs_.end(), 0.0);
// Copy attributes from source bank site
type_ = src->particle;
@ -129,258 +130,244 @@ Particle::from_source(const Bank* src)
}
void
Particle::transport()
Particle::event_calculate_xs()
{
// Display message if high verbosity or trace is on
if (settings::verbosity >= 9 || simulation::trace) {
write_message("Simulating Particle " + std::to_string(id_));
// Set the random number stream
if (type_ == Particle::Type::neutron) {
stream_ = STREAM_TRACKING;
} else {
stream_ = STREAM_PHOTON;
}
// Initialize number of events to zero
int n_event = 0;
// Store pre-collision particle properties
wgt_last_ = wgt_;
E_last_ = E_;
u_last_ = this->u();
r_last_ = this->r();
// Add paricle's starting weight to count for normalizing tallies later
#pragma omp atomic
simulation::total_weight += wgt_;
// Reset event variables
event_ = EVENT_KILL;
event_nuclide_ = NUCLIDE_NONE;
event_mt_ = REACTION_NONE;
// If the cell hasn't been determined based on the particle's location,
// initiate a search for the current cell. This generally happens at the
// beginning of the history and again for any secondary particles
if (coord_[n_coord_ - 1].cell == C_NONE) {
if (!find_cell(this, false)) {
this->mark_as_lost("Could not find the cell containing particle "
+ std::to_string(id_));
return;
}
// Set birth cell attribute
if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell;
}
// Write particle track.
if (write_track_) write_particle_track(*this);
if (settings::check_overlaps) check_cell_overlap(this);
// Calculate microscopic and macroscopic cross sections
if (material_ != MATERIAL_VOID) {
if (settings::run_CE) {
if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) {
// If the material is the same as the last material and the
// temperature hasn't changed, we don't need to lookup cross
// sections again.
model::materials[material_]->calculate_xs(*this);
}
} else {
// Get the MG data; unlike the CE case above, we have to re-calculate
// cross sections for every collision since the cross sections may
// be angle-dependent
data::mg.macro_xs_[material_].calculate_xs(*this);
// Update the particle's group while we know we are multi-group
g_last_ = g_;
}
} else {
macro_xs_.total = 0.0;
macro_xs_.absorption = 0.0;
macro_xs_.fission = 0.0;
macro_xs_.nu_fission = 0.0;
}
}
void
Particle::event_advance()
{
// Find the distance to the nearest boundary
boundary_ = distance_to_boundary(this);
// Sample a distance to collision
if (type_ == Particle::Type::electron ||
type_ == Particle::Type::positron) {
collision_distance_ = 0.0;
} else if (macro_xs_.total == 0.0) {
collision_distance_ = INFINITY;
} else {
collision_distance_ = -std::log(prn(this->current_seed())) / macro_xs_.total;
}
// Select smaller of the two distances
double distance = std::min(boundary_.distance, collision_distance_);
// Advance particle
for (int j = 0; j < n_coord_; ++j) {
coord_[j].r += distance * coord_[j].u;
}
// Score track-length tallies
if (!model::active_tracklength_tallies.empty()) {
score_tracklength_tally(this, distance);
}
// Score track-length estimate of k-eff
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
type_ == Particle::Type::neutron) {
keff_tally_tracklength_ += wgt_ * distance * macro_xs_.nu_fission;
}
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) {
score_track_derivative(this, distance);
}
}
void
Particle::event_cross_surface()
{
// 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) {
cell_last_[j] = coord_[j].cell;
}
n_coord_last_ = n_coord_;
if (boundary_.lattice_translation[0] != 0 ||
boundary_.lattice_translation[1] != 0 ||
boundary_.lattice_translation[2] != 0) {
// Particle crosses lattice boundary
cross_lattice(this, boundary_);
event_ = EVENT_LATTICE;
} else {
// Particle crosses surface
this->cross_surface();
event_ = EVENT_SURFACE;
}
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
score_surface_tally(this, model::active_surface_tallies);
}
}
void
Particle::event_collide()
{
// Score collision estimate of keff
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
type_ == Particle::Type::neutron) {
keff_tally_collision_ += wgt_ * macro_xs_.nu_fission
/ macro_xs_.total;
}
// Score surface current tallies -- this has to be done before the collision
// since the direction of the particle will change and we need to use the
// pre-collision direction to figure out what mesh surfaces were crossed
if (!model::active_meshsurf_tallies.empty())
score_surface_tally(this, model::active_meshsurf_tallies);
// Clear surface component
surface_ = 0;
// Force calculation of cross-sections by setting last energy to zero
if (settings::run_CE) {
for (auto& micro : neutron_xs_) micro.last_E = 0.0;
collision(this);
} else {
collision_mg(this);
}
// Prepare to write out particle track.
if (write_track_) add_particle_track();
// Every particle starts with no accumulated flux derivative.
if (!model::active_tallies.empty()) zero_flux_derivs();
while (true) {
// Set the random number stream
if (type_ == Particle::Type::neutron) {
stream_ = STREAM_TRACKING;
// Score collision estimator tallies -- this is done after a collision
// has occurred rather than before because we need information on the
// outgoing energy for any tallies with an outgoing energy filter
if (!model::active_collision_tallies.empty()) score_collision_tally(this);
if (!model::active_analog_tallies.empty()) {
if (settings::run_CE) {
score_analog_tally_ce(this);
} else {
stream_ = STREAM_PHOTON;
}
// Store pre-collision particle properties
wgt_last_ = wgt_;
E_last_ = E_;
u_last_ = this->u();
r_last_ = this->r();
// Reset event variables
event_ = TallyEvent::KILL;
event_nuclide_ = NUCLIDE_NONE;
event_mt_ = REACTION_NONE;
// If the cell hasn't been determined based on the particle's location,
// initiate a search for the current cell. This generally happens at the
// beginning of the history and again for any secondary particles
if (coord_[n_coord_ - 1].cell == C_NONE) {
if (!find_cell(this, false)) {
this->mark_as_lost("Could not find the cell containing particle "
+ std::to_string(id_));
return;
}
// Set birth cell attribute
if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell;
}
// Write particle track.
if (write_track_) write_particle_track(*this);
if (settings::check_overlaps) check_cell_overlap(this);
// Calculate microscopic and macroscopic cross sections
if (material_ != MATERIAL_VOID) {
if (settings::run_CE) {
if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) {
// If the material is the same as the last material and the
// temperature hasn't changed, we don't need to lookup cross
// sections again.
model::materials[material_]->calculate_xs(*this);
}
} else {
// Get the MG data; unlike the CE case above, we have to re-calculate
// cross sections for every collision since the cross sections may
// be angle-dependent
data::mg.macro_xs_[material_].calculate_xs(*this);
// Update the particle's group while we know we are multi-group
g_last_ = g_;
}
} else {
macro_xs_.total = 0.0;
macro_xs_.absorption = 0.0;
macro_xs_.fission = 0.0;
macro_xs_.nu_fission = 0.0;
}
// Find the distance to the nearest boundary
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 (macro_xs_.total == 0.0) {
d_collision = INFINITY;
} else {
d_collision = -std::log(prn(this->current_seed())) / macro_xs_.total;
}
// Select smaller of the two distances
double distance = std::min(boundary.distance, d_collision);
// Advance particle
for (int j = 0; j < n_coord_; ++j) {
coord_[j].r += distance * coord_[j].u;
}
// Score track-length tallies
if (!model::active_tracklength_tallies.empty()) {
score_tracklength_tally(this, distance);
}
// Score track-length estimate of k-eff
if (settings::run_mode == RunMode::EIGENVALUE &&
type_ == Particle::Type::neutron) {
global_tally_tracklength += wgt_ * distance * macro_xs_.nu_fission;
}
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) {
score_track_derivative(this, distance);
}
if (d_collision > boundary.distance) {
// ====================================================================
// PARTICLE CROSSES SURFACE
// 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) {
cell_last_[j] = coord_[j].cell;
}
n_coord_last_ = n_coord_;
if (boundary.lattice_translation[0] != 0 ||
boundary.lattice_translation[1] != 0 ||
boundary.lattice_translation[2] != 0) {
// Particle crosses lattice boundary
cross_lattice(this, boundary);
event_ = TallyEvent::LATTICE;
} else {
// Particle crosses surface
this->cross_surface();
event_ = TallyEvent::SURFACE;
}
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
score_surface_tally(this, model::active_surface_tallies);
}
} else {
// ====================================================================
// PARTICLE HAS COLLISION
// Score collision estimate of keff
if (settings::run_mode == RunMode::EIGENVALUE &&
type_ == Particle::Type::neutron) {
global_tally_collision += wgt_ * macro_xs_.nu_fission
/ macro_xs_.total;
}
// Score surface current tallies -- this has to be done before the collision
// since the direction of the particle will change and we need to use the
// pre-collision direction to figure out what mesh surfaces were crossed
if (!model::active_meshsurf_tallies.empty())
score_surface_tally(this, model::active_meshsurf_tallies);
// Clear surface component
surface_ = 0;
if (settings::run_CE) {
collision(this);
} else {
collision_mg(this);
}
// Score collision estimator tallies -- this is done after a collision
// has occurred rather than before because we need information on the
// outgoing energy for any tallies with an outgoing energy filter
if (!model::active_collision_tallies.empty()) score_collision_tally(this);
if (!model::active_analog_tallies.empty()) {
if (settings::run_CE) {
score_analog_tally_ce(this);
} else {
score_analog_tally_mg(this);
}
}
// Reset banked weight during collision
n_bank_ = 0;
n_bank_second_ = 0;
wgt_bank_ = 0.0;
for (int& v : n_delayed_bank_) v = 0;
// Reset fission logical
fission_ = false;
// Save coordinates for tallying purposes
r_last_current_ = this->r();
// Set last material to none since cross sections will need to be
// re-evaluated
material_last_ = C_NONE;
// Set all directions to base level -- right now, after a collision, only
// the base level directions are changed
for (int j = 0; j < n_coord_ - 1; ++j) {
if (coord_[j + 1].rotated) {
// If next level is rotated, apply rotation matrix
const auto& m {model::cells[coord_[j].cell]->rotation_};
const auto& u {coord_[j].u};
coord_[j + 1].u = u.rotate(m);
} else {
// Otherwise, copy this level's direction
coord_[j+1].u = coord_[j].u;
}
}
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) score_collision_derivative(this);
}
// If particle has too many events, display warning and kill it
++n_event;
if (n_event == MAX_EVENTS) {
warning("Particle " + std::to_string(id_) +
" underwent maximum number of events.");
alive_ = false;
}
// Check for secondary particles if this particle is dead
if (!alive_) {
// If no secondary particles, break out of event loop
if (simulation::secondary_bank.empty()) break;
this->from_source(&simulation::secondary_bank.back());
simulation::secondary_bank.pop_back();
n_event = 0;
// Enter new particle in particle track file
if (write_track_) add_particle_track();
score_analog_tally_mg(this);
}
}
// Reset banked weight during collision
n_bank_ = 0;
n_bank_second_ = 0;
wgt_bank_ = 0.0;
for (int& v : n_delayed_bank_) v = 0;
// Reset fission logical
fission_ = false;
// Save coordinates for tallying purposes
r_last_current_ = this->r();
// Set last material to none since cross sections will need to be
// re-evaluated
material_last_ = C_NONE;
// Set all directions to base level -- right now, after a collision, only
// the base level directions are changed
for (int j = 0; j < n_coord_ - 1; ++j) {
if (coord_[j + 1].rotated) {
// If next level is rotated, apply rotation matrix
const auto& m {model::cells[coord_[j].cell]->rotation_};
const auto& u {coord_[j].u};
coord_[j + 1].u = u.rotate(m);
} else {
// Otherwise, copy this level's direction
coord_[j+1].u = coord_[j].u;
}
}
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) score_collision_derivative(this);
}
void
Particle::event_revive_from_secondary()
{
// If particle has too many events, display warning and kill it
++n_event_;
if (n_event_ == MAX_EVENTS) {
warning("Particle " + std::to_string(id_) +
" underwent maximum number of events.");
alive_ = false;
}
// Check for secondary particles if this particle is dead
if (!alive_) {
// If no secondary particles, break out of event loop
if (secondary_bank_.empty()) return;
this->from_source(&secondary_bank_.back());
secondary_bank_.pop_back();
n_event_ = 0;
// Enter new particle in particle track file
if (write_track_) add_particle_track(*this);
}
}
void
Particle::event_death()
{
#ifdef DAGMC
if (settings::dagmc) simulation::history.reset();
if (settings::dagmc) history_.reset();
#endif
// Finish particle track output.
@ -388,15 +375,39 @@ Particle::transport()
write_particle_track(*this);
finalize_particle_track(*this);
}
// Contribute tally reduction variables to global accumulator
#pragma omp atomic
global_tally_absorption += keff_tally_absorption_;
#pragma omp atomic
global_tally_collision += keff_tally_collision_;
#pragma omp atomic
global_tally_tracklength += keff_tally_tracklength_;
#pragma omp atomic
global_tally_leakage += keff_tally_leakage_;
// Reset particle tallies once accumulated
keff_tally_absorption_ = 0.0;
keff_tally_collision_ = 0.0;
keff_tally_tracklength_ = 0.0;
keff_tally_leakage_ = 0.0;
// Record the number of progeny created by this particle.
// This data will be used to efficiently sort the fission bank.
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
int64_t offset = id_ - 1 - simulation::work_index[mpi::rank];
simulation::progeny_per_particle[offset] = n_progeny_;
}
}
void
Particle::cross_surface()
{
int i_surface = std::abs(surface_);
// TODO: off-by-one
const auto& surf {model::surfaces[i_surface - 1].get()};
if (settings::verbosity >= 10 || simulation::trace) {
if (settings::verbosity >= 10 || trace_) {
write_message(" Crossing surface " + std::to_string(surf->id_));
}
@ -420,10 +431,10 @@ Particle::cross_surface()
}
// Score to global leakage tally
global_tally_leakage += wgt_;
keff_tally_leakage_ += wgt_;
// Display message
if (settings::verbosity >= 10 || simulation::trace) {
if (settings::verbosity >= 10 || trace_) {
write_message(" Leaked out of surface " + std::to_string(surf->id_));
}
return;
@ -460,7 +471,7 @@ Particle::cross_surface()
}
Direction u = (surf->bc_ == Surface::BoundaryType::REFLECT) ?
surf->reflect(this->r(), this->u()) :
surf->reflect(this->r(), this->u(), this) :
surf->diffuse_reflect(this->r(), this->u(), this->current_seed());
// Make sure new particle direction is normalized
@ -487,7 +498,7 @@ Particle::cross_surface()
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || simulation::trace) {
if (settings::verbosity >= 10 || trace_) {
write_message(" Reflected from surface " + std::to_string(surf->id_));
}
return;
@ -541,7 +552,7 @@ Particle::cross_surface()
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || simulation::trace) {
if (settings::verbosity >= 10 || trace_) {
write_message(" Hit periodic boundary on surface " +
std::to_string(surf->id_));
}
@ -667,7 +678,7 @@ Particle::write_restart() const
write_dataset(file_id, "id", id_);
write_dataset(file_id, "type", static_cast<int>(type_));
int64_t i = simulation::current_work;
int64_t i = current_work_;
write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt);
write_dataset(file_id, "energy", simulation::source_bank[i-1].E);
write_dataset(file_id, "xyz", simulation::source_bank[i-1].r);

21
src/particle_restart.cpp
View file

@ -10,7 +10,9 @@
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/tally.h"
#include "openmc/track_output.h"
#include <algorithm> // for copy
#include <array>
@ -100,8 +102,25 @@ void run_particle_restart()
}
init_particle_seeds(particle_seed, p.seeds_);
// Force calculation of cross-sections by setting last energy to zero
if (settings::run_CE) {
for (auto& micro : p.neutron_xs_) micro.last_E = 0.0;
}
// Prepare to write out particle track.
if (p.write_track_) add_particle_track(p);
// Every particle starts with no accumulated flux derivative.
if (!model::active_tallies.empty()) {
p.flux_derivs_.resize(model::tally_derivs.size(), 0.0);
std::fill(p.flux_derivs_.begin(), p.flux_derivs_.end(), 0.0);
}
// Allocate space for tally filter matches
p.filter_matches_.resize(model::tally_filters.size());
// Transport neutron
p.transport();
transport_history_based_single_particle(p);
// Write output if particle made it
print_particle(&p);

89
src/physics.cpp
View file

@ -60,7 +60,7 @@ void collision(Particle* p)
}
// Display information about collision
if (settings::verbosity >= 10 || simulation::trace) {
if (settings::verbosity >= 10 || p->trace_) {
std::stringstream msg;
if (p->event_ == TallyEvent::KILL) {
msg << " Killed. Energy = " << p->E_ << " eV.";
@ -95,14 +95,14 @@ void sample_neutron_reaction(Particle* p)
if (nuc->fissionable_) {
Reaction* rx = sample_fission(i_nuclide, p);
if (settings::run_mode == RunMode::EIGENVALUE) {
create_fission_sites(p, i_nuclide, rx, simulation::fission_bank);
create_fission_sites(p, i_nuclide, rx);
} else if (settings::run_mode == RunMode::FIXED_SOURCE &&
settings::create_fission_neutrons) {
create_fission_sites(p, i_nuclide, rx, simulation::secondary_bank);
create_fission_sites(p, i_nuclide, rx);
// Make sure particle population doesn't grow out of control for
// subcritical multiplication problems.
if (simulation::secondary_bank.size() >= 10000) {
if (p->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.");
@ -146,8 +146,7 @@ void sample_neutron_reaction(Particle* p)
}
void
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
std::vector<Particle::Bank>& bank)
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
{
// If uniform fission source weighting is turned on, we increase or decrease
// the expected number of fission sites produced
@ -169,28 +168,74 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
// group.
double nu_d[MAX_DELAYED_GROUPS] = {0.};
p->fission_ = true;
for (int i = 0; i < nu; ++i) {
// Create new bank site and get reference to last element
bank.emplace_back();
auto& site {bank.back()};
// Clear out particle's nu fission bank
p->nu_bank_.clear();
// Bank source neutrons by copying the particle data
site.r = p->r();
site.particle = Particle::Type::neutron;
site.wgt = 1. / weight;
p->fission_ = true;
int skipped = 0;
// Determine whether to place fission sites into the shared fission bank
// or the secondary particle bank.
bool use_fission_bank = (settings::run_mode == RUN_MODE_EIGENVALUE);
for (int i = 0; i < nu; ++i) {
Particle::Bank* site;
if (use_fission_bank) {
int64_t idx;
#pragma omp atomic capture
idx = simulation::fission_bank_length++;
if (idx >= simulation::fission_bank_max) {
warning("The shared fission bank is full. Additional fission sites created "
"in this generation will not be banked.");
#pragma omp atomic write
simulation::fission_bank_length = simulation::fission_bank_max;
skipped++;
break;
}
site = &simulation::fission_bank[idx];
} else {
// Create new bank site and get reference to last element
auto& bank = p->secondary_bank_;
bank.emplace_back();
site = &bank.back();
}
site->r = p->r();
site->particle = Particle::Type::neutron;
site->wgt = 1. / weight;
site->parent_id = p->id_;
site->progeny_id = p->n_progeny_++;
// Sample delayed group and angle/energy for fission reaction
sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->current_seed());
sample_fission_neutron(i_nuclide, rx, p->E_, site, p->current_seed());
// Set the delayed group on the particle as well
p->delayed_group_ = site.delayed_group;
p->delayed_group_ = site->delayed_group;
// Increment the number of neutrons born delayed
if (p->delayed_group_ > 0) {
nu_d[p->delayed_group_-1]++;
}
// Write fission particles to nuBank
if (use_fission_bank) {
p->nu_bank_.emplace_back();
Particle::NuBank* nu_bank_entry = &p->nu_bank_.back();
nu_bank_entry->wgt = site->wgt;
nu_bank_entry->E = site->E;
nu_bank_entry->delayed_group = site->delayed_group;
}
}
// If shared fission bank was full, and no fissions could be added,
// set the particle fission flag to false.
if (nu == skipped) {
p->fission_ = false;
return;
}
// If shared fission bank was full, but some fissions could be added,
// reduce nu accordingly
nu -= skipped;
// Store the total weight banked for analog fission tallies
p->n_bank_ = nu;
@ -548,7 +593,7 @@ void absorption(Particle* p, int i_nuclide)
// Score implicit absorption estimate of keff
if (settings::run_mode == RunMode::EIGENVALUE) {
global_tally_absorption += p->wgt_absorb_ * p->neutron_xs_[
p->keff_tally_absorption += p->wgt_absorb_ * p->neutron_xs_[
i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
}
} else {
@ -557,7 +602,7 @@ void absorption(Particle* p, int i_nuclide)
prn(p->current_seed()) * p->neutron_xs_[i_nuclide].total) {
// Score absorption estimate of keff
if (settings::run_mode == RunMode::EIGENVALUE) {
global_tally_absorption += p->wgt_ * p->neutron_xs_[
p->keff_tally_absorption_ += p->wgt_ * p->neutron_xs_[
i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
}
@ -696,12 +741,12 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
// Find speed of neutron in CM
vel = v_n.norm();
// Sample scattering angle, checking if it is an ncorrelated angle-energy
// distribution
// Sample scattering angle, checking if angle distribution is present (assume
// isotropic otherwise)
double mu_cm;
auto& d = rx.products_[0].distribution_[0];
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (d_) {
if (!d_->angle().empty()) {
mu_cm = d_->angle().sample(p->E_, p->current_seed());
} else {
mu_cm = 2.0*prn(p->current_seed()) - 1.0;

89
src/physics_mg.cpp
View file

@ -31,7 +31,7 @@ collision_mg(Particle* p)
sample_reaction(p);
// Display information about collision
if ((settings::verbosity >= 10) || (simulation::trace)) {
if ((settings::verbosity >= 10) || (p->trace_)) {
std::stringstream msg;
msg << " Energy Group = " << p->g_;
write_message(msg, 1);
@ -47,11 +47,10 @@ sample_reaction(Particle* p)
// absorption (including fission)
if (model::materials[p->material_]->fissionable_) {
if (settings::run_mode == RunMode::EIGENVALUE) {
create_fission_sites(p, simulation::fission_bank);
} else if ((settings::run_mode == RunMode::FIXED_SOURCE) &&
(settings::create_fission_neutrons)) {
create_fission_sites(p, simulation::secondary_bank);
if (settings::run_mode == RunMode::EIGENVALUE ||
(settings::run_mode == RunMode::FIXED_SOURCE &&
settings::create_fission_neutrons)) {
create_fission_sites(p);
}
}
@ -91,7 +90,7 @@ scatter(Particle* p)
}
void
create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
create_fission_sites(Particle* p)
{
// If uniform fission source weighting is turned on, we increase or decrease
// the expected number of fission sites produced
@ -114,17 +113,45 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
// Initialize the counter of delayed neutrons encountered for each delayed
// group.
double nu_d[MAX_DELAYED_GROUPS] = {0.};
// Clear out particle's nu fission bank
p->nu_bank_.clear();
p->fission_ = true;
int skipped = 0;
// Determine whether to place fission sites into the shared fission bank
// or the secondary particle bank.
bool use_fission_bank = (settings::run_mode == RUN_MODE_EIGENVALUE);
for (int i = 0; i < nu; ++i) {
// Create new bank site and get reference to last element
bank.emplace_back();
auto& site {bank.back()};
Particle::Bank* site;
if (use_fission_bank) {
int64_t idx;
#pragma omp atomic capture
idx = simulation::fission_bank_length++;
if (idx >= simulation::fission_bank_max) {
warning("The shared fission bank is full. Additional fission sites created "
"in this generation will not be banked.");
#pragma omp atomic write
simulation::fission_bank_length = simulation::fission_bank_max;
skipped++;
break;
}
site = &simulation::fission_bank[idx];
} else {
// Create new bank site and get reference to last element
auto& bank = p->secondary_bank_;
bank.emplace_back();
site = &bank.back();
}
// Bank source neutrons by copying the particle data
site.r = p->r();
site.particle = Particle::Type::neutron;
site.wgt = 1. / weight;
site->r = p->r();
site->particle = Particle::Type::neutron;
site->wgt = 1. / weight;
site->parent_id = p->id_;
site->progeny_id = p->n_progeny_++;
// Sample the cosine of the angle, assuming fission neutrons are emitted
// isotropically
@ -132,9 +159,9 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
// Sample the azimuthal angle uniformly in [0, 2.pi)
double phi = 2. * PI * prn(p->current_seed() );
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);
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
int dg;
@ -142,10 +169,10 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout,
p->current_seed());
// Store the energy and delayed groups on the fission bank
site.E = gout;
site->E = gout;
// We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest
// of the code, 0 is prompt.
site.delayed_group = dg + 1;
site->delayed_group = dg + 1;
// Set the delayed group on the particle as well
p->delayed_group_ = dg + 1;
@ -154,7 +181,27 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
if (p->delayed_group_ > 0) {
nu_d[dg]++;
}
// Write fission particles to nuBank
if (use_fission_bank) {
p->nu_bank_.emplace_back();
Particle::NuBank* nu_bank_entry = &p->nu_bank_.back();
nu_bank_entry->wgt = site->wgt;
nu_bank_entry->E = site->E;
nu_bank_entry->delayed_group = site->delayed_group;
}
}
// If shared fission bank was full, and no fissions could be added,
// set the particle fission flag to false.
if (nu == skipped) {
p->fission_ = false;
return;
}
// If shared fission bank was full, but some fissions could be added,
// reduce nu accordingly
nu -= skipped;
// Store the total weight banked for analog fission tallies
p->n_bank_ = nu;
@ -176,13 +223,11 @@ absorption(Particle* p)
p->wgt_last_ = p->wgt_;
// Score implicit absorpion estimate of keff
#pragma omp atomic
global_tally_absorption += p->wgt_absorb_ * p->macro_xs_.nu_fission /
p->keff_tally_absorption_ += p->wgt_absorb_ * p->macro_xs_.nu_fission /
p->macro_xs_.absorption;
} else {
if (p->macro_xs_.absorption > prn(p->current_seed()) * p->macro_xs_.total) {
#pragma omp atomic
global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission /
p->keff_tally_absorption_ += p->wgt_ * p->macro_xs_.nu_fission /
p->macro_xs_.absorption;
p->alive_ = false;
p->event_ = TallyEvent::ABSORB;

176
src/simulation.cpp
View file

@ -17,18 +17,25 @@
#include "openmc/source.h"
#include "openmc/state_point.h"
#include "openmc/timer.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/trigger.h"
#include "openmc/track_output.h"
#ifdef _OPENMP
#include <omp.h>
#endif
#include "xtensor/xview.hpp"
#ifdef OPENMC_MPI
#include <mpi.h>
#endif
#include <algorithm>
#include <string>
//==============================================================================
// C API functions
//==============================================================================
@ -61,12 +68,6 @@ int openmc_simulation_init()
// Determine how much work each process should do
calculate_work();
// Allocate array for matching filter bins
#pragma omp parallel
{
simulation::filter_matches.resize(model::tally_filters.size());
}
// Allocate source bank, and for eigenvalue simulations also allocate the
// fission bank
allocate_banks();
@ -139,11 +140,6 @@ int openmc_simulation_finalize()
// Write tally results to tallies.out
if (settings::output_tallies && mpi::master) write_tallies();
#pragma omp parallel
{
simulation::filter_matches.clear();
}
// Deactivate all tallies
for (auto& t : model::tallies) {
t->active_ = false;
@ -186,20 +182,8 @@ int openmc_next_batch(int* status)
// Start timer for transport
simulation::time_transport.start();
// ====================================================================
// LOOP OVER PARTICLES
#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
Particle p;
initialize_history(&p, simulation::current_work);
// transport particle
p.transport();
}
// Transport loop
transport_history_based();
// Accumulate time for transport
simulation::time_transport.stop();
@ -242,7 +226,6 @@ namespace simulation {
int current_batch;
int current_gen;
int64_t current_work;
bool initialized {false};
double keff {1.0};
double keff_std;
@ -264,8 +247,6 @@ const RegularMesh* ufs_mesh {nullptr};
std::vector<double> k_generation;
std::vector<int64_t> work_index;
// Threadprivate variables
bool trace; //!< flag to show debug information
} // namespace simulation
@ -279,26 +260,9 @@ void allocate_banks()
simulation::source_bank.resize(simulation::work_per_rank);
if (settings::run_mode == RunMode::EIGENVALUE) {
#ifdef _OPENMP
// If OpenMP is being used, each thread needs its own private fission
// bank. Since the private fission banks need to be combined at the end of
// a generation, there is also a 'master_fission_bank' that is used to
// collect the sites from each thread.
#pragma omp parallel
{
if (omp_get_thread_num() == 0) {
simulation::fission_bank.reserve(3*simulation::work_per_rank);
} else {
int n_threads = omp_get_num_threads();
simulation::fission_bank.reserve(3*simulation::work_per_rank / n_threads);
}
}
simulation::master_fission_bank.reserve(3*simulation::work_per_rank);
#else
simulation::fission_bank.reserve(3*simulation::work_per_rank);
#endif
init_fission_bank(3*simulation::work_per_rank);
}
}
void initialize_batch()
@ -399,7 +363,7 @@ void initialize_generation()
{
if (settings::run_mode == RunMode::EIGENVALUE) {
// Clear out the fission bank
simulation::fission_bank.clear();
simulation::fission_bank_length = 0;
// Count source sites if using uniform fission source weighting
if (settings::ufs_on) ufs_count_sites();
@ -414,34 +378,27 @@ void finalize_generation()
{
auto& gt = simulation::global_tallies;
// Update global tallies with the omp private accumulation variables
#pragma omp parallel
{
#pragma omp critical(increment_global_tallies)
{
if (settings::run_mode == RunMode::EIGENVALUE) {
gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision;
gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += global_tally_absorption;
gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += global_tally_tracklength;
}
gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage;
}
// reset threadprivate tallies
if (settings::run_mode == RunMode::EIGENVALUE) {
global_tally_collision = 0.0;
global_tally_absorption = 0.0;
global_tally_tracklength = 0.0;
}
global_tally_leakage = 0.0;
}
// Update global tallies with the accumulation variables
if (settings::run_mode == RunMode::EIGENVALUE) {
#ifdef _OPENMP
// Join the fission bank from each thread into one global fission bank
join_bank_from_threads();
#endif
gt(K_COLLISION, RESULT_VALUE) += global_tally_collision;
gt(K_ABSORPTION, RESULT_VALUE) += global_tally_absorption;
gt(K_TRACKLENGTH, RESULT_VALUE) += global_tally_tracklength;
}
gt(LEAKAGE, RESULT_VALUE) += global_tally_leakage;
// reset tallies
if (settings::run_mode == RunMode::EIGENVALUE) {
global_tally_collision = 0.0;
global_tally_absorption = 0.0;
global_tally_tracklength = 0.0;
}
global_tally_leakage = 0.0;
if (settings::run_mode == RunMode::EIGENVALUE) {
// If using shared memory, stable sort the fission bank (by parent IDs)
// so as to allow for reproducibility regardless of which order particles
// are run in.
sort_fission_bank();
// Distribute fission bank across processors evenly
synchronize_bank();
@ -462,7 +419,9 @@ void finalize_generation()
} else if (settings::run_mode == RunMode::FIXED_SOURCE) {
// For fixed-source mode, we need to sample the external source
simulation::time_sample_source.start();
fill_source_bank_fixedsource();
simulation::time_sample_source.stop();
}
}
@ -470,20 +429,24 @@ void initialize_history(Particle* p, int64_t index_source)
{
// set defaults
p->from_source(&simulation::source_bank[index_source - 1]);
p->current_work_ = index_source;
// set identifier for particle
p->id_ = simulation::work_index[mpi::rank] + index_source;
// set progeny count to zero
p->n_progeny_ = 0;
// set random number seed
int64_t particle_seed = (simulation::total_gen + overall_generation() - 1)
* settings::n_particles + p->id_;
init_particle_seeds(particle_seed, p->seeds_);
// set particle trace
simulation::trace = false;
p->trace_ = false;
if (simulation::current_batch == settings::trace_batch &&
simulation::current_gen == settings::trace_gen &&
p->id_ == settings::trace_particle) simulation::trace = true;
p->id_ == settings::trace_particle) p->trace_ = true;
// Set particle track.
p->write_track_ = false;
@ -499,6 +462,38 @@ void initialize_history(Particle* p, int64_t index_source)
}
}
}
// Display message if high verbosity or trace is on
if (settings::verbosity >= 9 || p->trace_) {
write_message("Simulating Particle " + std::to_string(p->id_));
}
// Add paricle's starting weight to count for normalizing tallies later
#pragma omp atomic
simulation::total_weight += p->wgt_;
initialize_history_partial(p);
}
void initialize_history_partial(Particle* p)
{
// Force calculation of cross-sections by setting last energy to zero
if (settings::run_CE) {
for (auto& micro : p->neutron_xs_) micro.last_E = 0.0;
}
// Prepare to write out particle track.
if (p->write_track_) add_particle_track(*p);
// Every particle starts with no accumulated flux derivative.
if (!model::active_tallies.empty())
{
p->flux_derivs_.resize(model::tally_derivs.size(), 0.0);
std::fill(p->flux_derivs_.begin(), p->flux_derivs_.end(), 0.0);
}
// Allocate space for tally filter matches
p->filter_matches_.resize(model::tally_filters.size());
}
int overall_generation()
@ -571,4 +566,31 @@ void free_memory_simulation()
simulation::entropy.clear();
}
void transport_history_based_single_particle(Particle& p)
{
while (true) {
p.event_calculate_xs();
p.event_advance();
if (p.collision_distance_ > p.boundary_.distance) {
p.event_cross_surface();
} else {
p.event_collide();
}
p.event_revive_from_secondary();
if (!p.alive_)
break;
}
p.event_death();
}
void transport_history_based()
{
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
Particle p;
initialize_history(&p, i_work);
transport_history_based_single_particle(p);
}
}
} // namespace openmc

1
src/source.cpp
View file

@ -328,6 +328,7 @@ void free_memory_source()
void fill_source_bank_fixedsource()
{
if (settings::path_source.empty()) {
#pragma omp parallel for
for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
// initialize random number seed
int64_t id = (simulation::total_gen + overall_generation()) *

2
src/state_point.cpp
View file

@ -283,6 +283,8 @@ openmc_statepoint_write(const char* filename, bool* write_source)
write_dataset(runtime_group, "synchronizing fission bank", time_bank.elapsed());
write_dataset(runtime_group, "sampling source sites", time_bank_sample.elapsed());
write_dataset(runtime_group, "SEND-RECV source sites", time_bank_sendrecv.elapsed());
} else {
write_dataset(runtime_group, "sampling source sites", time_sample_source.elapsed());
}
write_dataset(runtime_group, "accumulating tallies", time_tallies.elapsed());
write_dataset(runtime_group, "total", time_total.elapsed());

10
src/surface.cpp
View file

@ -175,7 +175,7 @@ Surface::sense(Position r, Direction u) const
}
Direction
Surface::reflect(Position r, Direction u) const
Surface::reflect(Position r, Direction u, Particle* p) const
{
// Determine projection of direction onto normal and squared magnitude of
// normal.
@ -283,11 +283,11 @@ Direction DAGSurface::normal(Position r) const
return dir;
}
Direction DAGSurface::reflect(Position r, Direction u) const
Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const
{
simulation::history.reset_to_last_intersection();
simulation::last_dir = Surface::reflect(r, u);
return simulation::last_dir;
p->history_.reset_to_last_intersection();
p->last_dir_ = Surface::reflect(r, u, p);
return p->last_dir_;
}
void DAGSurface::to_hdf5(hid_t group_id) const {}

44
src/tallies/derivative.cpp
View file

@ -81,13 +81,9 @@ TallyDerivative::TallyDerivative(pugi::xml_node node)
void
read_tally_derivatives(pugi::xml_node node)
{
// Populate the derivatives array. This must be done in parallel because
// the derivatives are threadprivate.
#pragma omp parallel
{
for (auto deriv_node : node.children("derivative"))
model::tally_derivs.emplace_back(deriv_node);
}
// Populate the derivatives array.
for (auto deriv_node : node.children("derivative"))
model::tally_derivs.emplace_back(deriv_node);
// Fill the derivative map.
for (auto i = 0; i < model::tally_derivs.size(); ++i) {
@ -120,7 +116,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
// perturbated variable.
const auto& deriv {model::tally_derivs[tally.deriv_]};
auto flux_deriv = deriv.flux_deriv;
const auto flux_deriv = p->flux_derivs_[tally.deriv_];
// Handle special cases where we know that d_c/d_p must be zero.
if (score_bin == SCORE_FLUX) {
@ -567,13 +563,15 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
}
void
score_track_derivative(const Particle* p, double distance)
score_track_derivative(Particle* p, double distance)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material_ == MATERIAL_VOID) return;
const Material& material {*model::materials[p->material_]};
for (auto& deriv : model::tally_derivs) {
for (auto idx = 0; idx < model::tally_derivs.size(); idx++) {
const auto& deriv = model::tally_derivs[idx];
auto& flux_deriv = p->flux_derivs_[idx];
if (deriv.diff_material != material.id_) continue;
switch (deriv.variable) {
@ -582,7 +580,7 @@ score_track_derivative(const Particle* p, double distance)
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
// (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
deriv.flux_deriv -= distance * p->macro_xs_.total
flux_deriv -= distance * p->macro_xs_.total
/ material.density_gpcc_;
break;
@ -590,7 +588,7 @@ score_track_derivative(const Particle* p, double distance)
// phi is proportional to e^(-Sigma_tot * dist)
// (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
flux_deriv -= distance
* p->neutron_xs_[deriv.diff_nuclide].total;
break;
@ -604,7 +602,7 @@ score_track_derivative(const Particle* p, double distance)
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_, p->sqrtkT_);
deriv.flux_deriv -= distance * (dsig_s + dsig_a)
flux_deriv -= distance * (dsig_s + dsig_a)
* material.atom_density_(i);
}
}
@ -613,14 +611,17 @@ score_track_derivative(const Particle* p, double distance)
}
}
void score_collision_derivative(const Particle* p)
void score_collision_derivative(Particle* p)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material_ == MATERIAL_VOID) return;
const Material& material {*model::materials[p->material_]};
for (auto& deriv : model::tally_derivs) {
for (auto idx = 0; idx < model::tally_derivs.size(); idx++) {
const auto& deriv = model::tally_derivs[idx];
auto& flux_deriv = p->flux_derivs_[idx];
if (deriv.diff_material != material.id_) continue;
switch (deriv.variable) {
@ -629,7 +630,7 @@ void score_collision_derivative(const Particle* p)
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s
// (1 / phi) * (d_phi / d_rho) = 1 / rho
deriv.flux_deriv += 1. / material.density_gpcc_;
flux_deriv += 1. / material.density_gpcc_;
break;
case DerivativeVariable::NUCLIDE_DENSITY:
@ -650,7 +651,7 @@ void score_collision_derivative(const Particle* p)
// (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s
// (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s
// (1 / phi) * (d_phi / d_N) = 1 / N
deriv.flux_deriv += 1. / material.atom_density_(i);
flux_deriv += 1. / material.atom_density_(i);
break;
case DerivativeVariable::TEMPERATURE:
@ -665,7 +666,7 @@ void score_collision_derivative(const Particle* p)
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
// Note that this is an approximation! The real scattering cross
// section is
// Sigma_s(E'->E, u'->u) = Sigma_s(E') * P(E'->E, u'->u).
@ -680,9 +681,4 @@ void score_collision_derivative(const Particle* p)
}
}
void zero_flux_derivs()
{
for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.;
}
}// namespace openmc

4
src/tallies/filter.cpp
View file

@ -39,10 +39,6 @@ namespace openmc {
// Global variables
//==============================================================================
namespace simulation {
std::vector<FilterMatch> filter_matches;
}
namespace model {
std::vector<std::unique_ptr<Filter>> tally_filters;
std::unordered_map<int, int> filter_map;

5
src/tallies/tally.cpp
View file

@ -1044,10 +1044,7 @@ setup_active_tallies()
void
free_memory_tally()
{
#pragma omp parallel
{
model::tally_derivs.clear();
}
model::tally_derivs.clear();
model::tally_deriv_map.clear();
model::tally_filters.clear();

111
src/tallies/tally_scoring.cpp
View file

@ -26,13 +26,13 @@ namespace openmc {
// FilterBinIter implementation
//==============================================================================
FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p)
: tally_{tally}
FilterBinIter::FilterBinIter(const Tally& tally, Particle* p)
: tally_{tally}, filter_matches_{p->filter_matches_}
{
// Find all valid bins in each relevant filter if they have not already been
// found for this event.
for (auto i_filt : tally_.filters()) {
auto& match {simulation::filter_matches[i_filt]};
auto& match {filter_matches_[i_filt]};
if (!match.bins_present_) {
match.bins_.clear();
match.weights_.clear();
@ -55,8 +55,9 @@ FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p)
this->compute_index_weight();
}
FilterBinIter::FilterBinIter(const Tally& tally, bool end)
: tally_{tally}
FilterBinIter::FilterBinIter(const Tally& tally, bool end,
std::vector<FilterMatch>* particle_filter_matches)
: tally_{tally}, filter_matches_{*particle_filter_matches}
{
// Handle the special case for an iterator that points to the end.
if (end) {
@ -65,7 +66,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, bool end)
}
for (auto i_filt : tally_.filters()) {
auto& match {simulation::filter_matches[i_filt]};
auto& match {filter_matches_[i_filt]};
if (!match.bins_present_) {
match.bins_.clear();
match.weights_.clear();
@ -97,7 +98,7 @@ FilterBinIter::operator++()
bool done_looping = true;
for (int i = tally_.filters().size()-1; i >= 0; --i) {
auto i_filt = tally_.filters(i);
auto& match {simulation::filter_matches[i_filt]};
auto& match {filter_matches_[i_filt]};
if (match.i_bin_ < match.bins_.size()-1) {
// The bin for this filter can be incremented. Increment it and do not
// touch any of the remaining filters.
@ -130,7 +131,7 @@ FilterBinIter::compute_index_weight()
weight_ = 1.;
for (auto i = 0; i < tally_.filters().size(); ++i) {
auto i_filt = tally_.filters(i);
auto& match {simulation::filter_matches[i_filt]};
auto& match {filter_matches_[i_filt]};
auto i_bin = match.i_bin_;
index_ += match.bins_[i_bin] * tally_.strides(i);
weight_ *= match.weights_[i_bin];
@ -144,12 +145,13 @@ FilterBinIter::compute_index_weight()
//! Helper function used to increment tallies with a delayed group filter.
void
score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index)
score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index,
std::vector<FilterMatch>& filter_matches)
{
// Save the original delayed group bin
auto& tally {*model::tallies[i_tally]};
auto i_filt = tally.filters(tally.delayedgroup_filter_);
auto& dg_match {simulation::filter_matches[i_filt]};
auto& dg_match {filter_matches[i_filt]};
auto i_bin = dg_match.i_bin_;
auto original_bin = dg_match.bins_[i_bin];
dg_match.bins_[i_bin] = d_bin;
@ -158,7 +160,7 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index)
auto filter_index = 0;
for (auto i = 0; i < tally.filters().size(); ++i) {
auto i_filt = tally.filters(i);
auto& match {simulation::filter_matches[i_filt]};
auto& match {filter_matches[i_filt]};
auto i_bin = match.i_bin_;
filter_index += match.bins_[i_bin] * tally.strides(i);
}
@ -314,12 +316,12 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux,
//! neutrons produced with different energies.
void
score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin)
{
auto& tally {*model::tallies[i_tally]};
auto i_eout_filt = tally.filters()[tally.energyout_filter_];
auto i_bin = simulation::filter_matches[i_eout_filt].i_bin_;
auto bin_energyout = simulation::filter_matches[i_eout_filt].bins_[i_bin];
auto i_bin = p->filter_matches_[i_eout_filt].i_bin_;
auto bin_energyout = p->filter_matches_[i_eout_filt].bins_[i_bin];
const EnergyoutFilter& eo_filt
{*dynamic_cast<EnergyoutFilter*>(model::tally_filters[i_eout_filt].get())};
@ -331,8 +333,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::fission_bank.size() - p->n_bank_ + i;
const auto& bank = simulation::fission_bank[i_bank];
const auto& bank = p->nu_bank_[i];
// get the delayed group
auto g = bank.delayed_group;
@ -356,7 +357,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
g_out = eo_filt.n_bins() - g_out - 1;
// change outgoing energy bin
simulation::filter_matches[i_eout_filt].bins_[i_bin] = g_out;
p->filter_matches_[i_eout_filt].bins_[i_bin] = g_out;
} else {
@ -373,7 +374,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
} else {
auto i_match = lower_bound_index(eo_filt.bins().begin(),
eo_filt.bins().end(), E_out);
simulation::filter_matches[i_eout_filt].bins_[i_bin] = i_match;
p->filter_matches_[i_eout_filt].bins_[i_bin] = i_match;
}
}
@ -387,7 +388,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
int filter_index = 0;
for (auto j = 0; j < tally.filters().size(); ++j) {
auto i_filt = tally.filters(j);
auto& match {simulation::filter_matches[i_filt]};
auto& match {p->filter_matches_[i_filt]};
auto i_bin = match.i_bin_;
filter_index += match.bins_[i_bin] * tally.strides(j);
}
@ -415,13 +416,13 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
double filter_weight = 1.;
for (auto j = 0; j < tally.filters().size(); ++j) {
auto i_filt = tally.filters(j);
auto& match {simulation::filter_matches[i_filt]};
auto& match {p->filter_matches_[i_filt]};
auto i_bin = match.i_bin_;
filter_weight *= match.weights_[i_bin];
}
score_fission_delayed_dg(i_tally, d_bin, score*filter_weight,
i_score);
i_score, p->filter_matches_);
}
}
@ -433,7 +434,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
double filter_weight = 1.;
for (auto j = 0; j < tally.filters().size(); ++j) {
auto i_filt = tally.filters(j);
auto& match {simulation::filter_matches[i_filt]};
auto& match {p->filter_matches_[i_filt]};
auto i_bin = match.i_bin_;
filter_index += match.bins_[i_bin] * tally.strides(j);
filter_weight *= match.weights_[i_bin];
@ -447,7 +448,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
}
// Reset outgoing energy bin and score index
simulation::filter_matches[i_eout_filt].bins_[i_bin] = bin_energyout;
p->filter_matches_[i_eout_filt].bins_[i_bin] = bin_energyout;
}
//! Update tally results for continuous-energy tallies with any estimator.
@ -777,7 +778,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
* 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);
score_index, p->filter_matches_);
}
continue;
} else {
@ -811,7 +812,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
auto d = filt.groups()[d_bin];
score = simulation::keff * p->wgt_bank_ / p->n_bank_
* p->n_delayed_bank_[d-1] * flux;
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_);
}
continue;
} else {
@ -836,7 +837,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
->nu(E, ReactionProduct::EmissionMode::delayed, d);
score = p->neutron_xs_[i_nuclide].fission * yield
* atom_density * flux;
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_);
}
continue;
} else {
@ -866,7 +867,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
score = p->neutron_xs_[j_nuclide].fission * yield
* atom_density * flux;
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
score_index, p->filter_matches_);
}
}
}
@ -917,7 +918,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
/ p->neutron_xs_[p->event_nuclide_].absorption
* rate * flux;
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
score_index, p->filter_matches_);
}
continue;
} else {
@ -953,8 +954,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
// contribution to the fission bank to the score.
score = 0.;
for (auto i = 0; i < p->n_bank_; ++i) {
auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i;
const auto& bank = simulation::fission_bank[i_bank];
const auto& bank = p->nu_bank_[i];
auto g = bank.delayed_group;
if (g != 0) {
const auto& nuc {*data::nuclides[p->event_nuclide_]};
@ -971,7 +971,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
auto d = filt.groups()[d_bin];
if (d == g)
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
score_index, p->filter_matches_);
}
score = 0.;
}
@ -996,7 +996,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
auto rate = rxn.products_[d].decay_rate_;
score = p->neutron_xs_[i_nuclide].fission * yield * flux
* atom_density * rate;
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_);
}
continue;
} else {
@ -1037,7 +1037,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
score = p->neutron_xs_[j_nuclide].fission * yield
* flux * atom_density * rate;
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
score_index, p->filter_matches_);
}
}
}
@ -1225,8 +1225,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
// We need to substract the energy of the secondary particles since
// they will be transported individually later
for (auto i = 0; i < p->n_bank_second_; ++i) {
auto i_bank = simulation::secondary_bank.size() - p->n_bank_second_ + i;
const auto& bank = simulation::secondary_bank[i_bank];
auto i_bank = p->secondary_bank_.size() - p->n_bank_second_ + i;
const auto& bank = p->secondary_bank_[i_bank];
if (bank.particle == Particle::Type::photon ||
bank.particle == Particle::Type::neutron) {
score -= bank.E;
@ -1342,7 +1342,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
//! argument is really just used for filter weights.
void
score_general_mg(const Particle* p, int i_tally, int start_index,
score_general_mg(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]};
@ -1723,7 +1723,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
nullptr, nullptr, &d)
/ abs_xs;
}
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_);
}
continue;
} else {
@ -1764,7 +1764,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g)
/ macro_xs.get_xs(MgxsType::FISSION, p_g);
}
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_);
}
continue;
} else {
@ -1795,7 +1795,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
score = flux * macro_xs.get_xs(
MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
}
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_);
}
continue;
} else {
@ -1838,7 +1838,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
* macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d) / abs_xs;
}
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_);
}
continue;
} else {
@ -1876,8 +1876,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::fission_bank.size() - p->n_bank_ + i;
const auto& bank = simulation::fission_bank[i_bank];
const auto& bank = p->nu_bank_[i];
auto d = bank.delayed_group - 1;
if (d != -1) {
if (i_nuclide >= 0) {
@ -1899,7 +1898,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
auto dg = filt.groups()[d_bin];
if (dg == d + 1)
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
score_index, p->filter_matches_);
}
score = 0.;
}
@ -1929,7 +1928,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
* macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr,
nullptr, &d);
}
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_);
}
continue;
} else {
@ -2046,7 +2045,7 @@ void score_analog_tally_ce(Particle* p)
// no valid combinations, use a continue statement to ensure we skip the
// assume_separate break below.
auto filter_iter = FilterBinIter(tally, p);
auto end = FilterBinIter(tally, true);
auto end = FilterBinIter(tally, true, &p->filter_matches_);
if (filter_iter == end) continue;
// Loop over filter bins.
@ -2090,11 +2089,11 @@ void score_analog_tally_ce(Particle* p)
}
// Reset all the filter matches for the next tally event.
for (auto& match : simulation::filter_matches)
for (auto& match : p->filter_matches_)
match.bins_present_ = false;
}
void score_analog_tally_mg(const Particle* p)
void score_analog_tally_mg(Particle* p)
{
for (auto i_tally : model::active_analog_tallies) {
const Tally& tally {*model::tallies[i_tally]};
@ -2103,7 +2102,7 @@ void score_analog_tally_mg(const Particle* p)
// no valid combinations, use a continue statement to ensure we skip the
// assume_separate break below.
auto filter_iter = FilterBinIter(tally, p);
auto end = FilterBinIter(tally, true);
auto end = FilterBinIter(tally, true, &p->filter_matches_);
if (filter_iter == end) continue;
// Loop over filter bins.
@ -2135,7 +2134,7 @@ void score_analog_tally_mg(const Particle* p)
}
// Reset all the filter matches for the next tally event.
for (auto& match : simulation::filter_matches)
for (auto& match : p->filter_matches_)
match.bins_present_ = false;
}
@ -2152,7 +2151,7 @@ score_tracklength_tally(Particle* p, double distance)
// no valid combinations, use a continue statement to ensure we skip the
// assume_separate break below.
auto filter_iter = FilterBinIter(tally, p);
auto end = FilterBinIter(tally, true);
auto end = FilterBinIter(tally, true, &p->filter_matches_);
if (filter_iter == end) continue;
// Loop over filter bins.
@ -2199,7 +2198,7 @@ score_tracklength_tally(Particle* p, double distance)
}
// Reset all the filter matches for the next tally event.
for (auto& match : simulation::filter_matches)
for (auto& match : p->filter_matches_)
match.bins_present_ = false;
}
@ -2220,7 +2219,7 @@ void score_collision_tally(Particle* p)
// no valid combinations, use a continue statement to ensure we skip the
// assume_separate break below.
auto filter_iter = FilterBinIter(tally, p);
auto end = FilterBinIter(tally, true);
auto end = FilterBinIter(tally, true, &p->filter_matches_);
if (filter_iter == end) continue;
// Loop over filter bins.
@ -2263,12 +2262,12 @@ void score_collision_tally(Particle* p)
}
// Reset all the filter matches for the next tally event.
for (auto& match : simulation::filter_matches)
for (auto& match : p->filter_matches_)
match.bins_present_ = false;
}
void
score_surface_tally(const Particle* p, const std::vector<int>& tallies)
score_surface_tally(Particle* p, const std::vector<int>& tallies)
{
// No collision, so no weight change when survival biasing
double flux = p->wgt_;
@ -2280,7 +2279,7 @@ score_surface_tally(const Particle* p, const std::vector<int>& tallies)
// no valid combinations, use a continue statement to ensure we skip the
// assume_separate break below.
auto filter_iter = FilterBinIter(tally, p);
auto end = FilterBinIter(tally, true);
auto end = FilterBinIter(tally, true, &p->filter_matches_);
if (filter_iter == end) continue;
// Loop over filter bins.
@ -2307,7 +2306,7 @@ score_surface_tally(const Particle* p, const std::vector<int>& tallies)
}
// Reset all the filter matches for the next tally event.
for (auto& match : simulation::filter_matches)
for (auto& match : p->filter_matches_)
match.bins_present_ = false;
}

2
src/timer.cpp
View file

@ -16,6 +16,7 @@ Timer time_finalize;
Timer time_inactive;
Timer time_initialize;
Timer time_read_xs;
Timer time_sample_source;
Timer time_tallies;
Timer time_total;
Timer time_transport;
@ -68,6 +69,7 @@ void reset_timers()
simulation::time_inactive.reset();
simulation::time_initialize.reset();
simulation::time_read_xs.reset();
simulation::time_sample_source.reset();
simulation::time_tallies.reset();
simulation::time_total.reset();
simulation::time_transport.reset();

30
src/track_output.cpp
View file

@ -13,37 +13,27 @@
#include <string>
#include <vector>
// Explicit vector template specialization of threadprivate variable outside of
// the openmc namespace for the picky Intel compiler.
template class std::vector<std::vector<openmc::Position>>;
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
// Forward declaration needed in order to declare tracks as threadprivate
extern std::vector<std::vector<Position>> tracks;
#pragma omp threadprivate(tracks)
std::vector<std::vector<Position>> tracks;
//==============================================================================
// Non-member functions
//==============================================================================
void add_particle_track()
void add_particle_track(Particle& p)
{
tracks.emplace_back();
p.tracks_.emplace_back();
}
void write_particle_track(const Particle& p)
void write_particle_track(Particle& p)
{
tracks.back().push_back(p.r());
p.tracks_.back().push_back(p.r());
}
void finalize_particle_track(const Particle& p)
void finalize_particle_track(Particle& p)
{
std::stringstream filename;
filename << settings::path_output << "track_" << simulation::current_batch
@ -51,7 +41,7 @@ void finalize_particle_track(const Particle& p)
// Determine number of coordinates for each particle
std::vector<int> n_coords;
for (auto& coords : tracks) {
for (auto& coords : p.tracks_) {
n_coords.push_back(coords.size());
}
@ -60,10 +50,10 @@ void finalize_particle_track(const Particle& p)
hid_t file_id = file_open(filename.str().c_str(), 'w');
write_attribute(file_id, "filetype", "track");
write_attribute(file_id, "version", VERSION_TRACK);
write_attribute(file_id, "n_particles", tracks.size());
write_attribute(file_id, "n_particles", p.tracks_.size());
write_attribute(file_id, "n_coords", n_coords);
for (int i = 1; i <= tracks.size(); ++i) {
const auto& t {tracks[i-1]};
for (auto i = 1; i <= p.tracks_.size(); ++i) {
const auto& t {p.tracks_[i-1]};
size_t n = t.size();
xt::xtensor<double, 2> data({n,3});
for (int j = 0; j < n; ++j) {
@ -78,7 +68,7 @@ void finalize_particle_track(const Particle& p)
}
// Clear particle tracks
tracks.clear();
p.tracks_.clear();
}
} // namespace openmc

47
tests/unit_tests/test_geometry.py
View file

@ -301,3 +301,50 @@ def test_rotation_matrix():
assert geom.find((0.0, 0.5, 0.0))[-1] == c3
assert geom.find((0.0, -0.5, 0.0))[-1] == c1
assert geom.find((0.0, -1.5, 0.0))[-1] == c2
def test_remove_redundant_surfaces():
"""Test ability to remove redundant surfaces"""
m1 = openmc.Material()
m1.add_nuclide('U235', 1.0, 'wo')
m1.add_nuclide('O16', 2.0, 'wo')
m1.set_density('g/cm3', 10.0)
m2 = openmc.Material()
m2.add_element('Zr', 1.0)
m2.set_density('g/cm3', 2.0)
m3 = openmc.Material()
m3.add_nuclide('H1', 2.0)
m3.add_nuclide('O16', 1.0)
m3.set_density('g/cm3', 1.0)
def get_cyl_cell(r1, r2, z1, z2, fill):
"""Create a finite height cylindrical cell with a fill"""
cyl2 = openmc.ZCylinder(r=r2)
zplane1 = openmc.ZPlane(z1)
zplane2 = openmc.ZPlane(z2)
if np.isclose(r1, 0):
region = -cyl2 & +zplane1 & -zplane2
else:
cyl1 = openmc.ZCylinder(r=r1)
region = +cyl1 & -cyl2 & +zplane1 & -zplane2
return openmc.Cell(region=region, fill=fill)
r1, r2, r3 = 1., 2., 3.
z1, z2 = -2., 2.
fuel = get_cyl_cell(0, r1, z1, z2, m1)
clad = get_cyl_cell(r1, r2, z1, z2, m2)
water = get_cyl_cell(r2, r3, z1, z2, m3)
root = openmc.Universe(cells=[fuel, clad, water])
geom = openmc.Geometry(root)
# There should be 6 redundant surfaces in this geometry
n_redundant_surfs = len(geom.get_redundant_surfaces().keys())
assert n_redundant_surfs == 6
# Remove redundant surfaces
geom.remove_redundant_surfaces()
# There should be 0 remaining redundant surfaces
n_redundant_surfs = len(geom.get_redundant_surfaces().keys())
assert n_redundant_surfs == 0