• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 19518401907

19 Nov 2025 10:25PM UTC coverage: 82.121% (+0.05%) from 82.071%
19518401907

Pull #3591

github

web-flow
Merge 6cf5bfc86 into f544d02e4
Pull Request #3591: Generalize RotationalPeriodicBC for X-, Y-, or Z-axis

16992 of 23535 branches covered (72.2%)

Branch coverage included in aggregate %.

53 of 59 new or added lines in 2 files covered. (89.83%)

62 existing lines in 3 files now uncovered.

54937 of 64054 relevant lines covered (85.77%)

42160829.94 hits per line

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

79.67
/src/boundary_condition.cpp
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
57,022,606✔
19
{
20
  // Random ray and Monte Carlo need different treatments at vacuum BCs
21
  if (settings::solver_type == SolverType::RANDOM_RAY) {
57,022,606✔
22
    // Reflect ray off of the surface
23
    ReflectiveBC().handle_particle(p, surf);
23,673,129✔
24

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

29
  } else {
30
    p.cross_vacuum_bc(surf);
33,349,477✔
31
  }
32
}
57,022,606✔
33

34
//==============================================================================
35
// ReflectiveBC implementation
36
//==============================================================================
37

38
void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
612,887,246✔
39
{
40
  Direction u = surf.reflect(p.r(), p.u(), &p);
612,887,246✔
41
  u /= u.norm();
612,887,246✔
42

43
  // Handle the effects of the surface albedo on the particle's weight.
44
  BoundaryCondition::handle_albedo(p, surf);
612,887,246✔
45

46
  p.cross_reflective_bc(surf, u);
612,887,246✔
47
}
612,887,246✔
48

49
//==============================================================================
50
// WhiteBC implementation
51
//==============================================================================
52

53
void WhiteBC::handle_particle(Particle& p, const Surface& surf) const
1,007,149✔
54
{
55
  Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed());
1,007,149✔
56
  u /= u.norm();
1,007,149✔
57

58
  // Handle the effects of the surface albedo on the particle's weight.
59
  BoundaryCondition::handle_albedo(p, surf);
1,007,149✔
60

61
  p.cross_reflective_bc(surf, u);
1,007,149✔
62
}
1,007,149✔
63

64
//==============================================================================
65
// TranslationalPeriodicBC implementation
66
//==============================================================================
67

68
TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf)
234✔
69
  : PeriodicBC(i_surf, j_surf)
234✔
70
{
71
  Surface& surf1 {*model::surfaces[i_surf_]};
234✔
72
  Surface& surf2 {*model::surfaces[j_surf_]};
234✔
73

74
  // Make sure the first surface has an appropriate type.
75
  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
234!
76
  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
180!
77
  } else if (const auto* ptr = dynamic_cast<const SurfaceZPlane*>(&surf1)) {
158!
78
  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf1)) {
64!
79
  } else {
80
    throw std::invalid_argument(fmt::format(
×
81
      "Surface {} is an invalid type for "
82
      "translational periodic BCs. Only planes are supported for these BCs.",
83
      surf1.id_));
×
84
  }
85

86
  // Make sure the second surface has an appropriate type.
87
  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
234!
88
  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
180!
89
  } else if (const auto* ptr = dynamic_cast<const SurfaceZPlane*>(&surf2)) {
158!
90
  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf2)) {
96!
91
  } else {
92
    throw std::invalid_argument(fmt::format(
×
93
      "Surface {} is an invalid type for "
94
      "translational periodic BCs. Only planes are supported for these BCs.",
95
      surf2.id_));
×
96
  }
97

98
  // Compute the distance from the first surface to the origin.  Check the
99
  // surface evaluate function to decide if the distance is positive, negative,
100
  // or zero.
101
  Position origin {0, 0, 0};
234✔
102
  Direction u = surf1.normal(origin);
234✔
103
  double d1;
104
  double e1 = surf1.evaluate(origin);
234✔
105
  if (e1 > FP_COINCIDENT) {
234✔
106
    d1 = -surf1.distance(origin, -u, false);
138✔
107
  } else if (e1 < -FP_COINCIDENT) {
96!
108
    d1 = surf1.distance(origin, u, false);
96✔
109
  } else {
110
    d1 = 0.0;
×
111
  }
112

113
  // Compute the distance from the second surface to the origin.
114
  double d2;
115
  double e2 = surf2.evaluate(origin);
