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

Use references consistently in geometry functions

Browse source
This commit is contained in:
Paul Romano 2020-04-23 12:49:48 -05:00
parent 43eae72aeb
commit 9dab8bf507
15 changed files with 170 additions and 167 deletions
Download patch file Download diff file Expand all files Collapse all files

8
include/openmc/geometry.h
View file

@ -36,7 +36,7 @@ inline bool coincident(double d1, double d2) {
//! Check for overlapping cells at a particle's position.
//==============================================================================
bool check_cell_overlap(Particle* p, bool error=true);
bool check_cell_overlap(Particle& p, bool error=true);
//==============================================================================
//! Locate a particle in the geometry tree and set its geometry data fields.
@ -50,19 +50,19 @@ bool check_cell_overlap(Particle* p, bool error=true);
//! valid geometry coordinate stack.
//==============================================================================
bool find_cell(Particle* p, bool use_neighbor_lists);
bool find_cell(Particle& p, bool use_neighbor_lists);
//==============================================================================
//! Move a particle into a new lattice tile.
//==============================================================================
void cross_lattice(Particle* p, const BoundaryInfo& boundary);
void cross_lattice(Particle& p, const BoundaryInfo& boundary);
//==============================================================================
//! Find the next boundary a particle will intersect.
//==============================================================================
BoundaryInfo distance_to_boundary(Particle* p);
BoundaryInfo distance_to_boundary(Particle& p);
} // namespace openmc

2
include/openmc/output.h
View file

@ -26,7 +26,7 @@ void header(const char* msg, int level);
std::string time_stamp();
//! Display the attributes of a particle.
extern "C" void print_particle(Particle* p);
void print_particle(Particle& p);
//! Display plot information.
void print_plot();

4
include/openmc/plot.h
View file

@ -182,13 +182,13 @@ T PlotBase::get_map() const {
p.r()[in_i] = xyz[in_i] + in_pixel * x;
p.n_coord_ = 1;
// local variables
bool found_cell = find_cell(&p, 0);
bool found_cell = find_cell(p, 0);
j = p.n_coord_ - 1;
if (level >=0) {j = level + 1;}
if (found_cell) {
data.set_value(y, x, p, j);
}
if (color_overlaps_ && check_cell_overlap(&p, false)) {
if (color_overlaps_ && check_cell_overlap(p, false)) {
data.set_overlap(y, x);
}
} // inner for

4
include/openmc/simulation.h
View file

@ -64,10 +64,10 @@ void initialize_batch();
void initialize_generation();
//! Full initialization of a particle history
void initialize_history(Particle* p, int64_t index_source);
void initialize_history(Particle& p, int64_t index_source);
//! Helper function for initialize_history() that is called independently elsewhere
void initialize_history_partial(Particle* p);
void initialize_history_partial(Particle& p);
//! Finalize a batch
//!

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

@ -53,13 +53,13 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide,
//! further tallies are scored.
//
//! \param p The particle being tracked
void score_collision_derivative(Particle* p);
void score_collision_derivative(Particle& p);
//! Adjust diff tally flux derivatives for a particle tracking event.
//
//! \param p The particle being tracked
//! \param distance The distance in [cm] traveled by the particle
void score_track_derivative(Particle* p, double distance);
void score_track_derivative(Particle& p, double distance);
} // namespace openmc

1
src/cell.cpp
View file

@ -767,6 +767,7 @@ DAGCell::DAGCell() : Cell{} {};
std::pair<double, int32_t>
DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
{
Expects(p);
// if we've changed direction or we're not on a surface,
// reset the history and update last direction
if (u != p->last_dir_ || on_surface == 0) {

12
src/event.cpp
View file

@ -62,7 +62,7 @@ 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);
initialize_history(simulation::particles[i], source_offset + i + 1);
dispatch_xs_event(i);
}
simulation::time_event_init.stop();
@ -76,17 +76,17 @@ void process_calculate_xs_events(SharedArray<EventQueueItem>& 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.
// 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];
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.
@ -136,7 +136,7 @@ void process_surface_crossing_events()
}
simulation::surface_crossing_queue.resize(0);
simulation::time_event_surface_crossing.stop();
}

