30828111152

Committed 03 Aug 2026 03:34PM UTC coverage: 81.425% (-0.004%) from 81.429%

Build # 30828111152

Build Type

Pull #3933

github

Committed by

web-flow

Commit Message

Merge 53dbb570d into 5982acdf8

Pull Request Pull Request #3933: Reevaluate cell in a close collision after a surface crossing

Coverage Stats

18533 of 26846 branches covered (69.03%)

Branch coverage included in aggregate %.

15 of 18 new or added lines in 2 files covered. (83.33%)

60260 of 69922 relevant lines covered (86.18%)

50032313.18 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

82.54

/src/geometry.cpp

#include "openmc/geometry.h"

#include <fmt/core.h>
#include <fmt/ostream.h>

#include "openmc/array.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/lattice.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/surface.h"

namespace openmc {

//==============================================================================
// Global variables
//==============================================================================

namespace model {

int root_universe {-1};
int n_coord_levels;

vector<int64_t> overlap_check_count;

vector<OverlapKey> overlap_keys;
std::unordered_map<OverlapKey, int, OverlapKeyHash> overlap_key_index;

} // namespace model

//==============================================================================
// Non-member functions
//==============================================================================

int check_cell_overlap(GeometryState& p, bool error)
{
  int n_coord = p.n_coord();

  // If no overlap found, return a nonphysical index
  int overlap_index = -1;

  // Loop through each coordinate level
  for (int j = 0; j < n_coord; j++) {
    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())) {
#pragma omp atomic
        ++model::overlap_check_count[index_cell];
        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_));
          }

          // With no fatal error (plotter is calling), now adds overlaps and
          // ensures order does not matter when making overlap key
          int cell_a = model::cells[index_cell]->id_;
          int cell_b = model::cells[p.coord(j).cell()]->id_;
          int a = std::min(cell_a, cell_b);
          int b = std::max(cell_a, cell_b);
          OverlapKey key {univ.id_, a, b};
#pragma omp critical(overlap_key_update)
          {
            auto it = model::overlap_key_index.find(key);
            if (it != model::overlap_key_index.end()) {
              overlap_index = it->second; // already exists, reuse index
            } else {
              int idx = int(model::overlap_keys.size());
              model::overlap_keys.push_back(key);
              model::overlap_key_index[key] = idx;
              overlap_index = idx;
            }
          }
          break;
        }
      }
    }
  }
  return overlap_index;
}

//==============================================================================

int cell_instance_at_level(const GeometryState& p, int level)
{
  // throw error if the requested level is too deep for the geometry
  if (level > model::n_coord_levels) {
    fatal_error(fmt::format("Cell instance at level {} requested, but only {} "
                            "levels exist in the geometry.",
      level, p.n_coord()));
  }

  // determine the cell instance
  Cell& c {*model::cells[p.coord(level).cell()]};

  // quick exit if this cell doesn't have distribcell instances
  if (c.distribcell_index_ == C_NONE)
    return C_NONE;

  // compute the cell's instance
  int instance = 0;
  for (int i = 0; i < level; i++) {
    const auto& c_i {*model::cells[p.coord(i).cell()]};
    if (c_i.type_ == Fill::UNIVERSE) {
      instance += c_i.offset_[c.distribcell_index_];
    } else if (c_i.type_ == Fill::LATTICE) {
      instance += c_i.offset_[c.distribcell_index_];
      auto& lat {*model::lattices[p.coord(i + 1).lattice()]};
      const auto& i_xyz {p.coord(i + 1).lattice_index()};
      if (lat.are_valid_indices(i_xyz)) {
        instance += lat.offset(c.distribcell_index_, i_xyz);
      }
    }
  }
  return instance;
}

//==============================================================================

