22073868368

Committed 16 Feb 2026 06:32PM UTC coverage: 81.546% (-0.2%) from 81.721%

Build # 22073868368

Build Type

Pull #3810

github

Committed by

web-flow

Commit Message

Merge fe77e3fdb into c6ef84d1d

Pull Request Pull Request #3810: Implementation of migration area score

Run Details

17050 of 23727 branches covered (71.86%)

Branch coverage included in aggregate %.

36 of 55 new or added lines in 5 files covered. (65.45%)

300 existing lines in 24 files now uncovered.

56162 of 66053 relevant lines covered (85.03%)

43232737.44 hits per line

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

70.86

/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/simulation.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);
  u /= u.norm();

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

  // Handle phantom birth location if migration present
  if (simulation::migration_present) {
    auto r = p.r();
    auto n = surf.normal(r);
    n /= n.norm();
    auto r_born = p.r_born();
    p.r_born() = r_born - 2.0 * (r_born - r).dot(n) * n;
  }

  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());
  u /= u.norm();

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

  // Handle phantom birth location if migration present
  if (simulation::migration_present)
    fatal_error("Cannot tally migration area with white boundary conditions.");

  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);

  // Handle phantom birth location if migration present
  if (simulation::migration_present) {
    p.r_born() += translation_;
  }

  // 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_];

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

  // Handle phantom birth location if migration present
  if (simulation::migration_present) {
    auto r_born = p.r_born();
    Position new_r_born;
    new_r_born[zero_axis_idx_] = r_born[zero_axis_idx_];
    new_r_born[axis_1_idx_] =
      cos_theta * r_born[axis_1_idx_] - sin_theta * r_born[axis_2_idx_];
    new_r_born[axis_2_idx_] =
      sin_theta * r_born[axis_1_idx_] + cos_theta * r_born[axis_2_idx_];
    p.r_born() = new_r_born;
  }

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

openmc-dev / openmc / 22073868368

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