150
src/geometry.cpp
View file

@ -35,23 +35,23 @@ std::vector<int64_t> overlap_check_count;
// Non-member functions
//==============================================================================
bool check_cell_overlap(Particle* p, bool error)
bool check_cell_overlap(Particle& p, bool error)
{
int n_coord = p->n_coord_;
int n_coord = p.n_coord_;
// Loop through each coordinate level
for (int j = 0; j < n_coord; j++) {
Universe& univ = *model::universes[p->coord_[j].universe];
Universe& univ = *model::universes[p.coord_[j].universe];
// Loop through each cell on this level
for (auto index_cell : univ.cells_) {
Cell& c = *model::cells[index_cell];
if (c.contains(p->coord_[j].r, p->coord_[j].u, p->surface_)) {
if (index_cell != p->coord_[j].cell) {
if (c.contains(p.coord_[j].r, p.coord_[j].u, p.surface_)) {
if (index_cell != p.coord_[j].cell) {
if (error) {
fatal_error(fmt::format(
"Overlapping cells detected: {}, {} on universe {}",
c.id_, model::cells[p->coord_[j].cell]->id_, univ.id_));
c.id_, model::cells[p.coord_[j].cell]->id_, univ.id_));
}
return true;
}
@ -67,7 +67,7 @@ bool check_cell_overlap(Particle* p, bool error)
//==============================================================================
bool
find_cell_inner(Particle* p, const NeighborList* neighbor_list)
find_cell_inner(Particle& p, const NeighborList* neighbor_list)
{
// Find which cell of this universe the particle is in. Use the neighbor list
// to shorten the search if one was provided.
@ -78,41 +78,41 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
i_cell = *it;
// Make sure the search cell is in the same universe.
int i_universe = p->coord_[p->n_coord_-1].universe;
int i_universe = p.coord_[p.n_coord_-1].universe;
if (model::cells[i_cell]->universe_ != i_universe) continue;
// Check if this cell contains the particle.
Position r {p->r_local()};
Direction u {p->u_local()};
auto surf = p->surface_;
Position r {p.r_local()};
Direction u {p.u_local()};
auto surf = p.surface_;
if (model::cells[i_cell]->contains(r, u, surf)) {
p->coord_[p->n_coord_-1].cell = i_cell;
p.coord_[p.n_coord_-1].cell = i_cell;
found = true;
break;
}
}
} else {
int i_universe = p->coord_[p->n_coord_-1].universe;
int i_universe = p.coord_[p.n_coord_-1].universe;
const auto& univ {*model::universes[i_universe]};
const auto& cells {
!univ.partitioner_
? model::universes[i_universe]->cells_
: univ.partitioner_->get_cells(p->r_local(), p->u_local())
: univ.partitioner_->get_cells(p.r_local(), p.u_local())
};
for (auto it = cells.cbegin(); it != cells.cend(); it++) {
i_cell = *it;
// Make sure the search cell is in the same universe.
int i_universe = p->coord_[p->n_coord_-1].universe;
int i_universe = p.coord_[p.n_coord_-1].universe;
if (model::cells[i_cell]->universe_ != i_universe) continue;
// Check if this cell contains the particle.
Position r {p->r_local()};
Direction u {p->u_local()};
auto surf = p->surface_;
Position r {p.r_local()};
Direction u {p.u_local()};
auto surf = p.surface_;
if (model::cells[i_cell]->contains(r, u, surf)) {
p->coord_[p->n_coord_-1].cell = i_cell;
p.coord_[p.n_coord_-1].cell = i_cell;
found = true;
break;
}
@ -120,7 +120,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
}
// Announce the cell that the particle is entering.
if (found && (settings::verbosity >= 10 || p->trace_)) {
if (found && (settings::verbosity >= 10 || p.trace_)) {
auto msg = fmt::format(" Entering cell {}", model::cells[i_cell]->id_);
write_message(msg, 1);
}
@ -134,35 +134,35 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
// Find the distribcell instance number.
int offset = 0;
if (c.distribcell_index_ >= 0) {
for (int i = 0; i < p->n_coord_; i++) {
const auto& c_i {*model::cells[p->coord_[i].cell]};
for (int i = 0; i < p.n_coord_; i++) {
const auto& c_i {*model::cells[p.coord_[i].cell]};
if (c_i.type_ == Fill::UNIVERSE) {
offset += c_i.offset_[c.distribcell_index_];
} else if (c_i.type_ == Fill::LATTICE) {
auto& lat {*model::lattices[p->coord_[i+1].lattice]};
int i_xyz[3] {p->coord_[i+1].lattice_x,
p->coord_[i+1].lattice_y,
p->coord_[i+1].lattice_z};
auto& lat {*model::lattices[p.coord_[i+1].lattice]};
int i_xyz[3] {p.coord_[i+1].lattice_x,
p.coord_[i+1].lattice_y,
p.coord_[i+1].lattice_z};
if (lat.are_valid_indices(i_xyz)) {
offset += lat.offset(c.distribcell_index_, i_xyz);
}
}
}
}
p->cell_instance_ = offset;
p.cell_instance_ = offset;
// Set the material and temperature.
p->material_last_ = p->material_;
p.material_last_ = p.material_;
if (c.material_.size() > 1) {
p->material_ = c.material_[p->cell_instance_];
p.material_ = c.material_[p.cell_instance_];
} else {
p->material_ = c.material_[0];
p.material_ = c.material_[0];
}
p->sqrtkT_last_ = p->sqrtkT_;
p.sqrtkT_last_ = p.sqrtkT_;
if (c.sqrtkT_.size() > 1) {
p->sqrtkT_ = c.sqrtkT_[p->cell_instance_];
p.sqrtkT_ = c.sqrtkT_[p.cell_instance_];
} else {
p->sqrtkT_ = c.sqrtkT_[0];
p.sqrtkT_ = c.sqrtkT_[0];
}
return true;
@ -172,12 +172,12 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
//! Found a lower universe, update this coord level then search the next.
// Set the lower coordinate level universe.
auto& coord {p->coord_[p->n_coord_]};
auto& coord {p.coord_[p.n_coord_]};
coord.universe = c.fill_;
// Set the position and direction.
coord.r = p->r_local();
coord.u = p->u_local();
coord.r = p.r_local();
coord.u = p.u_local();
// Apply translation.
coord.r -= c.translation_;
@ -188,7 +188,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
}
// Update the coordinate level and recurse.
++p->n_coord_;
++p.n_coord_;
return find_cell_inner(p, nullptr);
} else if (c.type_ == Fill::LATTICE) {
@ -198,9 +198,9 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
Lattice& lat {*model::lattices[c.fill_]};
// Set the position and direction.
auto& coord {p->coord_[p->n_coord_]};
coord.r = p->r_local();
coord.u = p->u_local();
auto& coord {p.coord_[p.n_coord_]};
coord.r = p.r_local();
coord.u = p.u_local();
// Apply translation.
coord.r -= c.translation_;
@ -230,13 +230,13 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
coord.universe = lat.outer_;
} else {
warning(fmt::format("Particle {} is outside lattice {} but the "
"lattice has no defined outer universe.", p->id_, lat.id_));
"lattice has no defined outer universe.", p.id_, lat.id_));
return false;
}
}
// Update the coordinate level and recurse.
++p->n_coord_;
++p.n_coord_;
return find_cell_inner(p, nullptr);
}
}
@ -247,25 +247,25 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
//==============================================================================
bool
find_cell(Particle* p, bool use_neighbor_lists)
find_cell(Particle& p, bool use_neighbor_lists)
{
// Determine universe (if not yet set, use root universe).
int i_universe = p->coord_[p->n_coord_-1].universe;
int i_universe = p.coord_[p.n_coord_-1].universe;
if (i_universe == C_NONE) {
p->coord_[0].universe = model::root_universe;
p->n_coord_ = 1;
p.coord_[0].universe = model::root_universe;
p.n_coord_ = 1;
i_universe = model::root_universe;
}
// Reset all the deeper coordinate levels.
for (int i = p->n_coord_; i < p->coord_.size(); i++) {
p->coord_[i].reset();
for (int i = p.n_coord_; i < p.coord_.size(); i++) {
p.coord_[i].reset();
}
if (use_neighbor_lists) {
// Get the cell this particle was in previously.
auto coord_lvl = p->n_coord_ - 1;
auto i_cell = p->coord_[coord_lvl].cell;
auto coord_lvl = p.n_coord_ - 1;
auto i_cell = p.coord_[coord_lvl].cell;
Cell& c {*model::cells[i_cell]};
// Search for the particle in that cell's neighbor list. Return if we
@ -277,7 +277,7 @@ find_cell(Particle* p, bool use_neighbor_lists)
// cells in this universe, and update the neighbor list if we find a new
// neighboring cell.
found = find_cell_inner(p, nullptr);
if (found) c.neighbors_.push_back(p->coord_[coord_lvl].cell);
if (found) c.neighbors_.push_back(p.coord_[coord_lvl].cell);
return found;
} else {
@ -289,15 +289,15 @@ find_cell(Particle* p, bool use_neighbor_lists)
//==============================================================================
void
cross_lattice(Particle* p, const BoundaryInfo& boundary)
cross_lattice(Particle& p, const BoundaryInfo& boundary)
{
auto& coord {p->coord_[p->n_coord_ - 1]};
auto& coord {p.coord_[p.n_coord_ - 1]};
auto& lat {*model::lattices[coord.lattice]};
if (settings::verbosity >= 10 || p->trace_) {
if (settings::verbosity >= 10 || p.trace_) {
write_message(fmt::format(
" Crossing lattice {}. Current position ({},{},{}). r={}",
lat.id_, coord.lattice_x, coord.lattice_y, coord.lattice_z, p->r()), 1);
lat.id_, coord.lattice_x, coord.lattice_y, coord.lattice_z, p.r()), 1);
}
// Set the lattice indices.
@ -307,37 +307,37 @@ cross_lattice(Particle* p, const BoundaryInfo& boundary)
std::array<int, 3> i_xyz {coord.lattice_x, coord.lattice_y, coord.lattice_z};
// Set the new coordinate position.
const auto& upper_coord {p->coord_[p->n_coord_ - 2]};
const auto& upper_coord {p.coord_[p.n_coord_ - 2]};
const auto& cell {model::cells[upper_coord.cell]};
Position r = upper_coord.r;
r -= cell->translation_;
if (!cell->rotation_.empty()) {
r = r.rotate(cell->rotation_);
}
p->r_local() = lat.get_local_position(r, i_xyz);
p.r_local() = lat.get_local_position(r, i_xyz);
if (!lat.are_valid_indices(i_xyz)) {
// The particle is outside the lattice. Search for it from the base coords.
p->n_coord_ = 1;
p.n_coord_ = 1;
bool found = find_cell(p, 0);
if (!found && p->alive_) {
p->mark_as_lost(fmt::format("Could not locate particle {} after "
"crossing a lattice boundary", p->id_));
if (!found && p.alive_) {
p.mark_as_lost(fmt::format("Could not locate particle {} after "
"crossing a lattice boundary", p.id_));
}
} else {
// Find cell in next lattice element.
p->coord_[p->n_coord_-1].universe = lat[i_xyz];
p.coord_[p.n_coord_-1].universe = lat[i_xyz];
bool found = find_cell(p, 0);
if (!found) {
// A particle crossing the corner of a lattice tile may not be found. In
// this case, search for it from the base coords.
p->n_coord_ = 1;
p.n_coord_ = 1;
bool found = find_cell(p, 0);
if (!found && p->alive_) {
p->mark_as_lost(fmt::format("Could not locate particle {} after "
"crossing a lattice boundary", p->id_));
if (!found && p.alive_) {
p.mark_as_lost(fmt::format("Could not locate particle {} after "
"crossing a lattice boundary", p.id_));
}
}
}
@ -345,7 +345,7 @@ cross_lattice(Particle* p, const BoundaryInfo& boundary)
//==============================================================================
BoundaryInfo distance_to_boundary(Particle* p)
BoundaryInfo distance_to_boundary(Particle& p)
{
BoundaryInfo info;
double d_lat = INFINITY;
@ -354,14 +354,14 @@ BoundaryInfo distance_to_boundary(Particle* p)
std::array<int, 3> level_lat_trans {};
// Loop over each coordinate level.
for (int i = 0; i < p->n_coord_; i++) {
const auto& coord {p->coord_[i]};
for (int i = 0; i < p.n_coord_; i++) {
const auto& coord {p.coord_[i]};
Position r {coord.r};
Direction u {coord.u};
Cell& c {*model::cells[coord.cell]};
// Find the oncoming surface in this cell and the distance to it.
auto surface_distance = c.distance(r, u, p->surface_, p);
auto surface_distance = c.distance(r, u, p.surface_, &p);
d_surf = surface_distance.first;
level_surf_cross = surface_distance.second;
@ -377,8 +377,8 @@ BoundaryInfo distance_to_boundary(Particle* p)
lattice_distance = lat.distance(r, u, i_xyz);
break;
case LatticeType::hex:
auto& cell_above {model::cells[p->coord_[i-1].cell]};
Position r_hex {p->coord_[i-1].r};
auto& cell_above {model::cells[p.coord_[i-1].cell]};
Position r_hex {p.coord_[i-1].r};
r_hex -= cell_above->translation_;
if (coord.rotated) {
r_hex = r_hex.rotate(cell_above->rotation_);
@ -391,8 +391,8 @@ BoundaryInfo distance_to_boundary(Particle* p)
level_lat_trans = lattice_distance.second;
if (d_lat < 0) {
p->mark_as_lost(fmt::format(
"Particle {} had a negative distance to a lattice boundary", p->id_));
p.mark_as_lost(fmt::format(
"Particle {} had a negative distance to a lattice boundary", p.id_));
}
}
@ -450,7 +450,7 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance)
p.r() = Position{xyz};
p.u() = {0.0, 0.0, 1.0};
if (!find_cell(&p, false)) {
if (!find_cell(p, false)) {
set_errmsg(fmt::format("Could not find cell at position {}.", p.r()));
return OPENMC_E_GEOMETRY;
}

42
src/output.cpp
View file

@ -143,10 +143,10 @@ std::string time_stamp()
//==============================================================================
extern "C" void print_particle(Particle* p)
void print_particle(Particle& p)
{
// Display particle type and ID.
switch (p->type_) {
switch (p.type_) {
case Particle::Type::neutron:
fmt::print("Neutron ");
break;
@ -162,46 +162,46 @@ extern "C" void print_particle(Particle* p)
default:
fmt::print("Unknown Particle ");
}
fmt::print("{}\n", p->id_);
fmt::print("{}\n", p.id_);
// Display particle geometry hierarchy.
for (auto i = 0; i < p->n_coord_; i++) {
for (auto i = 0; i < p.n_coord_; i++) {
fmt::print(" Level {}\n", i);
if (p->coord_[i].cell != C_NONE) {
const Cell& c {*model::cells[p->coord_[i].cell]};
if (p.coord_[i].cell != C_NONE) {
const Cell& c {*model::cells[p.coord_[i].cell]};
fmt::print(" Cell = {}\n", c.id_);
}
if (p->coord_[i].universe != C_NONE) {
const Universe& u {*model::universes[p->coord_[i].universe]};
if (p.coord_[i].universe != C_NONE) {
const Universe& u {*model::universes[p.coord_[i].universe]};
fmt::print(" Universe = {}\n", u.id_);
}
if (p->coord_[i].lattice != C_NONE) {
const Lattice& lat {*model::lattices[p->coord_[i].lattice]};
if (p.coord_[i].lattice != C_NONE) {
const Lattice& lat {*model::lattices[p.coord_[i].lattice]};
fmt::print(" Lattice = {}\n", lat.id_);
fmt::print(" Lattice position = ({},{},{})\n", p->coord_[i].lattice_x,
p->coord_[i].lattice_y, p->coord_[i].lattice_z);
fmt::print(" Lattice position = ({},{},{})\n", p.coord_[i].lattice_x,
p.coord_[i].lattice_y, p.coord_[i].lattice_z);
}
fmt::print(" r = {}\n", p->coord_[i].r);
fmt::print(" u = {}\n", p->coord_[i].u);
fmt::print(" r = {}\n", p.coord_[i].r);
fmt::print(" u = {}\n", p.coord_[i].u);
}
// Display miscellaneous info.
if (p->surface_ != 0) {
if (p.surface_ != 0) {
// Surfaces identifiers are >= 1, but indices are >= 0 so we need -1
const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]};
fmt::print(" Surface = {}\n", (p->surface_ > 0) ? surf.id_ : -surf.id_);
const Surface& surf {*model::surfaces[std::abs(p.surface_)-1]};
fmt::print(" Surface = {}\n", (p.surface_ > 0) ? surf.id_ : -surf.id_);
}
fmt::print(" Weight = {}\n", p->wgt_);
fmt::print(" Weight = {}\n", p.wgt_);
if (settings::run_CE) {
fmt::print(" Energy = {}\n", p->E_);
fmt::print(" Energy = {}\n", p.E_);
} else {
fmt::print(" Energy Group = {}\n", p->g_);
fmt::print(" Energy Group = {}\n", p.g_);
}
fmt::print(" Delayed Group = {}\n\n", p->delayed_group_);
fmt::print(" Delayed Group = {}\n\n", p.delayed_group_);
}
//==============================================================================

22
src/particle.cpp
View file

@ -156,7 +156,7 @@ Particle::event_calculate_xs()
// initiate a search for the current cell. This generally happens at the
// beginning of the history and again for any secondary particles
if (coord_[n_coord_ - 1].cell == C_NONE) {
if (!find_cell(this, false)) {
if (!find_cell(*this, false)) {
this->mark_as_lost("Could not find the cell containing particle "
+ std::to_string(id_));
return;
@ -169,7 +169,7 @@ Particle::event_calculate_xs()
// Write particle track.
if (write_track_) write_particle_track(*this);
if (settings::check_overlaps) check_cell_overlap(this);
if (settings::check_overlaps) check_cell_overlap(*this);
// Calculate microscopic and macroscopic cross sections
if (material_ != MATERIAL_VOID) {
@ -201,7 +201,7 @@ 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 ||
@ -234,7 +234,7 @@ Particle::event_advance()
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) {
score_track_derivative(this, distance);
score_track_derivative(*this, distance);
}
}
@ -255,7 +255,7 @@ Particle::event_cross_surface()
boundary_.lattice_translation[1] != 0 ||
boundary_.lattice_translation[2] != 0) {
// Particle crosses lattice boundary
cross_lattice(this, boundary_);
cross_lattice(*this, boundary_);
event_ = TallyEvent::LATTICE;
} else {
// Particle crosses surface
@ -337,7 +337,7 @@ Particle::event_collide()
}
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) score_collision_derivative(this);
if (!model::active_tallies.empty()) score_collision_derivative(*this);
}
void
@ -490,7 +490,7 @@ Particle::cross_surface()
// (unless we're using a dagmc model, which has exactly one universe)
if (!settings::dagmc) {
n_coord_ = 1;
if (!find_cell(this, true)) {
if (!find_cell(*this, true)) {
this->mark_as_lost("Couldn't find particle after reflecting from surface "
+ std::to_string(surf->id_) + ".");
return;
@ -545,7 +545,7 @@ Particle::cross_surface()
// Figure out what cell particle is in now
n_coord_ = 1;
if (!find_cell(this, true)) {
if (!find_cell(*this, true)) {
this->mark_as_lost("Couldn't find particle after hitting periodic "
"boundary on surface " + std::to_string(surf->id_) + ".");
return;
@ -583,7 +583,7 @@ Particle::cross_surface()
}
#endif
if (find_cell(this, true)) return;
if (find_cell(*this, true)) return;
// ==========================================================================
// COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
@ -591,7 +591,7 @@ Particle::cross_surface()
// Remove lower coordinate levels and assignment of surface
surface_ = 0;
n_coord_ = 1;
bool found = find_cell(this, false);
bool found = find_cell(*this, false);
if (settings::run_mode != RunMode::PLOTTING && (!found)) {
// If a cell is still not found, there are two possible causes: 1) there is
@ -605,7 +605,7 @@ Particle::cross_surface()
// Couldn't find next cell anywhere! This probably means there is an actual
// undefined region in the geometry.
if (!find_cell(this, false)) {
if (!find_cell(*this, false)) {
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.");

4
src/particle_restart.cpp
View file

@ -115,7 +115,7 @@ void run_particle_restart()
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());
@ -123,7 +123,7 @@ void run_particle_restart()
transport_history_based_single_particle(p);
// Write output if particle made it
print_particle(&p);
print_particle(p);
}
} // namespace openmc

48
src/simulation.cpp
View file

@ -443,7 +443,7 @@ void finalize_generation()
}
}
void initialize_history(Particle* p, int64_t index_source)
void initialize_history(Particle& p, int64_t index_source)
{
// set defaults
if (settings::run_mode == RunMode::FIXED_SOURCE) {
@ -453,76 +453,76 @@ void initialize_history(Particle* p, int64_t index_source)
uint64_t seed = init_seed(id, STREAM_SOURCE);
// sample from external source distribution or custom library then set
auto site = sample_external_source(&seed);
p->from_source(&site);
p.from_source(&site);
} else if (settings::run_mode == RunMode::EIGENVALUE) {
// set defaults for eigenvalue simulations from primary bank
p->from_source(&simulation::source_bank[index_source - 1]);
p.from_source(&simulation::source_bank[index_source - 1]);
}
p->current_work_ = index_source;
p.current_work_ = index_source;
// set identifier for particle
p->id_ = simulation::work_index[mpi::rank] + index_source;
p.id_ = simulation::work_index[mpi::rank] + index_source;
// set progeny count to zero
p->n_progeny_ = 0;
p.n_progeny_ = 0;
// set random number seed
int64_t particle_seed = (simulation::total_gen + overall_generation() - 1)
* settings::n_particles + p->id_;
init_particle_seeds(particle_seed, p->seeds_);
* settings::n_particles + p.id_;
init_particle_seeds(particle_seed, p.seeds_);
// set particle trace
p->trace_ = false;
p.trace_ = false;
if (simulation::current_batch == settings::trace_batch &&
simulation::current_gen == settings::trace_gen &&
p->id_ == settings::trace_particle) p->trace_ = true;
p.id_ == settings::trace_particle) p.trace_ = true;
// Set particle track.
p->write_track_ = false;
p.write_track_ = false;
if (settings::write_all_tracks) {
p->write_track_ = true;
p.write_track_ = true;
} else if (settings::track_identifiers.size() > 0) {
for (const auto& t : settings::track_identifiers) {
if (simulation::current_batch == t[0] &&
simulation::current_gen == t[1] &&
p->id_ == t[2]) {
p->write_track_ = true;
p.id_ == t[2]) {
p.write_track_ = true;
break;
}
}
}
// Display message if high verbosity or trace is on
if (settings::verbosity >= 9 || p->trace_) {
write_message("Simulating Particle " + std::to_string(p->id_));
if (settings::verbosity >= 9 || p.trace_) {
write_message("Simulating Particle " + std::to_string(p.id_));
}
// Add paricle's starting weight to count for normalizing tallies later
#pragma omp atomic
simulation::total_weight += p->wgt_;
simulation::total_weight += p.wgt_;
initialize_history_partial(p);
}
void initialize_history_partial(Particle* p)
void initialize_history_partial(Particle& p)
{
// Force calculation of cross-sections by setting last energy to zero
if (settings::run_CE) {
for (auto& micro : p->neutron_xs_) micro.last_E = 0.0;
for (auto& micro : p.neutron_xs_) micro.last_E = 0.0;
}
// Prepare to write out particle track.
if (p->write_track_) add_particle_track(*p);
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);
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());
p.filter_matches_.resize(model::tally_filters.size());
}
int overall_generation()
@ -617,7 +617,7 @@ void transport_history_based()
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
Particle p;
initialize_history(&p, i_work);
initialize_history(p, i_work);
transport_history_based_single_particle(p);
}
}

2
src/surface.cpp
View file

@ -5,6 +5,7 @@
#include <utility>
#include <fmt/core.h>
#include <gsl/gsl>
#include "openmc/error.h"
#include "openmc/dagmc.h"
@ -276,6 +277,7 @@ Direction DAGSurface::normal(Position r) const
Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const
{
Expects(p);
p->history_.reset_to_last_intersection();
p->last_dir_ = Surface::reflect(r, u, p);
return p->last_dir_;

32
src/tallies/derivative.cpp
View file

@ -558,15 +558,15 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide,
}
void
score_track_derivative(Particle* p, double distance)
score_track_derivative(Particle& p, double distance)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material_ == MATERIAL_VOID) return;
const Material& material {*model::materials[p->material_]};
if (p.material_ == MATERIAL_VOID) return;
const Material& material {*model::materials[p.material_]};
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) {
@ -575,7 +575,7 @@ score_track_derivative(Particle* p, double distance)
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
// (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
flux_deriv -= distance * p->macro_xs_.total
flux_deriv -= distance * p.macro_xs_.total
/ material.density_gpcc_;
break;
@ -584,19 +584,19 @@ score_track_derivative(Particle* p, double distance)
// (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
// (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
flux_deriv -= distance
* p->neutron_xs_[deriv.diff_nuclide].total;
* p.neutron_xs_[deriv.diff_nuclide].total;
break;
case DerivativeVariable::TEMPERATURE:
for (auto i = 0; i < material.nuclide_.size(); ++i) {
const auto& nuc {*data::nuclides[material.nuclide_[i]]};
if (multipole_in_range(&nuc, p->E_last_)) {
if (multipole_in_range(&nuc, p.E_last_)) {
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
// (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_, p->sqrtkT_);
= nuc.multipole_->evaluate_deriv(p.E_, p.sqrtkT_);
flux_deriv -= distance * (dsig_s + dsig_a)
* material.atom_density_(i);
}
@ -606,16 +606,16 @@ score_track_derivative(Particle* p, double distance)
}
}
void score_collision_derivative(Particle* p)
void score_collision_derivative(Particle& p)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material_ == MATERIAL_VOID) return;
if (p.material_ == MATERIAL_VOID) return;
const Material& material {*model::materials[p->material_]};
const Material& material {*model::materials[p.material_]};
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;
@ -629,7 +629,7 @@ void score_collision_derivative(Particle* p)
break;
case DerivativeVariable::NUCLIDE_DENSITY:
if (p->event_nuclide_ != deriv.diff_nuclide) continue;
if (p.event_nuclide_ != deriv.diff_nuclide) continue;
// Find the index in this material for the diff_nuclide.
int i;
for (i = 0; i < material.nuclide_.size(); ++i)
@ -651,14 +651,14 @@ void score_collision_derivative(Particle* p)
// Loop over the material's nuclides until we find the event nuclide.
for (auto i_nuc : material.nuclide_) {
const auto& nuc {*data::nuclides[i_nuc]};
if (i_nuc == p->event_nuclide_ && multipole_in_range(&nuc, p->E_last_)) {
if (i_nuc == p.event_nuclide_ && multipole_in_range(&nuc, p.E_last_)) {
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
const auto& micro_xs {p->neutron_xs_[i_nuc]};
const auto& micro_xs {p.neutron_xs_[i_nuc]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
= nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_);
flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
// Note that this is an approximation! The real scattering cross
// section is

2
src/volume_calc.cpp
View file

@ -135,7 +135,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
p.u() = {0.5, 0.5, 0.5};
// If this location is not in the geometry at all, move on to next block
if (!find_cell(&p, false)) continue;
if (!find_cell(p, false)) continue;
if (domain_type_ == TallyDomain::MATERIAL) {
if (p.material_ != MATERIAL_VOID) {