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

Merge pull request #1461 from jtramm/event_based_update

Browse source 
Event-Based Mode
This commit is contained in:
Sterling Harper 2020-01-29 10:18:10 -05:00 • committed by GitHub
commit 328d9bbb34
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
30 changed files with 734 additions and 102 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.travis.yml
View file

@ -45,6 +45,8 @@ matrix:
env: OMP=n MPI=y PHDF5=y
- python: "3.7"
env: OMP=y MPI=y PHDF5=y DAGMC=y
- python: "3.7"
env: OMP=y MPI=n PHDF5=n EVENT=y
notifications:
webhooks: https://coveralls.io/webhook?repo_token=$COVERALLS_REPO_TOKEN
install:

1
CMakeLists.txt
View file

@ -223,6 +223,7 @@ list(APPEND libopenmc_SOURCES
src/eigenvalue.cpp
src/endf.cpp
src/error.cpp
src/event.cpp
src/initialize.cpp
src/finalize.cpp
src/geometry.cpp

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

@ -137,6 +137,15 @@ The ``<entropy_mesh>`` element indicates the ID of a mesh that is to be used for
calculating Shannon entropy. The mesh should cover all possible fissionable
materials in the problem and is specified using a :ref:`mesh_element`.
----------------------------
``<event_based>``
----------------------------
Determines whether to use event-based parallelism instead of the default
history-based parallelism.
*Default*: false
-----------------------------------
``<generations_per_batch>`` Element
-----------------------------------
@ -218,6 +227,16 @@ to false.
*Default*: true
----------------------------------------
``<max_particles_in_flight>`` Element
----------------------------------------
This element indicates the number of neutrons to run in flight concurrently
when using event-based parallelism. A higher value uses more memory, but
may be more efficient computationally.
*Default*: 100000
---------------------------
``<max_order>`` Element
---------------------------

1
docs/source/usersguide/scripts.rst
View file

@ -33,6 +33,7 @@ working directory which may be different from
flags:
-c, --volume Run in stochastic volume calculation mode
-e, --event Run using event-based parallelism
-g, --geometry-debug Run in geometry debugging mode, where cell overlaps are
checked for after each move of a particle
-n, --particles N Use *N* particles per generation or batch

6
include/openmc/bank.h
View file

@ -4,6 +4,7 @@
#include <cstdint>
#include <vector>
#include "openmc/shared_array.h"
#include "openmc/particle.h"
#include "openmc/position.h"
@ -17,9 +18,8 @@ namespace simulation {
extern std::vector<Particle::Bank> source_bank;
extern std::unique_ptr<Particle::Bank[]> fission_bank;
extern int64_t fission_bank_length;
extern int64_t fission_bank_max;
extern SharedArray<Particle::Bank> fission_bank;
extern std::vector<int64_t> progeny_per_particle;
} // namespace simulation

115
include/openmc/event.h Normal file
View file

@ -0,0 +1,115 @@
#ifndef OPENMC_EVENT_H
#define OPENMC_EVENT_H
//! \file event.h
//! \brief Event-based data structures and methods
#include "openmc/particle.h"
#include "openmc/shared_array.h"
namespace openmc {
//==============================================================================
// Structs
//==============================================================================
// In the event-based model, instead of moving or sorting the particles
// themselves based on which event they need, a queue is used to store the
// index (and other useful info) for each event type.
// The EventQueueItem struct holds the relevant information about a particle needed
// for sorting the queue. For very high particle counts, a sorted queue has the
// potential to result in greatly improved cache efficiency. However, sorting
// will introduce some overhead due to the sorting process itself, and may not
// result in any benefits if not enough particles are present for them to achieve
// consistent locality improvements.
struct EventQueueItem{
int64_t idx; //!< particle index in event-based particle buffer
Particle::Type type; //!< particle type
int64_t material; //!< material that particle is in
double E; //!< particle energy
// Constructors
EventQueueItem() = default;
EventQueueItem(const Particle& p, int64_t buffer_idx) :
idx(buffer_idx), type(p.type_), material(p.material_), E(p.E_) {}
// Compare by particle type, then by material type (4.5% fuel/7.0% fuel/cladding/etc),
// then by energy.
// TODO: Currently in OpenMC, the material ID corresponds not only to a general
// type, but also specific isotopic densities. Ideally we would
// like to be able to just sort by general material type, regardless of densities.
// A more general material type ID may be added in the future, in which case we
// can update the material field of this struct to contain the more general id.
bool operator<(const EventQueueItem& rhs) const
{
return std::tie(type, material, E) < std::tie(rhs.type, rhs.material, rhs.E);
}
};
//==============================================================================
// Global variable declarations
//==============================================================================
namespace simulation {
// Event queues. These use the special SharedArray type, rather than a normal
// vector, as they will be shared between threads and may be appended to at the
// same time. To facilitate this, the SharedArray thread_safe_append() method
// is provided which controls the append operations using atomics.
extern SharedArray<EventQueueItem> calculate_fuel_xs_queue;
extern SharedArray<EventQueueItem> calculate_nonfuel_xs_queue;
extern SharedArray<EventQueueItem> advance_particle_queue;
extern SharedArray<EventQueueItem> surface_crossing_queue;
extern SharedArray<EventQueueItem> collision_queue;
// Particle buffer
extern std::vector<Particle> particles;
} // namespace simulation
//==============================================================================
// Functions
//==============================================================================
//! Allocate space for the event queues and particle buffer
//
//! \param n_particles The number of particles in the particle buffer
void init_event_queues(int64_t n_particles);
//! Free the event queues and particle buffer
void free_event_queues(void);
//! Enqueue a particle based on if it is in fuel or a non-fuel material
//
//! \param buffer_idx The particle's actual index in the particle buffer
void dispatch_xs_event(int64_t buffer_idx);
//! Execute the initialization event for all particles
//
//! \param n_particles The number of particles in the particle buffer
//! \param source_offset The offset index in the source bank to use
void process_init_events(int64_t n_particles, int64_t source_offset);
//! Execute the calculate XS event for all particles in this event's buffer
//
//! \param queue A reference to the desired XS lookup queue
void process_calculate_xs_events(SharedArray<EventQueueItem>& queue);
//! Execute the advance particle event for all particles in this event's buffer
void process_advance_particle_events();
//! Execute the surface crossing event for all particles in this event's buffer
void process_surface_crossing_events();
//! Execute the collision event for all particles in this event's buffer
void process_collision_events();
//! Execute the death event for all particles
//
//! \param n_particles The number of particles in the particle buffer
void process_death_events(int64_t n_particles);
} // namespace openmc
#endif // OPENMC_EVENT_H

4
include/openmc/settings.h
View file

@ -31,6 +31,7 @@ extern "C" bool cmfd_run; //!< is a CMFD run?
extern "C" bool dagmc; //!< indicator of DAGMC geometry
extern bool delayed_photon_scaling; //!< Scale fission photon yield to include delayed
extern "C" bool entropy_on; //!< calculate Shannon entropy?
extern bool event_based; //!< use event-based mode (instead of history-based)
extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular?
extern bool material_cell_offsets; //!< create material cells offsets?
extern "C" bool output_summary; //!< write summary.h5?
@ -67,6 +68,9 @@ extern "C" int32_t n_inactive; //!< number of inactive batches
extern "C" int32_t gen_per_batch; //!< number of generations per batch
extern "C" int64_t n_particles; //!< number of particles per generation
extern int64_t max_particles_in_flight; //!< Max num. event-based particles in flight
extern ElectronTreatment electron_treatment; //!< how to treat secondary electrons
extern std::array<double, 4> energy_cutoff; //!< Energy cutoff in [eV] for each particle type
extern int legendre_to_tabular_points; //!< number of points to convert Legendres

131
include/openmc/shared_array.h Normal file
View file

@ -0,0 +1,131 @@
#ifndef OPENMC_SHARED_ARRAY_H
#define OPENMC_SHARED_ARRAY_H
//! \file shared_array.h
//! \brief Shared array data structure
#include <memory>
namespace openmc {
//==============================================================================
// Class declarations
//==============================================================================
// This container is an array that is capable of being appended to in an
// thread safe manner by use of atomics. It only provides protection for the
// use cases currently present in OpenMC. Namely, it covers the scenario where
// multiple threads are appending to an array, but no threads are reading from
// or operating on it in any other way at the same time. Multiple threads can
// call the thread_safe_append() function concurrently and store data to the
// object at the index returned from thread_safe_append() safely, but no other
// operations are protected.
template <typename T>
class SharedArray {
public:
//==========================================================================
// Constructors
//! Default constructor.
SharedArray() = default;
//! Construct a zero size container with space to hold capacity number of
//! elements.
//
//! \param capacity The number of elements for the container to allocate
//! space for
SharedArray(int64_t capacity) : capacity_(capacity)
{
data_ = std::make_unique<T[]>(capacity);
}
//==========================================================================
// Methods and Accessors
//! Return a reference to the element at specified location i. No bounds
//! checking is performed.
T& operator[](int64_t i) {return data_[i];}
const T& operator[](int64_t i) const { return data_[i]; }
//! Allocate space in the container for the specified number of elements.
//! reserve() does not change the size of the container.
//
//! \param capacity The number of elements to allocate in the container
void reserve(int64_t capacity)
{
data_ = std::make_unique<T[]>(capacity);
capacity_ = capacity;
}
//! Increase the size of the container by one and append value to the
//! array. Returns an index to the element of the array written to. Also
//! tests to enforce that the append operation does not read off the end
//! of the array. In the event that this does happen, set the size to be
//! equal to the capacity and return -1.
//
//! \value The value of the element to append
//! \return The index in the array written to. In the event that this
//! index would be greater than what was allocated for the container,
//! return -1.
int64_t thread_safe_append(const T& value)
{
// Atomically capture the index we want to write to
int64_t idx;
#pragma omp atomic capture
idx = size_++;
// Check that we haven't written off the end of the array
if (idx >= capacity_) {
#pragma omp atomic write
size_ = capacity_;
return -1;
}
// Copy element value to the array
data_[idx] = value;
return idx;
}
//! Free any space that was allocated for the container. Set the
//! container's size and capacity to 0.
void clear()
{
data_.reset();
size_ = 0;
capacity_ = 0;
}
//! Return the number of elements in the container
int64_t size() {return size_;}
//! Resize the container to contain a specified number of elements. This is
//! useful in cases where the container is written to in a non-thread safe manner,
//! where the internal size of the array needs to be manually updated.
//
//! \param size The new size of the container
void resize(int64_t size) {size_ = size;}
//! Return the number of elements that the container has currently allocated
//! space for.
int64_t capacity() {return capacity_;}
//! Return pointer to the underlying array serving as element storage.
T* data() {return data_.get();}
const T* data() const {return data_.get();}
private:
//==========================================================================
// Data members
std::unique_ptr<T[]> data_; //!< An RAII handle to the elements
int64_t size_ {0}; //!< The current number of elements
int64_t capacity_ {0}; //!< The total space allocated for elements
};
} // namespace openmc
#endif // OPENMC_SHARED_ARRAY_H

3
include/openmc/simulation.h
View file

@ -95,6 +95,9 @@ void transport_history_based_single_particle(Particle& p);
//! Simulate all particle histories using history-based parallelism
void transport_history_based();
//! Simulate all particle histories using event-based parallelism
void transport_event_based();
} // namespace openmc
#endif // OPENMC_SIMULATION_H