234✔
116
  if (e2 > FP_COINCIDENT) {
234✔
117
    d2 = -surf2.distance(origin, -u, false);
96✔
118
  } else if (e2 < -FP_COINCIDENT) {
138!
119
    d2 = surf2.distance(origin, u, false);
138✔
120
  } else {
121
    d2 = 0.0;
×
122
  }
123

124
  // Set the translation vector; it's length is the difference in the two
125
  // distances.
126
  translation_ = u * (d2 - d1);
234✔
127
}
234✔
128

129
void TranslationalPeriodicBC::handle_particle(
296,797✔
130
  Particle& p, const Surface& surf) const
131
{
132
  int i_particle_surf = p.surface_index();
296,797✔
133

134
  // Figure out which of the two BC surfaces were struck then find the
135
  // particle's new location and surface.
136
  Position new_r;
296,797✔
137
  int new_surface;
138
  if (i_particle_surf == i_surf_) {
296,797✔
139
    new_r = p.r() + translation_;
147,047✔
140
    new_surface = p.surface() > 0 ? j_surf_ + 1 : -(j_surf_ + 1);
147,047✔
141
  } else if (i_particle_surf == j_surf_) {
149,750!
142
    new_r = p.r() - translation_;
149,750✔
143
    new_surface = p.surface() > 0 ? i_surf_ + 1 : -(i_surf_ + 1);
149,750✔
144
  } else {
145
    throw std::runtime_error(
×
146
      "Called BoundaryCondition::handle_particle after "
147
      "hitting a surface, but that surface is not recognized by the BC.");
×
148
  }
149

150
  // Handle the effects of the surface albedo on the particle's weight.
151
  BoundaryCondition::handle_albedo(p, surf);
296,797✔
152

153
  // Pass the new location and surface to the particle.
154
  p.cross_periodic_bc(surf, new_r, p.u(), new_surface);
296,797✔
155
}
296,797✔
156

157
//==============================================================================
158
// RotationalPeriodicBC implementation
159
//==============================================================================
160

161
RotationalPeriodicBC::RotationalPeriodicBC(
128✔
162
  int i_surf, int j_surf, PeriodicAxis axis)
128✔
163
  : PeriodicBC(i_surf, j_surf)
128✔
164
{
165
  Surface& surf1 {*model::surfaces[i_surf_]};
128✔
166
  Surface& surf2 {*model::surfaces[j_surf_]};
128✔
167

168
  // below convention for right handed coordinate system
169
  switch (axis) {
128!
170
  case x:
32✔
171
    zero_axis_idx_ = 0; // x component of plane must be zero
32✔
172
    axis_1_idx_ = 1;    // y component independent
32✔
173
    axis_2_idx_ = 2;    // z component dependent
32✔
174
    break;
32✔
175
  case y:
32✔
176
    zero_axis_idx_ = 1; // y component of plane must be zero
32✔
177
    axis_1_idx_ = 2;    // z component independent
32✔
178
    axis_2_idx_ = 0;    // x component dependent
32✔
179
    break;
32✔
180
  case z:
64✔
181
    zero_axis_idx_ = 2; // z component of plane must be zero
64✔
182
    axis_1_idx_ = 0;    // x component independent
64✔
183
    axis_2_idx_ = 1;    // y component dependent
64✔
184
    break;
64✔
NEW
185
  default:
×
NEW
186
    throw std::invalid_argument(
×
NEW
187
      fmt::format("You've specified an axis that is not x, y, or z."));
×
188
  }
189

190
  // Compute the surface normal vectors and make sure they are perpendicular
191
  // to the correct axis
192
  Direction norm1 = surf1.normal({0, 0, 0});
128✔
193
  Direction norm2 = surf2.normal({0, 0, 0});
128✔
194
  // Make sure both surfaces intersect the origin
195
  if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {
128!
196
    throw std::invalid_argument(fmt::format(
×
197
      "Rotational periodic BCs are only "
198
      "supported for rotations about the origin, but surface {} does not "
199
      "intersect the origin.",
200
      surf1.id_));
×
201
  }
202
  if (std::abs(surf2.evaluate({0, 0, 0})) > FP_COINCIDENT) {
128!
UNCOV
203
    throw std::invalid_argument(fmt::format(
×
204
      "Rotational periodic BCs are only "
205
      "supported for rotations about the origin, but surface {} does not "
206
      "intersect the origin.",
UNCOV
207
      surf2.id_));
×
208
  }
209

210
  // reverse the angle computed if there is y-periodicity to account for
211
  // reversal of axes cross product direction
212
  angle_ =
128✔
213
    (zero_axis_idx_ != 1)
128✔
214
      ? compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_],
