From bcd691a807f86fc8a1e5b282ab93fb8ee9bdef50 Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Wed, 1 Oct 2025 22:20:42 -0500 Subject: [PATCH 1/2] allow RotationalPeriodicBC around x,y, or z axis --- include/openmc/boundary_condition.h | 9 +++- src/boundary_condition.cpp | 84 ++++++++++++++--------------- 2 files changed, 47 insertions(+), 46 deletions(-) diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index af40131f1c8..0e7cc717679 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -138,18 +138,23 @@ class TranslationalPeriodicBC : public PeriodicBC { //============================================================================== //! 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 diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 7216ac89649..b3cee5673e8 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -158,56 +158,46 @@ void TranslationalPeriodicBC::handle_particle( // 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(&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_)); + // 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(&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 + // 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 " @@ -231,15 +221,7 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) 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)); @@ -251,6 +233,20 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) } } +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 { @@ -279,9 +275,9 @@ void RotationalPeriodicBC::handle_particle( 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); From 5dea22ba4b0fb7a04a1d98d87c410dbba2ad099d Mon Sep 17 00:00:00 2001 From: lewisgross1296 Date: Sat, 1 Nov 2025 17:36:22 -0500 Subject: [PATCH 2/2] added two periodic cylinders with 1/6th symmetry to regression tests. they should generate the same answers, need to run transport to get the results_true.dat --- .../periodic_cyls/__init__.py | 0 tests/regression_tests/periodic_cyls/test.py | 85 +++++++++++++++++++ .../periodic_cyls/xcyl_model/inputs_true.dat | 29 +++++++ .../periodic_cyls/xcyl_model/results_true.dat | 0 .../periodic_cyls/ycyl_model/inputs_true.dat | 29 +++++++ .../periodic_cyls/ycyl_model/results_true.dat | 0 6 files changed, 143 insertions(+) create mode 100644 tests/regression_tests/periodic_cyls/__init__.py create mode 100644 tests/regression_tests/periodic_cyls/test.py create mode 100644 tests/regression_tests/periodic_cyls/xcyl_model/inputs_true.dat create mode 100644 tests/regression_tests/periodic_cyls/xcyl_model/results_true.dat create mode 100644 tests/regression_tests/periodic_cyls/ycyl_model/inputs_true.dat create mode 100644 tests/regression_tests/periodic_cyls/ycyl_model/results_true.dat diff --git a/tests/regression_tests/periodic_cyls/__init__.py b/tests/regression_tests/periodic_cyls/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/periodic_cyls/test.py b/tests/regression_tests/periodic_cyls/test.py new file mode 100644 index 00000000000..e3bc293a320 --- /dev/null +++ b/tests/regression_tests/periodic_cyls/test.py @@ -0,0 +1,85 @@ +import openmc +import numpy +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def xcyl_model(): + model = openmc.Model() + # Define materials + fuel = openmc.Material() + fuel.add_nuclide('U235', 0.2) + fuel.add_nuclide('U238', 0.8) + fuel.set_density('g/cc', 19.1) + model.materials = openmc.Materials([fuel]) + + # Define geometry + # finite cylinder + x_min = openmc.XPlane(x0=0.0, boundary_type='reflective') + x_max = openmc.XPlane(x0=20.0, boundary_type='reflective') + x_cyl = openmc.XCylinder(r=20.0,boundary_type='vacuum') + # slice cylinder for periodic BC + periodic_bounding_yplane = openmc.YPlane(y0=0, boundary_type='periodic') + periodic_bounding_plane = openmc.Plane( + a=0.0, b=-np.sqrt(3) / 3, c=1, boundary_type='periodic', + ) + sixth_cyl_cell = openmc.Cell(1, fill=fuel, region = + +x_min &- x_max & -x_cyl & +periodic_bounding_yplane & +periodic_bounding_plane) + periodic_bounding_yplane.periodic_surface = periodic_bounding_plane + periodic_bounding_plane.periodic_surface = periodic_bounding_yplane + + model.geometry = openmc.Geometry([sixth_cyl_cell]) + + + # Define settings + model.settings.particles = 1000 + model.settings.batches = 4 + model.settings.inactive = 0 + model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( + (0, 0, 0), (20, 20, 20)) + ) + +@pytest.fixture +def ycyl_model(): + model = openmc.Model() + # Define materials + fuel = openmc.Material() + fuel.add_nuclide('U235', 0.2) + fuel.add_nuclide('U238', 0.8) + fuel.set_density('g/cc', 19.1) + model.materials = openmc.Materials([fuel]) + + # Define geometry + # finite cylinder + y_min = openmc.YPlane(y0=0.0, boundary_type='reflective') + y_max = openmc.YPlane(y0=20.0, boundary_type='reflective') + y_cyl = openmc.YCylinder(r=20.0,boundary_type='vacuum') + # slice cylinder for periodic BC + periodic_bounding_xplane = openmc.XPlane(x0=0, boundary_type='periodic') + periodic_bounding_plane = openmc.Plane( + a=-np.sqrt(3) / 3, b=0.0, c=1, boundary_type='periodic', + ) + sixth_cyl_cell = openmc.Cell(1, fill=fuel, region = + +y_min &- y_max & -y_cyl & +periodic_bounding_xplane & +periodic_bounding_plane) + periodic_bounding_xplane.periodic_surface = periodic_bounding_plane + periodic_bounding_plane.periodic_surface = periodic_bounding_xplane + model.geometry = openmc.Geometry([sixth_cyl_cell]) + + + # Define settings + model.settings.particles = 1000 + model.settings.batches = 4 + model.settings.inactive = 0 + model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( + (0, 0, 0), (20, 20, 20)) + ) + + +@pytest.mark.parametrize('cyl_model', ['xcyl_model','ycyl_model']) +def test_periodic(cyl_model): + with change_directory(cyl_model): + harness = PyAPITestHarness('statepoint.4.h5', cyl_model) + harness.main() + assert \ No newline at end of file diff --git a/tests/regression_tests/periodic_cyls/xcyl_model/inputs_true.dat b/tests/regression_tests/periodic_cyls/xcyl_model/inputs_true.dat new file mode 100644 index 00000000000..7d1ecf426a8 --- /dev/null +++ b/tests/regression_tests/periodic_cyls/xcyl_model/inputs_true.dat @@ -0,0 +1,29 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 4 + 0 + + + 0 0 0 20 20 20 + + + + diff --git a/tests/regression_tests/periodic_cyls/xcyl_model/results_true.dat b/tests/regression_tests/periodic_cyls/xcyl_model/results_true.dat new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/periodic_cyls/ycyl_model/inputs_true.dat b/tests/regression_tests/periodic_cyls/ycyl_model/inputs_true.dat new file mode 100644 index 00000000000..8ce785ed749 --- /dev/null +++ b/tests/regression_tests/periodic_cyls/ycyl_model/inputs_true.dat @@ -0,0 +1,29 @@ + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 4 + 0 + + + 0 0 0 20 20 20 + + + + diff --git a/tests/regression_tests/periodic_cyls/ycyl_model/results_true.dat b/tests/regression_tests/periodic_cyls/ycyl_model/results_true.dat new file mode 100644 index 00000000000..e69de29bb2d