From 9dab8bf507731da358edb2111e298b36f2900bf5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Apr 2020 12:49:48 -0500 Subject: [PATCH] Use references consistently in geometry functions --- include/openmc/geometry.h | 8 +- include/openmc/output.h | 2 +- include/openmc/plot.h | 4 +- include/openmc/simulation.h | 4 +- include/openmc/tallies/derivative.h | 4 +- src/cell.cpp | 1 + src/event.cpp | 12 +-- src/geometry.cpp | 150 ++++++++++++++-------------- src/output.cpp | 42 ++++---- src/particle.cpp | 22 ++-- src/particle_restart.cpp | 4 +- src/simulation.cpp | 48 ++++----- src/surface.cpp | 2 + src/tallies/derivative.cpp | 32 +++--- src/volume_calc.cpp | 2 +- 15 files changed, 170 insertions(+), 167 deletions(-) diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index ef226db66..6fd0d1ad9 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -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 diff --git a/include/openmc/output.h b/include/openmc/output.h index 439820bd3..2ebe2711f 100644 --- a/include/openmc/output.h +++ b/include/openmc/output.h @@ -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(); diff --git a/include/openmc/plot.h b/include/openmc/plot.h index e05d6f2e2..2468af1bf 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -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 diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index d02308fe1..a3032ee01 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -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 //! diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index 07ce03c94..421094c99 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -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 diff --git a/src/cell.cpp b/src/cell.cpp index 1d175f05d..324257022 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -767,6 +767,7 @@ DAGCell::DAGCell() : Cell{} {}; std::pair 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) { diff --git a/src/event.cpp b/src/event.cpp index 2f22d0954..d1b4a36a7 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -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& 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(); } diff --git a/src/geometry.cpp b/src/geometry.cpp index 7af1da95c..4b5c9b47c 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -35,23 +35,23 @@ std::vector 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 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 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; } diff --git a/src/output.cpp b/src/output.cpp index c429b2f59..338dda634 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -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_); } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index 8008310fd..481d7755a 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -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."); diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 97b82c5c6..3347030ba 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -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 diff --git a/src/simulation.cpp b/src/simulation.cpp index 4da79ce53..d3e0f48c7 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -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); } } diff --git a/src/surface.cpp b/src/surface.cpp index a2c54260a..8133a6a7e 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -5,6 +5,7 @@ #include #include +#include #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_; diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index d8d9d8820..1e3ca4884 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -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 diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 97594c3b5..d1b14e329 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -135,7 +135,7 @@ std::vector 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) {