20368472998

Committed 19 Dec 2025 11:20AM UTC coverage: 82.143% (+0.004%) from 82.139%

Build # 20368472998

Build Type

Pull #3566

github

Committed by

web-flow

Commit Message

Merge a657bcb6c into a230b8612

Pull Request Pull Request #3566: Resolve merge conflicts for PR #3516 - Implement Virtual Lattice Method for Efficient Simulation of Dispersed Fuels

Run Details

17200 of 23810 branches covered (72.24%)

Branch coverage included in aggregate %.

266 of 322 new or added lines in 9 files covered. (82.61%)

2 existing lines in 1 file now uncovered.

55371 of 64537 relevant lines covered (85.8%)

43403485.97 hits per line

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

81.65

/src/universe.cpp

#include "openmc/universe.h"

#include <set>

#include "openmc/hdf5_interface.h"
#include "openmc/particle.h"

namespace openmc {

namespace model {

std::unordered_map<int32_t, int32_t> universe_map;
vector<unique_ptr<Universe>> universes;

} // namespace model

//==============================================================================
// Universe implementation
//==============================================================================

void Universe::to_hdf5(hid_t universes_group) const
{
  // Create a group for this universe.
  auto group = create_group(universes_group, fmt::format("universe {}", id_));

  // Write the geometry representation type.
  write_string(group, "geom_type", "csg", false);

  // Write the contained cells.
  if (cells_.size() > 0) {
    vector<int32_t> cell_ids;
    for (auto i_cell : cells_)
      cell_ids.push_back(model::cells[i_cell]->id_);
    write_dataset(group, "cells", cell_ids);
  }

  close_group(group);
}

bool Universe::find_cell(GeometryState& p) const
{
  if (filled_with_triso_base_ != -1) {
    bool found = find_cell_in_virtual_lattice(p);
    if (found) {
      return found;
    }
  }
  const auto& cells {
    !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};

  Position r {p.r_local()};
  Position u {p.u_local()};
  auto surf = p.surface();
  int32_t i_univ = p.lowest_coord().universe();

  for (auto i_cell : cells) {
    if (model::cells[i_cell]->universe_ != i_univ)
      continue;
    // Check if this cell contains the particle
    if (model::cells[i_cell]->contains(r, u, surf)) {
      p.lowest_coord().cell() = i_cell;
      return true;
    }
  }
  return false;
}
bool Universe::find_cell_in_virtual_lattice(GeometryState& p) const
{
  Cell& c {*model::cells[model::cell_map[filled_with_triso_base_]]};
  vector<int> lat_ind(3);
  Position r {p.r_local()};
  lat_ind[0] = std::max<int>(
    std::min<int>(
      floor((r.x - c.vl_lower_left_[0]) / c.vl_pitch_[0]), c.vl_shape_[0] - 1),
    0);
  lat_ind[1] = std::max<int>(
    std::min<int>(
      floor((r.y - c.vl_lower_left_[1]) / c.vl_pitch_[1]), c.vl_shape_[1] - 1),
    0);
  lat_ind[2] = std::max<int>(
    std::min<int>(
      floor((r.z - c.vl_lower_left_[2]) / c.vl_pitch_[2]), c.vl_shape_[2] - 1),
    0);

  int32_t i_univ = p.lowest_coord().universe();
  for (int token :
    c.vl_triso_distribution_[lat_ind[0] + lat_ind[1] * c.vl_shape_[0] +
                             lat_ind[2] * c.vl_shape_[0] * c.vl_shape_[1]]) {
    vector<double> triso_center = model::surfaces[abs(token) - 1]->get_center();
    double triso_radius = model::surfaces[abs(token) - 1]->get_radius();
    if (model::cells
          [model::cell_map[model::surfaces[abs(token) - 1]->triso_base_index_]]
            ->universe_ != i_univ)
      continue;
    if (abs(token) == abs(p.surface())) {
      if (p.surface() < 0) {
        p.lowest_coord().cell() =
          model::cell_map[model::surfaces[abs(token) - 1]
                            ->triso_particle_index_];
        return true;
      } else {
        p.lowest_coord().cell() = model::cell_map[filled_with_triso_base_];
        return true;
      }
    }
    if (pow(r.x - triso_center[0], 2) + pow(r.y - triso_center[1], 2) +
          pow(r.z - triso_center[2], 2) <
        pow(triso_radius, 2)) {
      p.lowest_coord().cell() =
        model::cell_map[model::surfaces[abs(token) - 1]->triso_particle_index_];
      return true;
    }
  }
  if (model::cells[model::cell_map[filled_with_triso_base_]]->universe_ ==
      i_univ) {
    p.lowest_coord().cell() = model::cell_map[filled_with_triso_base_];
    return true;
  }
  return false;
}

BoundingBox Universe::bounding_box() const
{
  BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
  if (cells_.size() == 0) {
    return {};
  } else {
    for (const auto& cell : cells_) {
      auto& c = model::cells[cell];
      bbox |= c->bounding_box();
    }
  }
  return bbox;
}

//==============================================================================
// UniversePartitioner implementation
//==============================================================================

UniversePartitioner::UniversePartitioner(const Universe& univ)
{
  // Define an ordered set of surface indices that point to z-planes.  Use a
  // functor to to order the set by the z0_ values of the corresponding planes.
  struct compare_surfs {
    bool operator()(const int32_t& i_surf, const int32_t& j_surf) const
    {
      const auto* surf = model::surfaces[i_surf].get();
      const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);
      double zi = zplane->z0_;
      surf = model::surfaces[j_surf].get();
      zplane = dynamic_cast<const SurfaceZPlane*>(surf);
      double zj = zplane->z0_;
      return zi < zj;
    }
  };
  std::set<int32_t, compare_surfs> surf_set;

