Skip to content

Generalize RotationalPeriodicBC for X-, Y-, or Z-axis #3591

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 1 commit into
base: develop
Choose a base branch
from
Draft

Generalize RotationalPeriodicBC for X-, Y-, or Z-axis #3591

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 7 additions & 2 deletions include/openmc/boundary_condition.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
#ifndef OPENMC_BOUNDARY_CONDITION_H
#define OPENMC_BOUNDARY_CONDITION_H

Expand Down Expand Up @@ -138,18 +138,23 @@
//==============================================================================
//! A BC that rotates particles about a global axis.
//
//! Currently only rotations about the z-axis are supported.
//! Only rotations about the x, y, and z axes are supported.
//==============================================================================

class RotationalPeriodicBC : public PeriodicBC {
public:
RotationalPeriodicBC(int i_surf, int j_surf);

float compute_periodic_rotation(float rise_1, float run_1, float rise_2, float run_2) const;
void handle_particle(Particle& p, const Surface& surf) const override;

protected:
//! Angle about the axis by which particle coordinates will be rotated
double angle_;
enum PeriodicAxis {x, y, z};
//! Ensure that choice of axes is right handed. axis_1_idx_ corresponds to the independent axis and axis_2_idx_ corresponds to the dependent axis in the 2D plane perpendicular to the planes' axis of rotation
int zero_axis_idx;
int axis_1_idx_;
int axis_2_idx_;
};

} // namespace openmc
Expand Down
84 changes: 40 additions & 44 deletions src/boundary_condition.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
#include "openmc/boundary_condition.h"

#include <exception>
Expand Down Expand Up @@ -158,56 +158,46 @@
// RotationalPeriodicBC implementation
//==============================================================================

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

// Check the type of the first surface
bool surf1_is_xyplane;
if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
surf1_is_xyplane = true;
} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
surf1_is_xyplane = true;
} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&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_));
// 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."),axis)
}

// Check the type of the second surface
bool surf2_is_xyplane;
if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
surf2_is_xyplane = true;
} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
surf2_is_xyplane = true;
} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&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
// to the correct axis
Direction norm1 = surf1.normal({0, 0, 0});
Direction norm2 = surf2.normal({0, 0, 0});
if (std::abs(norm1.z) > FP_PRECISION) {
if (std::abs(norm1[zero_axis_idx]) > 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) {
if (std::abs(norm2[zero_axis_idx]) > FP_PRECISION) {
throw std::invalid_argument(fmt::format(
"Rotational periodic BCs are only "
"supported for rotations about the z-axis, but surface {} is not "
Expand All @@ -231,15 +221,7 @@
surf2.id_));
}

// Compute the BC rotation angle. Here it is assumed that both surface
// normal vectors point inwards---towards the valid geometry region.
// Consequently, the rotation angle is not the difference between the two
// normals, but is instead the difference between one normal and one
// anti-normal. (An incident ray on one surface must be an outgoing ray on
// the other surface after rotation hence the anti-normal.)
double theta1 = std::atan2(norm1.y, norm1.x);
double theta2 = std::atan2(norm2.y, norm2.x) + PI;
angle_ = theta2 - theta1;
angle_ = compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_], norm2[axis_2_idx_], norm2[axis_1_idx_])

// Warn the user if the angle does not evenly divide a circle
double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
Expand All @@ -251,6 +233,20 @@
}
}

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

void RotationalPeriodicBC::handle_particle(
Particle& p, const Surface& surf) const
{
Expand Down Expand Up @@ -279,9 +275,9 @@
double cos_theta = std::cos(theta);
double sin_theta = std::sin(theta);
Position new_r = {
cos_theta * r.x - sin_theta * r.y, sin_theta * r.x + cos_theta * r.y, r.z};
cos_theta * r[axis_1_idx_] - sin_theta * r[axis_2_idx_], sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_], r[zero_axis_idx]};
Direction new_u = {
cos_theta * u.x - sin_theta * u.y, sin_theta * u.x + cos_theta * u.y, u.z};
cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_], sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_], u[zero_axis_idx]};

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