bool find_cell_inner(
  GeometryState& p, const NeighborList* neighbor_list, bool verbose)
{
  // Find which cell of this universe the particle is in.  Use the neighbor list
  // to shorten the search if one was provided.
  bool found = false;
  int32_t i_cell = C_NONE;
  if (neighbor_list) {
    for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) {
      i_cell = *it;

      // Make sure the search cell is in the same universe.
      int i_universe = p.lowest_coord().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();
      if (model::cells[i_cell]->contains(r, u, surf)) {
        p.lowest_coord().cell() = i_cell;
        found = true;
        break;
      }
    }

    // If we're attempting a neighbor list search and fail, we
    // now know we should return false. This will trigger an
    // exhaustive search from neighbor_list_find_cell and make
    // the result from that be appended to the neighbor list.
    if (!found) {
      return found;
    }
  }

  // Check successively lower coordinate levels until finding material fill
  for (;; ++p.n_coord()) {
    // If we did not attempt to use neighbor lists, i_cell is still C_NONE.  In
    // that case, we should now do an exhaustive search to find the right value
    // of i_cell.
    //
    // Alternatively, neighbor list searches could have succeeded, but we found
    // that the fill of the neighbor cell was another universe. As such, in the
    // code below this conditional, we set i_cell back to C_NONE to indicate
    // that.
    if (i_cell == C_NONE) {
      int i_universe = p.lowest_coord().universe();
      const auto& univ {model::universes[i_universe]};
      found = univ->find_cell(p);
    }

    if (!found) {
      return found;
    }
    i_cell = p.lowest_coord().cell();

    // Announce the cell that the particle is entering.
    if (found && verbose) {
      auto msg = fmt::format("    Entering cell {}", model::cells[i_cell]->id_);
      write_message(msg, 1);
    }

    Cell& c {*model::cells[i_cell]};
    if (c.type_ == Fill::MATERIAL) {
      // Found a material cell which means this is the lowest coord level.

      p.cell_instance() = 0;
      // Find the distribcell instance number.
      if (c.distribcell_index_ >= 0) {
        p.cell_instance() = cell_instance_at_level(p, p.n_coord() - 1);
      }

      // Set the material, temperature and density multiplier.
      p.material_last() = p.material();
      p.material() = c.material(p.cell_instance());
      p.sqrtkT_last() = p.sqrtkT();
      p.sqrtkT() = c.sqrtkT(p.cell_instance());
      p.density_mult_last() = p.density_mult();
      p.density_mult() = c.density_mult(p.cell_instance());

      return true;

    } else if (c.type_ == Fill::UNIVERSE) {
      //========================================================================
      //! 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())};
      coord.universe() = c.fill_;

      // Set the position and direction.
      coord.r() = p.r_local();
      coord.u() = p.u_local();

      // Apply translation.
      coord.r() -= c.translation_;

      // Apply rotation.
      if (!c.rotation_.empty()) {
        coord.rotate(c.rotation_);
      }

    } else if (c.type_ == Fill::LATTICE) {
      //========================================================================
      //! Found a lower lattice, update this coord level then search the next.

      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();

      // Apply translation.
      coord.r() -= c.translation_;

      // Apply rotation.
      if (!c.rotation_.empty()) {
        coord.rotate(c.rotation_);
      }

      // Determine lattice indices.
      auto& i_xyz {coord.lattice_index()};
      lat.get_indices(coord.r(), coord.u(), i_xyz);

      // Get local position in appropriate lattice cell
      coord.r() = lat.get_local_position(coord.r(), i_xyz);

      // Set lattice indices.
      coord.lattice() = c.fill_;

      // Set the lower coordinate level universe.
      if (lat.are_valid_indices(i_xyz)) {
        coord.universe() = lat[i_xyz];
      } else {
        if (lat.outer_ != NO_OUTER_UNIVERSE) {
          coord.universe() = lat.outer_;
        } else {
          p.mark_as_lost(fmt::format(
            "Particle {} left lattice {}, but it has no outer definition.",
            p.id(), lat.id_));
        }
      }
    }
    i_cell = C_NONE; // trip non-neighbor cell search at next iteration
    found = false;
  }

  return found;
}

//==============================================================================

bool neighbor_list_find_cell(GeometryState& p, bool verbose)
{

  // Reset all the deeper coordinate levels.
  for (int i = p.n_coord(); i < model::n_coord_levels; i++) {
    p.coord(i).reset();
  }

  // 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_, verbose);
  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, verbose);
  if (found)
    c.neighbors_.push_back(p.coord(coord_lvl).cell());
  return found;
}

void reconcile_cell_after_collision(GeometryState& p)
{
  // Find the first coordinate level whose current cell is inconsistent with
  // the particle's post-collision direction.
  int invalid_level = C_NONE;
  for (int level = 0; level < p.n_coord(); ++level) {
    const auto& coord {p.coord(level)};
    if (coord.cell() == C_NONE || !model::cells[coord.cell()]->contains(
                                    coord.r(), coord.u(), SURFACE_NONE)) {
      invalid_level = level;
      break;
    }
  }

  if (invalid_level == C_NONE)
    return;

  // Search from the first inconsistent level so that parent cells and
  // transformed lower universes are both handled correctly.
  if (p.coord(invalid_level).cell() != C_NONE) {
    p.n_coord() = invalid_level + 1;
    if (neighbor_list_find_cell(p))
      return;
  }

  // The current cell may not have a complete neighbor list yet. Fall back to
  // the normal exhaustive search used after a failed surface crossing.
  p.n_coord() = 1;
  if (!exhaustive_find_cell(p)) {
    p.mark_as_lost("Could not find particle after a collision near a surface.");
  }
}

bool exhaustive_find_cell(GeometryState& p, bool verbose)
{
  int i_universe = p.lowest_coord().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 < model::n_coord_levels; i++) {
    p.coord(i).reset();
  }
  return find_cell_inner(p, nullptr, verbose);
}

//==============================================================================

