13894606366

Committed 17 Mar 2025 08:15AM UTC coverage: 84.924% (-0.1%) from 85.042%

Build # 13894606366

Build Type

Pull #3305

github

Committed by

web-flow

Commit Message

Merge 042852a3c into 58f2a2177

Pull Request Pull Request #3305: New Feature Development: The Implementation of chord length sampling (CLS) method for stochastic media( #3286 )

Run Details

20 of 161 new or added lines in 8 files covered. (12.42%)

1043 existing lines in 45 files now uncovered.

51567 of 60721 relevant lines covered (84.92%)

36823581.23 hits per line

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

89.87

/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"
#include <openmc/stochastic_media.h>

namespace openmc {

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

namespace model {

int root_universe {-1};
int n_coord_levels;

vector<int64_t> overlap_check_count;

} // namespace model

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

bool check_cell_overlap(GeometryState& p, bool error)
{
  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];

    // 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 (error) {
            fatal_error(
              fmt::format("Overlapping cells detected: {}, {} on universe {}",
                c.id_, model::cells[p.coord(j).cell]->id_, univ.id_));
          }
          return true;
        }
#pragma omp atomic
        ++model::overlap_check_count[index_cell];
      }
    }
  }

  return false;
}

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

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_i};
      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 and temperature.
      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.status() = ParticleStatus::OUTSIDE;

      return true;

    } else if (c.type_ == Fill::STOCHASTIC_MEDIA) {
      //========================================================================
      //! Found stochastic media cell, means this is the lowest coord level.
      model::stochastic_media[c.fill_]->sample_material(p);
      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_i};
      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;
}

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_i[0], coord.lattice_i[1], coord.lattice_i[2],
        p.r()),
      1);
  }

  // Set the lattice indices.
  coord.lattice_i[0] += boundary.lattice_translation[0];
  coord.lattice_i[1] += boundary.lattice_translation[1];
  coord.lattice_i[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_i);

  if (!lat.are_valid_indices(coord.lattice_i)) {
    // 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_i];
    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_i);
        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_i);
        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.xmin;
  llc[1] = bbox.ymin;
  llc[2] = bbox.zmin;

  // set upper right corner values
  urc[0] = bbox.xmax;
  urc[1] = bbox.ymax;
  urc[2] = bbox.zmax;

  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	#include <openmc/stochastic_media.h>