  // Find all of the z-planes in this universe.  A set is used here for the
  // O(log(n)) insertions that will ensure entries are not repeated.
  for (auto i_cell : univ.cells_) {
    for (auto token : model::cells[i_cell]->surfaces()) {
      auto i_surf = std::abs(token) - 1;
      const auto* surf = model::surfaces[i_surf].get();
      if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf))
        surf_set.insert(i_surf);
    }
  }

  // Populate the surfs_ vector from the ordered set.
  surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end());

  // Populate the partition lists.
  partitions_.resize(surfs_.size() + 1);
  for (auto i_cell : univ.cells_) {
    // It is difficult to determine the bounds of a complex cell, so add complex
    // cells to all partitions.
    if (!model::cells[i_cell]->is_simple()) {
      for (auto& p : partitions_)
        p.push_back(i_cell);
      continue;
    }

    // Find the tokens for bounding z-planes.
    int32_t lower_token = 0, upper_token = 0;
    double min_z, max_z;
    for (auto token : model::cells[i_cell]->surfaces()) {
      const auto* surf = model::surfaces[std::abs(token) - 1].get();
      if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf)) {
        if (lower_token == 0 || zplane->z0_ < min_z) {
          lower_token = token;
          min_z = zplane->z0_;
        }
        if (upper_token == 0 || zplane->z0_ > max_z) {
          upper_token = token;
          max_z = zplane->z0_;
        }
      }
    }

    // If there are no bounding z-planes, add this cell to all partitions.
    if (lower_token == 0) {
      for (auto& p : partitions_)
        p.push_back(i_cell);
      continue;
    }

    // Find the first partition this cell lies in.  If the lower_token indicates
    // a negative halfspace, then the cell is unbounded in the lower direction
    // and it lies in the first partition onward.  Otherwise, it is bounded by
    // the positive halfspace given by the lower_token.
    int first_partition = 0;
    if (lower_token > 0) {
      for (int i = 0; i < surfs_.size(); ++i) {
        if (lower_token == surfs_[i] + 1) {
          first_partition = i + 1;
          break;
        }
      }
    }

    // Find the last partition this cell lies in.  The logic is analogous to the
    // logic for first_partition.
    int last_partition = surfs_.size();
    if (upper_token < 0) {
      for (int i = first_partition; i < surfs_.size(); ++i) {
        if (upper_token == -(surfs_[i] + 1)) {
          last_partition = i;
          break;
        }
      }
    }

    // Add the cell to all relevant partitions.
    for (int i = first_partition; i <= last_partition; ++i) {
      partitions_[i].push_back(i_cell);
    }
  }
}