224✔
215
          norm2[axis_2_idx_], norm2[axis_1_idx_])
96✔
216
      : -1 * compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_],
32✔
217
               norm2[axis_2_idx_], norm2[axis_1_idx_]);
32✔
218

219
  // Warn the user if the angle does not evenly divide a circle
220
  double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
128✔
221
  if (rem > FP_REL_PRECISION && rem < 1 - FP_REL_PRECISION) {
128!
222
    warning(fmt::format(
64✔
223
      "Rotational periodic BC specified with a rotation "
224
      "angle of {} degrees which does not evenly divide 360 degrees.",
225
      angle_ * 180 / PI));
128✔
226
  }
227
}
128✔
228

229
double RotationalPeriodicBC::compute_periodic_rotation(
128✔
230
  double rise_1, double run_1, double rise_2, double run_2) const
231
{
232
  // Compute the BC rotation angle.  Here it is assumed that both surface
233
  // normal vectors point inwards---towards the valid geometry region.
234
  // Consequently, the rotation angle is not the difference between the two
235
  // normals, but is instead the difference between one normal and one
236
  // anti-normal.  (An incident ray on one surface must be an outgoing ray on
237
  // the other surface after rotation hence the anti-normal.)
238
  double theta1 = std::atan2(rise_1, run_1);
128✔
239
  double theta2 = std::atan2(rise_2, run_2) + PI;
128✔
240
  return theta2 - theta1;
128✔
241
}
242

243
void RotationalPeriodicBC::handle_particle(
510,136✔
244
  Particle& p, const Surface& surf) const
245
{
246
  int i_particle_surf = p.surface_index();
510,136✔
247

248
  // Figure out which of the two BC surfaces were struck to figure out if a
249
  // forward or backward rotation is required.  Specify the other surface as
250
  // the particle's new surface.
251
  double theta;
252
  int new_surface;
253
  if (i_particle_surf == i_surf_) {
510,136✔
254
    theta = angle_;
255,860✔
255
    new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1;
255,860!
256
  } else if (i_particle_surf == j_surf_) {
254,276!
257
    theta = -angle_;
254,276✔
258
    new_surface = p.surface() > 0 ? -(i_surf_ + 1) : i_surf_ + 1;
254,276!
259
  } else {
UNCOV
260
    throw std::runtime_error(
×
261
      "Called BoundaryCondition::handle_particle after "
UNCOV
262
      "hitting a surface, but that surface is not recognized by the BC.");
×
263
  }
264

265
  // Rotate the particle's position and direction about the z-axis.
266
  Position r = p.r();
510,136✔
267
  Direction u = p.u();
510,136✔
268
  double cos_theta = std::cos(theta);
510,136✔
269
  double sin_theta = std::sin(theta);
510,136✔
270

271
  // rotations around the y-axis have sign flipped for the sin_theta terms
272
  int flip;
273
  if (zero_axis_idx_ == 1) {
510,136✔
274
    flip = -1;
75,119✔
275
  } else {
276
    flip = 1;
435,017✔
277
  }
278

279
  Position new_r;
510,136✔
280
  new_r[zero_axis_idx_] = r[zero_axis_idx_];
510,136✔
281
  new_r[axis_1_idx_] =
1,020,272✔
282
    cos_theta * r[axis_1_idx_] - flip * sin_theta * r[axis_2_idx_];
510,136✔
283
  new_r[axis_2_idx_] =
1,020,272✔
284
    flip * sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_];
510,136✔
285

286
  Direction new_u;
510,136✔
287
  new_u[zero_axis_idx_] = u[zero_axis_idx_];
510,136✔
288
  new_u[axis_1_idx_] =
1,020,272✔
289
    cos_theta * u[axis_1_idx_] - flip * sin_theta * u[axis_2_idx_];
510,136✔
290
  new_u[axis_2_idx_] =
1,020,272✔
291
    flip * sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];
510,136✔
292

293
  // Handle the effects of the surface albedo on the particle's weight.
294
  BoundaryCondition::handle_albedo(p, surf);
510,136✔
295

296
  // Pass the new location, direction, and surface to the particle.
297
  p.cross_periodic_bc(surf, new_r, new_u, new_surface);
510,136✔
298
}
510,136✔
299

300
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc