29642010356

Committed 18 Jul 2026 11:07AM UTC coverage: 81.335% (+0.03%) from 81.301%

Build # 29642010356

Build Type

Pull #4014

github

Committed by

web-flow

Commit Message

Merge 32f249b12 into db673b9ac

Pull Request Pull Request #4014: Add analytic tests for ray-traced intersection distances

Coverage Stats

18352 of 26615 branches covered (68.95%)

Branch coverage included in aggregate %.

59867 of 69554 relevant lines covered (86.07%)

47903540.77 hits per line

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

68.83

/src/boundary_condition.cpp

#include "openmc/boundary_condition.h"

#include <exception>

#include <fmt/core.h>

#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/random_ray/random_ray.h"
#include "openmc/surface.h"

namespace openmc {

//==============================================================================
// VacuumBC implementation
//==============================================================================

void VacuumBC::handle_particle(Particle& p, const Surface& surf) const
{
  // Random ray and Monte Carlo need different treatments at vacuum BCs
  if (settings::solver_type == SolverType::RANDOM_RAY) {
    // Reflect ray off of the surface
    ReflectiveBC().handle_particle(p, surf);

    // Set ray's angular flux spectrum to vacuum conditions (zero)
    RandomRay* r = static_cast<RandomRay*>(&p);
    std::fill(r->angular_flux_.begin(), r->angular_flux_.end(), 0.0);

  } else {
    p.cross_vacuum_bc(surf);
  }
}

//==============================================================================
// ReflectiveBC implementation
//==============================================================================

void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
{
  Direction u = surf.reflect(p.r(), p.u(), &p);

  // normalize reflected u to ensure no floating point error leads to
  // unnormalized directions
  u /= u.norm();

  // Handle the effects of the surface albedo on the particle's weight.
  BoundaryCondition::handle_albedo(p, surf);

  p.cross_reflective_bc(surf, u);
}

//==============================================================================
// WhiteBC implementation
//==============================================================================

void WhiteBC::handle_particle(Particle& p, const Surface& surf) const
{
  Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed());

  // normalize outgoing u to ensure no floating point error leads to
  // unnormalized directions
  u /= u.norm();

  // Handle the effects of the surface albedo on the particle's weight.
  BoundaryCondition::handle_albedo(p, surf);

  p.cross_reflective_bc(surf, u);
}

//==============================================================================
// TranslationalPeriodicBC implementation
//==============================================================================

TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf)
  : PeriodicBC(i_surf, j_surf)
{
  Surface& surf1 {*model::surfaces[i_surf_]};
  Surface& surf2 {*model::surfaces[j_surf_]};

  // Make sure the first surface has an appropriate type.
  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
  } else if (const auto* ptr = dynamic_cast<const SurfaceZPlane*>(&surf1)) {
  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf1)) {
  } else {
    throw std::invalid_argument(fmt::format(
      "Surface {} is an invalid type for "
      "translational periodic BCs. Only planes are supported for these BCs.",
      surf1.id_));
  }

  // Make sure the second surface has an appropriate type.
  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
  } else if (const auto* ptr = dynamic_cast<const SurfaceZPlane*>(&surf2)) {
  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf2)) {
  } else {
    throw std::invalid_argument(fmt::format(
      "Surface {} is an invalid type for "
      "translational periodic BCs. Only planes are supported for these BCs.",
      surf2.id_));
  }

  // Compute the distance from the first surface to the origin.  Check the
  // surface evaluate function to decide if the distance is positive, negative,
  // or zero.
  Position origin {0, 0, 0};
  Direction u = surf1.normal(origin);
  double d1;
  double e1 = surf1.evaluate(origin);
  if (e1 > FP_COINCIDENT) {
    d1 = -surf1.distance(origin, -u, false);
  } else if (e1 < -FP_COINCIDENT) {
    d1 = surf1.distance(origin, u, false);
  } else {
    d1 = 0.0;
  }

  // Compute the distance from the second surface to the origin.
  double d2;
  double e2 = surf2.evaluate(origin);
  if (e2 > FP_COINCIDENT) {
    d2 = -surf2.distance(origin, -u, false);
  } else if (e2 < -FP_COINCIDENT) {
    d2 = surf2.distance(origin, u, false);
  } else {
    d2 = 0.0;
  }

  // Set the translation vector; it's length is the difference in the two
  // distances.
  translation_ = u * (d2 - d1);
}