16
17	namespace openmc {
18
19	//==============================================================================
20	// Global variables
21	//==============================================================================
22
23	namespace model {
24
25	int root_universe {-1};
26	int n_coord_levels;
27
28	vector<int64_t> overlap_check_count;
29
30	} // namespace model
31
32	//==============================================================================
33	// Non-member functions
34	//==============================================================================
35
36	bool check_cell_overlap(GeometryState& p, bool error)	1,320,000✔
37	{
38	int n_coord = p.n_coord();	1,320,000✔
39
40	// Loop through each coordinate level
41	for (int j = 0; j < n_coord; j++) {	2,266,044✔
42	Universe& univ = *model::universes[p.coord(j).universe];	1,320,000✔
43
44	// Loop through each cell on this level
45	for (auto index_cell : univ.cells_) {	4,466,044✔
46	Cell& c = *model::cells[index_cell];	3,520,000✔
47	if (c.contains(p.coord(j).r, p.coord(j).u, p.surface())) {	3,520,000✔
48	if (index_cell != p.coord(j).cell) {	1,249,622✔
49	if (error) {	373,956✔
50	fatal_error(	×
51	fmt::format("Overlapping cells detected: {}, {} on universe {}",	×
52	c.id_, model::cells[p.coord(j).cell]->id_, univ.id_));	×
53	}
54	return true;	373,956✔
55	}
56	#pragma omp atomic	477,636✔
57	++model::overlap_check_count[index_cell];	875,666✔
58	}
59	}
60	}
61
62	return false;	946,044✔
63	}
64
65	//==============================================================================
66
67	int cell_instance_at_level(const GeometryState& p, int level)	2,147,483,647✔
68	{
69	// throw error if the requested level is too deep for the geometry
70	if (level > model::n_coord_levels) {	2,147,483,647✔
71	fatal_error(fmt::format("Cell instance at level {} requested, but only {} "	×
72	"levels exist in the geometry.",
73	level, p.n_coord()));
74	}
75
76	// determine the cell instance
77	Cell& c {*model::cells[p.coord(level).cell]};	2,147,483,647✔
78
79	// quick exit if this cell doesn't have distribcell instances
80	if (c.distribcell_index_ == C_NONE)	2,147,483,647✔
81	return C_NONE;	×
82
83	// compute the cell's instance
84	int instance = 0;	2,147,483,647✔
85	for (int i = 0; i < level; i++) {	2,147,483,647✔
86	const auto& c_i {*model::cells[p.coord(i).cell]};	1,956,762,317✔
87	if (c_i.type_ == Fill::UNIVERSE) {	1,956,762,317✔
88	instance += c_i.offset_[c.distribcell_index_];	866,022,912✔
89	} else if (c_i.type_ == Fill::LATTICE) {	1,090,739,405✔
90	instance += c_i.offset_[c.distribcell_index_];	1,090,739,405✔
91	auto& lat {*model::lattices[p.coord(i + 1).lattice]};	1,090,739,405✔
92	const auto& i_xyz {p.coord(i + 1).lattice_i};	1,090,739,405✔
93	if (lat.are_valid_indices(i_xyz)) {	1,090,739,405✔
94	instance += lat.offset(c.distribcell_index_, i_xyz);	1,085,853,480✔
95	}
96	}
97	}
98	return instance;	2,147,483,647✔
99	}
100
101	//==============================================================================
102
103	bool find_cell_inner(	2,147,483,647✔
104	GeometryState& p, const NeighborList* neighbor_list, bool verbose)
105	{
106	// Find which cell of this universe the particle is in. Use the neighbor list
107	// to shorten the search if one was provided.
108	bool found = false;	2,147,483,647✔
109	int32_t i_cell = C_NONE;	2,147,483,647✔
110	if (neighbor_list) {	2,147,483,647✔
111	for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) {	1,576,654,721✔
112	i_cell = *it;	1,575,059,146✔
113
114	// Make sure the search cell is in the same universe.
115	int i_universe = p.lowest_coord().universe;	1,575,059,146✔
116	if (model::cells[i_cell]->universe_ != i_universe)	1,575,059,146✔
117	continue;	×
118
119	// Check if this cell contains the particle.
120	Position r {p.r_local()};	1,575,059,146✔
121	Direction u {p.u_local()};	1,575,059,146✔
122	auto surf = p.surface();	1,575,059,146✔
123	if (model::cells[i_cell]->contains(r, u, surf)) {	1,575,059,146✔
124	p.lowest_coord().cell = i_cell;	1,280,498,179✔
125	found = true;	1,280,498,179✔
126	break;	1,280,498,179✔
127	}
128	}
129
130	// If we're attempting a neighbor list search and fail, we
131	// now know we should return false. This will trigger an
132	// exhaustive search from neighbor_list_find_cell and make
133	// the result from that be appended to the neighbor list.
134	if (!found) {	1,282,093,754✔
135	return found;	1,595,575✔
136	}
137	}
138
139	// Check successively lower coordinate levels until finding material fill
140	for (;; ++p.n_coord()) {	651,349,434✔
141	// If we did not attempt to use neighbor lists, i_cell is still C_NONE. In
142	// that case, we should now do an exhaustive search to find the right value
143	// of i_cell.
144	//
145	// Alternatively, neighbor list searches could have succeeded, but we found
146	// that the fill of the neighbor cell was another universe. As such, in the
147	// code below this conditional, we set i_cell back to C_NONE to indicate
148	// that.
149	if (i_cell == C_NONE) {	2,147,483,647✔
150	int i_universe = p.lowest_coord().universe;	1,478,721,424✔
151	const auto& univ {model::universes[i_universe]};	1,478,721,424✔
152	found = univ->find_cell(p);	1,478,721,424✔
153	}
154
155	if (!found) {	2,147,483,647✔
156	return found;	173,970,759✔
157	}
158	i_cell = p.lowest_coord().cell;	2,147,483,647✔
159
160	// Announce the cell that the particle is entering.
161	if (found && verbose) {	2,147,483,647✔
162	auto msg = fmt::format(" Entering cell {}", model::cells[i_cell]->id_);	×
163	write_message(msg, 1);	×
164	}
165
166	Cell& c {*model::cells[i_cell]};	2,147,483,647✔
167	if (c.type_ == Fill::MATERIAL) {	2,147,483,647✔
168	// Found a material cell which means this is the lowest coord level.
169
170	p.cell_instance() = 0;	2,147,483,647✔
171	// Find the distribcell instance number.
172	if (c.distribcell_index_ >= 0) {	2,147,483,647✔
173	p.cell_instance() = cell_instance_at_level(p, p.n_coord() - 1);	2,147,483,647✔
174	}
175
176	// Set the material and temperature.
177	p.material_last() = p.material();	2,147,483,647✔
178	p.material() = c.material(p.cell_instance());	2,147,483,647✔
179	p.sqrtkT_last() = p.sqrtkT();	2,147,483,647✔
180	p.sqrtkT() = c.sqrtkT(p.cell_instance());	2,147,483,647✔
181	p.status() = ParticleStatus::OUTSIDE;	2,147,483,647✔
182
183	return true;	2,147,483,647✔
184
185	} else if (c.type_ == Fill::STOCHASTIC_MEDIA) {	325,674,717✔
186	//========================================================================
187	//! Found stochastic media cell, means this is the lowest coord level.
NEW 188	model::stochastic_media[c.fill_]->sample_material(p);	×
NEW 189	return true;	×
190
191	} else if (c.type_ == Fill::UNIVERSE) {	325,674,717✔
192	//========================================================================
193	//! Found a lower universe, update this coord level then search the next.
194
195	// Set the lower coordinate level universe.
196	auto& coord {p.coord(p.n_coord())};	203,108,157✔
197	coord.universe = c.fill_;	203,108,157✔
198
199	// Set the position and direction.
200	coord.r = p.r_local();	203,108,157✔
201	coord.u = p.u_local();	203,108,157✔
202
203	// Apply translation.
204	coord.r -= c.translation_;	203,108,157✔
205
206	// Apply rotation.
207	if (!c.rotation_.empty()) {	203,108,157✔
208	coord.rotate(c.rotation_);	2,189,110✔
209	}
210
211	} else if (c.type_ == Fill::LATTICE) {	122,566,560✔
212	//========================================================================
213	//! Found a lower lattice, update this coord level then search the next.
214
215	Lattice& lat {*model::lattices[c.fill_]};	122,566,560✔
216
217	// Set the position and direction.
218	auto& coord {p.coord(p.n_coord())};	122,566,560✔
219	coord.r = p.r_local();	122,566,560✔
220	coord.u = p.u_local();	122,566,560✔
221
222	// Apply translation.
223	coord.r -= c.translation_;	122,566,560✔
224
225	// Apply rotation.
226	if (!c.rotation_.empty()) {	122,566,560✔
227	coord.rotate(c.rotation_);	354,519✔
228	}
229
230	// Determine lattice indices.
231	auto& i_xyz {coord.lattice_i};	122,566,560✔
232	lat.get_indices(coord.r, coord.u, i_xyz);	122,566,560✔
233
234	// Get local position in appropriate lattice cell
235	coord.r = lat.get_local_position(coord.r, i_xyz);	122,566,560✔
236
237	// Set lattice indices.
238	coord.lattice = c.fill_;	122,566,560✔
239
240	// Set the lower coordinate level universe.
241	if (lat.are_valid_indices(i_xyz)) {	122,566,560✔
242	coord.universe = lat[i_xyz];	117,680,635✔
243	} else {
244	if (lat.outer_ != NO_OUTER_UNIVERSE) {	4,885,925✔
245	coord.universe = lat.outer_;	4,885,925✔
246	} else {
247	p.mark_as_lost(fmt::format(	×
248	"Particle {} left lattice {}, but it has no outer definition.",
249	p.id(), lat.id_));	×
250	}
251	}
252	}
253	i_cell = C_NONE; // trip non-neighbor cell search at next iteration	325,674,717✔
254	found = false;	325,674,717✔
255	}	325,674,717✔
256
257	return found;
258	}
259
260	//==============================================================================
261
262	bool neighbor_list_find_cell(GeometryState& p, bool verbose)	1,282,093,754✔
263	{
264
265	// Reset all the deeper coordinate levels.
266	for (int i = p.n_coord(); i < model::n_coord_levels; i++) {	1,967,802,112✔
267	p.coord(i).reset();	685,708,358✔
268	}
269
270	// Get the cell this particle was in previously.
271	auto coord_lvl = p.n_coord() - 1;	1,282,093,754✔
272	auto i_cell = p.coord(coord_lvl).cell;	1,282,093,754✔
273	Cell& c {*model::cells[i_cell]};	1,282,093,754✔
274
275	// Search for the particle in that cell's neighbor list. Return if we
276	// found the particle.
277	bool found = find_cell_inner(p, &c.neighbors_, verbose);	1,282,093,754✔
278	if (found)	1,282,093,754✔
279	return found;	1,280,498,179✔
280
281	// The particle could not be found in the neighbor list. Try searching all
282	// cells in this universe, and update the neighbor list if we find a new
283	// neighboring cell.
284	found = find_cell_inner(p, nullptr, verbose);	1,595,575✔
285	if (found)	1,595,575✔
286	c.neighbors_.push_back(p.coord(coord_lvl).cell);	26,544✔
287	return found;	1,595,575✔
288	}
289
290	bool exhaustive_find_cell(GeometryState& p, bool verbose)	1,151,451,132✔
291	{
292	int i_universe = p.lowest_coord().universe;	1,151,451,132✔
293	if (i_universe == C_NONE) {	1,151,451,132✔
294	p.coord(0).universe = model::root_universe;	257,325,753✔
295	p.n_coord() = 1;	257,325,753✔
296	i_universe = model::root_universe;	257,325,753✔
297	}
298	// Reset all the deeper coordinate levels.
299	for (int i = p.n_coord(); i < model::n_coord_levels; i++) {	1,228,241,607✔
300	p.coord(i).reset();	76,790,475✔
301	}
302	return find_cell_inner(p, nullptr, verbose);	1,151,451,132✔
303	}
304
305	//==============================================================================
306
307	void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)	671,876,182✔
308	{
309	auto& coord {p.lowest_coord()};	671,876,182✔
310	auto& lat {*model::lattices[coord.lattice]};	671,876,182✔
311
312	if (verbose) {	671,876,182✔
313	write_message(	×
314	fmt::format(" Crossing lattice {}. Current position ({},{},{}). r={}",	×
315	lat.id_, coord.lattice_i[0], coord.lattice_i[1], coord.lattice_i[2],	×
316	p.r()),
317	1);
318	}
319
320	// Set the lattice indices.
321	coord.lattice_i[0] += boundary.lattice_translation[0];	671,876,182✔
322	coord.lattice_i[1] += boundary.lattice_translation[1];	671,876,182✔
323	coord.lattice_i[2] += boundary.lattice_translation[2];	671,876,182✔
324
325	// Set the new coordinate position.
326	const auto& upper_coord {p.coord(p.n_coord() - 2)};	671,876,182✔
327	const auto& cell {model::cells[upper_coord.cell]};	671,876,182✔
328	Position r = upper_coord.r;	671,876,182✔
329	r -= cell->translation_;	671,876,182✔
330	if (!cell->rotation_.empty()) {	671,876,182✔
331	r = r.rotate(cell->rotation_);	470,965✔
332	}
333	p.r_local() = lat.get_local_position(r, coord.lattice_i);	671,876,182✔
334
335	if (!lat.are_valid_indices(coord.lattice_i)) {	671,876,182✔
336	// The particle is outside the lattice. Search for it from the base coords.
337	p.n_coord() = 1;	3,726,107✔
338	bool found = exhaustive_find_cell(p);	3,726,107✔
339
340	if (!found) {	3,726,107✔
341	p.mark_as_lost(fmt::format("Particle {} could not be located after "	×
342	"crossing a boundary of lattice {}",
343	p.id(), lat.id_));	×
344	}
345
346	} else {
347	// Find cell in next lattice element.
348	p.lowest_coord().universe = lat[coord.lattice_i];	668,150,075✔
349	bool found = exhaustive_find_cell(p);	668,150,075✔
350
351	if (!found) {	668,150,075✔
352	// A particle crossing the corner of a lattice tile may not be found. In
353	// this case, search for it from the base coords.
354	p.n_coord() = 1;	×
355	bool found = exhaustive_find_cell(p);	×
356	if (!found) {	×
357	p.mark_as_lost(fmt::format("Particle {} could not be located after "	×
358	"crossing a boundary of lattice {}",
359	p.id(), lat.id_));	×
360	}
361	}
362	}
363	}	671,876,182✔
364
365	//==============================================================================
366
367	BoundaryInfo distance_to_boundary(GeometryState& p)	2,147,483,647✔
368	{
369	BoundaryInfo info;	2,147,483,647✔
370	double d_lat = INFINITY;	2,147,483,647✔
371	double d_surf = INFINITY;	2,147,483,647✔
372	int32_t level_surf_cross;
373	array<int, 3> level_lat_trans {};	2,147,483,647✔
374
375	// Loop over each coordinate level.
376	for (int i = 0; i < p.n_coord(); i++) {	2,147,483,647✔
377	const auto& coord {p.coord(i)};	2,147,483,647✔
378	const Position& r {coord.r};	2,147,483,647✔
379	const Direction& u {coord.u};	2,147,483,647✔
380	Cell& c {*model::cells[coord.cell]};	2,147,483,647✔
381
382	// Find the oncoming surface in this cell and the distance to it.
383	auto surface_distance = c.distance(r, u, p.surface(), &p);	2,147,483,647✔
384	d_surf = surface_distance.first;	2,147,483,647✔
385	level_surf_cross = surface_distance.second;	2,147,483,647✔
386
387	// Find the distance to the next lattice tile crossing.
388	if (coord.lattice != C_NONE) {	2,147,483,647✔
389	auto& lat {*model::lattices[coord.lattice]};	1,124,791,135✔
390	// TODO: refactor so both lattice use the same position argument (which
391	// also means the lat.type attribute can be removed)
392	std::pair<double, array<int, 3>> lattice_distance;	1,124,791,135✔
393	switch (lat.type_) {	1,124,791,135✔
394	case LatticeType::rect:	1,039,072,645✔
395	lattice_distance = lat.distance(r, u, coord.lattice_i);	1,039,072,645✔
396	break;	1,039,072,645✔
397	case LatticeType::hex:	85,718,490✔
398	auto& cell_above {model::cells[p.coord(i - 1).cell]};	85,718,490✔
399	Position r_hex {p.coord(i - 1).r};	85,718,490✔
400	r_hex -= cell_above->translation_;	85,718,490✔
401	if (coord.rotated) {	85,718,490✔
402	r_hex = r_hex.rotate(cell_above->rotation_);	681,417✔
403	}
404	r_hex.z = coord.r.z;	85,718,490✔
405	lattice_distance = lat.distance(r_hex, u, coord.lattice_i);	85,718,490✔
406	break;	85,718,490✔
407	}
408	d_lat = lattice_distance.first;	1,124,791,135✔
409	level_lat_trans = lattice_distance.second;	1,124,791,135✔
410
411	if (d_lat < 0) {	1,124,791,135✔
412	p.mark_as_lost(fmt::format("Particle {} had a negative distance "	×
413	"to a lattice boundary.",
414	p.id()));	×
415	}
416	}
417
418	// If the boundary on this coordinate level is coincident with a boundary on
419	// a higher level then we need to make sure that the higher level boundary
420	// is selected. This logic must consider floating point precision.
421	double& d = info.distance;	2,147,483,647✔
422	if (d_surf < d_lat - FP_COINCIDENT) {	2,147,483,647✔
423	if (d == INFINITY \|\| (d - d_surf) / d >= FP_REL_PRECISION) {	2,147,483,647✔
424	// Update closest distance
425	d = d_surf;	2,147,483,647✔
426
427	// If the cell is not simple, it is possible that both the negative and
428	// positive half-space were given in the region specification. Thus, we
429	// have to explicitly check which half-space the particle would be
430	// traveling into if the surface is crossed
431	if (c.is_simple() \|\| d == INFTY) {	2,147,483,647✔
432	info.surface = level_surf_cross;	2,147,483,647✔
433	} else {
434	Position r_hit = r + d_surf * u;	12,743,337✔
435	Surface& surf {*model::surfaces[std::abs(level_surf_cross) - 1]};	12,743,337✔
436	Direction norm = surf.normal(r_hit);	12,743,337✔
437	if (u.dot(norm) > 0) {	12,743,337✔
438	info.surface = std::abs(level_surf_cross);	6,317,408✔
439	} else {
440	info.surface = -std::abs(level_surf_cross);	6,425,929✔
441	}
442	}
443
444	info.lattice_translation[0] = 0;	2,147,483,647✔
445	info.lattice_translation[1] = 0;	2,147,483,647✔
446	info.lattice_translation[2] = 0;	2,147,483,647✔
447	info.coord_level = i + 1;	2,147,483,647✔
448	}
449	} else {
450	if (d == INFINITY \|\| (d - d_lat) / d >= FP_REL_PRECISION) {	932,586,638✔
451	d = d_lat;	758,347,769✔
452	info.surface = SURFACE_NONE;	758,347,769✔
453	info.lattice_translation = level_lat_trans;	758,347,769✔
454	info.coord_level = i + 1;	758,347,769✔
455	}
456	}
457	}
458	return info;	2,147,483,647✔
459	}
460
461	//==============================================================================
462	// C API
463	//==============================================================================
464
465	extern "C" int openmc_find_cell(	600,308✔
466	const double* xyz, int32_t* index, int32_t* instance)
467	{
468	GeometryState geom_state;	600,308✔
469
470	geom_state.r() = Position {xyz};	600,308✔
471	geom_state.u() = {0.0, 0.0, 1.0};	600,308✔
472
473	if (!exhaustive_find_cell(geom_state)) {	600,308✔
474	set_errmsg(	11✔
475	fmt::format("Could not find cell at position {}.", geom_state.r()));	31✔
476	return OPENMC_E_GEOMETRY;	11✔
477	}
478
479	*index = geom_state.lowest_coord().cell;	600,297✔
480	*instance = geom_state.cell_instance();	600,297✔
481	return 0;	600,297✔
482	}	600,308✔
483
484	extern "C" int openmc_global_bounding_box(double* llc, double* urc)	11✔
485	{
486	auto bbox = model::universes.at(model::root_universe)->bounding_box();	11✔
487
488	// set lower left corner values
489	llc[0] = bbox.xmin;	11✔
490	llc[1] = bbox.ymin;	11✔
491	llc[2] = bbox.zmin;	11✔
492
493	// set upper right corner values
494	urc[0] = bbox.xmax;	11✔
495	urc[1] = bbox.ymax;	11✔
496	urc[2] = bbox.zmax;	11✔
497
498	return 0;	11✔
499	}
500
501	} // namespace openmc

openmc-dev / openmc / 13894606366

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