34703966194

Committed 12 Sep 2026 04:01PM UTC coverage: 81.582% (+0.05%) from 81.533%

Build # 34703966194

Build Type

Pull #3757

github

Committed by

web-flow

Commit Message

Merge c7f3903d7 into 0888e2ff3

Pull Request Pull Request #3757: Implementation of point detectors

Coverage Stats

19334 of 27917 branches covered (69.26%)

Branch coverage included in aggregate %.

487 of 575 new or added lines in 24 files covered. (84.7%)

197 existing lines in 6 files now uncovered.

61598 of 71286 relevant lines covered (86.41%)

49750709.47 hits per line

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

66.83

/src/ray.cpp

#include "openmc/ray.h"

#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/material.h"
#include "openmc/mgxs_interface.h"
#include "openmc/settings.h"

namespace openmc {

namespace {
// Max intersections before we assume ray tracing is caught in an infinite loop
constexpr int MAX_INTERSECTIONS = 1000000;
} // namespace

template<typename RayT>
void trace_ray(RayT& ray, double max_distance)
{
  // Clear anything left over from a previous trace on this object.
  ray.reset_trace_state();

  // To trace the ray from its origin all the way through the model, we have
  // to proceed in two phases. In the first, the ray may or may not be found
  // inside the model. If the ray is already in the model, phase one can be
  // skipped. Otherwise, the ray has to be advanced to the boundary of the
  // model where all the cells are defined. Importantly, this is assuming that
  // the model is convex, which is a very reasonable assumption for any
  // radiation transport model.
  //
  // After phase one is done, we can starting tracing from cell to cell within
  // the model. This step can use neighbor lists to accelerate the ray tracing.

  // Remaining distance budget
  double max = max_distance;

  bool inside_cell;
  // Check for location if the particle is already known
  if (ray.lowest_coord().cell() == C_NONE) {
    // The geometry position of the particle is either unknown or outside of the
    // edge of the model.
    if (ray.lowest_coord().universe() == C_NONE) {
      // Attempt to initialize the particle. We may have to
      // enter a loop to move it up to the edge of the model.
      inside_cell = exhaustive_find_cell(ray, settings::verbosity >= 10);
    } else {
      // It has been already calculated that the current position is outside of
      // the edge of the model.
      inside_cell = false;
    }
  } else {
    // Availability of the cell means that the particle is located inside the
    // edge.
    inside_cell = true;
  }

  // Advance to the boundary of the model
  while (!inside_cell) {
    ray.advance_to_boundary_from_void();

    // Flight through void is real flight: it has to be charged against the
    // distance budget, otherwise max_distance ends up being measured from the
    // point where the ray entered the model rather than from its origin. It
    // is accumulated into total_distance_ only -- traversal_distance_ stays
    // zero until the model is reached.
    if (ray.boundary().surface() != SURFACE_NONE &&
        ray.boundary().distance() < INFTY) {
      // advance_to_boundary_from_void() has already moved the ray to the
      // boundary plus the TINY_BIT of padding, so that is the distance to
      // account for.
      double advance = ray.boundary().distance() + TINY_BIT;

      if (advance >= max) {
        // The budget runs out before the ray even reaches the model. Back it
        // up so the net movement is exactly max.
        ray.move_distance(max - advance);
        ray.update_distance(max);
        ray.completed_ = true;
        return;
      }
      ray.update_distance(advance);
      max -= advance;
    }

    inside_cell = exhaustive_find_cell(ray, settings::verbosity >= 10);

    // If true this means no surface was intersected. See cell.cpp and search
    // for numeric_limits to see where we return it.
    if (ray.surface() == std::numeric_limits<int>::max()) {
      warning(fmt::format("Lost a ray, r = {}, u = {}", ray.r(), ray.u()));
      return;
    }

    // Exit this loop and enter into cell-to-cell ray tracing (which uses
    // neighbor lists)
    if (inside_cell)
      break;

    // if there is no intersection with the model, we're done
    if (ray.boundary().surface() == SURFACE_NONE)
      return;

    ray.event_counter_++;
    if (ray.event_counter_ > MAX_INTERSECTIONS) {
      warning("Likely infinite loop while ray tracing");
      return;
    }
  }

  // From here on the ray is inside the model, so its flight counts toward
  // traversal_distance_ as well as total_distance_.
  ray.in_model_ = true;

  // Call the specialized logic for this type of ray. This is for the
  // intersection for the first intersection if we had one.
  if (ray.boundary().surface() != SURFACE_NONE) {
    // set the geometry state's surface attribute to be used for
    // surface normal computation
    ray.surface() = ray.boundary().surface();
    ray.on_intersection();
    if (ray.stop_)
      return;
  }

  // reset surface attribute to zero after the first intersection so that it
  // doesn't perturb surface crossing logic from here on out
  ray.surface() = 0;

  // This is the ray tracing loop within the model. It exits after exiting
  // the model, which is equivalent to assuming that the model is convex.
  // It would be nice to factor out the on_intersection at the end of this
  // loop and then do "while (inside_cell)", but we can't guarantee it's
  // on a surface in that case. There might be some other way to set it
  // up that is perhaps a little more elegant, but this is what works just
  // fine.
  while (true) {

    ray.boundary() = distance_to_boundary(ray);

    // There are no more intersections to process
    // if we hit the edge of the model, so stop
    // the particle in that case. Also, just exit
    // if a negative distance was somehow computed.
    if (ray.boundary().distance() == INFTY ||
        ray.boundary().distance() == INFINITY ||
        ray.boundary().distance() < 0) {
      return;
    }

    // Distance from the ray's current position to the next surface.
    const double surface_distance = ray.boundary().distance();

    // See below comment where call_on_intersection is checked in an
    // if statement for an explanation of this.
    bool call_on_intersection {surface_distance >= 10 * TINY_BIT};

    // DAGMC surfaces expect us to go a little bit further than the advance
    // distance to properly check cell inclusion.
    double advance = surface_distance + TINY_BIT;

    if (advance >= max) {
      // The ray runs out of budget inside this cell, so no surface is
      // crossed. Only the truncated distance was actually travelled.
      ray.move_distance(max);
      ray.update_distance(max);
      ray.completed_ = true;
      return;
    }

    ray.move_distance(advance);

    max -= advance;

    ray.surface() = ray.boundary().surface();
    // Initialize last cells from the current cell, because the cell() variable
    // does not contain the data for the case of a single-segment ray
    for (int j = 0; j < ray.n_coord(); ++j) {
      ray.cell_last(j) = ray.coord(j).cell();
    }
    ray.n_coord_last() = ray.n_coord();
    ray.n_coord() = ray.boundary().coord_level();
    if (ray.boundary().lattice_translation()[0] != 0 ||
        ray.boundary().lattice_translation()[1] != 0 ||
        ray.boundary().lattice_translation()[2] != 0) {
      cross_lattice(ray, ray.boundary(), settings::verbosity >= 10);
    }

    // Accumulate before the cell search, while material() still refers to the
    // cell the ray just crossed.
    ray.update_distance(advance);

    inside_cell = neighbor_list_find_cell(ray, settings::verbosity >= 10);

    // Call the specialized logic for this type of ray. Note that we do not
    // call this if the advance distance is very small. Unfortunately, it seems
    // darn near impossible to get the particle advanced to the model boundary
    // and through it without sometimes accidentally calling on_intersection
    // twice. This incorrectly shades the region as occluded when it might not
    // actually be. By screening out intersection distances smaller than a
    // threshold 10x larger than the scoot distance used to advance up to the
    // model boundary, we can avoid that situation.
    if (call_on_intersection) {
      ray.on_intersection();
      if (ray.stop_)
        return;
    }

    if (!inside_cell)
      return;

    ray.event_counter_++;
    if (ray.event_counter_ > MAX_INTERSECTIONS) {
      warning("Likely infinite loop while ray tracing");
      return;
    }
  }
}

void RayState::accumulate_distance(double distance)
{
  total_distance_ += distance;
  if (in_model_) {
    traversal_distance_ += distance;
  }
}

void Ray::trace(double max_distance)
{
  trace_ray(*this, max_distance);
}

void ParticleRay::trace(double max_distance)
{
  trace_ray(*this, max_distance);
}

void ParticleRay::init_physics(ParticleType type_, double time_, double E_)
{
  type() = type_;
  time() = time_;
  time_start_ = time_;

  E() = E_;
  E_last() = E_;

  // In multigroup mode the physics is driven by the group index rather than
  // the energy. Particle::from_source() takes the group straight from the
  // source site; here it has to be derived from the energy, otherwise every
  // ray silently uses group 0.
  if (!settings::run_CE) {
    g() = data::mg.get_group_index(E_);
    g_last() = g();
    E() = data::mg.energy_bin_avg_[g()];
    E_last() = E();
  }

  // MacroXS has no default member initializers, and the void branch of
  // update_distance() only writes four of its fields.
  macro_xs() = {};
}

void ParticleRay::mark_as_lost(const char* message)
{
  if (settings::verbosity >= 10) {
    warning(message);
  }
  stop();
}

void ParticleRay::update_distance(double distance)
{
  accumulate_distance(distance);

  time() += distance / speed();

  // Calculate microscopic and macroscopic cross sections
  if (material() != MATERIAL_VOID) {
    if (settings::run_CE) {
      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
          density_mult() != density_mult_last()) {
        // If the material is the same as the last material and the
        // temperature hasn't changed, we don't need to lookup cross
        // sections again.
        model::materials[material()]->calculate_xs(*this);
      }
    } else {
      // Get the MG data; unlike the CE case above, we have to re-calculate
      // cross sections for every collision since the cross sections may
      // be angle-dependent
      data::mg.macro_xs_[material()].calculate_xs(*this);

      // Update the particle's group while we know we are multi-group
      g_last() = g();
    }
  } else {
    macro_xs() = {};
  }

  traversal_mfp_ += macro_xs().total * distance;
}

// Explicit instantiations: the two kinds of ray that share the tracing loop
template void trace_ray<Ray>(Ray&, double);
template void trace_ray<ParticleRay>(ParticleRay&, double);

} // namespace openmc

1	#include "openmc/ray.h"
2
3	#include "openmc/error.h"
4	#include "openmc/geometry.h"
5	#include "openmc/material.h"
6	#include "openmc/mgxs_interface.h"
7	#include "openmc/settings.h"
8
9	namespace openmc {
10
11	namespace {
12	// Max intersections before we assume ray tracing is caught in an infinite loop
13	constexpr int MAX_INTERSECTIONS = 1000000;
14	} // namespace
15
16	template<typename RayT>
17	void trace_ray(RayT& ray, double max_distance)	13,112,011 !
18	{
19	// Clear anything left over from a previous trace on this object.
20	ray.reset_trace_state();	13,112,011 ✔
21
22	// To trace the ray from its origin all the way through the model, we have
23	// to proceed in two phases. In the first, the ray may or may not be found
24	// inside the model. If the ray is already in the model, phase one can be
25	// skipped. Otherwise, the ray has to be advanced to the boundary of the
26	// model where all the cells are defined. Importantly, this is assuming that
27	// the model is convex, which is a very reasonable assumption for any
28	// radiation transport model.
29	//
30	// After phase one is done, we can starting tracing from cell to cell within
31	// the model. This step can use neighbor lists to accelerate the ray tracing.
32
33	// Remaining distance budget
34	double max = max_distance;	13,112,011 ✔
35
36	bool inside_cell;
37	// Check for location if the particle is already known
38	if (ray.lowest_coord().cell() == C_NONE) {	13,112,011 !
39	// The geometry position of the particle is either unknown or outside of the
40	// edge of the model.
41	if (ray.lowest_coord().universe() == C_NONE) {	13,112,011 !
42	// Attempt to initialize the particle. We may have to
43	// enter a loop to move it up to the edge of the model.
44	inside_cell = exhaustive_find_cell(ray, settings::verbosity >= 10);	13,112,011 ✔
45	} else {
46	// It has been already calculated that the current position is outside of
47	// the edge of the model.
48	inside_cell = false;
49	}
50	} else {
51	// Availability of the cell means that the particle is located inside the
52	// edge.
53	inside_cell = true;
54	}
55
56	// Advance to the boundary of the model
57	while (!inside_cell) {	25,209,393 ✔
58	ray.advance_to_boundary_from_void();	15,618,460 ✔
59
60	// Flight through void is real flight: it has to be charged against the
61	// distance budget, otherwise max_distance ends up being measured from the
62	// point where the ray entered the model rather than from its origin. It
63	// is accumulated into total_distance_ only -- traversal_distance_ stays
64	// zero until the model is reached.
65	if (ray.boundary().surface() != SURFACE_NONE &&	15,618,460 !
66	ray.boundary().distance() < INFTY) {	13,647,238 !
67	// advance_to_boundary_from_void() has already moved the ray to the
68	// boundary plus the TINY_BIT of padding, so that is the distance to
69	// account for.
70	double advance = ray.boundary().distance() + TINY_BIT;	13,647,238 !
71
72	if (advance >= max) {	13,647,238 !
73	// The budget runs out before the ray even reaches the model. Back it
74	// up so the net movement is exactly max.
NEW 75	ray.move_distance(max - advance);	×
NEW 76	ray.update_distance(max);	×
NEW 77	ray.completed_ = true;	×
NEW 78	return;	×
79	}
80	ray.update_distance(advance);	13,647,238 ✔
81	max -= advance;	13,647,238 ✔
82	}
83
84	inside_cell = exhaustive_find_cell(ray, settings::verbosity >= 10);	15,618,460 ✔
85
86	// If true this means no surface was intersected. See cell.cpp and search
87	// for numeric_limits to see where we return it.
88	if (ray.surface() == std::numeric_limits<int>::max()) {	15,618,460 !
NEW 89	warning(fmt::format("Lost a ray, r = {}, u = {}", ray.r(), ray.u()));	×
UNCOV 90	return;	×
91	}
92
93	// Exit this loop and enter into cell-to-cell ray tracing (which uses
94	// neighbor lists)
95	if (inside_cell)	15,618,460 ✔
96	break;
97
98	// if there is no intersection with the model, we're done
99	if (ray.boundary().surface() == SURFACE_NONE)	14,068,604 ✔
100	return;
101
102	ray.event_counter_++;	12,097,382 ✔
103	if (ray.event_counter_ > MAX_INTERSECTIONS) {	12,097,382 !
NEW 104	warning("Likely infinite loop while ray tracing");	×
UNCOV 105	return;	×
106	}
107	}
108
109	// From here on the ray is inside the model, so its flight counts toward
110	// traversal_distance_ as well as total_distance_.
111	ray.in_model_ = true;	11,140,789 ✔
112
113	// Call the specialized logic for this type of ray. This is for the
114	// intersection for the first intersection if we had one.
115	if (ray.boundary().surface() != SURFACE_NONE) {	11,140,789 ✔
116	// set the geometry state's surface attribute to be used for
117	// surface normal computation
118	ray.surface() = ray.boundary().surface();	1,549,856 !
119	ray.on_intersection();	1,549,856 ✔
120	if (ray.stop_)	1,549,856 !
121	return;
122	}
123
124	// reset surface attribute to zero after the first intersection so that it
125	// doesn't perturb surface crossing logic from here on out
126	ray.surface() = 0;	11,140,789 ✔
127
128	// This is the ray tracing loop within the model. It exits after exiting
129	// the model, which is equivalent to assuming that the model is convex.
130	// It would be nice to factor out the on_intersection at the end of this
131	// loop and then do "while (inside_cell)", but we can't guarantee it's
132	// on a surface in that case. There might be some other way to set it
133	// up that is perhaps a little more elegant, but this is what works just
134	// fine.
135	while (true) {
136
137	ray.boundary() = distance_to_boundary(ray);	16,961,901 ✔
138
139	// There are no more intersections to process
140	// if we hit the edge of the model, so stop
141	// the particle in that case. Also, just exit
142	// if a negative distance was somehow computed.
143	if (ray.boundary().distance() == INFTY \|\|	16,961,901 !
144	ray.boundary().distance() == INFINITY \|\|	16,961,901 !
145	ray.boundary().distance() < 0) {	16,961,901 !
146	return;
147	}
148
149	// Distance from the ray's current position to the next surface.
150	const double surface_distance = ray.boundary().distance();	16,961,901 ✔
151
152	// See below comment where call_on_intersection is checked in an
153	// if statement for an explanation of this.
154	bool call_on_intersection {surface_distance >= 10 * TINY_BIT};	16,961,901 ✔
155
156	// DAGMC surfaces expect us to go a little bit further than the advance
157	// distance to properly check cell inclusion.
158	double advance = surface_distance + TINY_BIT;	16,961,901 ✔
159
160	if (advance >= max) {	16,961,901 ✔
161	// The ray runs out of budget inside this cell, so no surface is
162	// crossed. Only the truncated distance was actually travelled.
163	ray.move_distance(max);	9,590,933 ✔
164	ray.update_distance(max);	9,590,933 ✔
165	ray.completed_ = true;	9,590,933 ✔
166	return;	9,590,933 ✔
167	}
168
169	ray.move_distance(advance);	7,370,968 ✔
170
171	max -= advance;	7,370,968 ✔
172
173	ray.surface() = ray.boundary().surface();	7,370,968 ✔
174	// Initialize last cells from the current cell, because the cell() variable
175	// does not contain the data for the case of a single-segment ray
176	for (int j = 0; j < ray.n_coord(); ++j) {	14,741,936 ✔
177	ray.cell_last(j) = ray.coord(j).cell();	7,370,968 ✔
178	}
179	ray.n_coord_last() = ray.n_coord();	7,370,968 ✔
180	ray.n_coord() = ray.boundary().coord_level();	7,370,968 !
181	if (ray.boundary().lattice_translation()[0] != 0 \|\|	7,370,968 !
182	ray.boundary().lattice_translation()[1] != 0 \|\|	7,370,968 !
183	ray.boundary().lattice_translation()[2] != 0) {	7,370,968 !
NEW 184	cross_lattice(ray, ray.boundary(), settings::verbosity >= 10);	×
185	}
186
187	// Accumulate before the cell search, while material() still refers to the
188	// cell the ray just crossed.
189	ray.update_distance(advance);	7,370,968 ✔
190
191	inside_cell = neighbor_list_find_cell(ray, settings::verbosity >= 10);	7,370,968 ✔
192
193	// Call the specialized logic for this type of ray. Note that we do not
194	// call this if the advance distance is very small. Unfortunately, it seems
195	// darn near impossible to get the particle advanced to the model boundary
196	// and through it without sometimes accidentally calling on_intersection
197	// twice. This incorrectly shades the region as occluded when it might not
198	// actually be. By screening out intersection distances smaller than a
199	// threshold 10x larger than the scoot distance used to advance up to the
200	// model boundary, we can avoid that situation.
201	if (call_on_intersection) {	7,370,968 ✔
202	ray.on_intersection();	1,715,329 ✔
203	if (ray.stop_)	6,777,683 ✔
204	return;
205	}
206
207	if (!inside_cell)	7,335,405 ✔
208	return;
209
210	ray.event_counter_++;	5,821,112 ✔
211	if (ray.event_counter_ > MAX_INTERSECTIONS) {	5,821,112 !
NEW 212	warning("Likely infinite loop while ray tracing");	×
UNCOV 213	return;	×
214	}
215	}
216	}
217
218	void RayState::accumulate_distance(double distance)	30,609,139 ✔
219	{
220	total_distance_ += distance;	30,609,139 ✔
221	if (in_model_) {	30,609,139 ✔
222	traversal_distance_ += distance;	16,961,901 ✔
223	}
224	}	30,609,139 ✔
225
226	void Ray::trace(double max_distance)	3,521,078 ✔
227	{
228	trace_ray(*this, max_distance);	3,521,078 ✔
229	}	3,521,078 ✔
230
231	void ParticleRay::trace(double max_distance)	9,590,933 ✔
232	{
233	trace_ray(*this, max_distance);	9,590,933 ✔
234	}	9,590,933 ✔
235
236	void ParticleRay::init_physics(ParticleType type_, double time_, double E_)	9,590,933 ✔
237	{
238	type() = type_;	9,590,933 ✔
239	time() = time_;	9,590,933 ✔
240	time_start_ = time_;	9,590,933 ✔
241
242	E() = E_;	9,590,933 ✔
243	E_last() = E_;	9,590,933 ✔
244
245	// In multigroup mode the physics is driven by the group index rather than
246	// the energy. Particle::from_source() takes the group straight from the
247	// source site; here it has to be derived from the energy, otherwise every
248	// ray silently uses group 0.
249	if (!settings::run_CE) {	9,590,933 !
NEW 250	g() = data::mg.get_group_index(E_);	×
NEW 251	g_last() = g();	×
NEW 252	E() = data::mg.energy_bin_avg_[g()];	×
NEW 253	E_last() = E();	×
254	}
255
256	// MacroXS has no default member initializers, and the void branch of
257	// update_distance() only writes four of its fields.
258	macro_xs() = {};	9,590,933 ✔
259	}	9,590,933 ✔
260
NEW 261	void ParticleRay::mark_as_lost(const char* message)	×
262	{
NEW 263	if (settings::verbosity >= 10) {	×
NEW 264	warning(message);	×
265	}
NEW 266	stop();	×
NEW 267	}	×
268
269	void ParticleRay::update_distance(double distance)	14,653,287 ✔
270	{
271	accumulate_distance(distance);	14,653,287 ✔
272
273	time() += distance / speed();	14,653,287 ✔
274
275	// Calculate microscopic and macroscopic cross sections
276	if (material() != MATERIAL_VOID) {	14,653,287 !
277	if (settings::run_CE) {	14,653,287 !
278	if (material() != material_last() \|\| sqrtkT() != sqrtkT_last() \|\|	14,653,287 !
NEW 279	density_mult() != density_mult_last()) {	×
280	// If the material is the same as the last material and the
281	// temperature hasn't changed, we don't need to lookup cross
282	// sections again.
283	model::materials[material()]->calculate_xs(*this);	14,653,287 ✔
284	}
285	} else {
286	// Get the MG data; unlike the CE case above, we have to re-calculate
287	// cross sections for every collision since the cross sections may
288	// be angle-dependent
NEW 289	data::mg.macro_xs_[material()].calculate_xs(*this);	×
290
291	// Update the particle's group while we know we are multi-group
NEW 292	g_last() = g();	×
293	}
294	} else {
NEW 295	macro_xs() = {};	×
296	}
297
298	traversal_mfp_ += macro_xs().total * distance;	14,653,287 ✔
299	}	14,653,287 ✔
300
301	// Explicit instantiations: the two kinds of ray that share the tracing loop
302	template void trace_ray<Ray>(Ray&, double);
303	template void trace_ray<ParticleRay>(ParticleRay&, double);
304
305	} // namespace openmc

openmc-dev / openmc / 34703966194

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