void TranslationalPeriodicBC::handle_particle(
  Particle& p, const Surface& surf) const
{
  auto new_r = p.r() + translation_;
  int new_surface = p.surface() > 0 ? j_surf_ + 1 : -(j_surf_ + 1);

  // Handle the effects of the surface albedo on the particle's weight.
  BoundaryCondition::handle_albedo(p, surf);

  // Pass the new location and surface to the particle.
  p.cross_periodic_bc(surf, new_r, p.u(), new_surface);
}

//==============================================================================
// RotationalPeriodicBC implementation
//==============================================================================

RotationalPeriodicBC::RotationalPeriodicBC(
  int i_surf, int j_surf, PeriodicAxis axis)
  : PeriodicBC(std::abs(i_surf) - 1, std::abs(j_surf) - 1)
{
  Surface& surf1 {*model::surfaces[i_surf_]};
  Surface& surf2 {*model::surfaces[j_surf_]};

  // below convention for right handed coordinate system
  switch (axis) {
  case x:
    zero_axis_idx_ = 0; // x component of plane must be zero
    axis_1_idx_ = 1;    // y component independent
    axis_2_idx_ = 2;    // z component dependent
    break;
  case y:
    zero_axis_idx_ = 1; // y component of plane must be zero
    axis_1_idx_ = 2;    // z component independent
    axis_2_idx_ = 0;    // x component dependent
    break;
  case z:
    zero_axis_idx_ = 2; // z component of plane must be zero
    axis_1_idx_ = 0;    // x component independent
    axis_2_idx_ = 1;    // y component dependent
    break;
  default:
    throw std::invalid_argument(
      fmt::format("You've specified an axis that is not x, y, or z."));
  }

  Direction ax = {0.0, 0.0, 0.0};
  ax[zero_axis_idx_] = 1.0;

  auto i_sign = std::copysign(1, i_surf);
  auto j_sign = -std::copysign(1, j_surf);

  // Compute the surface normal vectors and make sure they are perpendicular
  // to the correct axis
  Direction norm1 = i_sign * surf1.normal({0, 0, 0});
  Direction norm2 = j_sign * surf2.normal({0, 0, 0});
  // Make sure both surfaces intersect the origin
  if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {
    throw std::invalid_argument(fmt::format(
      "Rotational periodic BCs are only "
      "supported for rotations about the origin, but surface {} does not "
      "intersect the origin.",
      surf1.id_));
  }
  if (std::abs(surf2.evaluate({0, 0, 0})) > FP_COINCIDENT) {
    throw std::invalid_argument(fmt::format(
      "Rotational periodic BCs are only "
      "supported for rotations about the origin, but surface {} does not "
      "intersect the origin.",
      surf2.id_));
  }

  // Compute the signed rotation angle about the periodic axis. Note that
  // (n1×n2)·a = |n1||n2|sin(θ) and n1·n2 = |n1||n2|cos(θ), where a is the axis
  // of rotation.
  auto c = norm1.cross(norm2);
  angle_ = std::atan2(c.dot(ax), norm1.dot(norm2));

  // If the normals point in the same general direction, the surface sense
  // should change when crossing the boundary
  flip_sense_ = (i_sign * j_sign > 0.0);

  // Warn the user if the angle does not evenly divide a circle
  double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
  if (rem > FP_REL_PRECISION && rem < 1 - FP_REL_PRECISION) {
    warning(fmt::format(
      "Rotational periodic BC specified with a rotation "
      "angle of {} degrees which does not evenly divide 360 degrees.",
      angle_ * 180 / PI));
  }
}