6
include/openmc/timer.h
View file

@ -25,6 +25,12 @@ extern Timer time_sample_source;
extern Timer time_tallies;
extern Timer time_total;
extern Timer time_transport;
extern Timer time_event_init;
extern Timer time_event_calculate_xs;
extern Timer time_event_advance_particle;
extern Timer time_event_surface_crossing;
extern Timer time_event_collision;
extern Timer time_event_death;
} // namespace simulation

3
man/man1/openmc.1
View file

@ -16,6 +16,9 @@ is specified, the XML input files are present in the current directory.
.B "\-c\fR, \fP\-\-volume"
Run in stochastic volume calculation mode
.TP
.B "\-e\fR, \fP\-\-event"
Run using event-based parallelism
.TP
.B "\-g\fR, \fP\-\-geometry-debug"
Run with geometry debugging turned on, where cell overlaps are checked for after
each move of a particle

7
openmc/executor.py
View file

@ -154,7 +154,7 @@ def calculate_volumes(threads=None, output=True, cwd='.',
def run(particles=None, threads=None, geometry_debug=False,
restart_file=None, tracks=False, output=True, cwd='.',
openmc_exec='openmc', mpi_args=None):
openmc_exec='openmc', mpi_args=None, event_based=False):
"""Run an OpenMC simulation.
Parameters
@ -182,6 +182,8 @@ def run(particles=None, threads=None, geometry_debug=False,
mpi_args : list of str, optional
MPI execute command and any additional MPI arguments to pass,
e.g. ['mpiexec', '-n', '8'].
event_based : bool, optional
Turns on event-based parallelism, instead of default history-based
Raises
------
@ -199,6 +201,9 @@ def run(particles=None, threads=None, geometry_debug=False,
if geometry_debug:
args.append('-g')
if event_based:
args.append('-e')
if isinstance(restart_file, str):
args += ['-r', restart_file]

52
openmc/settings.py
View file

@ -53,6 +53,9 @@ class Settings(object):
Mesh to be used to calculate Shannon entropy. If the mesh dimensions are
not specified. OpenMC assigns a mesh such that 20 source sites per mesh
cell are to be expected on average.
event_based : bool
Indicate whether to use event-based parallelism instead of the default
history-based parallelism.
generations_per_batch : int
Number of generations per batch
inactive : int
@ -68,6 +71,9 @@ class Settings(object):
material_cell_offsets : bool
Generate an "offset table" for material cells by default. These tables
are necessary when a particular instance of a cell needs to be tallied.
max_particles_in_flight : int
Number of neutrons to run concurrently when using event-based
parallelism.
max_order : None or int
Maximum scattering order to apply globally when in multi-group mode.
no_reduce : bool
@ -229,6 +235,9 @@ class Settings(object):
self._dagmc = False
self._event_based = None
self._max_particles_in_flight = None
@property
def run_mode(self):
return self._run_mode
@ -376,6 +385,14 @@ class Settings(object):
@property
def dagmc(self):
return self._dagmc
@property
def event_based(self):
return self._event_based
@property
def max_particles_in_flight(self):
return self._max_particles_in_flight
@run_mode.setter
def run_mode(self, run_mode):
@ -704,6 +721,17 @@ class Settings(object):
def delayed_photon_scaling(self, value):
cv.check_type('delayed photon scaling', value, bool)
self._delayed_photon_scaling = value
@event_based.setter
def event_based(self, value):
cv.check_type('event based', value, bool)
self._event_based = value
@max_particles_in_flight.setter
def max_particles_in_flight(self, value):
cv.check_type('max particles in flight', value, Integral)
cv.check_greater_than('max particles in flight', value, 0)
self._max_particles_in_flight = value
@material_cell_offsets.setter
def material_cell_offsets(self, value):
@ -948,6 +976,16 @@ class Settings(object):
if self._delayed_photon_scaling is not None:
elem = ET.SubElement(root, "delayed_photon_scaling")
elem.text = str(self._delayed_photon_scaling).lower()
def _create_event_based_subelement(self, root):
if self._event_based is not None:
elem = ET.SubElement(root, "event_based")
elem.text = str(self._event_based).lower()
def _create_max_particles_in_flight_subelement(self, root):
if self._max_particles_in_flight is not None:
elem = ET.SubElement(root, "max_particles_in_flight")
elem.text = str(self._max_particles_in_flight).lower()
def _create_material_cell_offsets_subelement(self, root):
if self._material_cell_offsets is not None:
@ -1189,6 +1227,16 @@ class Settings(object):
text = get_text(root, 'delayed_photon_scaling')
if text is not None:
self.delayed_photon_scaling = text in ('true', '1')
def _event_based_from_xml_element(self, root):
text = get_text(root, 'event_based')
if text is not None:
self.event_based = text in ('true', '1')
def _max_particles_in_flight_from_xml_element(self, root):
text = get_text(root, 'max_particles_in_flight')
if text is not None:
self.max_particles_in_flight = int(text)
def _material_cell_offsets_from_xml_element(self, root):
text = get_text(root, 'material_cell_offsets')
@ -1250,6 +1298,8 @@ class Settings(object):
self._create_volume_calcs_subelement(root_element)
self._create_create_fission_neutrons_subelement(root_element)
self._create_delayed_photon_scaling_subelement(root_element)
self._create_event_based_subelement(root_element)
self._create_max_particles_in_flight_subelement(root_element)
self._create_material_cell_offsets_subelement(root_element)
self._create_log_grid_bins_subelement(root_element)
self._create_dagmc_subelement(root_element)
@ -1317,6 +1367,8 @@ class Settings(object):
settings._resonance_scattering_from_xml_element(root)
settings._create_fission_neutrons_from_xml_element(root)
settings._delayed_photon_scaling_from_xml_element(root)
settings._event_based_from_xml_element(root)
settings._max_particles_in_flight_from_xml_element(root)
settings._material_cell_offsets_from_xml_element(root)
settings._log_grid_bins_from_xml_element(root)
settings._dagmc_from_xml_element(root)

35
src/bank.cpp
View file

@ -17,15 +17,11 @@ namespace simulation {
std::vector<Particle::Bank> source_bank;
// The fission bank is allocated as a pointer, rather than a vector, as it will
// The fission bank is allocated as a SharedArray, 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;
// added to it by using SharedArray's special thread_safe_append() function.
SharedArray<Particle::Bank> fission_bank;
// Each entry in this vector corresponds to the number of progeny produced
// this generation for the particle located at that index. This vector is
@ -41,16 +37,13 @@ std::vector<int64_t> progeny_per_particle;
void free_memory_bank()
{
simulation::source_bank.clear();
simulation::fission_bank.reset();
simulation::fission_bank_length = 0;
simulation::fission_bank.clear();
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::fission_bank.reserve(max);
simulation::progeny_per_particle.resize(simulation::work_per_rank);
}
@ -87,15 +80,15 @@ void sort_fission_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);
if (simulation::fission_bank.size() > simulation::fission_bank.capacity() / 2) {
sorted_bank_holder.resize(simulation::fission_bank.size());
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];
sorted_bank = &simulation::fission_bank[simulation::fission_bank.size()];
}
// Use parent and progeny indices to sort fission bank
for (int64_t i = 0; i < simulation::fission_bank_length; i++) {
for (int64_t i = 0; i < simulation::fission_bank.size(); 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;
@ -103,8 +96,8 @@ void sort_fission_bank()
}
// Copy sorted bank into the fission bank
std::copy(sorted_bank, sorted_bank + simulation::fission_bank_length,
simulation::fission_bank.get());
std::copy(sorted_bank, sorted_bank + simulation::fission_bank.size(),
simulation::fission_bank.data());
}
//==============================================================================
@ -125,12 +118,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_length == 0) {
if (simulation::fission_bank.size() == 0) {
set_errmsg("Fission bank has not been allocated.");
return OPENMC_E_ALLOCATE;
} else {
*ptr = simulation::fission_bank.get();
*n = simulation::fission_bank_length;
*ptr = simulation::fission_bank.data();
*n = simulation::fission_bank.size();
return 0;
}
}

14
src/eigenvalue.cpp
View file

@ -83,7 +83,7 @@ void synchronize_bank()
#ifdef OPENMC_MPI
int64_t start = 0;
int64_t n_bank = simulation::fission_bank_length;
int64_t n_bank = simulation::fission_bank.size();
MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
// While we would expect the value of start on rank 0 to be 0, the MPI
@ -91,13 +91,13 @@ void synchronize_bank()
// significant
if (mpi::rank == 0) start = 0;
int64_t finish = start + simulation::fission_bank_length;
int64_t finish = start + simulation::fission_bank.size();
int64_t total = finish;
MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
#else
int64_t start = 0;
int64_t finish = simulation::fission_bank_length;
int64_t finish = simulation::fission_bank.size();
int64_t total = finish;
#endif
@ -106,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_length == 0) {
if (simulation::fission_bank.size() == 0) {
fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
}
@ -138,7 +138,7 @@ void synchronize_bank()
int64_t index_temp = 0;
std::vector<Particle::Bank> temp_sites(3*simulation::work_per_rank);
for (int64_t i = 0; i < simulation::fission_bank_length; i++ ) {
for (int64_t i = 0; i < simulation::fission_bank.size(); i++ ) {
const auto& site = simulation::fission_bank[i];
// If there are less than n_particles particles banked, automatically add
@ -195,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_length - sites_needed + i;
int i_bank = simulation::fission_bank.size() - sites_needed + i;
temp_sites[index_temp] = simulation::fission_bank[i_bank];
++index_temp;
}
@ -506,7 +506,7 @@ 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.get(), simulation::fission_bank_length,
simulation::fission_bank.data(), simulation::fission_bank.size(),
&sites_outside);
// display warning message if there were sites outside entropy box

173
src/event.cpp Normal file
View file

@ -0,0 +1,173 @@
#include "openmc/event.h"
#include "openmc/material.h"
#include "openmc/simulation.h"
#include "openmc/timer.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace simulation {
SharedArray<EventQueueItem> calculate_fuel_xs_queue;
SharedArray<EventQueueItem> calculate_nonfuel_xs_queue;
SharedArray<EventQueueItem> advance_particle_queue;
SharedArray<EventQueueItem> surface_crossing_queue;
SharedArray<EventQueueItem> collision_queue;
std::vector<Particle> particles;
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================
void init_event_queues(int64_t n_particles)
{
simulation::calculate_fuel_xs_queue.reserve(n_particles);
simulation::calculate_nonfuel_xs_queue.reserve(n_particles);
simulation::advance_particle_queue.reserve(n_particles);
simulation::surface_crossing_queue.reserve(n_particles);
simulation::collision_queue.reserve(n_particles);
simulation::particles.resize(n_particles);
}
void free_event_queues(void)
{
simulation::calculate_fuel_xs_queue.clear();
simulation::calculate_nonfuel_xs_queue.clear();
simulation::advance_particle_queue.clear();
simulation::surface_crossing_queue.clear();
simulation::collision_queue.clear();
simulation::particles.clear();
}
void dispatch_xs_event(int64_t buffer_idx)
{
Particle& p = simulation::particles[buffer_idx];
if (p.material_ == MATERIAL_VOID || !model::materials[p.material_]->fissionable_) {
simulation::calculate_nonfuel_xs_queue.thread_safe_append({p, buffer_idx});
} else {
simulation::calculate_fuel_xs_queue.thread_safe_append({p, buffer_idx});
}
}
void process_init_events(int64_t n_particles, int64_t source_offset)
{
simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < n_particles; i++) {
initialize_history(&simulation::particles[i], source_offset + i + 1);
dispatch_xs_event(i);
}
simulation::time_event_init.stop();
}
void process_calculate_xs_events(SharedArray<EventQueueItem>& queue)
{
simulation::time_event_calculate_xs.start();
// TODO: If using C++17, perform a parallel sort of the queue
// by particle type, material type, and then energy, in order to
// improve cache locality and reduce thread divergence on GPU. Prior
// to C++17, std::sort is a serial only operation, which in this case
// makes it too slow to be practical for most test problems.
//
// std::sort(std::execution::par_unseq, queue.data(), queue.data() + queue.size());
int64_t offset = simulation::advance_particle_queue.size();;
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < queue.size(); i++) {
Particle* p = &simulation::particles[queue[i].idx];
p->event_calculate_xs();
// After executing a calculate_xs event, particles will
// always require an advance event. Therefore, we don't need to use
// the protected enqueuing function.
simulation::advance_particle_queue[offset + i] = queue[i];
}
simulation::advance_particle_queue.resize(offset + queue.size());
queue.resize(0);
simulation::time_event_calculate_xs.stop();
}
void process_advance_particle_events()
{
simulation::time_event_advance_particle.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
Particle& p = simulation::particles[buffer_idx];
p.event_advance();
if (p.collision_distance_ > p.boundary_.distance) {
simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});
} else {
simulation::collision_queue.thread_safe_append({p, buffer_idx});
}
}
simulation::advance_particle_queue.resize(0);
simulation::time_event_advance_particle.stop();
}
void process_surface_crossing_events()
{
simulation::time_event_surface_crossing.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {
int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
Particle& p = simulation::particles[buffer_idx];
p.event_cross_surface();
p.event_revive_from_secondary();
if (p.alive_)
dispatch_xs_event(buffer_idx);
}
simulation::surface_crossing_queue.resize(0);
simulation::time_event_surface_crossing.stop();
}
void process_collision_events()
{
simulation::time_event_collision.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {
int64_t buffer_idx = simulation::collision_queue[i].idx;
Particle& p = simulation::particles[buffer_idx];
p.event_collide();
p.event_revive_from_secondary();
if (p.alive_)
dispatch_xs_event(buffer_idx);
}
simulation::collision_queue.resize(0);
simulation::time_event_collision.stop();
}
void process_death_events(int64_t n_particles)
{
simulation::time_event_death.start();
#pragma omp parallel for schedule(runtime)
for (int64_t i = 0; i < n_particles; i++) {
Particle& p = simulation::particles[i];
p.event_death();
}
simulation::time_event_death.stop();
}
} // namespace openmc

6
src/finalize.cpp
View file

@ -7,6 +7,7 @@
#include "openmc/cross_sections.h"
#include "openmc/dagmc.h"
#include "openmc/eigenvalue.h"
#include "openmc/event.h"
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/material.h"
@ -47,6 +48,9 @@ void free_memory()
if (mpi::master) {
free_memory_cmfd();
}
if (settings::event_based) {
free_event_queues();
}
#ifdef DAGMC
free_memory_dagmc();
#endif
@ -76,7 +80,9 @@ int openmc_finalize()
settings::gen_per_batch = 1;
settings::legendre_to_tabular = true;
settings::legendre_to_tabular_points = -1;
settings::event_based = false;
settings::material_cell_offsets = true;
settings::max_particles_in_flight = 100000;
settings::n_particles = -1;
settings::output_summary = true;
settings::output_tallies = true;

3
src/initialize.cpp
View file

@ -134,6 +134,9 @@ parse_command_line(int argc, char* argv[])
} else if (arg == "-n" || arg == "--particles") {
i += 1;
settings::n_particles = std::stoll(argv[i]);
} else if (arg == "-e" || arg == "--event") {
settings::event_based = true;
} else if (arg == "-r" || arg == "--restart") {
i += 1;

9
src/output.cpp
View file

@ -327,6 +327,7 @@ void print_usage()
" or a particle restart file\n"
" -s, --threads Number of OpenMP threads\n"
" -t, --track Write tracks for all particles\n"
" -e, --event Run using event-based parallelism\n"
" -v, --version Show version information\n"
" -h, --help Show this message\n";
}
@ -466,6 +467,14 @@ void print_runtime()
show_time("Total time in simulation", time_inactive.elapsed() +
time_active.elapsed());
show_time("Time in transport only", time_transport.elapsed(), 1);
if (settings::event_based) {
show_time("Particle initialization", time_event_init.elapsed(), 2);
show_time("XS lookups", time_event_calculate_xs.elapsed(), 2);
show_time("Advancing", time_event_advance_particle.elapsed(), 2);
show_time("Surface crossings", time_event_surface_crossing.elapsed(), 2);
show_time("Collisions", time_event_collision.elapsed(), 2);
show_time("Particle death", time_event_death.elapsed(), 2);
}
if (settings::run_mode == RunMode::EIGENVALUE) {
show_time("Time in inactive batches", time_inactive.elapsed(), 1);
}

43
src/physics.cpp
View file

@ -179,37 +179,32 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
for (int i = 0; i < nu; ++i) {
Particle::Bank* site;
// Initialize fission site object with particle data
Particle::Bank site;
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());
// Store fission site in bank
if (use_fission_bank) {
int64_t idx;
#pragma omp atomic capture
idx = simulation::fission_bank_length++;
if (idx >= simulation::fission_bank_max) {
int64_t idx = simulation::fission_bank.thread_safe_append(site);
if (idx == -1) {
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();
p->secondary_bank_.push_back(site);
}
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());
// 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) {
@ -220,9 +215,9 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
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;
nu_bank_entry->wgt = site.wgt;
nu_bank_entry->E = site.E;
nu_bank_entry->delayed_group = site.delayed_group;
}
}

65
src/physics_mg.cpp
View file

@ -125,33 +125,13 @@ create_fission_sites(Particle* p)
bool use_fission_bank = (settings::run_mode == RunMode::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();
}
// Bank source neutrons by copying the particle data
site->r = p->r();
site->particle = Particle::Type::neutron;
site->wgt = 1. / weight;
site->parent_id = p->id_;
site->progeny_id = p->n_progeny_++;
// Initialize fission site object with particle data
Particle::Bank site;
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
@ -159,20 +139,35 @@ create_fission_sites(Particle* p)
// 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;
int gout;
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;
// Store fission site in bank
if (use_fission_bank) {
int64_t idx = simulation::fission_bank.thread_safe_append(site);
if (idx == -1) {
warning("The shared fission bank is full. Additional fission sites created "
"in this generation will not be banked.");
skipped++;
break;
}
} else {
p->secondary_bank_.push_back(site);
}
// Set the delayed group on the particle as well
p->delayed_group_ = dg + 1;
@ -186,9 +181,9 @@ create_fission_sites(Particle* p)
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;
nu_bank_entry->wgt = site.wgt;
nu_bank_entry->E = site.E;
nu_bank_entry->delayed_group = site.delayed_group;
}
}

6
src/relaxng/settings.rnc
View file

@ -23,7 +23,9 @@ element settings {
element energy_mode { ( "continuous-energy" | "ce" | "CE" | "multi-group" | "mg" | "MG" ) }? &
element entropy_mesh { xsd:positiveInteger }? &
element event_based { xsd:boolean }? &
element generations_per_batch { xsd:positiveInteger }? &
element inactive { xsd:nonNegativeInteger }? &
@ -36,6 +38,8 @@ element settings {
element log_grid_bins { xsd:positiveInteger }? &
element material_cell_offsets { xsd:boolean }? &
element max_particles_in_flight { xsd:positiveInteger }? &
element max_order { xsd:nonNegativeInteger }? &

10
src/relaxng/settings.rng
View file

@ -87,6 +87,16 @@
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="event_based">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="max_particles_in_flight">
<data type="positiveInteger"/>
</element>
</optional>
<optional>
<element name="electron_treatment">
<choice>

14
src/settings.cpp
View file

@ -46,6 +46,7 @@ bool create_fission_neutrons {true};
bool dagmc {false};
bool delayed_photon_scaling {true};
bool entropy_on {false};
bool event_based {false};
bool legendre_to_tabular {true};
bool material_cell_offsets {true};
bool output_summary {true};
@ -81,6 +82,8 @@ int32_t n_inactive {0};
int32_t gen_per_batch {1};
int64_t n_particles {-1};
int64_t max_particles_in_flight {100000};
ElectronTreatment electron_treatment {ElectronTreatment::TTB};
std::array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
int legendre_to_tabular_points {C_NONE};
@ -127,6 +130,12 @@ void get_run_parameters(pugi::xml_node node_base)
if (n_particles == -1) {
n_particles = std::stoll(get_node_value(node_base, "particles"));
}
// Get maximum number of in flight particles for event-based mode
if (check_for_node(node_base, "max_particles_in_flight")) {
max_particles_in_flight = std::stoll(get_node_value(node_base,
"max_particles_in_flight"));
}
// Get number of basic batches
if (check_for_node(node_base, "batches")) {
@ -783,6 +792,11 @@ void read_settings_xml()
if (check_for_node(root, "delayed_photon_scaling")) {
delayed_photon_scaling = get_node_value_bool(root, "delayed_photon_scaling");
}
// Check whether to use event-based parallelism
if (check_for_node(root, "event_based")) {
event_based = get_node_value_bool(root, "event_based");
}
// Check whether material cell offsets should be generated
if (check_for_node(root, "material_cell_offsets")) {

70
src/simulation.cpp
View file

@ -5,6 +5,7 @@
#include "openmc/container_util.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/event.h"
#include "openmc/geometry_aux.h"
#include "openmc/material.h"
#include "openmc/message_passing.h"
@ -72,6 +73,14 @@ int openmc_simulation_init()
// fission bank
allocate_banks();
// If doing an event-based simulation, intialize the particle buffer
// and event queues
if (settings::event_based) {
int64_t event_buffer_length = std::min(simulation::work_per_rank,
settings::max_particles_in_flight);
init_event_queues(event_buffer_length);
}
// Allocate tally results arrays if they're not allocated yet
for (auto& t : model::tallies) {
t->init_results();
@ -183,7 +192,11 @@ int openmc_next_batch(int* status)
simulation::time_transport.start();
// Transport loop
transport_history_based();
if (settings::event_based) {
transport_event_based();
} else {
transport_history_based();
}
// Accumulate time for transport
simulation::time_transport.stop();
@ -262,7 +275,6 @@ void allocate_banks()
if (settings::run_mode == RunMode::EIGENVALUE) {
init_fission_bank(3*simulation::work_per_rank);
}
}
void initialize_batch()
@ -363,7 +375,7 @@ void initialize_generation()
{
if (settings::run_mode == RunMode::EIGENVALUE) {
// Clear out the fission bank
simulation::fission_bank_length = 0;
simulation::fission_bank.resize(0);
// Count source sites if using uniform fission source weighting
if (settings::ufs_on) ufs_count_sites();
@ -593,4 +605,56 @@ void transport_history_based()
}
}
void transport_event_based()
{
int64_t remaining_work = simulation::work_per_rank;
int64_t source_offset = 0;
// To cap the total amount of memory used to store particle object data, the
// number of particles in flight at any point in time can bet set. In the case
// that the maximum in flight particle count is lower than the total number
// of particles that need to be run this iteration, the event-based transport
// loop is executed multiple times until all particles have been completed.
while (remaining_work > 0) {
// Figure out # of particles to run for this subiteration
int64_t n_particles = std::min(remaining_work, settings::max_particles_in_flight);
// Initialize all particle histories for this subiteration
process_init_events(n_particles, source_offset);
// Event-based transport loop
while (true) {
// Determine which event kernel has the longest queue
int64_t max = std::max({
simulation::calculate_fuel_xs_queue.size(),
simulation::calculate_nonfuel_xs_queue.size(),
simulation::advance_particle_queue.size(),
simulation::surface_crossing_queue.size(),
simulation::collision_queue.size()});
// Execute event with the longest queue
if (max == 0) {
break;
} else if (max == simulation::calculate_fuel_xs_queue.size()) {
process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
} else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
} else if (max == simulation::advance_particle_queue.size()) {
process_advance_particle_events();
} else if (max == simulation::surface_crossing_queue.size()) {
process_surface_crossing_events();
} else if (max == simulation::collision_queue.size()) {
process_collision_events();
}
}
// Execute death event for all particles
process_death_events(n_particles);
// Adjust remaining work and source offset variables
remaining_work -= n_particles;
source_offset += n_particles;
}
}
} // namespace openmc

12
src/timer.cpp
View file

@ -20,6 +20,12 @@ Timer time_sample_source;
Timer time_tallies;
Timer time_total;
Timer time_transport;
Timer time_event_init;
Timer time_event_calculate_xs;
Timer time_event_advance_particle;
Timer time_event_surface_crossing;
Timer time_event_collision;
Timer time_event_death;
} // namespace simulation
@ -73,6 +79,12 @@ void reset_timers()
simulation::time_tallies.reset();
simulation::time_total.reset();
simulation::time_transport.reset();
simulation::time_event_init.reset();
simulation::time_event_calculate_xs.reset();
simulation::time_event_advance_particle.reset();
simulation::time_event_surface_crossing.reset();
simulation::time_event_collision.reset();
simulation::time_event_death.reset();
}
} // namespace openmc

3
tests/conftest.py
View file

@ -10,10 +10,11 @@ def pytest_addoption(parser):
parser.addoption('--mpi-np')
parser.addoption('--update', action='store_true')
parser.addoption('--build-inputs', action='store_true')
parser.addoption('--event', action='store_true')
def pytest_configure(config):
opts = ['exe', 'mpi', 'mpiexec', 'mpi_np', 'update', 'build_inputs']
opts = ['exe', 'mpi', 'mpiexec', 'mpi_np', 'update', 'build_inputs', 'event']
for opt in opts:
if config.getoption(opt) is not None:
regression_config[opt] = config.getoption(opt)

1
tests/regression_tests/__init__.py
View file

@ -1,5 +1,6 @@
# Test configuration options for regression tests
config = {
'event' : False,
'exe': 'openmc',
'mpi': False,
'mpiexec': 'mpiexec',

5
tests/testing_harness.py
View file

@ -68,9 +68,10 @@ class TestHarness(object):
def _run_openmc(self):
if config['mpi']:
mpi_args = [config['mpiexec'], '-n', config['mpi_np']]
openmc.run(openmc_exec=config['exe'], mpi_args=mpi_args)
openmc.run(openmc_exec=config['exe'], mpi_args=mpi_args,
event_based=config['event'])
else:
openmc.run(openmc_exec=config['exe'])
openmc.run(openmc_exec=config['exe'], event_based=config['event'])
def _test_output_created(self):
"""Make sure statepoint.* and tallies.out have been created."""

17
tools/ci/travis-script.sh
View file

@ -1,9 +1,18 @@
#!/bin/bash
set -ex
# Run regression and unit tests
# Argument List
args=" "
# Check for MPI
if [[ $MPI == 'y' ]]; then
pytest --cov=openmc -v --mpi tests
else
pytest --cov=openmc -v tests
args="${args} --mpi "
fi
# Check for event-based
if [[ $EVENT == 'y' ]]; then
args="${args} --event "
fi
# Run regression and unit tests
pytest --cov=openmc -v $args tests