void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)
{
  auto& coord {p.lowest_coord()};
  auto& lat {*model::lattices[coord.lattice()]};

  if (verbose) {
    write_message(
      fmt::format("    Crossing lattice {}. Current position ({},{},{}). r={}",
        lat.id_, coord.lattice_index()[0], coord.lattice_index()[1],
        coord.lattice_index()[2], p.r()),
      1);
  }

  // Set the lattice indices.
  coord.lattice_index()[0] += boundary.lattice_translation()[0];
  coord.lattice_index()[1] += boundary.lattice_translation()[1];
  coord.lattice_index()[2] += boundary.lattice_translation()[2];

  // Set the new coordinate position.
  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, coord.lattice_index());

  if (!lat.are_valid_indices(coord.lattice_index())) {
    // The particle is outside the lattice.  Search for it from the base coords.
    p.n_coord() = 1;
    bool found = exhaustive_find_cell(p);

    if (!found) {
      p.mark_as_lost(fmt::format("Particle {} could not be located after "
                                 "crossing a boundary of lattice {}",
        p.id(), lat.id_));
    }

  } else {
    // Find cell in next lattice element.
    p.lowest_coord().universe() = lat[coord.lattice_index()];
    bool found = exhaustive_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 = exhaustive_find_cell(p);
      if (!found) {
        p.mark_as_lost(fmt::format("Particle {} could not be located after "
                                   "crossing a boundary of lattice {}",
          p.id(), lat.id_));
      }
    }
  }
}

//==============================================================================

BoundaryInfo distance_to_boundary(GeometryState& p)
{
  BoundaryInfo info;
  double d_lat = INFINITY;
  double d_surf = INFINITY;
  int32_t level_surf_cross;
  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)};
    const Position& r {coord.r()};
    const 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);
    d_surf = surface_distance.first;
    level_surf_cross = surface_distance.second;

    // Find the distance to the next lattice tile crossing.
    if (coord.lattice() != C_NONE) {
      auto& lat {*model::lattices[coord.lattice()]};
      // TODO: refactor so both lattice use the same position argument (which
      // also means the lat.type attribute can be removed)
      std::pair<double, array<int, 3>> lattice_distance;
      switch (lat.type_) {
      case LatticeType::rect:
        lattice_distance = lat.distance(r, u, coord.lattice_index());
        break;
      case LatticeType::hex:
        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_);
        }
        r_hex.z = coord.r().z;
        lattice_distance = lat.distance(r_hex, u, coord.lattice_index());
        break;
      }
      d_lat = lattice_distance.first;
      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()));
      }
    }

    // If the boundary on this coordinate level is coincident with a boundary on
    // a higher level then we need to make sure that the higher level boundary
    // is selected.  This logic must consider floating point precision.
    double& d = info.distance();
    if (d_surf < d_lat - FP_COINCIDENT) {
      if (d == INFINITY || (d - d_surf) / d >= FP_REL_PRECISION) {
        // Update closest distance
        d = d_surf;

        // If the cell is not simple, it is possible that both the negative and
        // positive half-space were given in the region specification. Thus, we
        // have to explicitly check which half-space the particle would be
        // traveling into if the surface is crossed
        if (c.is_simple() || d == INFTY) {
          info.surface() = level_surf_cross;
        } else {
          Position r_hit = r + d_surf * u;
          Surface& surf {*model::surfaces[std::abs(level_surf_cross) - 1]};
          Direction norm = surf.normal(r_hit);
          if (u.dot(norm) > 0) {
            info.surface() = std::abs(level_surf_cross);
          } else {
            info.surface() = -std::abs(level_surf_cross);
          }
        }

        info.lattice_translation()[0] = 0;
        info.lattice_translation()[1] = 0;
        info.lattice_translation()[2] = 0;
        info.coord_level() = i + 1;
      }
    } else {
      if (d == INFINITY || (d - d_lat) / d >= FP_REL_PRECISION) {
        d = d_lat;
        info.surface() = SURFACE_NONE;
        info.lattice_translation() = level_lat_trans;
        info.coord_level() = i + 1;
      }
    }
  }
  return info;
}

//==============================================================================
// C API
//==============================================================================

extern "C" int openmc_find_cell(
  const double* xyz, int32_t* index, int32_t* instance)
{
  GeometryState geom_state;

  geom_state.r() = Position {xyz};
  geom_state.u() = {0.0, 0.0, 1.0};

  if (!exhaustive_find_cell(geom_state)) {
    set_errmsg(
      fmt::format("Could not find cell at position {}.", geom_state.r()));
    return OPENMC_E_GEOMETRY;
  }

  *index = geom_state.lowest_coord().cell();
  *instance = geom_state.cell_instance();
  return 0;
}

extern "C" int openmc_global_bounding_box(double* llc, double* urc)
{
  auto bbox = model::universes.at(model::root_universe)->bounding_box();

  // set lower left corner values
  llc[0] = bbox.min.x;
  llc[1] = bbox.min.y;
  llc[2] = bbox.min.z;

  // set upper right corner values
  urc[0] = bbox.max.x;
  urc[1] = bbox.max.y;
  urc[2] = bbox.max.z;

  return 0;
}

} // namespace openmc