void RotationalPeriodicBC::handle_particle(
  Particle& p, const Surface& surf) const
{
  int new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1;
  if (flip_sense_)
    new_surface = -new_surface;

  // Rotate the particle's position and direction.
  Position r = p.r();
  Direction u = p.u();
  double cos_theta = std::cos(angle_);
  double sin_theta = std::sin(angle_);

  Position new_r;
  new_r[zero_axis_idx_] = r[zero_axis_idx_];
  new_r[axis_1_idx_] = cos_theta * r[axis_1_idx_] - sin_theta * r[axis_2_idx_];
  new_r[axis_2_idx_] = sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_];

  Direction new_u;
  new_u[zero_axis_idx_] = u[zero_axis_idx_];
  new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_];
  new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];

  // normalize new_u to ensure no floating point error leads to unnormalized
  // directions
  new_u /= new_u.norm();

  // Handle the effects of the surface albedo on the particle's weight.
  BoundaryCondition::handle_albedo(p, surf);

  // Pass the new location, direction, and surface to the particle.
  p.cross_periodic_bc(surf, new_r, new_u, new_surface);
}

} // namespace openmc

1	#include "openmc/boundary_condition.h"
2
3	#include <exception>
4
5	#include <fmt/core.h>
6
7	#include "openmc/constants.h"
8	#include "openmc/error.h"
9	#include "openmc/random_ray/random_ray.h"
10	#include "openmc/surface.h"
11
12	namespace openmc {
13
14	//==============================================================================
15	// VacuumBC implementation
16	//==============================================================================
17
18	void VacuumBC::handle_particle(Particle& p, const Surface& surf) const	63,442,873✔
19	{
20	// Random ray and Monte Carlo need different treatments at vacuum BCs
21	if (settings::solver_type == SolverType::RANDOM_RAY) {	63,442,873✔
22	// Reflect ray off of the surface
23	ReflectiveBC().handle_particle(p, surf);	27,421,522✔
24
25	// Set ray's angular flux spectrum to vacuum conditions (zero)
26	RandomRay* r = static_cast<RandomRay*>(&p);	27,421,522✔
27	std::fill(r->angular_flux_.begin(), r->angular_flux_.end(), 0.0);	27,421,522✔
28
29	} else {
30	p.cross_vacuum_bc(surf);	36,021,351✔
31	}
32	}	63,442,873✔
33
34	//==============================================================================
35	// ReflectiveBC implementation
36	//==============================================================================
37
38	void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const	969,428,442✔
39	{
40	Direction u = surf.reflect(p.r(), p.u(), &p);	969,428,442✔
41
42	// normalize reflected u to ensure no floating point error leads to
43	// unnormalized directions
44	u /= u.norm();	969,428,442✔
45
46	// Handle the effects of the surface albedo on the particle's weight.
47	BoundaryCondition::handle_albedo(p, surf);	969,428,442✔
48
49	p.cross_reflective_bc(surf, u);	969,428,442✔
50	}	969,428,442✔
51
52	//==============================================================================
53	// WhiteBC implementation
54	//==============================================================================
55
56	void WhiteBC::handle_particle(Particle& p, const Surface& surf) const	1,007,149✔
57	{
58	Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed());	1,007,149✔
59
60	// normalize outgoing u to ensure no floating point error leads to
61	// unnormalized directions
62	u /= u.norm();	1,007,149✔
63
64	// Handle the effects of the surface albedo on the particle's weight.
65	BoundaryCondition::handle_albedo(p, surf);	1,007,149✔
66
67	p.cross_reflective_bc(surf, u);	1,007,149✔
68	}	1,007,149✔
69
70	//==============================================================================
71	// TranslationalPeriodicBC implementation
72	//==============================================================================
73
74	TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf)	226✔
75	: PeriodicBC(i_surf, j_surf)	226!
76	{
77	Surface& surf1 {*model::surfaces[i_surf_]};	226!
78	Surface& surf2 {*model::surfaces[j_surf_]};	226✔
79
80	// Make sure the first surface has an appropriate type.
81	if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {	226!
82	} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {	174!
83	} else if (const auto* ptr = dynamic_cast<const SurfaceZPlane*>(&surf1)) {	152!
84	} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf1)) {	75!
85	} else {
86	throw std::invalid_argument(fmt::format(	×
87	"Surface {} is an invalid type for "
88	"translational periodic BCs. Only planes are supported for these BCs.",
89	surf1.id_));	×
90	}
91
92	// Make sure the second surface has an appropriate type.
93	if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {	226!
94	} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {	174!
95	} else if (const auto* ptr = dynamic_cast<const SurfaceZPlane*>(&surf2)) {	152!
96	} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf2)) {	75!
97	} else {
98	throw std::invalid_argument(fmt::format(	×
99	"Surface {} is an invalid type for "
100	"translational periodic BCs. Only planes are supported for these BCs.",
101	surf2.id_));	×
102	}
103
104	// Compute the distance from the first surface to the origin. Check the
105	// surface evaluate function to decide if the distance is positive, negative,
106	// or zero.
107	Position origin {0, 0, 0};	226✔
108	Direction u = surf1.normal(origin);	226✔
109	double d1;	226✔
110	double e1 = surf1.evaluate(origin);	226✔
111	if (e1 > FP_COINCIDENT) {	226✔
112	d1 = -surf1.distance(origin, -u, false);	113✔
113	} else if (e1 < -FP_COINCIDENT) {	113!
114	d1 = surf1.distance(origin, u, false);	113✔
115	} else {
116	d1 = 0.0;
117	}
118
119	// Compute the distance from the second surface to the origin.
120	double d2;	226✔
121	double e2 = surf2.evaluate(origin);	226✔
122	if (e2 > FP_COINCIDENT) {	226✔
123	d2 = -surf2.distance(origin, -u, false);	113✔
124	} else if (e2 < -FP_COINCIDENT) {	113!
125	d2 = surf2.distance(origin, u, false);	113✔
126	} else {
127	d2 = 0.0;
128	}
129
130	// Set the translation vector; it's length is the difference in the two
131	// distances.
132	translation_ = u * (d2 - d1);	226✔
133	}	226✔
134
135	void TranslationalPeriodicBC::handle_particle(	296,797✔
136	Particle& p, const Surface& surf) const
137	{
138	auto new_r = p.r() + translation_;	296,797✔
139	int new_surface = p.surface() > 0 ? j_surf_ + 1 : -(j_surf_ + 1);	296,797✔
140
141	// Handle the effects of the surface albedo on the particle's weight.
142	BoundaryCondition::handle_albedo(p, surf);	296,797✔
143
144	// Pass the new location and surface to the particle.
145	p.cross_periodic_bc(surf, new_r, p.u(), new_surface);	296,797✔
146	}	296,797✔
147
148	//==============================================================================
149	// RotationalPeriodicBC implementation
150	//==============================================================================
151
152	RotationalPeriodicBC::RotationalPeriodicBC(	320✔
153	int i_surf, int j_surf, PeriodicAxis axis)	320✔
154	: PeriodicBC(std::abs(i_surf) - 1, std::abs(j_surf) - 1)	320!
155	{
156	Surface& surf1 {*model::surfaces[i_surf_]};	320!
157	Surface& surf2 {*model::surfaces[j_surf_]};	320✔
158
159	// below convention for right handed coordinate system
160	switch (axis) {	320!
161	case x:	30✔
162	zero_axis_idx_ = 0; // x component of plane must be zero	30✔
163	axis_1_idx_ = 1; // y component independent	30✔
164	axis_2_idx_ = 2; // z component dependent	30✔
165	break;	30✔
166	case y:	30✔
167	zero_axis_idx_ = 1; // y component of plane must be zero	30✔
168	axis_1_idx_ = 2; // z component independent	30✔
169	axis_2_idx_ = 0; // x component dependent	30✔
170	break;	30✔
171	case z:	260✔
172	zero_axis_idx_ = 2; // z component of plane must be zero	260✔
173	axis_1_idx_ = 0; // x component independent	260✔
174	axis_2_idx_ = 1; // y component dependent	260✔
175	break;	260✔
176	default:	×
177	throw std::invalid_argument(	×
178	fmt::format("You've specified an axis that is not x, y, or z."));	×
179	}
180
181	Direction ax = {0.0, 0.0, 0.0};	320✔
182	ax[zero_axis_idx_] = 1.0;	320✔
183
184	auto i_sign = std::copysign(1, i_surf);	320✔
185	auto j_sign = -std::copysign(1, j_surf);	320✔
186
187	// Compute the surface normal vectors and make sure they are perpendicular
188	// to the correct axis
189	Direction norm1 = i_sign * surf1.normal({0, 0, 0});	320✔
190	Direction norm2 = j_sign * surf2.normal({0, 0, 0});	320✔
191	// Make sure both surfaces intersect the origin
192	if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {	320!
193	throw std::invalid_argument(fmt::format(	×
194	"Rotational periodic BCs are only "
195	"supported for rotations about the origin, but surface {} does not "
196	"intersect the origin.",
197	surf1.id_));	×
198	}
199	if (std::abs(surf2.evaluate({0, 0, 0})) > FP_COINCIDENT) {	320!
200	throw std::invalid_argument(fmt::format(	×
201	"Rotational periodic BCs are only "
202	"supported for rotations about the origin, but surface {} does not "
203	"intersect the origin.",
204	surf2.id_));	×
205	}
206
207	// Compute the signed rotation angle about the periodic axis. Note that
208	// (n1×n2)·a = \|n1\|\|n2\|sin(θ) and n1·n2 = \|n1\|\|n2\|cos(θ), where a is the axis
209	// of rotation.
210	auto c = norm1.cross(norm2);	320✔
211	angle_ = std::atan2(c.dot(ax), norm1.dot(norm2));	320✔
212
213	// If the normals point in the same general direction, the surface sense
214	// should change when crossing the boundary
215	flip_sense_ = (i_sign * j_sign > 0.0);	320✔
216
217	// Warn the user if the angle does not evenly divide a circle
218	double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));	320!
219	if (rem > FP_REL_PRECISION && rem < 1 - FP_REL_PRECISION) {	320!
220	warning(fmt::format(	×
221	"Rotational periodic BC specified with a rotation "
222	"angle of {} degrees which does not evenly divide 360 degrees.",
223	angle_ * 180 / PI));	×
224	}
225	}	320✔
226
227	void RotationalPeriodicBC::handle_particle(	1,947,640✔
228	Particle& p, const Surface& surf) const
229	{
230	int new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1;	1,947,640✔
231	if (flip_sense_)	1,947,640✔
232	new_surface = -new_surface;	851,269✔
233
234	// Rotate the particle's position and direction.
235	Position r = p.r();	1,947,640✔
236	Direction u = p.u();	1,947,640✔
237	double cos_theta = std::cos(angle_);	1,947,640✔
238	double sin_theta = std::sin(angle_);	1,947,640✔
239
240	Position new_r;	1,947,640✔
241	new_r[zero_axis_idx_] = r[zero_axis_idx_];	1,947,640✔
242	new_r[axis_1_idx_] = cos_theta * r[axis_1_idx_] - sin_theta * r[axis_2_idx_];	1,947,640✔
243	new_r[axis_2_idx_] = sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_];	1,947,640✔
244
245	Direction new_u;	1,947,640✔
246	new_u[zero_axis_idx_] = u[zero_axis_idx_];	1,947,640✔
247	new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_];	1,947,640✔
248	new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];	1,947,640✔
249
250	// normalize new_u to ensure no floating point error leads to unnormalized
251	// directions
252	new_u /= new_u.norm();	1,947,640✔
253
254	// Handle the effects of the surface albedo on the particle's weight.
255	BoundaryCondition::handle_albedo(p, surf);	1,947,640✔
256
257	// Pass the new location, direction, and surface to the particle.
258	p.cross_periodic_bc(surf, new_r, new_u, new_surface);	1,947,640✔
259	}	1,947,640✔
260
261	} // namespace openmc

openmc-dev / openmc / 29642010356

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