From 1fa2c0e8893f993f44f78a3e0abdc5e8a54f8a3d Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 4 Mar 2021 17:56:34 -0500 Subject: [PATCH] separate find_cell into two functions --- include/openmc/geometry.h | 7 ++--- include/openmc/plot.h | 2 +- src/geometry.cpp | 63 +++++++++++++++++---------------------- src/particle.cpp | 13 ++++---- src/volume_calc.cpp | 3 +- 5 files changed, 40 insertions(+), 48 deletions(-) diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 0ad891269..258153502 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -43,14 +43,11 @@ bool check_cell_overlap(Particle& p, bool error=true); //! //! \param p A particle to be located. This function will populate the //! geometry-dependent data fields of the particle. -//! \param use_neighbor_lists If true, neighbor lists will be used to accelerate -//! the geometry search, but this only works if the cell attribute of the -//! particle's lowest coordinate level is valid and meaningful. //! \return True if the particle's location could be found and ascribed to a //! valid geometry coordinate stack. //============================================================================== - -bool find_cell(Particle& p, bool use_neighbor_lists); +bool brute_force_find_cell(Particle& p); +bool neighbor_list_find_cell(Particle& p); // Only usable on surface crossings //============================================================================== //! Move a particle into a new lattice tile. diff --git a/include/openmc/plot.h b/include/openmc/plot.h index a15bc7e38..295a7b874 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -182,7 +182,7 @@ 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 = brute_force_find_cell(p); j = p.n_coord_ - 1; if (level >= 0) { j = level; } if (found_cell) { diff --git a/src/geometry.cpp b/src/geometry.cpp index 4b5c9b47c..819e069af 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -246,44 +246,37 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) //============================================================================== -bool -find_cell(Particle& p, bool use_neighbor_lists) +bool neighbor_list_find_cell(Particle& p) +{ + // Get the cell this particle was in previously. + 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 + // found the particle. + bool found = find_cell_inner(p, &c.neighbors_); + if (found) + return found; + + // The particle could not be found in the neighbor list. Try searching all + // 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); + return found; +} + +bool brute_force_find_cell(Particle& p) { - // Determine universe (if not yet set, use root 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; 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(); - } - - 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; - Cell& c {*model::cells[i_cell]}; - - // Search for the particle in that cell's neighbor list. Return if we - // found the particle. - bool found = find_cell_inner(p, &c.neighbors_); - if (found) return found; - - // The particle could not be found in the neighbor list. Try searching all - // 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); - return found; - - } else { - // Search all cells in this universe for the particle. - return find_cell_inner(p, nullptr); - } + return find_cell_inner(p, nullptr); } //============================================================================== @@ -319,7 +312,7 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary) if (!lat.are_valid_indices(i_xyz)) { // The particle is outside the lattice. Search for it from the base coords. p.n_coord_ = 1; - bool found = find_cell(p, 0); + bool found = brute_force_find_cell(p); if (!found && p.alive_) { p.mark_as_lost(fmt::format("Could not locate particle {} after " "crossing a lattice boundary", p.id_)); @@ -328,13 +321,13 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary) } else { // Find cell in next lattice element. p.coord_[p.n_coord_-1].universe = lat[i_xyz]; - bool found = find_cell(p, 0); + bool found = brute_force_find_cell(p); 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; - bool found = find_cell(p, 0); + bool found = brute_force_find_cell(p); if (!found && p.alive_) { p.mark_as_lost(fmt::format("Could not locate particle {} after " "crossing a lattice boundary", p.id_)); @@ -450,7 +443,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 (!brute_force_find_cell(p)) { set_errmsg(fmt::format("Could not find cell at position {}.", p.r())); return OPENMC_E_GEOMETRY; } diff --git a/src/particle.cpp b/src/particle.cpp index d03511e67..3273ef4cc 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 (!brute_force_find_cell(*this)) { this->mark_as_lost("Could not find the cell containing particle " + std::to_string(id_)); return; @@ -454,7 +454,8 @@ Particle::cross_surface() } #endif - if (find_cell(*this, true)) return; + if (neighbor_list_find_cell(*this)) + return; // ========================================================================== // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS @@ -462,7 +463,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 = brute_force_find_cell(*this); if (settings::run_mode != RunMode::PLOTTING && (!found)) { // If a cell is still not found, there are two possible causes: 1) there is @@ -476,7 +477,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 (!brute_force_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."); @@ -553,7 +554,7 @@ Particle::cross_reflective_bc(const Surface& surf, Direction new_u) // (unless we're using a dagmc model, which has exactly one universe) if (!settings::dagmc) { n_coord_ = 1; - if (!find_cell(*this, true)) { + if (!neighbor_list_find_cell(*this)) { this->mark_as_lost("Couldn't find particle after reflecting from surface " + std::to_string(surf.id_) + "."); return; @@ -601,7 +602,7 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r, // Figure out what cell particle is in now n_coord_ = 1; - if (!find_cell(*this, true)) { + 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."); diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 96ef55234..46fc8352a 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -135,7 +135,8 @@ 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 (!brute_force_find_cell(p)) + continue; if (domain_type_ == TallyDomain::MATERIAL) { if (p.material_ != MATERIAL_VOID) {