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

use base class to handle layout of particle data

Browse source
This commit is contained in:
Gavin Ridley 2021-04-16 15:35:33 -04:00
parent 5cf8482d9a
commit 2af4c9cd92
39 changed files with 834 additions and 795 deletions
Download patch file Download diff file Expand all files Collapse all files

10
src/bank.cpp
View file

@ -15,15 +15,15 @@ namespace openmc {
namespace simulation {
std::vector<Particle::Bank> source_bank;
std::vector<ParticleBank> source_bank;
SharedArray<Particle::Bank> surf_source_bank;
SharedArray<ParticleBank> surf_source_bank;
// 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 using SharedArray's special thread_safe_append() function.
SharedArray<Particle::Bank> fission_bank;
SharedArray<ParticleBank> 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
@ -79,8 +79,8 @@ void sort_fission_bank()
// 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;
ParticleBank* sorted_bank;
std::vector<ParticleBank> sorted_bank_holder;
// If there is not enough space, allocate a temporary vector and point to it
if (simulation::fission_bank.size() > simulation::fission_bank.capacity() / 2) {

6
src/bremsstrahlung.cpp
View file

@ -31,13 +31,13 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
if (p.material() == MATERIAL_VOID)
return;
int photon = static_cast<int>(Particle::Type::photon);
int photon = static_cast<int>(ParticleType::photon);
if (p.E() < settings::energy_cutoff[photon])
return;
// Get bremsstrahlung data for this material and particle type
BremsstrahlungData* mat;
if (p.type() == Particle::Type::positron) {
if (p.type() == ParticleType::positron) {
mat = &model::materials[p.material()]->ttb_->positron;
} else {
mat = &model::materials[p.material()]->ttb_->electron;
@ -110,7 +110,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
if (w > settings::energy_cutoff[photon]) {
// Create secondary photon
p.create_secondary(p.wgt(), p.u(), w, Particle::Type::photon);
p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon);
*E_lost += w;
}
}

2
src/eigenvalue.cpp
View file

@ -137,7 +137,7 @@ void synchronize_bank()
// Allocate temporary source bank -- we don't really know how many fission
// sites were created, so overallocate by a factor of 3
int64_t index_temp = 0;
std::vector<Particle::Bank> temp_sites(3*simulation::work_per_rank);
std::vector<ParticleBank> temp_sites(3 * simulation::work_per_rank);
for (int64_t i = 0; i < simulation::fission_bank.size(); i++ ) {
const auto& site = simulation::fission_bank[i];

2
src/initialize.cpp
View file

@ -132,7 +132,7 @@ void initialize_mpi(MPI_Comm intracomm)
mpi::master = (mpi::rank == 0);
// Create bank datatype
Particle::Bank b;
ParticleBank b;
MPI_Aint disp[9];
MPI_Get_address(&b.r, &disp[0]);
MPI_Get_address(&b.u, &disp[1]);

6
src/material.cpp
View file

@ -753,9 +753,9 @@ void Material::calculate_xs(Particle& p) const
p.macro_xs().fission = 0.0;
p.macro_xs().nu_fission = 0.0;
if (p.type() == Particle::Type::neutron) {
if (p.type() == ParticleType::neutron) {
this->calculate_neutron_xs(p);
} else if (p.type() == Particle::Type::photon) {
} else if (p.type() == ParticleType::photon) {
this->calculate_photon_xs(p);
}
}
@ -763,7 +763,7 @@ void Material::calculate_xs(Particle& p) const
void Material::calculate_neutron_xs(Particle& p) const
{
// Find energy index on energy grid
int neutron = static_cast<int>(Particle::Type::neutron);
int neutron = static_cast<int>(ParticleType::neutron);
int i_grid =
std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing;

12
src/mesh.cpp
View file

@ -283,10 +283,8 @@ int StructuredMesh::n_surface_bins() const
return 4 * n_dimension_ * n_bins();
}
xt::xtensor<double, 1>
StructuredMesh::count_sites(const Particle::Bank* bank,
int64_t length,
bool* outside) const
xt::xtensor<double, 1> StructuredMesh::count_sites(
const ParticleBank* bank, int64_t length, bool* outside) const
{
// Determine shape of array for counts
std::size_t m = this->n_bins();
@ -989,10 +987,8 @@ void RegularMesh::to_hdf5(hid_t group) const
close_group(mesh_group);
}
xt::xtensor<double, 1>
RegularMesh::count_sites(const Particle::Bank* bank,
int64_t length,
bool* outside) const
xt::xtensor<double, 1> RegularMesh::count_sites(
const ParticleBank* bank, int64_t length, bool* outside) const
{
// Determine shape of array for counts
std::size_t m = this->n_bins();

2
src/mgxs_interface.cpp
View file

@ -230,7 +230,7 @@ void MgxsInterface::read_header(const std::string& path_cross_sections)
void put_mgxs_header_data_to_globals()
{
// Get the minimum and maximum energies
int neutron = static_cast<int>(Particle::Type::neutron);
int neutron = static_cast<int>(ParticleType::neutron);
data::energy_min[neutron] = data::mg.energy_bins_.back();
data::energy_max[neutron] = data::mg.energy_bins_.front();

7
src/nuclide.cpp
View file

@ -327,7 +327,7 @@ void Nuclide::create_derived(const Function1D* prompt_photons, const Function1D*
auto xs = xt::adapt(rx->xs_[t].value);
for (const auto& p : rx->products_) {
if (p.particle_ == Particle::Type::photon) {
if (p.particle_ == ParticleType::photon) {
auto pprod = xt::view(xs_[t], xt::range(j, j+n), XS_PHOTON_PROD);
for (int k = 0; k < n; ++k) {
double E = grid_[t].energy[k+j];
@ -445,7 +445,7 @@ void Nuclide::create_derived(const Function1D* prompt_photons, const Function1D*
void Nuclide::init_grid()
{
int neutron = static_cast<int>(Particle::Type::neutron);
int neutron = static_cast<int>(ParticleType::neutron);
double E_min = data::energy_min[neutron];
double E_max = data::energy_max[neutron];
int M = settings::n_log_bins;
@ -494,7 +494,8 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const
for (int i = 1; i < rx->products_.size(); ++i) {
// Skip any non-neutron products
const auto& product = rx->products_[i];
if (product.particle_ != Particle::Type::neutron) continue;
if (product.particle_ != ParticleType::neutron)
continue;
// Evaluate yield
if (product.emission_mode_ == EmissionMode::delayed) {

8
src/output.cpp
View file

@ -147,16 +147,16 @@ void print_particle(Particle& p)
{
// Display particle type and ID.
switch (p.type()) {
case Particle::Type::neutron:
case ParticleType::neutron:
fmt::print("Neutron ");
break;
case Particle::Type::photon:
case ParticleType::photon:
fmt::print("Photon ");
break;
case Particle::Type::electron:
case ParticleType::electron:
fmt::print("Electron ");
break;
case Particle::Type::positron:
case ParticleType::positron:
fmt::print("Positron ");
break;
default:

423
src/particle.cpp
View file

@ -38,9 +38,9 @@ namespace openmc {
void
LocalCoord::rotate(const std::vector<double>& rotation)
{
this->r = this->r.rotate(rotation);
this->u = this->u.rotate(rotation);
this->rotated = true;
r = r.rotate(rotation);
u = u.rotate(rotation);
rotated = true;
}
void
@ -59,141 +59,118 @@ LocalCoord::reset()
// Particle implementation
//==============================================================================
Particle::Particle()
void Particle::create_secondary(
double wgt, Direction u, double E, ParticleType type)
{
// Create and clear coordinate levels
coord_.resize(model::n_coord_levels);
cell_last_.resize(model::n_coord_levels);
clear();
secondary_bank().emplace_back();
for (int& n : n_delayed_bank_) {
n = 0;
}
// Create microscopic cross section caches
neutron_xs_.resize(data::nuclides.size());
photon_xs_.resize(data::elements.size());
}
void
Particle::clear()
{
// Reset any coordinate levels
for (auto& level : coord_) level.reset();
n_coord_ = 1;
}
void
Particle::create_secondary(double wgt, Direction u, double E, Type type)
{
secondary_bank_.emplace_back();
auto& bank {secondary_bank_.back()};
auto& bank {secondary_bank().back()};
bank.particle = type;
bank.wgt = wgt;
bank.r = this->r();
bank.r = r();
bank.u = u;
bank.E = settings::run_CE ? E : g_;
bank.E = settings::run_CE ? E : g();
n_bank_second_ += 1;
n_bank_second() += 1;
}
void
Particle::from_source(const Bank* src)
void Particle::from_source(const ParticleBank* src)
{
// Reset some attributes
this->clear();
alive_ = true;
surface_ = 0;
cell_born_ = C_NONE;
material_ = C_NONE;
n_collision_ = 0;
fission_ = false;
std::fill(flux_derivs_.begin(), flux_derivs_.end(), 0.0);
clear();
alive() = true;
surface() = 0;
cell_born() = C_NONE;
material() = C_NONE;
n_collision() = 0;
fission() = false;
zero_flux_derivs();
// Copy attributes from source bank site
type_ = src->particle;
wgt_ = src->wgt;
wgt_last_ = src->wgt;
this->r() = src->r;
this->u() = src->u;
r_last_current_ = src->r;
r_last_ = src->r;
u_last_ = src->u;
type() = src->particle;
wgt() = src->wgt;
wgt_last() = src->wgt;
r() = src->r;
u() = src->u;
r_last_current() = src->r;
r_last() = src->r;
u_last() = src->u;
if (settings::run_CE) {
E_ = src->E;
g_ = 0;
E() = src->E;
g() = 0;
} else {
g_ = static_cast<int>(src->E);
g_last_ = static_cast<int>(src->E);
E_ = data::mg.energy_bin_avg_[g_];
g() = static_cast<int>(src->E);
g_last() = static_cast<int>(src->E);
E() = data::mg.energy_bin_avg_[g()];
}
E_last_ = E_;
E_last() = E();
}
void
Particle::event_calculate_xs()
{
// Set the random number stream
if (type_ == Particle::Type::neutron) {
stream_ = STREAM_TRACKING;
if (type() == ParticleType::neutron) {
stream() = STREAM_TRACKING;
} else {
stream_ = STREAM_PHOTON;
stream() = STREAM_PHOTON;
}
// Store pre-collision particle properties
wgt_last_ = wgt_;
E_last_ = E_;
u_last_ = this->u();
r_last_ = this->r();
wgt_last() = wgt();
E_last() = E();
u_last() = u();
r_last() = r();
// Reset event variables
event_ = TallyEvent::KILL;
event_nuclide_ = NUCLIDE_NONE;
event_mt_ = REACTION_NONE;
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 (coord(n_coord() - 1).cell == C_NONE) {
if (!exhaustive_find_cell(*this)) {
this->mark_as_lost("Could not find the cell containing particle "
+ std::to_string(id_));
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;
if (cell_born() == C_NONE)
cell_born() = coord(n_coord() - 1).cell;
}
// Write particle track.
if (write_track_) write_particle_track(*this);
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 (material() != MATERIAL_VOID) {
if (settings::run_CE) {
if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) {
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);
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);
data::mg.macro_xs_[material()].calculate_xs(*this);
// Update the particle's group while we know we are multi-group
g_last_ = g_;
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;
macro_xs().total = 0.0;
macro_xs().absorption = 0.0;
macro_xs().fission = 0.0;
macro_xs().nu_fission = 0.0;
}
}
@ -201,24 +178,23 @@ void
Particle::event_advance()
{
// Find the distance to the nearest boundary
boundary_ = distance_to_boundary(*this);
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;
if (type() == ParticleType::electron || type() == ParticleType::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;
collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
}
// Select smaller of the two distances
double distance = std::min(boundary_.distance, collision_distance_);
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;
for (int j = 0; j < n_coord(); ++j) {
coord(j).r += distance * coord(j).u;
}
// Score track-length tallies
@ -228,8 +204,8 @@ Particle::event_advance()
// Score track-length estimate of k-eff
if (settings::run_mode == RunMode::EIGENVALUE &&
type_ == Particle::Type::neutron) {
keff_tally_tracklength_ += wgt_ * distance * macro_xs_.nu_fission;
type() == ParticleType::neutron) {
keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
}
// Score flux derivative accumulators for differential tallies.
@ -242,25 +218,25 @@ void
Particle::event_cross_surface()
{
// Set surface that particle is on and adjust coordinate levels
surface_ = boundary_.surface_index;
n_coord_ = boundary_.coord_level;
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;
for (int j = 0; j < n_coord(); ++j) {
cell_last(j) = coord(j).cell;
}
n_coord_last_ = n_coord_;
n_coord_last() = n_coord();
if (boundary_.lattice_translation[0] != 0 ||
boundary_.lattice_translation[1] != 0 ||
boundary_.lattice_translation[2] != 0) {
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;
cross_lattice(*this, boundary());
event() = TallyEvent::LATTICE;
} else {
// Particle crosses surface
this->cross_surface();
event_ = TallyEvent::SURFACE;
cross_surface();
event() = TallyEvent::SURFACE;
}
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
@ -273,9 +249,8 @@ Particle::event_collide()
{
// Score collision estimate of keff
if (settings::run_mode == RunMode::EIGENVALUE &&
type_ == Particle::Type::neutron) {
keff_tally_collision_ += wgt_ * macro_xs_.nu_fission
/ macro_xs_.total;
type() == ParticleType::neutron) {
keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
}
// Score surface current tallies -- this has to be done before the collision
@ -286,7 +261,7 @@ Particle::event_collide()
score_surface_tally(*this, model::active_meshsurf_tallies);
// Clear surface component
surface_ = 0;
surface() = 0;
if (settings::run_CE) {
collision(*this);
@ -307,32 +282,32 @@ Particle::event_collide()
}
// Reset banked weight during collision
n_bank_ = 0;
n_bank_second_ = 0;
wgt_bank_ = 0.0;
for (int& v : n_delayed_bank_) v = 0;
n_bank() = 0;
n_bank_second() = 0;
wgt_bank() = 0.0;
zero_delayed_bank();
// Reset fission logical
fission_ = false;
fission() = false;
// Save coordinates for tallying purposes
r_last_current_ = this->r();
r_last_current() = r();
// Set last material to none since cross sections will need to be
// re-evaluated
material_last_ = C_NONE;
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) {
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);
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;
coord(j + 1).u = coord(j).u;
}
}
@ -344,24 +319,26 @@ 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;
++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 (!alive()) {
// If no secondary particles, break out of event loop
if (secondary_bank_.empty()) return;
if (secondary_bank().empty())
return;
this->from_source(&secondary_bank_.back());
secondary_bank_.pop_back();
n_event_ = 0;
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);
if (write_track())
add_particle_track(*this);
}
}
@ -373,32 +350,32 @@ Particle::event_death()
#endif
// Finish particle track output.
if (write_track_) {
if (write_track()) {
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_;
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;
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 == RunMode::EIGENVALUE) {
int64_t offset = id_ - 1 - simulation::work_index[mpi::rank];
simulation::progeny_per_particle[offset] = n_progeny_;
int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
simulation::progeny_per_particle[offset] = n_progeny();
}
}
@ -406,24 +383,24 @@ Particle::event_death()
void
Particle::cross_surface()
{
int i_surface = std::abs(surface_);
int i_surface = std::abs(surface());
// TODO: off-by-one
const auto& surf {model::surfaces[i_surface - 1].get()};
if (settings::verbosity >= 10 || trace_) {
if (settings::verbosity >= 10 || trace()) {
write_message(1, " Crossing surface {}", surf->id_);
}
if (surf->surf_source_ && simulation::current_batch == settings::n_batches) {
Particle::Bank site;
site.r = this->r();
site.u = this->u();
site.E = this->E_;
site.wgt = this->wgt_;
site.delayed_group = this->delayed_group_;
ParticleBank site;
site.r = r();
site.u = u();
site.E = E();
site.wgt = wgt();
site.delayed_group = delayed_group();
site.surf_id = surf->id_;
site.particle = this->type_;
site.parent_id = this->id_;
site.progeny_id = this->n_progeny_;
site.particle = type();
site.parent_id = id();
site.progeny_id = n_progeny();
int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
}
@ -438,18 +415,19 @@ Particle::cross_surface()
#ifdef DAGMC
if (settings::dagmc) {
auto cellp = dynamic_cast<DAGCell*>(model::cells[cell_last_[0]].get());
auto cellp = dynamic_cast<DAGCell*>(model::cells[cell_last(0)].get());
// TODO: off-by-one
auto surfp = dynamic_cast<DAGSurface*>(model::surfaces[std::abs(surface_) - 1].get());
auto surfp =
dynamic_cast<DAGSurface*>(model::surfaces[std::abs(surface()) - 1].get());
int32_t i_cell = next_cell(cellp, surfp) - 1;
// save material and temp
material_last_ = material_;
sqrtkT_last_ = sqrtkT_;
material_last() = material();
sqrtkT_last() = sqrtkT();
// set new cell value
coord_[0].cell = i_cell;
cell_instance_ = 0;
material_ = model::cells[i_cell]->material_[0];
sqrtkT_ = model::cells[i_cell]->sqrtkT_[0];
coord(0).cell = i_cell;
cell_instance() = 0;
material() = model::cells[i_cell]->material_[0];
sqrtkT() = model::cells[i_cell]->sqrtkT_[0];
return;
}
#endif
@ -461,8 +439,8 @@ Particle::cross_surface()
// COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
// Remove lower coordinate levels and assignment of surface
surface_ = 0;
n_coord_ = 1;
surface() = 0;
n_coord() = 1;
bool found = exhaustive_find_cell(*this);
if (settings::run_mode != RunMode::PLOTTING && (!found)) {
@ -471,16 +449,16 @@ Particle::cross_surface()
// the particle is really traveling tangent to a surface, if we move it
// forward a tiny bit it should fix the problem.
n_coord_ = 1;
this->r() += TINY_BIT * this->u();
n_coord() = 1;
r() += TINY_BIT * u();
// Couldn't find next cell anywhere! This probably means there is an actual
// undefined region in the geometry.
if (!exhaustive_find_cell(*this)) {
this->mark_as_lost("After particle " + std::to_string(id_) +
" crossed surface " + std::to_string(surf->id_) +
" it could not be located in any cell and it did not leak.");
mark_as_lost("After particle " + std::to_string(id()) +
" crossed surface " + std::to_string(surf->id_) +
" it could not be located in any cell and it did not leak.");
return;
}
}
@ -490,7 +468,7 @@ void
Particle::cross_vacuum_bc(const Surface& surf)
{
// Kill the particle
alive_ = false;
alive() = false;
// Score any surface current tallies -- note that the particle is moved
// forward slightly so that if the mesh boundary is on the surface, it is
@ -500,15 +478,15 @@ Particle::cross_vacuum_bc(const Surface& surf)
// TODO: Find a better solution to score surface currents than
// physically moving the particle forward slightly
this->r() += TINY_BIT * this->u();
r() += TINY_BIT * u();
score_surface_tally(*this, model::active_meshsurf_tallies);
}
// Score to global leakage tally
keff_tally_leakage_ += wgt_;
keff_tally_leakage() += wgt();
// Display message
if (settings::verbosity >= 10 || trace_) {
if (settings::verbosity >= 10 || trace()) {
write_message(1, " Leaked out of surface {}", surf.id_);
}
}
@ -517,9 +495,9 @@ void
Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
{
// Do not handle reflective boundary conditions on lower universes
if (n_coord_ != 1) {
this->mark_as_lost("Cannot reflect particle " + std::to_string(id_) +
" off surface in a lower universe.");
if (n_coord() != 1) {
mark_as_lost("Cannot reflect particle " + std::to_string(id()) +
" off surface in a lower universe.");
return;
}
@ -536,36 +514,36 @@ Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
if (!model::active_meshsurf_tallies.empty()) {
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
this->r() -= TINY_BIT * u();
score_surface_tally(*this, model::active_meshsurf_tallies);
this->r() = r;
}
// Set the new particle direction
this->u() = new_u;
u() = new_u;
// Reassign particle's cell and surface
coord_[0].cell = cell_last_[n_coord_last_ - 1];
surface_ = -surface_;
coord(0).cell = cell_last(n_coord_last() - 1);
surface() = -surface();
// If a reflective surface is coincident with a lattice or universe
// boundary, it is necessary to redetermine the particle's coordinates in
// the lower universes.
// (unless we're using a dagmc model, which has exactly one universe)
if (!settings::dagmc) {
n_coord_ = 1;
n_coord() = 1;
if (!neighbor_list_find_cell(*this)) {
this->mark_as_lost("Couldn't find particle after reflecting from surface "
+ std::to_string(surf.id_) + ".");
mark_as_lost("Couldn't find particle after reflecting from surface " +
std::to_string(surf.id_) + ".");
return;
}
}
// Set previous coordinate going slightly past surface crossing
r_last_current_ = this->r() + TINY_BIT*this->u();
r_last_current() = r() + TINY_BIT * u();
// Diagnostic message
if (settings::verbosity >= 10 || trace_) {
if (settings::verbosity >= 10 || trace()) {
write_message(1, " Reflected from surface {}", surf.id_);
}
}
@ -575,8 +553,9 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r,
Direction new_u, int new_surface)
{
// Do not handle periodic boundary conditions on lower universes
if (n_coord_ != 1) {
this->mark_as_lost("Cannot transfer particle " + std::to_string(id_) +
if (n_coord() != 1) {
mark_as_lost(
"Cannot transfer particle " + std::to_string(id()) +
" across surface in a lower universe. Boundary conditions must be "
"applied to root universe.");
return;
@ -587,7 +566,7 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r,
// case the surface crossing is coincident with a mesh boundary
if (!model::active_meshsurf_tallies.empty()) {
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
this->r() -= TINY_BIT * u();
score_surface_tally(*this, model::active_meshsurf_tallies);
this->r() = r;
}
@ -597,23 +576,25 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r,
u() = new_u;
// Reassign particle's surface
surface_ = new_surface;
surface() = new_surface;
// Figure out what cell particle is in now
n_coord_ = 1;
n_coord() = 1;
if (!neighbor_list_find_cell(*this)) {
this->mark_as_lost("Couldn't find particle after hitting periodic "
"boundary on surface " + std::to_string(surf.id_) + ". The normal vector "
"of one periodic surface may need to be reversed.");
mark_as_lost("Couldn't find particle after hitting periodic "
"boundary on surface " +
std::to_string(surf.id_) +
". The normal vector "
"of one periodic surface may need to be reversed.");
return;
}
// Set previous coordinate going slightly past surface crossing
r_last_current_ = this->r() + TINY_BIT*this->u();
r_last_current() = r() + TINY_BIT * u();
// Diagnostic message
if (settings::verbosity >= 10 || trace_) {
if (settings::verbosity >= 10 || trace()) {
write_message(1, " Hit periodic boundary on surface {}", surf.id_);
}
}
@ -626,8 +607,8 @@ Particle::mark_as_lost(const char* message)
write_restart();
// Increment number of lost particles
alive_ = false;
#pragma omp atomic
alive() = false;
#pragma omp atomic
simulation::n_lost_particles += 1;
// Count the total number of simulated particles (on this processor)
@ -650,9 +631,9 @@ Particle::write_restart() const
// Set up file name
auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,
simulation::current_batch, id_);
simulation::current_batch, id());
#pragma omp critical (WriteParticleRestart)
#pragma omp critical (WriteParticleRestart)
{
// Create file
hid_t file_id = file_open(filename, 'w');
@ -683,10 +664,10 @@ Particle::write_restart() const
default:
break;
}
write_dataset(file_id, "id", id_);
write_dataset(file_id, "type", static_cast<int>(type_));
write_dataset(file_id, "id", id());
write_dataset(file_id, "type", static_cast<int>(type()));
int64_t i = current_work_;
int64_t i = current_work();
if (settings::run_mode == RunMode::EIGENVALUE) {
//take source data from primary bank for eigenvalue simulation
write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt);
@ -711,31 +692,31 @@ Particle::write_restart() const
} // #pragma omp critical
}
std::string particle_type_to_str(Particle::Type type)
std::string particle_type_to_str(ParticleType type)
{
switch (type) {
case Particle::Type::neutron:
return "neutron";
case Particle::Type::photon:
return "photon";
case Particle::Type::electron:
return "electron";
case Particle::Type::positron:
return "positron";
case ParticleType::neutron:
return "neutron";
case ParticleType::photon:
return "photon";
case ParticleType::electron:
return "electron";
case ParticleType::positron:
return "positron";
}
UNREACHABLE();
}
Particle::Type str_to_particle_type(std::string str)
ParticleType str_to_particle_type(std::string str)
{
if (str == "neutron") {
return Particle::Type::neutron;
return ParticleType::neutron;
} else if (str == "photon") {
return Particle::Type::photon;
return ParticleType::photon;
} else if (str == "electron") {
return Particle::Type::electron;
return ParticleType::electron;
} else if (str == "positron") {
return Particle::Type::positron;
return ParticleType::positron;
} else {
throw std::invalid_argument{fmt::format("Invalid particle name: {}", str)};
}

35
src/particle_data.cpp Normal file
View file

@ -0,0 +1,35 @@
#include "openmc/particle_data.h"
#include "openmc/geometry.h"
#include "openmc/nuclide.h"
#include "openmc/photon.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
namespace openmc {
ParticleData::ParticleData()
{
// Create and clear coordinate levels
coord_.resize(model::n_coord_levels);
cell_last_.resize(model::n_coord_levels);
clear();
zero_delayed_bank();
// Every particle starts with no accumulated flux derivative.
if (!model::active_tallies.empty()) {
flux_derivs_.resize(model::tally_derivs.size());
zero_flux_derivs();
}
// Allocate space for tally filter matches
filter_matches_.resize(model::tally_filters.size());
// Create microscopic cross section caches
neutron_xs_.resize(data::nuclides.size());
photon_xs_.resize(data::elements.size());
}
} // namespace openmc

12
src/particle_restart.cpp
View file

@ -44,7 +44,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
read_dataset(file_id, "id", p.id());
int type;
read_dataset(file_id, "type", type);
p.type() = static_cast<Particle::Type>(type);
p.type() = static_cast<ParticleType>(type);
read_dataset(file_id, "weight", p.wgt());
read_dataset(file_id, "energy", p.E());
read_dataset(file_id, "xyz", p.r());
@ -116,16 +116,6 @@ void run_particle_restart()
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 (TODO shouldn't this be in the
// particle constructor, instead?)
p.filter_matches().resize(model::tally_filters.size());
// Transport neutron
transport_history_based_single_particle(p);

8
src/photon.cpp
View file

@ -231,7 +231,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
close_group(rgroup);
// Truncate the bremsstrahlung data at the cutoff energy
int photon = static_cast<int>(Particle::Type::photon);
int photon = static_cast<int>(ParticleType::photon);
const auto& E {electron_energy};
double cutoff = settings::energy_cutoff[photon];
if (cutoff > E(0)) {
@ -669,7 +669,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
double E = shell.binding_energy;
p.create_secondary(p.wgt(), u, E, Particle::Type::photon);
p.create_secondary(p.wgt(), u, E, ParticleType::photon);
return;
}
@ -701,7 +701,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Non-radiative transition -- Auger/Coster-Kronig effect
// Create auger electron
p.create_secondary(p.wgt(), u, E, Particle::Type::electron);
p.create_secondary(p.wgt(), u, E, ParticleType::electron);
// Fill hole left by emitted auger electron
int i_hole = shell_map_.at(secondary);
@ -711,7 +711,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Radiative transition -- get X-ray energy
// Create fluorescent photon
p.create_secondary(p.wgt(), u, E, Particle::Type::photon);
p.create_secondary(p.wgt(), u, E, ParticleType::photon);
}
// Fill hole created by electron transitioning to the photoelectron hole

48
src/physics.cpp
View file

@ -40,16 +40,16 @@ void collision(Particle& p)
// Sample reaction for the material the particle is in
switch (p.type()) {
case Particle::Type::neutron:
case ParticleType::neutron:
sample_neutron_reaction(p);
break;
case Particle::Type::photon:
case ParticleType::photon:
sample_photon_reaction(p);
break;
case Particle::Type::electron:
case ParticleType::electron:
sample_electron_reaction(p);
break;
case Particle::Type::positron:
case ParticleType::positron:
sample_positron_reaction(p);
break;
}
@ -66,11 +66,11 @@ void collision(Particle& p)
std::string msg;
if (p.event() == TallyEvent::KILL) {
msg = fmt::format(" Killed. Energy = {} eV.", p.E());
} else if (p.type() == Particle::Type::neutron) {
} else if (p.type() == ParticleType::neutron) {
msg = fmt::format(" {} with {}. Energy = {} eV.",
reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_,
p.E());
} else if (p.type() == Particle::Type::photon) {
} else if (p.type() == ParticleType::photon) {
msg = fmt::format(" {} with {}. Energy = {} eV.",
reaction_name(p.event_mt()),
to_element(data::nuclides[p.event_nuclide()]->name_), p.E());
@ -187,9 +187,9 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
for (int i = 0; i < nu; ++i) {
// Initialize fission site object with particle data
Particle::Bank site;
ParticleBank site;
site.r = p.r();
site.particle = Particle::Type::neutron;
site.particle = ParticleType::neutron;
site.wgt = 1. / weight;
site.parent_id = p.id();
site.progeny_id = p.n_progeny()++;
@ -221,7 +221,7 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
// Write fission particles to nuBank
p.nu_bank().emplace_back();
Particle::NuBank* nu_bank_entry = &p.nu_bank().back();
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;
@ -251,7 +251,7 @@ void sample_photon_reaction(Particle& p)
// Kill photon if below energy cutoff -- an extra check is made here because
// photons with energy below the cutoff may have been produced by neutrons
// reactions or atomic relaxation
int photon = static_cast<int>(Particle::Type::photon);
int photon = static_cast<int>(ParticleType::photon);
if (p.E() < settings::energy_cutoff[photon]) {
p.E() = 0.0;
p.alive() = false;
@ -300,12 +300,12 @@ void sample_photon_reaction(Particle& p)
// Create Compton electron
double phi = 2.0*PI*prn(p.current_seed());
double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b;
int electron = static_cast<int>(Particle::Type::electron);
int electron = static_cast<int>(ParticleType::electron);
if (E_electron >= settings::energy_cutoff[electron]) {
double mu_electron = (alpha - alpha_out*mu)
/ std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed());
p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron);
p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
}
// TODO: Compton subshell data does not match atomic relaxation data
@ -366,7 +366,7 @@ void sample_photon_reaction(Particle& p)
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
// Create secondary electron
p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron);
p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
// Allow electrons to fill orbital and produce auger electrons
// and fluorescent photons
@ -391,11 +391,11 @@ void sample_photon_reaction(Particle& p)
// Create secondary electron
Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed());
p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron);
p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
// Create secondary positron
u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed());
p.create_secondary(p.wgt(), u, E_positron, Particle::Type::positron);
p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron);
p.event() = TallyEvent::ABSORB;
p.event_mt() = PAIR_PROD;
@ -436,8 +436,8 @@ void sample_positron_reaction(Particle& p)
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
// Create annihilation photon pair traveling in opposite directions
p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, Particle::Type::photon);
p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, Particle::Type::photon);
p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon);
p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon);
p.E() = 0.0;
p.alive() = false;
@ -569,7 +569,7 @@ void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product
+ f*(rx->xs_[i_temp].value[i_grid - threshold + 1]));
for (int j = 0; j < rx->products_.size(); ++j) {
if (rx->products_[j].particle_ == Particle::Type::photon) {
if (rx->products_[j].particle_ == ParticleType::photon) {
// For fission, artificially increase the photon yield to account
// for delayed photons
double f = 1.0;
@ -1014,7 +1014,8 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
return vt * rotate_angle(u, mu, nullptr, seed);
}
void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Particle::Bank* site, uint64_t* seed)
void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in,
ParticleBank* site, uint64_t* seed)
{
// Sample cosine of angle -- fission neutrons are always emitted
// isotropically. Sometimes in ACE data, fission reactions actually have
@ -1066,7 +1067,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part
rx.products_[group].sample(E_in, site->E, mu, seed);
// resample if energy is greater than maximum neutron energy
constexpr int neutron = static_cast<int>(Particle::Type::neutron);
constexpr int neutron = static_cast<int>(ParticleType::neutron);
if (site->E < data::energy_max[neutron]) break;
// check for large number of resamples
@ -1091,7 +1092,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part
rx.products_[0].sample(E_in, site->E, mu, seed);
// resample if energy is greater than maximum neutron energy
constexpr int neutron = static_cast<int>(Particle::Type::neutron);
constexpr int neutron = static_cast<int>(ParticleType::neutron);
if (site->E < data::energy_max[neutron]) break;
// check for large number of resamples
@ -1146,7 +1147,7 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
if (std::floor(yield) == yield) {
// If yield is integral, create exactly that many secondary particles
for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
p.create_secondary(p.wgt(), p.u(), p.E(), Particle::Type::neutron);
p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron);
}
} else {
// Otherwise, change weight of particle based on yield
@ -1191,8 +1192,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
}
// Create the secondary photon
p.create_secondary(wgt, u, E, Particle::Type::photon);
p.create_secondary(wgt, u, E, ParticleType::photon);
}
}

7
src/physics_mg.cpp
View file

@ -13,6 +13,7 @@
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/particle.h"
#include "openmc/physics_common.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
@ -126,9 +127,9 @@ create_fission_sites(Particle& p)
for (int i = 0; i < nu; ++i) {
// Initialize fission site object with particle data
Particle::Bank site;
ParticleBank site;
site.r = p.r();
site.particle = Particle::Type::neutron;
site.particle = ParticleType::neutron;
site.wgt = 1. / weight;
site.parent_id = p.id();
site.progeny_id = p.n_progeny()++;
@ -179,7 +180,7 @@ create_fission_sites(Particle& p)
// Write fission particles to nuBank
p.nu_bank().emplace_back();
Particle::NuBank* nu_bank_entry = &p.nu_bank().back();
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;

2
src/reaction.cpp
View file

@ -73,7 +73,7 @@ Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
// mark fission reactions so that we avoid the angle sampling.
if (is_fission(mt_)) {
for (auto& p : products_) {
if (p.particle_ == Particle::Type::neutron) {
if (p.particle_ == ParticleType::neutron) {
for (auto& d : p.distribution_) {
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (d_) d_->fission() = true;

2
src/reaction_product.cpp
View file

@ -42,7 +42,7 @@ ReactionProduct::ReactionProduct(hid_t group)
if (emission_mode_ == EmissionMode::delayed) {
if (attribute_exists(group, "decay_rate")) {
read_attribute(group, "decay_rate", decay_rate_);
} else if (particle_ == Particle::Type::neutron) {
} else if (particle_ == ParticleType::neutron) {
warning(fmt::format("Decay rate doesn't exist for delayed neutron "
"emission ({}).", object_name(group)));
}

25
src/simulation.cpp
View file

@ -511,11 +511,6 @@ void initialize_history(Particle& p, int64_t index_source)
#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) {
p.invalidate_neutron_xs();
@ -524,16 +519,6 @@ void initialize_history_partial(Particle& p)
// 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()
@ -573,7 +558,7 @@ void initialize_data()
data::energy_min = {0.0, 0.0};
for (const auto& nuc : data::nuclides) {
if (nuc->grid_.size() >= 1) {
int neutron = static_cast<int>(Particle::Type::neutron);
int neutron = static_cast<int>(ParticleType::neutron);
data::energy_min[neutron] = std::max(data::energy_min[neutron],
nuc->grid_[0].energy.front());
data::energy_max[neutron] = std::min(data::energy_max[neutron],
@ -584,7 +569,7 @@ void initialize_data()
if (settings::photon_transport) {
for (const auto& elem : data::elements) {
if (elem->energy_.size() >= 1) {
int photon = static_cast<int>(Particle::Type::photon);
int photon = static_cast<int>(ParticleType::photon);
int n = elem->energy_.size();
data::energy_min[photon] = std::max(data::energy_min[photon],
std::exp(elem->energy_(1)));
@ -597,7 +582,7 @@ void initialize_data()
// Determine if minimum/maximum energy for bremsstrahlung is greater/less
// than the current minimum/maximum
if (data::ttb_e_grid.size() >= 1) {
int photon = static_cast<int>(Particle::Type::photon);
int photon = static_cast<int>(ParticleType::photon);
int n_e = data::ttb_e_grid.size();
data::energy_min[photon] = std::max(data::energy_min[photon],
std::exp(data::ttb_e_grid(1)));
@ -613,7 +598,7 @@ void initialize_data()
// grid has not been allocated
if (nuc->grid_.size() > 0) {
double max_E = nuc->grid_[0].energy.back();
int neutron = static_cast<int>(Particle::Type::neutron);
int neutron = static_cast<int>(ParticleType::neutron);
if (max_E == data::energy_max[neutron]) {
write_message(7, "Maximum neutron transport energy: {} eV for {}",
data::energy_max[neutron], nuc->name_);
@ -630,7 +615,7 @@ void initialize_data()
for (auto& nuc : data::nuclides) {
nuc->init_grid();
}
int neutron = static_cast<int>(Particle::Type::neutron);
int neutron = static_cast<int>(ParticleType::neutron);
simulation::log_spacing = std::log(data::energy_max[neutron] /
data::energy_min[neutron]) / settings::n_log_bins;
}

14
src/source.cpp
View file

@ -56,9 +56,9 @@ IndependentSource::IndependentSource(pugi::xml_node node)
if (check_for_node(node, "particle")) {
auto temp_str = get_node_value(node, "particle", true, true);
if (temp_str == "neutron") {
particle_ = Particle::Type::neutron;
particle_ = ParticleType::neutron;
} else if (temp_str == "photon") {
particle_ = Particle::Type::photon;
particle_ = ParticleType::photon;
settings::photon_transport = true;
} else {
fatal_error(std::string("Unknown source particle type: ") + temp_str);
@ -141,9 +141,9 @@ IndependentSource::IndependentSource(pugi::xml_node node)
}
}
Particle::Bank IndependentSource::sample(uint64_t* seed) const
ParticleBank IndependentSource::sample(uint64_t* seed) const
{
Particle::Bank site;
ParticleBank site;
// Set weight to one by default
site.wgt = 1.0;
@ -263,7 +263,7 @@ FileSource::FileSource(std::string path)
file_close(file_id);
}
Particle::Bank FileSource::sample(uint64_t* seed) const
ParticleBank FileSource::sample(uint64_t* seed) const
{
size_t i_site = sites_.size()*prn(seed);
return sites_[i_site];
@ -349,7 +349,7 @@ void initialize_source()
}
}
Particle::Bank sample_external_source(uint64_t* seed)
ParticleBank sample_external_source(uint64_t* seed)
{
// Determine total source strength
double total_strength = 0.0;
@ -368,7 +368,7 @@ Particle::Bank sample_external_source(uint64_t* seed)
}
// Sample source site from i-th source distribution
Particle::Bank site {model::external_sources[i]->sample(seed)};
ParticleBank site {model::external_sources[i]->sample(seed)};
// If running in MG, convert site.E to group
if (!settings::run_CE) {

29
src/state_point.cpp
View file

@ -512,14 +512,17 @@ hid_t h5banktype() {
// - openmc/statepoint.py
// - docs/source/io_formats/statepoint.rst
// - docs/source/io_formats/source.rst
hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct Particle::Bank));
H5Tinsert(banktype, "r", HOFFSET(Particle::Bank, r), postype);
H5Tinsert(banktype, "u", HOFFSET(Particle::Bank, u), postype);
H5Tinsert(banktype, "E", HOFFSET(Particle::Bank, E), H5T_NATIVE_DOUBLE);
H5Tinsert(banktype, "wgt", HOFFSET(Particle::Bank, wgt), H5T_NATIVE_DOUBLE);
H5Tinsert(banktype, "delayed_group", HOFFSET(Particle::Bank, delayed_group), H5T_NATIVE_INT);
H5Tinsert(banktype, "surf_id", HOFFSET(Particle::Bank, surf_id), H5T_NATIVE_INT);
H5Tinsert(banktype, "particle", HOFFSET(Particle::Bank, particle), H5T_NATIVE_INT);
hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct ParticleBank));
H5Tinsert(banktype, "r", HOFFSET(ParticleBank, r), postype);
H5Tinsert(banktype, "u", HOFFSET(ParticleBank, u), postype);
H5Tinsert(banktype, "E", HOFFSET(ParticleBank, E), H5T_NATIVE_DOUBLE);
H5Tinsert(banktype, "wgt", HOFFSET(ParticleBank, wgt), H5T_NATIVE_DOUBLE);
H5Tinsert(banktype, "delayed_group", HOFFSET(ParticleBank, delayed_group),
H5T_NATIVE_INT);
H5Tinsert(
banktype, "surf_id", HOFFSET(ParticleBank, surf_id), H5T_NATIVE_INT);
H5Tinsert(
banktype, "particle", HOFFSET(ParticleBank, particle), H5T_NATIVE_INT);
H5Tclose(postype);
return banktype;
@ -595,9 +598,9 @@ write_source_bank(hid_t group_id, bool surf_source_bank)
// Set vectors for source bank and starting bank index of each process
std::vector<int64_t>* bank_index = &simulation::work_index;
std::vector<Particle::Bank>* source_bank = &simulation::source_bank;
std::vector<ParticleBank>* source_bank = &simulation::source_bank;
std::vector<int64_t> surf_source_index_vector;
std::vector<Particle::Bank> surf_source_bank_vector;
std::vector<ParticleBank> surf_source_bank_vector;
// Reset dataspace sizes and vectors for surface source bank
if (surf_source_bank) {
@ -653,7 +656,8 @@ write_source_bank(hid_t group_id, bool surf_source_bank)
// Save source bank sites since the array is overwritten below
#ifdef OPENMC_MPI
std::vector<Particle::Bank> temp_source {source_bank->begin(), source_bank->end()};
std::vector<ParticleBank> temp_source {
source_bank->begin(), source_bank->end()};
#endif
for (int i = 0; i < mpi::n_procs; ++i) {
@ -710,7 +714,8 @@ std::string dtype_member_names(hid_t dtype_id)
return names;
}
void read_source_bank(hid_t group_id, std::vector<Particle::Bank>& sites, bool distribute)
void read_source_bank(
hid_t group_id, std::vector<ParticleBank>& sites, bool distribute)
{
hid_t banktype = h5banktype();

6
src/tallies/derivative.cpp
View file

@ -109,7 +109,7 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide,
// perturbated variable.
const auto& deriv {model::tally_derivs[tally.deriv_]};
const auto flux_deriv = p.flux_derivs()[tally.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) {
@ -558,7 +558,7 @@ score_track_derivative(Particle& p, double distance)
for (auto idx = 0; idx < model::tally_derivs.size(); idx++) {
const auto& deriv = model::tally_derivs[idx];
auto& flux_deriv = p.flux_derivs()[idx];
auto& flux_deriv = p.flux_derivs(idx);
if (deriv.diff_material != material.id_) continue;
switch (deriv.variable) {
@ -606,7 +606,7 @@ void score_collision_derivative(Particle& p)
for (auto idx = 0; idx < model::tally_derivs.size(); idx++) {
const auto& deriv = model::tally_derivs[idx];
auto& flux_deriv = p.flux_derivs()[idx];
auto& flux_deriv = p.flux_derivs(idx);
if (deriv.diff_material != material.id_) continue;

7
src/tallies/filter_particle.cpp
View file

@ -11,16 +11,15 @@ ParticleFilter::from_xml(pugi::xml_node node)
{
auto particles = get_node_array<std::string>(node, "bins");
// Convert to vector of Particle::Type
std::vector<Particle::Type> types;
// Convert to vector of ParticleType
std::vector<ParticleType> types;
for (auto& p : particles) {
types.push_back(str_to_particle_type(p));
}
this->set_particles(types);
}
void
ParticleFilter::set_particles(gsl::span<Particle::Type> particles)
void ParticleFilter::set_particles(gsl::span<ParticleType> particles)
{
// Clear existing particles
particles_.clear();

2
src/tallies/tally.cpp
View file

@ -196,7 +196,7 @@ Tally::Tally(pugi::xml_node node)
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (auto p : pf->particles()) {
if (p != Particle::Type::neutron) {
if (p != ParticleType::neutron) {
warning(fmt::format("Particle filter other than NEUTRON used with "
"photon transport turned off. All tallies for particle type {}"
" will have no scores", static_cast<int>(p)));

7
src/tallies/tally_scoring.cpp
View file

@ -547,7 +547,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
// Get the pre-collision energy of the particle.
auto E = p.E_last();
using Type = Particle::Type;
using Type = ParticleType;
for (auto i = 0; i < tally.scores_.size(); ++i) {
auto score_bin = tally.scores_[i];
@ -2181,7 +2181,7 @@ void score_analog_tally_ce(Particle& p)
// Note that the heating score does NOT use the flux and will be non-zero for
// electrons/positrons.
double flux =
(p.type() == Particle::Type::neutron || p.type() == Particle::Type::photon)
(p.type() == ParticleType::neutron || p.type() == ParticleType::photon)
? 1.0
: 0.0;
@ -2356,8 +2356,7 @@ void score_collision_tally(Particle& p)
{
// Determine the collision estimate of the flux
double flux = 0.0;
if (p.type() == Particle::Type::neutron ||
p.type() == Particle::Type::photon) {
if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) {
if (!settings::survival_biasing) {
flux = p.wgt_last() / p.macro_xs().total;
} else {