1	#include "openmc/geometry.h"
2
3	#include <fmt/core.h>
4	#include <fmt/ostream.h>
5
6	#include "openmc/array.h"
7	#include "openmc/cell.h"
8	#include "openmc/constants.h"
9	#include "openmc/error.h"
10	#include "openmc/lattice.h"
11	#include "openmc/settings.h"
12	#include "openmc/simulation.h"
13	#include "openmc/string_utils.h"
14	#include "openmc/surface.h"
15
16	namespace openmc {
17
18	//==============================================================================
19	// Global variables
20	//==============================================================================
21
22	namespace model {
23
24	int root_universe {-1};
25	int n_coord_levels;
26
27	vector<int64_t> overlap_check_count;
28
29	vector<OverlapKey> overlap_keys;
30	std::unordered_map<OverlapKey, int, OverlapKeyHash> overlap_key_index;
31
32	} // namespace model
33
34	//==============================================================================
35	// Non-member functions
36	//==============================================================================
37
38	int check_cell_overlap(GeometryState& p, bool error)	1,442,100✔
39	{
40	int n_coord = p.n_coord();	1,442,100✔
41
42	// If no overlap found, return a nonphysical index
43	int overlap_index = -1;	1,442,100✔
44
45	// Loop through each coordinate level
46	for (int j = 0; j < n_coord; j++) {	2,884,200✔
47	Universe& univ = *model::universes[p.coord(j).universe()];	1,442,100✔
48
49	// Loop through each cell on this level
50	for (auto index_cell : univ.cells_) {	4,968,568✔
51	Cell& c = *model::cells[index_cell];	3,920,488✔
52	if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) {	3,920,488✔
53	#pragma omp atomic	759,036✔
54	++model::overlap_check_count[index_cell];	1,391,566✔
55	if (index_cell != p.coord(j).cell()) {	1,391,566✔
56	if (error) {	394,020!
57	fatal_error(	×
58	fmt::format("Overlapping cells detected: {}, {} on universe {}",	×
59	c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_));	×
60	}
61
62	// With no fatal error (plotter is calling), now adds overlaps and
63	// ensures order does not matter when making overlap key
64	int cell_a = model::cells[index_cell]->id_;	394,020!
65	int cell_b = model::cells[p.coord(j).cell()]->id_;	394,020!
66	int a = std::min(cell_a, cell_b);	394,020!
67	int b = std::max(cell_a, cell_b);	394,020!
68	OverlapKey key {univ.id_, a, b};	394,020✔
69	#pragma omp critical(overlap_key_update)	429,840✔
70	{	394,020✔
71	auto it = model::overlap_key_index.find(key);	394,020✔
72	if (it != model::overlap_key_index.end()) {	394,020✔
73	overlap_index = it->second; // already exists, reuse index	393,965✔
74	} else {
75	int idx = int(model::overlap_keys.size());	55✔
76	model::overlap_keys.push_back(key);	55✔
77	model::overlap_key_index[key] = idx;	55✔
78	overlap_index = idx;	55✔
79	}
80	}
81	break;	394,020✔
82	}
83	}
84	}
85	}
86	return overlap_index;	1,442,100✔
87	}
88
89	//==============================================================================
90
91	int cell_instance_at_level(const GeometryState& p, int level)	2,147,483,647✔
92	{
93	// throw error if the requested level is too deep for the geometry
94	if (level > model::n_coord_levels) {	2,147,483,647!
95	fatal_error(fmt::format("Cell instance at level {} requested, but only {} "	×
96	"levels exist in the geometry.",
97	level, p.n_coord()));
98	}
99
100	// determine the cell instance
101	Cell& c {*model::cells[p.coord(level).cell()]};	2,147,483,647!
102
103	// quick exit if this cell doesn't have distribcell instances
104	if (c.distribcell_index_ == C_NONE)	2,147,483,647!
105	return C_NONE;
106
107	// compute the cell's instance
108	int instance = 0;
109	for (int i = 0; i < level; i++) {	2,147,483,647✔
110	const auto& c_i {*model::cells[p.coord(i).cell()]};	2,147,483,647✔
111	if (c_i.type_ == Fill::UNIVERSE) {	2,147,483,647✔
112	instance += c_i.offset_[c.distribcell_index_];	1,892,151,861✔
113	} else if (c_i.type_ == Fill::LATTICE) {	1,312,500,812!
114	instance += c_i.offset_[c.distribcell_index_];	1,312,500,812✔
115	auto& lat {*model::lattices[p.coord(i + 1).lattice()]};	1,312,500,812✔
116	const auto& i_xyz {p.coord(i + 1).lattice_index()};	1,312,500,812✔
117	if (lat.are_valid_indices(i_xyz)) {	1,312,500,812✔
118	instance += lat.offset(c.distribcell_index_, i_xyz);	1,307,616,845✔
119	}
120	}
121	}
122	return instance;
123	}
124
125	//==============================================================================
126
127	bool find_cell_inner(	2,147,483,647✔
128	GeometryState& p, const NeighborList* neighbor_list, bool verbose)
129	{
130	// Find which cell of this universe the particle is in. Use the neighbor list
131	// to shorten the search if one was provided.
132	bool found = false;	2,147,483,647✔
133	int32_t i_cell = C_NONE;	2,147,483,647✔
134	if (neighbor_list) {	2,147,483,647✔
135	for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) {	2,147,483,647✔
136	i_cell = *it;	2,147,483,647!
137
138	// Make sure the search cell is in the same universe.
139	int i_universe = p.lowest_coord().universe();	2,147,483,647!
140	if (model::cells[i_cell]->universe_ != i_universe)	2,147,483,647!
141	continue;	×
142
143	// Check if this cell contains the particle.
144	Position r {p.r_local()};	2,147,483,647✔
145	Direction u {p.u_local()};	2,147,483,647✔
146	auto surf = p.surface();	2,147,483,647✔
147	if (model::cells[i_cell]->contains(r, u, surf)) {	2,147,483,647✔
148	p.lowest_coord().cell() = i_cell;	2,147,483,647✔
149	found = true;	2,147,483,647✔
150	break;	2,147,483,647✔
151	}
152	}
153
154	// If we're attempting a neighbor list search and fail, we
155	// now know we should return false. This will trigger an
156	// exhaustive search from neighbor_list_find_cell and make
157	// the result from that be appended to the neighbor list.
158	if (!found) {	2,147,483,647✔
159	return found;
160	}
161	}
162
163	// Check successively lower coordinate levels until finding material fill
164	for (;; ++p.n_coord()) {	2,147,483,647✔
165	// If we did not attempt to use neighbor lists, i_cell is still C_NONE. In
166	// that case, we should now do an exhaustive search to find the right value
167	// of i_cell.
168	//
169	// Alternatively, neighbor list searches could have succeeded, but we found
170	// that the fill of the neighbor cell was another universe. As such, in the
171	// code below this conditional, we set i_cell back to C_NONE to indicate
172	// that.
173	if (i_cell == C_NONE) {	2,147,483,647✔
174	int i_universe = p.lowest_coord().universe();	2,147,483,647✔
175	const auto& univ {model::universes[i_universe]};	2,147,483,647✔
176	found = univ->find_cell(p);	2,147,483,647✔
177	}
178
179	if (!found) {	2,147,483,647✔
180	return found;
181	}
182	i_cell = p.lowest_coord().cell();	2,147,483,647!
183
184	// Announce the cell that the particle is entering.
185	if (found && verbose) {	2,147,483,647!
186	auto msg = fmt::format(" Entering cell {}", model::cells[i_cell]->id_);	×
187	write_message(msg, 1);	×
188	}	×
189
190	Cell& c {*model::cells[i_cell]};	2,147,483,647✔
191	if (c.type_ == Fill::MATERIAL) {	2,147,483,647✔
192	// Found a material cell which means this is the lowest coord level.
193
194	p.cell_instance() = 0;	2,147,483,647✔
195	// Find the distribcell instance number.
196	if (c.distribcell_index_ >= 0) {	2,147,483,647✔
197	p.cell_instance() = cell_instance_at_level(p, p.n_coord() - 1);	2,147,483,647✔
198	}
199
200	// Set the material, temperature and density multiplier.
201	p.material_last() = p.material();	2,147,483,647✔
202	p.material() = c.material(p.cell_instance());	2,147,483,647✔
203	p.sqrtkT_last() = p.sqrtkT();	2,147,483,647✔
204	p.sqrtkT() = c.sqrtkT(p.cell_instance());	2,147,483,647✔
205	p.density_mult_last() = p.density_mult();	2,147,483,647✔
206	p.density_mult() = c.density_mult(p.cell_instance());	2,147,483,647✔
207
208	return true;	2,147,483,647✔
209
210	} else if (c.type_ == Fill::UNIVERSE) {	1,156,349,926✔
211	//========================================================================
212	//! Found a lower universe, update this coord level then search the next.
213
214	// Set the lower coordinate level universe.
215	auto& coord {p.coord(p.n_coord())};	1,003,639,390✔
216	coord.universe() = c.fill_;	1,003,639,390✔
217
218	// Set the position and direction.
219	coord.r() = p.r_local();	1,003,639,390✔
220	coord.u() = p.u_local();	1,003,639,390✔
221
222	// Apply translation.
223	coord.r() -= c.translation_;	1,003,639,390✔
224
225	// Apply rotation.
226	if (!c.rotation_.empty()) {	1,003,639,390✔
227	coord.rotate(c.rotation_);	66,157,597✔
228	}
229
230	} else if (c.type_ == Fill::LATTICE) {	152,710,536!
231	//========================================================================
232	//! Found a lower lattice, update this coord level then search the next.
233
234	Lattice& lat {*model::lattices[c.fill_]};	152,710,536✔
235
236	// Set the position and direction.
237	auto& coord {p.coord(p.n_coord())};	152,710,536✔
238	coord.r() = p.r_local();	152,710,536✔
239	coord.u() = p.u_local();	152,710,536✔
240
241	// Apply translation.
242	coord.r() -= c.translation_;	152,710,536✔
243
244	// Apply rotation.
245	if (!c.rotation_.empty()) {	152,710,536✔
246	coord.rotate(c.rotation_);	358,336✔
247	}
248
249	// Determine lattice indices.
250	auto& i_xyz {coord.lattice_index()};	152,710,536✔
251	lat.get_indices(coord.r(), coord.u(), i_xyz);	152,710,536✔
252
253	// Get local position in appropriate lattice cell
254	coord.r() = lat.get_local_position(coord.r(), i_xyz);	152,710,536✔
255
256	// Set lattice indices.
257	coord.lattice() = c.fill_;	152,710,536✔
258
259	// Set the lower coordinate level universe.
260	if (lat.are_valid_indices(i_xyz)) {	152,710,536✔
261	coord.universe() = lat[i_xyz];	147,826,569✔
262	} else {
263	if (lat.outer_ != NO_OUTER_UNIVERSE) {	4,883,967!
264	coord.universe() = lat.outer_;	4,883,967✔
265	} else {
266	p.mark_as_lost(fmt::format(	×
267	"Particle {} left lattice {}, but it has no outer definition.",
268	p.id(), lat.id_));	×
269	}
270	}
271	}
272	i_cell = C_NONE; // trip non-neighbor cell search at next iteration	1,156,349,926✔
273	found = false;	1,156,349,926✔
274	}
275
276	return found;
277	}
278
279	//==============================================================================
280
281	bool neighbor_list_find_cell(GeometryState& p, bool verbose)	2,147,483,647✔
282	{
283
284	// Reset all the deeper coordinate levels.
285	for (int i = p.n_coord(); i < model::n_coord_levels; i++) {	2,147,483,647✔
286	p.coord(i).reset();	1,917,322,550✔
287	}
288
289	// Get the cell this particle was in previously.
290	auto coord_lvl = p.n_coord() - 1;	2,147,483,647✔
291	auto i_cell = p.coord(coord_lvl).cell();	2,147,483,647✔
292	Cell& c {*model::cells[i_cell]};	2,147,483,647✔
293
294	// Search for the particle in that cell's neighbor list. Return if we
295	// found the particle.
296	bool found = find_cell_inner(p, &c.neighbors_, verbose);	2,147,483,647✔
297	if (found)	2,147,483,647✔
298	return found;
299
300	// The particle could not be found in the neighbor list. Try searching all
301	// cells in this universe, and update the neighbor list if we find a new
302	// neighboring cell.
303	found = find_cell_inner(p, nullptr, verbose);	1,601,781✔
304	if (found)	1,601,781✔
305	c.neighbors_.push_back(p.coord(coord_lvl).cell());	30,022✔
306	return found;
307	}
308
309	void reconcile_cell_after_collision(GeometryState& p)	66✔
310	{
311	// Find the first coordinate level whose current cell is inconsistent with
312	// the particle's post-collision direction.
313	int invalid_level = C_NONE;	66✔
314	for (int level = 0; level < p.n_coord(); ++level) {	154✔
315	const auto& coord {p.coord(level)};	99!
316	if (coord.cell() == C_NONE \|\| !model::cells[coord.cell()]->contains(	99!
317	coord.r(), coord.u(), SURFACE_NONE)) {	99✔
318	invalid_level = level;
319	break;
320	}
321	}
322
323	if (invalid_level == C_NONE)	66✔
324	return;
325
326	// Search from the first inconsistent level so that parent cells and
327	// transformed lower universes are both handled correctly.
328	if (p.coord(invalid_level).cell() != C_NONE) {	11!
329	p.n_coord() = invalid_level + 1;	11✔
330	if (neighbor_list_find_cell(p))	11!
331	return;
332	}
333
334	// The current cell may not have a complete neighbor list yet. Fall back to
335	// the normal exhaustive search used after a failed surface crossing.
NEW 336	p.n_coord() = 1;	×
NEW 337	if (!exhaustive_find_cell(p)) {	×
NEW 338	p.mark_as_lost("Could not find particle after a collision near a surface.");	×
339	}
340	}
341
342	bool exhaustive_find_cell(GeometryState& p, bool verbose)	1,303,474,823✔
343	{
344	int i_universe = p.lowest_coord().universe();	1,303,474,823✔
345	if (i_universe == C_NONE) {	1,303,474,823✔
346	p.coord(0).universe() = model::root_universe;	267,620,950✔
347	p.n_coord() = 1;	267,620,950✔
348	i_universe = model::root_universe;	267,620,950✔
349	}
350	// Reset all the deeper coordinate levels.
351	for (int i = p.n_coord(); i < model::n_coord_levels; i++) {	1,393,744,413✔
352	p.coord(i).reset();	90,269,590✔
353	}
354	return find_cell_inner(p, nullptr, verbose);	1,303,474,823✔
355	}
356
357	//==============================================================================
358
359	void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)	813,765,152✔
360	{
361	auto& coord {p.lowest_coord()};	813,765,152!
362	auto& lat {*model::lattices[coord.lattice()]};	813,765,152!
363
364	if (verbose) {	813,765,152!
365	write_message(	×
366	fmt::format(" Crossing lattice {}. Current position ({},{},{}). r={}",	×
367	lat.id_, coord.lattice_index()[0], coord.lattice_index()[1],	×
368	coord.lattice_index()[2], p.r()),	×
369	1);
370	}
371
372	// Set the lattice indices.
373	coord.lattice_index()[0] += boundary.lattice_translation()[0];	813,765,152✔
374	coord.lattice_index()[1] += boundary.lattice_translation()[1];	813,765,152✔
375	coord.lattice_index()[2] += boundary.lattice_translation()[2];	813,765,152✔
376
377	// Set the new coordinate position.
378	const auto& upper_coord {p.coord(p.n_coord() - 2)};	813,765,152✔
379	const auto& cell {model::cells[upper_coord.cell()]};	813,765,152✔
380	Position r = upper_coord.r();	813,765,152✔
381	r -= cell->translation_;	813,765,152✔
382	if (!cell->rotation_.empty()) {	813,765,152✔
383	r = r.rotate(cell->rotation_);	471,647✔
384	}
385	p.r_local() = lat.get_local_position(r, coord.lattice_index());	813,765,152✔
386
387	if (!lat.are_valid_indices(coord.lattice_index())) {	813,765,152✔
388	// The particle is outside the lattice. Search for it from the base coords.
389	p.n_coord() = 1;	3,714,425✔
390	bool found = exhaustive_find_cell(p);	3,714,425✔
391
392	if (!found) {	3,714,425!
393	p.mark_as_lost(fmt::format("Particle {} could not be located after "	×
394	"crossing a boundary of lattice {}",
395	p.id(), lat.id_));	×
396	}
397
398	} else {
399	// Find cell in next lattice element.
400	p.lowest_coord().universe() = lat[coord.lattice_index()];	810,050,727✔
401	bool found = exhaustive_find_cell(p);	810,050,727✔
402
403	if (!found) {	810,050,727!
404	// A particle crossing the corner of a lattice tile may not be found. In
405	// this case, search for it from the base coords.
406	p.n_coord() = 1;	×
407	bool found = exhaustive_find_cell(p);	×
408	if (!found) {	×
409	p.mark_as_lost(fmt::format("Particle {} could not be located after "	×
410	"crossing a boundary of lattice {}",
411	p.id(), lat.id_));	×
412	}
413	}
414	}
415	}	813,765,152✔
416
417	//==============================================================================
418
419	BoundaryInfo distance_to_boundary(GeometryState& p)	2,147,483,647✔
420	{
421	BoundaryInfo info;	2,147,483,647✔
422	double d_lat = INFINITY;	2,147,483,647✔
423	double d_surf = INFINITY;	2,147,483,647✔
424	int32_t level_surf_cross;	2,147,483,647✔
425	array<int, 3> level_lat_trans {};	2,147,483,647✔
426
427	// Loop over each coordinate level.
428	for (int i = 0; i < p.n_coord(); i++) {	2,147,483,647✔
429	const auto& coord {p.coord(i)};	2,147,483,647✔
430	const Position& r {coord.r()};	2,147,483,647✔
431	const Direction& u {coord.u()};	2,147,483,647✔
432	Cell& c {*model::cells[coord.cell()]};	2,147,483,647✔
433
434	// Find the oncoming surface in this cell and the distance to it.
435	auto surface_distance = c.distance(r, u, p.surface(), &p);	2,147,483,647✔
436	d_surf = surface_distance.first;	2,147,483,647✔
437	level_surf_cross = surface_distance.second;	2,147,483,647✔
438
439	// Find the distance to the next lattice tile crossing.
440	if (coord.lattice() != C_NONE) {	2,147,483,647✔
441	auto& lat {*model::lattices[coord.lattice()]};	1,406,189,327!
442	// TODO: refactor so both lattice use the same position argument (which
443	// also means the lat.type attribute can be removed)
444	std::pair<double, array<int, 3>> lattice_distance;	1,406,189,327✔
445	switch (lat.type_) {	1,406,189,327!
446	case LatticeType::rect:	1,320,686,030✔
447	lattice_distance = lat.distance(r, u, coord.lattice_index());	1,320,686,030✔
448	break;	1,320,686,030✔
449	case LatticeType::hex:	85,503,297✔
450	auto& cell_above {model::cells[p.coord(i - 1).cell()]};	85,503,297✔
451	Position r_hex {p.coord(i - 1).r()};	85,503,297✔
452	r_hex -= cell_above->translation_;	85,503,297✔
453	if (coord.rotated()) {	85,503,297✔
454	r_hex = r_hex.rotate(cell_above->rotation_);	701,954✔
455	}
456	r_hex.z = coord.r().z;	85,503,297✔
457	lattice_distance = lat.distance(r_hex, u, coord.lattice_index());	85,503,297✔
458	break;	85,503,297✔
459	}
460	d_lat = lattice_distance.first;	1,406,189,327✔
461	level_lat_trans = lattice_distance.second;	1,406,189,327✔
462
463	if (d_lat < 0) {	1,406,189,327!
464	p.mark_as_lost(fmt::format("Particle {} had a negative distance "	×
465	"to a lattice boundary.",
466	p.id()));	×
467	}
468	}
469
470	// If the boundary on this coordinate level is coincident with a boundary on
471	// a higher level then we need to make sure that the higher level boundary
472	// is selected. This logic must consider floating point precision.
473	double& d = info.distance();	2,147,483,647✔
474	if (d_surf < d_lat - FP_COINCIDENT) {	2,147,483,647✔
475	if (d == INFINITY \|\| (d - d_surf) / d >= FP_REL_PRECISION) {	2,147,483,647✔
476	// Update closest distance
477	d = d_surf;	2,147,483,647✔
478
479	// If the cell is not simple, it is possible that both the negative and
480	// positive half-space were given in the region specification. Thus, we
481	// have to explicitly check which half-space the particle would be
482	// traveling into if the surface is crossed
483	if (c.is_simple() \|\| d == INFTY) {	2,147,483,647✔
484	info.surface() = level_surf_cross;	2,147,483,647✔
485	} else {
486	Position r_hit = r + d_surf * u;	186,569,834✔
487	Surface& surf {*model::surfaces[std::abs(level_surf_cross) - 1]};	186,569,834✔
488	Direction norm = surf.normal(r_hit);	186,569,834✔
489	if (u.dot(norm) > 0) {	186,569,834✔
490	info.surface() = std::abs(level_surf_cross);	176,702,200✔
491	} else {
492	info.surface() = -std::abs(level_surf_cross);	9,867,634✔
493	}
494	}
495
496	info.lattice_translation()[0] = 0;	2,147,483,647✔
497	info.lattice_translation()[1] = 0;	2,147,483,647✔
498	info.lattice_translation()[2] = 0;	2,147,483,647✔
499	info.coord_level() = i + 1;	2,147,483,647✔
500	}
501	} else {
502	if (d == INFINITY \|\| (d - d_lat) / d >= FP_REL_PRECISION) {	1,171,653,347!
503	d = d_lat;	962,488,272✔
504	info.surface() = SURFACE_NONE;	962,488,272✔
505	info.lattice_translation() = level_lat_trans;	962,488,272✔
506	info.coord_level() = i + 1;	962,488,272✔
507	}
508	}
509	}
510	return info;	2,147,483,647✔
511	}
512
513	//==============================================================================
514	// C API
515	//==============================================================================
516
517	extern "C" int openmc_find_cell(	600,308✔
518	const double* xyz, int32_t* index, int32_t* instance)
519	{
520	GeometryState geom_state;	600,308✔
521
522	geom_state.r() = Position {xyz};	600,308✔
523	geom_state.u() = {0.0, 0.0, 1.0};	600,308✔
524
525	if (!exhaustive_find_cell(geom_state)) {	600,308✔
526	set_errmsg(	11✔
527	fmt::format("Could not find cell at position {}.", geom_state.r()));	11✔
528	return OPENMC_E_GEOMETRY;	11✔
529	}
530
531	*index = geom_state.lowest_coord().cell();	600,297✔
532	*instance = geom_state.cell_instance();	600,297✔
533	return 0;	600,297✔
534	}	600,308✔
535
536	extern "C" int openmc_global_bounding_box(double* llc, double* urc)	11✔
537	{
538	auto bbox = model::universes.at(model::root_universe)->bounding_box();	11✔
539
540	// set lower left corner values
541	llc[0] = bbox.min.x;	11✔
542	llc[1] = bbox.min.y;	11✔
543	llc[2] = bbox.min.z;	11✔
544
545	// set upper right corner values
546	urc[0] = bbox.max.x;	11✔
547	urc[1] = bbox.max.y;	11✔
548	urc[2] = bbox.max.z;	11✔
549
550	return 0;	11✔
551	}
552
553	} // namespace openmc

openmc-dev / openmc / 30828111152

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous