30387632100

Committed 28 Jul 2026 06:29PM UTC coverage: 81.066% (-0.4%) from 81.454%

Build # 30387632100

Build Type

Pull #3757

github

Committed by

web-flow

Commit Message

Merge f0f0b9157 into f3a8ba29f

Pull Request Pull Request #3757: Implementation of point detectors

Coverage Stats

18401 of 26807 branches covered (68.64%)

Branch coverage included in aggregate %.

52 of 384 new or added lines in 25 files covered. (13.54%)

3 existing lines in 2 files now uncovered.

60130 of 70066 relevant lines covered (85.82%)

48889380.27 hits per line

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

59.62

/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 {

void Ray::compute_distance()
{
  boundary() = distance_to_boundary(*this);
}

void Ray::trace(double max_distance)
{
  // 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.

  double max = max_distance;

  bool inside_cell;
  // Check for location if the particle is already known
  if (lowest_coord().cell() == C_NONE) {
    // The geometry position of the particle is either unknown or outside of the
    // edge of the model.
    if (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(*this, 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) {
    advance_to_boundary_from_void();
    inside_cell = exhaustive_find_cell(*this, 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 (surface() == std::numeric_limits<int>::max()) {
      warning(fmt::format("Lost a ray, r = {}, u = {}", r(), 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 (boundary().surface() == SURFACE_NONE)
      return;

    event_counter_++;
    if (event_counter_ > MAX_INTERSECTIONS) {
      warning("Likely infinite loop in ray traced plot");
      return;
    }
  }

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

  // reset surface attribute to zero after the first intersection so that it
  // doesn't perturb surface crossing logic from here on out
  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) {

    compute_distance();

    // 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 (boundary().distance() == INFTY || boundary().distance() == INFINITY ||
        boundary().distance() < 0) {
      return;
    }

    // See below comment where call_on_intersection is checked in an
    // if statement for an explanation of this.
    bool call_on_intersection {true};
    if (boundary().distance() < 10 * TINY_BIT) {
      call_on_intersection = false;
    }

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

    double distance = std::min(boundary().distance(), max);

    // Advance particle, prepare for next intersection
    for (int lev = 0; lev < n_coord(); ++lev) {
      coord(lev).r() += distance * coord(lev).u();
    }

    max -= distance;

    if (max == 0.0) {
      update_distance();
      break;
    }

    surface() = 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 < n_coord(); ++j) {
      cell_last(j) = coord(j).cell();
    }
    n_coord_last() = n_coord();
    n_coord() = boundary().coord_level();
    if (boundary().lattice_translation()[0] != 0 ||
        boundary().lattice_translation()[1] != 0 ||
        boundary().lattice_translation()[2] != 0) {
      cross_lattice(*this, boundary(), settings::verbosity >= 10);
    }

    update_distance();

    inside_cell = neighbor_list_find_cell(*this, 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) {
      on_intersection();
      if (stop_)
        return;
    }

    if (!inside_cell)
      return;

    event_counter_++;
    if (event_counter_ > MAX_INTERSECTIONS) {
      warning("Likely infinite loop in ray traced plot");
      return;
    }
  }
}

void Ray::update_distance()
{
  // Record how far the ray has traveled
  traversal_distance_ += boundary().distance();
}

void ParticleRay::on_intersection() {}

void ParticleRay::update_distance()
{
  Ray::update_distance();

  time() += boundary().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().total = 0.0;
    macro_xs().absorption = 0.0;
    macro_xs().fission = 0.0;
    macro_xs().nu_fission = 0.0;
  }

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

} // 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	void Ray::compute_distance()	3,014,682✔
12	{
13	boundary() = distance_to_boundary(*this);	3,014,682✔
14	}	3,014,682✔
15
16	void Ray::trace(double max_distance)	3,521,078✔
17	{
18	// To trace the ray from its origin all the way through the model, we have
19	// to proceed in two phases. In the first, the ray may or may not be found
20	// inside the model. If the ray is already in the model, phase one can be
21	// skipped. Otherwise, the ray has to be advanced to the boundary of the
22	// model where all the cells are defined. Importantly, this is assuming that
23	// the model is convex, which is a very reasonable assumption for any
24	// radiation transport model.
25	//
26	// After phase one is done, we can starting tracing from cell to cell within
27	// the model. This step can use neighbor lists to accelerate the ray tracing.
28
29	double max = max_distance;	3,521,078✔
30
31	bool inside_cell;	3,521,078✔
32	// Check for location if the particle is already known
33	if (lowest_coord().cell() == C_NONE) {	3,521,078!
34	// The geometry position of the particle is either unknown or outside of the
35	// edge of the model.
36	if (lowest_coord().universe() == C_NONE) {	3,521,078!
37	// Attempt to initialize the particle. We may have to
38	// enter a loop to move it up to the edge of the model.
39	inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);	3,521,078✔
40	} else {
41	// It has been already calculated that the current position is outside of
42	// the edge of the model.
43	inside_cell = false;
44	}
45	} else {
46	// Availability of the cell means that the particle is located inside the
47	// edge.
48	inside_cell = true;
49	}
50
51	// Advance to the boundary of the model
52	while (!inside_cell) {	15,618,460!
53	advance_to_boundary_from_void();	15,618,460✔
54	inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);	15,618,460✔
55
56	// If true this means no surface was intersected. See cell.cpp and search
57	// for numeric_limits to see where we return it.
58	if (surface() == std::numeric_limits<int>::max()) {	15,618,460!
59	warning(fmt::format("Lost a ray, r = {}, u = {}", r(), u()));	×
60	return;	3,521,078✔
61	}
62
63	// Exit this loop and enter into cell-to-cell ray tracing (which uses
64	// neighbor lists)
65	if (inside_cell)	15,618,460✔
66	break;
67
68	// if there is no intersection with the model, we're done
69	if (boundary().surface() == SURFACE_NONE)	14,068,604✔
70	return;
71
72	event_counter_++;	12,097,382✔
73	if (event_counter_ > MAX_INTERSECTIONS) {	12,097,382!
74	warning("Likely infinite loop in ray traced plot");	×
75	return;	×
76	}
77	}
78
79	// Call the specialized logic for this type of ray. This is for the
80	// intersection for the first intersection if we had one.
81	if (boundary().surface() != SURFACE_NONE) {	1,549,856!
82	// set the geometry state's surface attribute to be used for
83	// surface normal computation
84	surface() = boundary().surface();	1,549,856✔
85	on_intersection();	1,549,856✔
86	if (stop_)	1,549,856!
87	return;
88	}
89
90	// reset surface attribute to zero after the first intersection so that it
91	// doesn't perturb surface crossing logic from here on out
92	surface() = 0;	1,549,856✔
93
94	// This is the ray tracing loop within the model. It exits after exiting
95	// the model, which is equivalent to assuming that the model is convex.
96	// It would be nice to factor out the on_intersection at the end of this
97	// loop and then do "while (inside_cell)", but we can't guarantee it's
98	// on a surface in that case. There might be some other way to set it
99	// up that is perhaps a little more elegant, but this is what works just
100	// fine.
101	while (true) {	2,308,614✔
102
103	compute_distance();	2,308,614✔
104
105	// There are no more intersections to process
106	// if we hit the edge of the model, so stop
107	// the particle in that case. Also, just exit
108	// if a negative distance was somehow computed.
109	if (boundary().distance() == INFTY \|\| boundary().distance() == INFINITY \|\|	2,308,614!
110	boundary().distance() < 0) {	2,308,614!
111	return;
112	}
113
114	// See below comment where call_on_intersection is checked in an
115	// if statement for an explanation of this.
116	bool call_on_intersection {true};	2,308,614✔
117	if (boundary().distance() < 10 * TINY_BIT) {	2,308,614✔
118	call_on_intersection = false;	593,285✔
119	}
120
121	// DAGMC surfaces expect us to go a little bit further than the advance
122	// distance to properly check cell inclusion.
123	boundary().distance() += TINY_BIT;	2,308,614!
124
125	double distance = std::min(boundary().distance(), max);	2,308,614!
126
127	// Advance particle, prepare for next intersection
128	for (int lev = 0; lev < n_coord(); ++lev) {	4,617,228✔
129	coord(lev).r() += distance * coord(lev).u();	2,308,614✔
130	}
131
132	max -= distance;	2,308,614✔
133
134	if (max == 0.0) {	2,308,614!
NEW 135	update_distance();	×
NEW 136	break;	×
137	}
138
139	surface() = boundary().surface();	2,308,614✔
140	// Initialize last cells from the current cell, because the cell() variable
141	// does not contain the data for the case of a single-segment ray
142	for (int j = 0; j < n_coord(); ++j) {	4,617,228✔
143	cell_last(j) = coord(j).cell();	2,308,614✔
144	}
145	n_coord_last() = n_coord();	2,308,614✔
146	n_coord() = boundary().coord_level();	2,308,614!
147	if (boundary().lattice_translation()[0] != 0 \|\|	2,308,614!
148	boundary().lattice_translation()[1] != 0 \|\|	2,308,614!
149	boundary().lattice_translation()[2] != 0) {	2,308,614!
150	cross_lattice(*this, boundary(), settings::verbosity >= 10);	×
151	}
152
153	update_distance();	2,308,614✔
154
155	inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);	2,308,614✔
156
157	// Call the specialized logic for this type of ray. Note that we do not
158	// call this if the advance distance is very small. Unfortunately, it seems
159	// darn near impossible to get the particle advanced to the model boundary
160	// and through it without sometimes accidentally calling on_intersection
161	// twice. This incorrectly shades the region as occluded when it might not
162	// actually be. By screening out intersection distances smaller than a
163	// threshold 10x larger than the scoot distance used to advance up to the
164	// model boundary, we can avoid that situation.
165	if (call_on_intersection) {	2,308,614✔
166	on_intersection();	1,715,329✔
167	if (stop_)	1,715,329✔
168	return;
169	}
170
171	if (!inside_cell)	2,273,051✔
172	return;
173
174	event_counter_++;	758,758✔
175	if (event_counter_ > MAX_INTERSECTIONS) {	758,758!
176	warning("Likely infinite loop in ray traced plot");	×
177	return;	×
178	}
179	}
180	}
181
182	void Ray::update_distance()	2,308,614✔
183	{
184	// Record how far the ray has traveled
185	traversal_distance_ += boundary().distance();	2,308,614✔
186	}	2,308,614✔
187
NEW 188	void ParticleRay::on_intersection() {}	×
189
NEW 190	void ParticleRay::update_distance()	×
191	{
NEW 192	Ray::update_distance();	×
193
NEW 194	time() += boundary().distance() / speed();	×
195
196	// Calculate microscopic and macroscopic cross sections
NEW 197	if (material() != MATERIAL_VOID) {	×
NEW 198	if (settings::run_CE) {	×
NEW 199	if (material() != material_last() \|\| sqrtkT() != sqrtkT_last() \|\|	×
NEW 200	density_mult() != density_mult_last()) {	×
201	// If the material is the same as the last material and the
202	// temperature hasn't changed, we don't need to lookup cross
203	// sections again.
NEW 204	model::materials[material()]->calculate_xs(*this);	×
205	}
206	} else {
207	// Get the MG data; unlike the CE case above, we have to re-calculate
208	// cross sections for every collision since the cross sections may
209	// be angle-dependent
NEW 210	data::mg.macro_xs_[material()].calculate_xs(*this);	×
211
212	// Update the particle's group while we know we are multi-group
NEW 213	g_last() = g();	×
214	}
215	} else {
NEW 216	macro_xs().total = 0.0;	×
NEW 217	macro_xs().absorption = 0.0;	×
NEW 218	macro_xs().fission = 0.0;	×
NEW 219	macro_xs().nu_fission = 0.0;	×
220	}
221
NEW 222	traversal_mfp_ += macro_xs().total * boundary().distance();	×
NEW 223	}	×
224
225	} // namespace openmc

openmc-dev / openmc / 30387632100

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