const vector<int32_t>& UniversePartitioner::get_cells(
  Position r, Direction u) const
{
  // Perform a binary search for the partition containing the given coordinates.
  int left = 0;
  int middle = (surfs_.size() - 1) / 2;
  int right = surfs_.size() - 1;
  while (true) {
    // Check the sense of the coordinates for the current surface.
    const auto& surf = *model::surfaces[surfs_[middle]];
    if (surf.sense(r, u)) {
      // The coordinates lie in the positive halfspace.  Recurse if there are
      // more surfaces to check.  Otherwise, return the cells on the positive
      // side of this surface.
      int right_leaf = right - (right - middle) / 2;
      if (right_leaf != middle) {
        left = middle + 1;
        middle = right_leaf;
      } else {
        return partitions_[middle + 1];
      }

    } else {
      // The coordinates lie in the negative halfspace.  Recurse if there are
      // more surfaces to check.  Otherwise, return the cells on the negative
      // side of this surface.
      int left_leaf = left + (middle - left) / 2;
      if (left_leaf != middle) {
        right = middle - 1;
        middle = left_leaf;
      } else {
        return partitions_[middle];
      }
    }
  }
}

} // namespace openmc

1	#include "openmc/universe.h"
2
3	#include <set>
4
5	#include "openmc/hdf5_interface.h"
6	#include "openmc/particle.h"
7
8	namespace openmc {
9
10	namespace model {
11
12	std::unordered_map<int32_t, int32_t> universe_map;
13	vector<unique_ptr<Universe>> universes;
14
15	} // namespace model
16
17	//==============================================================================
18	// Universe implementation
19	//==============================================================================
20
21	void Universe::to_hdf5(hid_t universes_group) const	16,096✔
22	{
23	// Create a group for this universe.
24	auto group = create_group(universes_group, fmt::format("universe {}", id_));	32,192✔
25
26	// Write the geometry representation type.
27	write_string(group, "geom_type", "csg", false);	16,096✔
28
29	// Write the contained cells.
30	if (cells_.size() > 0) {	16,096!
31	vector<int32_t> cell_ids;	16,096✔
32	for (auto i_cell : cells_)	43,874✔
33	cell_ids.push_back(model::cells[i_cell]->id_);	27,778✔
34	write_dataset(group, "cells", cell_ids);	16,096✔
35	}	16,096✔
36
37	close_group(group);	16,096✔
38	}	16,096✔
39
40	bool Universe::find_cell(GeometryState& p) const	1,560,728,172✔
41	{
42	if (filled_with_triso_base_ != -1) {	1,560,728,172✔
43	bool found = find_cell_in_virtual_lattice(p);	3,180,298✔
44	if (found) {	3,180,298!
45	return found;	3,180,298✔
46	}
47	}
48	const auto& cells {
49	!partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};	1,557,547,874✔
50
51	Position r {p.r_local()};	1,557,547,874✔
52	Position u {p.u_local()};	1,557,547,874✔
53	auto surf = p.surface();	1,557,547,874✔
54	int32_t i_univ = p.lowest_coord().universe();	1,557,547,874✔
55
56	for (auto i_cell : cells) {	2,147,483,647✔
57	if (model::cells[i_cell]->universe_ != i_univ)	2,042,446,998!
58	continue;	×
59	// Check if this cell contains the particle
60	if (model::cells[i_cell]->contains(r, u, surf)) {	2,042,446,998✔
61	p.lowest_coord().cell() = i_cell;	1,384,416,507✔
62	return true;	1,384,416,507✔
63	}
64	}
65	return false;	173,131,367✔
66	}
67	bool Universe::find_cell_in_virtual_lattice(GeometryState& p) const	3,180,298✔
68	{
69	Cell& c {*model::cells[model::cell_map[filled_with_triso_base_]]};	3,180,298✔
70	vector<int> lat_ind(3);	3,180,298✔
71	Position r {p.r_local()};	3,180,298✔
72	lat_ind[0] = std::max<int>(	3,180,298✔
73	std::min<int>(
74	floor((r.x - c.vl_lower_left_[0]) / c.vl_pitch_[0]), c.vl_shape_[0] - 1),	3,180,298✔
75	0);	3,180,298✔
76	lat_ind[1] = std::max<int>(	3,180,298✔
77	std::min<int>(
78	floor((r.y - c.vl_lower_left_[1]) / c.vl_pitch_[1]), c.vl_shape_[1] - 1),	3,180,298✔
79	0);	3,180,298✔
80	lat_ind[2] = std::max<int>(	3,180,298✔
81	std::min<int>(
82	floor((r.z - c.vl_lower_left_[2]) / c.vl_pitch_[2]), c.vl_shape_[2] - 1),	3,180,298✔
83	0);	3,180,298✔
84
85	int32_t i_univ = p.lowest_coord().universe();	3,180,298✔
86	for (int token :	3,180,298✔
87	c.vl_triso_distribution_[lat_ind[0] + lat_ind[1] * c.vl_shape_[0] +	3,180,298✔
88	lat_ind[2] * c.vl_shape_[0] * c.vl_shape_[1]]) {	27,550,545✔
89	vector<double> triso_center = model::surfaces[abs(token) - 1]->get_center();	18,012,951✔
90	double triso_radius = model::surfaces[abs(token) - 1]->get_radius();	18,012,951✔
91	if (model::cells	18,012,951✔
92	[model::cell_map[model::surfaces[abs(token) - 1]->triso_base_index_]]	18,012,951✔
93	->universe_ != i_univ)	18,012,951!
NEW 94	continue;	×
95	if (abs(token) == abs(p.surface())) {	18,012,951!
NEW 96	if (p.surface() < 0) {	×
NEW 97	p.lowest_coord().cell() =	×
NEW 98	model::cell_map[model::surfaces[abs(token) - 1]	×
NEW 99	->triso_particle_index_];	×
NEW 100	return true;	×
101	} else {
NEW 102	p.lowest_coord().cell() = model::cell_map[filled_with_triso_base_];	×
NEW 103	return true;	×
104	}
105	}
106	if (pow(r.x - triso_center[0], 2) + pow(r.y - triso_center[1], 2) +	18,012,951✔
107	pow(r.z - triso_center[2], 2) <	18,012,951✔
108	pow(triso_radius, 2)) {	18,012,951✔
109	p.lowest_coord().cell() =	3,300✔
110	model::cell_map[model::surfaces[abs(token) - 1]->triso_particle_index_];	3,300✔
111	return true;	3,300✔
112	}
113	}	18,012,951!
114	if (model::cells[model::cell_map[filled_with_triso_base_]]->universe_ ==	3,176,998!
115	i_univ) {
116	p.lowest_coord().cell() = model::cell_map[filled_with_triso_base_];	3,176,998✔
117	return true;	3,176,998✔
118	}
NEW 119	return false;	×
120	}	3,180,298✔
121
122	BoundingBox Universe::bounding_box() const	11✔
123	{
124	BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};	11✔
125	if (cells_.size() == 0) {	11!
126	return {};	×
127	} else {
128	for (const auto& cell : cells_) {	44✔
129	auto& c = model::cells[cell];	33✔
130	bbox \|= c->bounding_box();	33✔
131	}
132	}
133	return bbox;	11✔
134	}
135
136	//==============================================================================
137	// UniversePartitioner implementation
138	//==============================================================================
139
140	UniversePartitioner::UniversePartitioner(const Universe& univ)	96✔
141	{
142	// Define an ordered set of surface indices that point to z-planes. Use a
143	// functor to to order the set by the z0_ values of the corresponding planes.
144	struct compare_surfs {
145	bool operator()(const int32_t& i_surf, const int32_t& j_surf) const	9,872✔
146	{
147	const auto* surf = model::surfaces[i_surf].get();	9,872✔
148	const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);	9,872!
149	double zi = zplane->z0_;	9,872✔
150	surf = model::surfaces[j_surf].get();	9,872✔
151	zplane = dynamic_cast<const SurfaceZPlane*>(surf);	9,872!
152	double zj = zplane->z0_;	9,872✔
153	return zi < zj;	9,872✔
154	}
155	};
156	std::set<int32_t, compare_surfs> surf_set;	96✔
157
158	// Find all of the z-planes in this universe. A set is used here for the
159	// O(log(n)) insertions that will ensure entries are not repeated.
160	for (auto i_cell : univ.cells_) {	1,344✔
161	for (auto token : model::cells[i_cell]->surfaces()) {	5,520✔
162	auto i_surf = std::abs(token) - 1;	4,272✔
163	const auto* surf = model::surfaces[i_surf].get();	4,272✔
164	if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf))	4,272!
165	surf_set.insert(i_surf);	2,432✔
166	}	1,248✔
167	}
168
169	// Populate the surfs_ vector from the ordered set.
170	surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end());	96✔
171
172	// Populate the partition lists.
173	partitions_.resize(surfs_.size() + 1);	96✔
174	for (auto i_cell : univ.cells_) {	1,344✔
175	// It is difficult to determine the bounds of a complex cell, so add complex
176	// cells to all partitions.
177	if (!model::cells[i_cell]->is_simple()) {	1,248!
178	for (auto& p : partitions_)	×
179	p.push_back(i_cell);	×
180	continue;	×
181	}	×
182
183	// Find the tokens for bounding z-planes.
184	int32_t lower_token = 0, upper_token = 0;	1,248✔
185	double min_z, max_z;
186	for (auto token : model::cells[i_cell]->surfaces()) {	5,520✔
187	const auto* surf = model::surfaces[std::abs(token) - 1].get();	4,272✔
188	if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf)) {	4,272!
189	if (lower_token == 0 \|\| zplane->z0_ < min_z) {	2,432!
190	lower_token = token;	1,248✔
191	min_z = zplane->z0_;	1,248✔
192	}
193	if (upper_token == 0 \|\| zplane->z0_ > max_z) {	2,432!
194	upper_token = token;	2,432✔
195	max_z = zplane->z0_;	2,432✔
196	}
197	}
198	}	1,248✔
199
200	// If there are no bounding z-planes, add this cell to all partitions.
201	if (lower_token == 0) {	1,248!
202	for (auto& p : partitions_)	×
203	p.push_back(i_cell);	×
204	continue;	×
205	}	×
206
207	// Find the first partition this cell lies in. If the lower_token indicates
208	// a negative halfspace, then the cell is unbounded in the lower direction
209	// and it lies in the first partition onward. Otherwise, it is bounded by
210	// the positive halfspace given by the lower_token.
211	int first_partition = 0;	1,248✔
212	if (lower_token > 0) {	1,248✔
213	for (int i = 0; i < surfs_.size(); ++i) {	4,624!
214	if (lower_token == surfs_[i] + 1) {	4,624✔
215	first_partition = i + 1;	1,216✔
216	break;	1,216✔
217	}
218	}
219	}
220
221	// Find the last partition this cell lies in. The logic is analogous to the
222	// logic for first_partition.
223	int last_partition = surfs_.size();	1,248✔
224	if (upper_token < 0) {	1,248✔
225	for (int i = first_partition; i < surfs_.size(); ++i) {	2,624!
226	if (upper_token == -(surfs_[i] + 1)) {	2,624✔
227	last_partition = i;	1,216✔
228	break;	1,216✔
229	}
230	}
231	}
232
233	// Add the cell to all relevant partitions.
234	for (int i = first_partition; i <= last_partition; ++i) {	3,920✔
235	partitions_[i].push_back(i_cell);	2,672✔
236	}
237	}
238	}	96✔
239
240	const vector<int32_t>& UniversePartitioner::get_cells(	247,447✔
241	Position r, Direction u) const
242	{
243	// Perform a binary search for the partition containing the given coordinates.
244	int left = 0;	247,447✔
245	int middle = (surfs_.size() - 1) / 2;	247,447✔
246	int right = surfs_.size() - 1;	247,447✔
247	while (true) {
248	// Check the sense of the coordinates for the current surface.
249	const auto& surf = *model::surfaces[surfs_[middle]];	742,368✔
250	if (surf.sense(r, u)) {	742,368✔
251	// The coordinates lie in the positive halfspace. Recurse if there are
252	// more surfaces to check. Otherwise, return the cells on the positive
253	// side of this surface.
254	int right_leaf = right - (right - middle) / 2;	286,219✔
255	if (right_leaf != middle) {	286,219✔
256	left = middle + 1;	248,266✔
257	middle = right_leaf;	248,266✔
258	} else {
259	return partitions_[middle + 1];	37,953✔
260	}
261
262	} else {
263	// The coordinates lie in the negative halfspace. Recurse if there are
264	// more surfaces to check. Otherwise, return the cells on the negative
265	// side of this surface.
266	int left_leaf = left + (middle - left) / 2;	456,149✔
267	if (left_leaf != middle) {	456,149✔
268	right = middle - 1;	246,655✔
269	middle = left_leaf;	246,655✔
270	} else {
271	return partitions_[middle];	209,494✔
272	}
273	}
274	}	494,921✔
275	}
276
277	} // namespace openmc

openmc-dev / openmc / 20368472998

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