diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index 7524856c46..caf83b6117 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -84,6 +84,7 @@ public: protected: double angle_; + bool aligned_normals_; }; } // namespace openmc diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index e671ad554c..67496187f6 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -5,6 +5,7 @@ #include #include "openmc/constants.h" +#include "openmc/error.h" #include "openmc/surface.h" namespace openmc { @@ -55,9 +56,6 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) Surface& surf1 {*model::surfaces[i_surf_]}; Surface& surf2 {*model::surfaces[j_surf_]}; - // The following blocks will resolve the type of each surface and compute the - // appropriate translation vector - // Make sure the first surface has an appropriate type. if (const auto* ptr = dynamic_cast(&surf1)) { } else if (const auto* ptr = dynamic_cast(&surf1)) { @@ -146,18 +144,81 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) Surface& surf1 {*model::surfaces[i_surf_]}; Surface& surf2 {*model::surfaces[j_surf_]}; - // The following blocks will resolve the type of each surface and compute the - // appropriate rotation angle - - if (const auto* xplane = dynamic_cast(&surf1)) { - if (const auto* yplane = dynamic_cast(&surf2)) { - angle_ = 0.5 * PI; - } - } else if (const auto* yplane = dynamic_cast(&surf1)) { - if (const auto* xplane = dynamic_cast(&surf2)) { - angle_ = -0.5 * PI; - } + // Check the type of the first surface + bool surf1_is_xyplane; + if (const auto* ptr = dynamic_cast(&surf1)) { + surf1_is_xyplane = true; + } else if (const auto* ptr = dynamic_cast(&surf1)) { + surf1_is_xyplane = true; + } else if (const auto* ptr = dynamic_cast(&surf1)) { + surf1_is_xyplane = false; + } else { + throw std::invalid_argument(fmt::format("Surface {} is an invalid type for " + "rotational periodic BCs. Only x-planes, y-planes, or general planes " + "(that are perpendicular to z) are supported for these BCs.", surf1.id_)); } + + // Check the type of the second surface + bool surf2_is_xyplane; + if (const auto* ptr = dynamic_cast(&surf2)) { + surf2_is_xyplane = true; + } else if (const auto* ptr = dynamic_cast(&surf2)) { + surf2_is_xyplane = true; + } else if (const auto* ptr = dynamic_cast(&surf2)) { + surf2_is_xyplane = false; + } else { + throw std::invalid_argument(fmt::format("Surface {} is an invalid type for " + "rotational periodic BCs. Only x-planes, y-planes, or general planes " + "(that are perpendicular to z) are supported for these BCs.", surf2.id_)); + } + + // Compute the surface normal vectors and make sure they are perpendicular + // to the z-axis + Direction norm1 = surf1.normal({0, 0, 0}); + Direction norm2 = surf2.normal({0, 0, 0}); + if (std::abs(norm1.z) > FP_PRECISION) { + throw std::invalid_argument(fmt::format("Rotational periodic BCs are only " + "supported for rotations about the z-axis, but surface {} is not " + "perpendicular to the z-axis.", surf1.id_)); + } + if (std::abs(norm2.z) > FP_PRECISION) { + throw std::invalid_argument(fmt::format("Rotational periodic BCs are only " + "supported for rotations about the z-axis, but surface {} is not " + "perpendicular to the z-axis.", surf2.id_)); + } + + // 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 angle between the two surfaces; this is the BC rotation angle + double theta1 = std::atan2(norm1.y, norm1.x); + double theta2 = std::atan2(norm2.y, norm2.x); + angle_ = theta2 - theta1; + + // 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)); + } + + // Guess whether or not the normal vectors of the two planes are aligned, i.e. + // if an arc passing from one surface through the geometry to the other + // surface will pass through the same "side" of both surfaces. If the user + // specified an x-plane and a y-plane then the geometry likely lies in the + // first quadrant which means the normals are not aligned. Otherwise, assume + // the opposite. + aligned_normals_ = !(surf1_is_xyplane && surf2_is_xyplane); } void @@ -167,20 +228,28 @@ RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const int i_particle_surf = std::abs(p.surface_) - 1; // Figure out which of the two BC surfaces were struck to figure out if a - // forward or backward rotation is required. Also specify the new surface. + // forward or backward rotation is required. Specify the other surface as + // the particle's new surface. double theta; int new_surface; if (i_particle_surf == i_surf_) { - theta = -angle_; - new_surface = p.surface_ > 0 ? -(j_surf_ + 1) : j_surf_ + 1; - } else if (i_particle_surf == j_surf_) { theta = angle_; - new_surface = p.surface_ > 0 ? -(i_surf_ + 1) : i_surf_ + 1; + new_surface = p.surface_ > 0 ? j_surf_ + 1 : -(j_surf_ + 1); + } else if (i_particle_surf == j_surf_) { + theta = -angle_; + new_surface = p.surface_ > 0 ? i_surf_ + 1 : -(i_surf_ + 1); } else { throw std::runtime_error("Called BoundaryCondition::handle_particle after " "hitting a surface, but that surface is not recognized by the BC."); } + // If the normal vectors of the two surfaces are not aligned, then the logic + // must be reversed for rotation and picking a new surface halfspace. + if (not aligned_normals_) { + theta = -theta; + new_surface = -new_surface; + } + // Rotate the particle's position and direction about the z-axis. Position r = p.r(); Direction u = p.u();