From 4c6dd2c78bfdee6c905c922b08cc458641134d74 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 18 Dec 2025 17:48:40 +0200 Subject: [PATCH 01/21] simplify rotational periodic boundary conditions --- src/boundary_condition.cpp | 25 ++++--------------------- src/surface.cpp | 2 +- 2 files changed, 5 insertions(+), 22 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 2840b3c7d52..313b4de7d10 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -242,30 +242,13 @@ double RotationalPeriodicBC::compute_periodic_rotation( void RotationalPeriodicBC::handle_particle( Particle& p, const Surface& surf) const { - int i_particle_surf = p.surface_index(); - - // Figure out which of the two BC surfaces were struck to figure out if a - // 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; - } else { - throw std::runtime_error( - "Called BoundaryCondition::handle_particle after " - "hitting a surface, but that surface is not recognized by the BC."); - } + int new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1; - // Rotate the particle's position and direction about the z-axis. + // Rotate the particle's position and direction. Position r = p.r(); Direction u = p.u(); - double cos_theta = std::cos(theta); - double sin_theta = std::sin(theta); + double cos_theta = std::cos(angle_); + double sin_theta = std::sin(angle_); Position new_r; new_r[zero_axis_idx_] = r[zero_axis_idx_]; diff --git a/src/surface.cpp b/src/surface.cpp index 09fc0f24f09..19c0915b32d 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1371,7 +1371,7 @@ void read_surfaces(pugi::xml_node node) "axis, which is not supported.")); } surf1.bc_ = make_unique(i_surf, j_surf, axis); - surf2.bc_ = make_unique(i_surf, j_surf, axis); + surf2.bc_ = make_unique(j_surf, i_surf, axis); } // If albedo data is present in albedo map, set the boundary albedo. From 1216d9da93fd2eaa1aedea3456d8ef0fed1c5862 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 19 Dec 2025 01:56:02 +0200 Subject: [PATCH 02/21] fix a bug in periodic rotational --- src/boundary_condition.cpp | 18 +++++-- tests/regression_tests/periodic_6fold/test.py | 51 ++++++++++++++----- 2 files changed, 52 insertions(+), 17 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 313b4de7d10..1ff4607b20d 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -192,6 +192,9 @@ RotationalPeriodicBC::RotationalPeriodicBC( 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; + // Compute the surface normal vectors and make sure they are perpendicular // to the correct axis Direction norm1 = surf1.normal({0, 0, 0}); @@ -212,9 +215,10 @@ RotationalPeriodicBC::RotationalPeriodicBC( surf2.id_)); } - angle_ = compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_], - norm2[axis_2_idx_], norm2[axis_1_idx_]); - + auto c = norm1.cross(norm2); + angle_ = std::atan2(c.dot(ax),norm1.dot(norm2)); + if (norm1.dot(norm2)<0.0) + angle_ += PI; // 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) { @@ -243,6 +247,14 @@ void RotationalPeriodicBC::handle_particle( Particle& p, const Surface& surf) const { int new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1; + + Surface& surf1 {*model::surfaces[i_surf_]}; + Surface& surf2 {*model::surfaces[j_surf_]}; + Direction norm1 = surf1.normal({0, 0, 0}); + Direction norm2 = surf2.normal({0, 0, 0}); + + if (norm1.dot(norm2)>0.0) + new_surface = -new_surface; // Rotate the particle's position and direction. Position r = p.r(); diff --git a/tests/regression_tests/periodic_6fold/test.py b/tests/regression_tests/periodic_6fold/test.py index 272c65df57f..7a72b3fad20 100644 --- a/tests/regression_tests/periodic_6fold/test.py +++ b/tests/regression_tests/periodic_6fold/test.py @@ -1,13 +1,15 @@ from math import sin, cos, pi import openmc +from openmc.utility_funcs import change_directory import pytest from tests.testing_harness import PyAPITestHarness -@pytest.fixture -def model(): +@pytest.mark.parametrize("flip1", [False, True]) +@pytest.mark.parametrize("flip2", [False, True]) +def test_periodic(flip1, flip2): model = openmc.model.Model() # Define materials @@ -27,17 +29,40 @@ def model(): # answers. theta1 = (-1/6 + 1/2) * pi theta2 = (1/6 - 1/2) * pi - plane1 = openmc.Plane(a=cos(theta1), b=sin(theta1), boundary_type='periodic') - plane2 = openmc.Plane(a=cos(theta2), b=sin(theta2), boundary_type='periodic') + if flip1: + plane1 = openmc.Plane(a=-cos(theta1), b=-sin(theta1), boundary_type='periodic') + else: + plane1 = openmc.Plane(a=cos(theta1), b=sin(theta1), boundary_type='periodic') + if flip2: + plane2 = openmc.Plane(a=-cos(theta2), b=-sin(theta2), boundary_type='periodic') + else: + plane2 = openmc.Plane(a=cos(theta2), b=sin(theta2), boundary_type='periodic') x_max = openmc.XPlane(5., boundary_type='reflective') z_cyl = openmc.ZCylinder(x0=3*cos(pi/6), y0=3*sin(pi/6), r=2.0) - - outside_cyl = openmc.Cell(1, fill=water, region=( - +plane1 & +plane2 & -x_max & +z_cyl)) - inside_cyl = openmc.Cell(2, fill=fuel, region=( - +plane1 & +plane2 & -z_cyl)) + + match (flip1,flip2): + case (False,False): + outside_cyl = openmc.Cell(1, fill=water, region=( + +plane1 & +plane2 & -x_max & +z_cyl)) + inside_cyl = openmc.Cell(2, fill=fuel, region=( + +plane1 & +plane2 & -z_cyl)) + case (False,True): + outside_cyl = openmc.Cell(1, fill=water, region=( + +plane1 & -plane2 & -x_max & +z_cyl)) + inside_cyl = openmc.Cell(2, fill=fuel, region=( + +plane1 & -plane2 & -z_cyl)) + case (True,False): + outside_cyl = openmc.Cell(1, fill=water, region=( + -plane1 & +plane2 & -x_max & +z_cyl)) + inside_cyl = openmc.Cell(2, fill=fuel, region=( + -plane1 & +plane2 & -z_cyl)) + case (True,True): + outside_cyl = openmc.Cell(1, fill=water, region=( + -plane1 & -plane2 & -x_max & +z_cyl)) + inside_cyl = openmc.Cell(2, fill=fuel, region=( + -plane1 & -plane2 & -z_cyl)) root_universe = openmc.Universe(0, cells=(outside_cyl, inside_cyl)) model.geometry = openmc.Geometry(root_universe) @@ -49,9 +74,7 @@ def model(): model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( (0, 0, 0), (5, 5, 0)) ) - return model - + with change_directory(f'{flip1}-{flip2}'): + harness = PyAPITestHarness('statepoint.4.h5', model) + harness.main() -def test_periodic(model): - harness = PyAPITestHarness('statepoint.4.h5', model) - harness.main() From 4bcfba2c10a2bde808806a34ffc9770792ef4fcd Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 19 Dec 2025 02:42:49 +0200 Subject: [PATCH 03/21] simplify further --- include/openmc/boundary_condition.h | 2 ++ src/boundary_condition.cpp | 25 +++---------------------- 2 files changed, 5 insertions(+), 22 deletions(-) diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index 5a14239e8ba..89bf7ca8bb3 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -152,6 +152,8 @@ class RotationalPeriodicBC : public PeriodicBC { protected: //! Angle about the axis by which particle coordinates will be rotated double angle_; + //! Do we need to flip surfaces senses when applying the transformation? + bool flip_sense_; //! 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 diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 1ff4607b20d..0b8b8b24f45 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -214,10 +214,11 @@ RotationalPeriodicBC::RotationalPeriodicBC( "intersect the origin.", surf2.id_)); } + flip_sense_ = (norm1.dot(norm2)>0.0); auto c = norm1.cross(norm2); angle_ = std::atan2(c.dot(ax),norm1.dot(norm2)); - if (norm1.dot(norm2)<0.0) + if (!flip_sense_) angle_ += PI; // Warn the user if the angle does not evenly divide a circle double rem = std::abs(std::remainder((2 * PI / angle_), 1.0)); @@ -229,31 +230,11 @@ RotationalPeriodicBC::RotationalPeriodicBC( } } -double RotationalPeriodicBC::compute_periodic_rotation( - double rise_1, double run_1, double rise_2, double 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 { int new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1; - - Surface& surf1 {*model::surfaces[i_surf_]}; - Surface& surf2 {*model::surfaces[j_surf_]}; - Direction norm1 = surf1.normal({0, 0, 0}); - Direction norm2 = surf2.normal({0, 0, 0}); - - if (norm1.dot(norm2)>0.0) + if (flip_sense_) new_surface = -new_surface; // Rotate the particle's position and direction. From 75d681020617ed05ee465657b369cec805097b1a Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 19 Dec 2025 02:44:45 +0200 Subject: [PATCH 04/21] added tests --- .../False-False/inputs_true.dat | 34 +++++++++++++++++++ .../False-False/results_true.dat | 2 ++ .../periodic_6fold/False-True/inputs_true.dat | 34 +++++++++++++++++++ .../False-True/results_true.dat | 2 ++ .../periodic_6fold/True-False/inputs_true.dat | 34 +++++++++++++++++++ .../True-False/results_true.dat | 2 ++ .../periodic_6fold/True-True/inputs_true.dat | 34 +++++++++++++++++++ .../periodic_6fold/True-True/results_true.dat | 2 ++ 8 files changed, 144 insertions(+) create mode 100644 tests/regression_tests/periodic_6fold/False-False/inputs_true.dat create mode 100644 tests/regression_tests/periodic_6fold/False-False/results_true.dat create mode 100644 tests/regression_tests/periodic_6fold/False-True/inputs_true.dat create mode 100644 tests/regression_tests/periodic_6fold/False-True/results_true.dat create mode 100644 tests/regression_tests/periodic_6fold/True-False/inputs_true.dat create mode 100644 tests/regression_tests/periodic_6fold/True-False/results_true.dat create mode 100644 tests/regression_tests/periodic_6fold/True-True/inputs_true.dat create mode 100644 tests/regression_tests/periodic_6fold/True-True/results_true.dat diff --git a/tests/regression_tests/periodic_6fold/False-False/inputs_true.dat b/tests/regression_tests/periodic_6fold/False-False/inputs_true.dat new file mode 100644 index 00000000000..075cffc123b --- /dev/null +++ b/tests/regression_tests/periodic_6fold/False-False/inputs_true.dat @@ -0,0 +1,34 @@ + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 4 + 0 + + + 0 0 0 5 5 0 + + + + diff --git a/tests/regression_tests/periodic_6fold/False-False/results_true.dat b/tests/regression_tests/periodic_6fold/False-False/results_true.dat new file mode 100644 index 00000000000..04aa308781d --- /dev/null +++ b/tests/regression_tests/periodic_6fold/False-False/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.848492E+00 2.933785E-03 diff --git a/tests/regression_tests/periodic_6fold/False-True/inputs_true.dat b/tests/regression_tests/periodic_6fold/False-True/inputs_true.dat new file mode 100644 index 00000000000..cbc1e91003e --- /dev/null +++ b/tests/regression_tests/periodic_6fold/False-True/inputs_true.dat @@ -0,0 +1,34 @@ + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 4 + 0 + + + 0 0 0 5 5 0 + + + + diff --git a/tests/regression_tests/periodic_6fold/False-True/results_true.dat b/tests/regression_tests/periodic_6fold/False-True/results_true.dat new file mode 100644 index 00000000000..04aa308781d --- /dev/null +++ b/tests/regression_tests/periodic_6fold/False-True/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.848492E+00 2.933785E-03 diff --git a/tests/regression_tests/periodic_6fold/True-False/inputs_true.dat b/tests/regression_tests/periodic_6fold/True-False/inputs_true.dat new file mode 100644 index 00000000000..675a8b3000f --- /dev/null +++ b/tests/regression_tests/periodic_6fold/True-False/inputs_true.dat @@ -0,0 +1,34 @@ + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 4 + 0 + + + 0 0 0 5 5 0 + + + + diff --git a/tests/regression_tests/periodic_6fold/True-False/results_true.dat b/tests/regression_tests/periodic_6fold/True-False/results_true.dat new file mode 100644 index 00000000000..04aa308781d --- /dev/null +++ b/tests/regression_tests/periodic_6fold/True-False/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.848492E+00 2.933785E-03 diff --git a/tests/regression_tests/periodic_6fold/True-True/inputs_true.dat b/tests/regression_tests/periodic_6fold/True-True/inputs_true.dat new file mode 100644 index 00000000000..e5ac03b4834 --- /dev/null +++ b/tests/regression_tests/periodic_6fold/True-True/inputs_true.dat @@ -0,0 +1,34 @@ + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 4 + 0 + + + 0 0 0 5 5 0 + + + + diff --git a/tests/regression_tests/periodic_6fold/True-True/results_true.dat b/tests/regression_tests/periodic_6fold/True-True/results_true.dat new file mode 100644 index 00000000000..04aa308781d --- /dev/null +++ b/tests/regression_tests/periodic_6fold/True-True/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.848492E+00 2.933785E-03 From f3e58f774be6448062c127903fdc93bab0e10067 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 19 Dec 2025 02:45:16 +0200 Subject: [PATCH 05/21] removed uneeded files --- .../periodic_6fold/inputs_true.dat | 34 ------------------- .../periodic_6fold/results_true.dat | 2 -- 2 files changed, 36 deletions(-) delete mode 100644 tests/regression_tests/periodic_6fold/inputs_true.dat delete mode 100644 tests/regression_tests/periodic_6fold/results_true.dat diff --git a/tests/regression_tests/periodic_6fold/inputs_true.dat b/tests/regression_tests/periodic_6fold/inputs_true.dat deleted file mode 100644 index 075cffc123b..00000000000 --- a/tests/regression_tests/periodic_6fold/inputs_true.dat +++ /dev/null @@ -1,34 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - eigenvalue - 1000 - 4 - 0 - - - 0 0 0 5 5 0 - - - - diff --git a/tests/regression_tests/periodic_6fold/results_true.dat b/tests/regression_tests/periodic_6fold/results_true.dat deleted file mode 100644 index 04aa308781d..00000000000 --- a/tests/regression_tests/periodic_6fold/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -1.848492E+00 2.933785E-03 From 6fb2e587e8bc204ca2ba8aec32a1925bc72f4d13 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 19 Dec 2025 02:49:21 +0200 Subject: [PATCH 06/21] ran clang-format --- src/boundary_condition.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 0b8b8b24f45..fd5a5f4b926 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -192,7 +192,7 @@ RotationalPeriodicBC::RotationalPeriodicBC( fmt::format("You've specified an axis that is not x, y, or z.")); } - Direction ax = {0.0,0.0,0.0}; + Direction ax = {0.0, 0.0, 0.0}; ax[zero_axis_idx_] = 1.0; // Compute the surface normal vectors and make sure they are perpendicular @@ -214,10 +214,10 @@ RotationalPeriodicBC::RotationalPeriodicBC( "intersect the origin.", surf2.id_)); } - flip_sense_ = (norm1.dot(norm2)>0.0); + flip_sense_ = (norm1.dot(norm2) > 0.0); auto c = norm1.cross(norm2); - angle_ = std::atan2(c.dot(ax),norm1.dot(norm2)); + angle_ = std::atan2(c.dot(ax), norm1.dot(norm2)); if (!flip_sense_) angle_ += PI; // Warn the user if the angle does not evenly divide a circle From 72c73a44a0f3f15f1c4a820e2239e3c03f76f59f Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 19 Dec 2025 10:44:50 +0200 Subject: [PATCH 07/21] fix y rotation case --- src/boundary_condition.cpp | 9 ++------- src/surface.cpp | 2 +- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index fd5a5f4b926..cf1e14614f1 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -173,14 +173,9 @@ RotationalPeriodicBC::RotationalPeriodicBC( axis_2_idx_ = 2; // z component dependent break; case y: - // for a right handed coordinate system, z should be the independent axis - // but this would cause the y-rotation case to be different than the other - // two. using a left handed coordinate system and a negative rotation the - // compute angle and rotation matrix behavior mimics that of the x and z - // cases zero_axis_idx_ = 1; // y component of plane must be zero - axis_1_idx_ = 0; // x component independent - axis_2_idx_ = 2; // z component dependent + axis_1_idx_ = 2; // x component independent + axis_2_idx_ = 0; // z component dependent break; case z: zero_axis_idx_ = 2; // z component of plane must be zero diff --git a/src/surface.cpp b/src/surface.cpp index 19c0915b32d..a7406caca1f 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1186,7 +1186,7 @@ void read_surfaces(pugi::xml_node node) pugi::xml_node surf_node; int i_surf; for (surf_node = node.child("surface"), i_surf = 0; surf_node; - surf_node = surf_node.next_sibling("surface"), i_surf++) { + surf_node = surf_node.next_sibling("surface"), i_surf++) { std::string surf_type = get_node_value(surf_node, "type", true, true); // Allocate and initialize the new surface From 8111c941c99ac71ed7fe66f98b9562284441910e Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 19 Dec 2025 10:56:56 +0200 Subject: [PATCH 08/21] ran clang-format --- src/surface.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/surface.cpp b/src/surface.cpp index a7406caca1f..19c0915b32d 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1186,7 +1186,7 @@ void read_surfaces(pugi::xml_node node) pugi::xml_node surf_node; int i_surf; for (surf_node = node.child("surface"), i_surf = 0; surf_node; - surf_node = surf_node.next_sibling("surface"), i_surf++) { + surf_node = surf_node.next_sibling("surface"), i_surf++) { std::string surf_type = get_node_value(surf_node, "type", true, true); // Allocate and initialize the new surface From 9404f8888d9d896aeda09ef16144e5232e871a81 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 23 Dec 2025 13:29:20 -0600 Subject: [PATCH 09/21] Fix one comment and add others --- src/boundary_condition.cpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index cf1e14614f1..c274025cd52 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -174,8 +174,8 @@ RotationalPeriodicBC::RotationalPeriodicBC( break; case y: zero_axis_idx_ = 1; // y component of plane must be zero - axis_1_idx_ = 2; // x component independent - axis_2_idx_ = 0; // z component dependent + 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 @@ -209,12 +209,19 @@ RotationalPeriodicBC::RotationalPeriodicBC( "intersect the origin.", surf2.id_)); } - flip_sense_ = (norm1.dot(norm2) > 0.0); + // 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_ = (norm1.dot(norm2) > 0.0); if (!flip_sense_) angle_ += PI; + // 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) { From 064b0c3b2694ef342079637334cd000dae53cb73 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 23 Dec 2025 23:56:01 +0200 Subject: [PATCH 10/21] normlize angle range --- src/boundary_condition.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index c274025cd52..dc0ae25cd5d 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -222,6 +222,9 @@ RotationalPeriodicBC::RotationalPeriodicBC( if (!flip_sense_) angle_ += PI; + // Normalize range of angle to [-PI,PI]. + angle_ = std::remainder(angle_, 2 * PI); + // 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) { From 25fd91ef808da43cdb8ff05acd38bdde333bde9e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 24 Dec 2025 02:28:13 -0600 Subject: [PATCH 11/21] Add test for rotationally periodic BCs --- tests/unit_tests/test_periodic_bc.py | 47 ++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 tests/unit_tests/test_periodic_bc.py diff --git a/tests/unit_tests/test_periodic_bc.py b/tests/unit_tests/test_periodic_bc.py new file mode 100644 index 00000000000..4cc126abc8e --- /dev/null +++ b/tests/unit_tests/test_periodic_bc.py @@ -0,0 +1,47 @@ +from math import cos, sin, radians +import random + +import openmc +import pytest + + +@pytest.mark.parametrize("angle", [30., 45., 60., 90., 120.]) +def test_rotational_periodic_bc(angle): + # Pick random starting angle + start = random.uniform(0., 360.) + degrees = angle + ang1 = radians(start) + ang2 = radians(start + degrees) + + # Define three points on each plane and then randomly shuffle them + p1_points = [(0., 0., 0.), (cos(ang1), sin(ang1), 0.), (0., 0., 1.)] + p2_points = [(0., 0., 0.), (cos(ang2), sin(ang2), 0.), (0., 0., 1.)] + random.shuffle(p1_points) + random.shuffle(p2_points) + + # Create periodic planes and a cylinder + p1 = openmc.Plane.from_points(*p1_points, boundary_type='periodic') + p2 = openmc.Plane.from_points(*p2_points, boundary_type='periodic') + p1.periodic_surface = p2 + zcyl = openmc.ZCylinder(r=5., boundary_type='vacuum') + + # Figure out which side of planes to use based on a point in the middle + ang_mid = radians(start + degrees/2.) + mid_point = (cos(ang_mid), sin(ang_mid), 0.) + r1 = -p1 if mid_point in -p1 else +p1 + r2 = -p2 if mid_point in -p2 else +p2 + + # Create one cell bounded by the two planes and the cylinder + mat = openmc.Material(density=1.0, density_units='g/cm3', components={'U235': 1.0}) + cell = openmc.Cell(fill=mat, region=r1 & r2 & -zcyl) + + # Make the model complete + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.settings.source = openmc.IndependentSource(space=openmc.stats.Point(mid_point)) + model.settings.particles = 1000 + model.settings.batches = 10 + model.settings.inactive = 5 + + # Run the model + model.run() From fc35e80183f08e3cc420fc8556363753922cef17 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 25 Dec 2025 09:38:15 +0200 Subject: [PATCH 12/21] fix code --- src/boundary_condition.cpp | 16 +++++++--------- src/surface.cpp | 32 ++++++++++++++++++++++++++++++-- 2 files changed, 37 insertions(+), 11 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index dc0ae25cd5d..eb095d8c587 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -160,7 +160,7 @@ void TranslationalPeriodicBC::handle_particle( RotationalPeriodicBC::RotationalPeriodicBC( int i_surf, int j_surf, PeriodicAxis axis) - : PeriodicBC(i_surf, j_surf) + : PeriodicBC(std::abs(i_surf) - 1, std::abs(j_surf) - 1) { Surface& surf1 {*model::surfaces[i_surf_]}; Surface& surf2 {*model::surfaces[j_surf_]}; @@ -190,10 +190,13 @@ RotationalPeriodicBC::RotationalPeriodicBC( 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 = surf1.normal({0, 0, 0}); - Direction norm2 = surf2.normal({0, 0, 0}); + 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( @@ -218,12 +221,7 @@ RotationalPeriodicBC::RotationalPeriodicBC( // If the normals point in the same general direction, the surface sense // should change when crossing the boundary - flip_sense_ = (norm1.dot(norm2) > 0.0); - if (!flip_sense_) - angle_ += PI; - - // Normalize range of angle to [-PI,PI]. - angle_ = std::remainder(angle_, 2 * PI); + 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)); diff --git a/src/surface.cpp b/src/surface.cpp index 19c0915b32d..646041014ca 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -9,6 +9,7 @@ #include #include "openmc/array.h" +#include "openmc/cell.h" #include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/external/quartic_solver.h" @@ -1182,6 +1183,7 @@ void read_surfaces(pugi::xml_node node) model::surfaces.reserve(n_surfaces); std::set> periodic_pairs; std::unordered_map albedo_map; + std::unordered_map periodic_sense_map; { pugi::xml_node surf_node; int i_surf; @@ -1244,6 +1246,7 @@ void read_surfaces(pugi::xml_node node) if (check_for_node(surf_node, "boundary")) { std::string surf_bc = get_node_value(surf_node, "boundary", true, true); if (surf_bc == "periodic") { + periodic_sense_map[model::surfaces.back()->id_] = 0; // Check for surface albedo. Skip sanity check as it is already done // in the Surface class's constructor. if (check_for_node(surf_node, "albedo")) { @@ -1276,6 +1279,27 @@ void read_surfaces(pugi::xml_node node) } } + // Fill the senses map for periodic surfaces + for (pugi::xml_node cell_node : node.children("cell")) { + + // Read the region specification. + std::string region_spec; + if (check_for_node(cell_node, "region")) { + region_spec = get_node_value(cell_node, "region"); + + // Get a tokenized representation of the region specification and apply De + // Morgans law + Region region(region_spec, 0); + + for (auto s : region.surfaces()) { + auto id = model::surfaces[std::abs(s) - 1]->id_; + if (periodic_sense_map.find(id) != periodic_sense_map.end()) { + periodic_sense_map[id] = std::copysign(1, s); + } + } + } + } + // Resolve unpaired periodic surfaces. A lambda function is used with // std::find_if to identify the unpaired surfaces. auto is_unresolved_pair = [](const std::pair p) { @@ -1370,8 +1394,12 @@ void read_surfaces(pugi::xml_node node) "indicates that the two planes are not periodic about the X, Y, or Z " "axis, which is not supported.")); } - surf1.bc_ = make_unique(i_surf, j_surf, axis); - surf2.bc_ = make_unique(j_surf, i_surf, axis); + auto i_sign = periodic_sense_map[periodic_pair.first]; + auto j_sign = periodic_sense_map[periodic_pair.second]; + surf1.bc_ = make_unique( + i_sign * (i_surf + 1), j_sign * (j_surf + 1), axis); + surf2.bc_ = make_unique( + j_sign * (j_surf + 1), i_sign * (i_surf + 1), axis); } // If albedo data is present in albedo map, set the boundary albedo. From a236ec3de3a093cc44f8fa2a794279e115570f20 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 25 Dec 2025 10:47:01 +0200 Subject: [PATCH 13/21] improved documentation --- src/particle.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 2d70d715e22..b9f9c8f86b6 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -744,9 +744,7 @@ void Particle::cross_periodic_bc( if (!neighbor_list_find_cell(*this)) { mark_as_lost("Couldn't find particle after hitting periodic " "boundary on surface " + - std::to_string(surf.id_) + - ". The normal vector " - "of one periodic surface may need to be reversed."); + std::to_string(surf.id_) + "."); return; } From ad6d723a094a186e421779c400966445402c7cb3 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 25 Dec 2025 11:07:01 +0200 Subject: [PATCH 14/21] add fast path for common case --- src/surface.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/surface.cpp b/src/surface.cpp index 646041014ca..8c5b4b8b1e4 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1183,6 +1183,7 @@ void read_surfaces(pugi::xml_node node) model::surfaces.reserve(n_surfaces); std::set> periodic_pairs; std::unordered_map albedo_map; + int n_periodic = 0; std::unordered_map periodic_sense_map; { pugi::xml_node surf_node; @@ -1280,21 +1281,22 @@ void read_surfaces(pugi::xml_node node) } // Fill the senses map for periodic surfaces - for (pugi::xml_node cell_node : node.children("cell")) { - + auto v = node.children("cell"); + auto n_periodic = periodic_sense_map.size(); + for (auto it = v.begin(); (it != v.end()) && (n_periodic > 0); ++it) { + pugi::xml_node cell_node = *it; // Read the region specification. std::string region_spec; if (check_for_node(cell_node, "region")) { region_spec = get_node_value(cell_node, "region"); - - // Get a tokenized representation of the region specification and apply De - // Morgans law + // Get a tokenized representation of the region specification and apply + // De Morgans law Region region(region_spec, 0); - for (auto s : region.surfaces()) { auto id = model::surfaces[std::abs(s) - 1]->id_; if (periodic_sense_map.find(id) != periodic_sense_map.end()) { periodic_sense_map[id] = std::copysign(1, s); + --n_periodic; } } } From 26ca7aee66bbf356116a0aac241a8cb418700559 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 25 Dec 2025 11:09:27 +0200 Subject: [PATCH 15/21] fix typo --- src/surface.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/surface.cpp b/src/surface.cpp index 8c5b4b8b1e4..a75407251bb 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1183,7 +1183,6 @@ void read_surfaces(pugi::xml_node node) model::surfaces.reserve(n_surfaces); std::set> periodic_pairs; std::unordered_map albedo_map; - int n_periodic = 0; std::unordered_map periodic_sense_map; { pugi::xml_node surf_node; From 86ed68a7139712f552aca610598ed74fbc8ec69b Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 31 Dec 2025 01:04:11 +0200 Subject: [PATCH 16/21] restructured code to read cell definitions once --- include/openmc/surface.h | 9 ++++++++- src/geometry_aux.cpp | 7 ++++++- src/surface.cpp | 37 +++++++++++++++++-------------------- 3 files changed, 31 insertions(+), 22 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 549e8f6abac..31c94bbe4b8 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -378,7 +378,14 @@ class SurfaceZTorus : public Surface { // Non-member functions //============================================================================== -void read_surfaces(pugi::xml_node node); +void read_surfaces(pugi::xml_node node, + std::set>& periodic_pairs, + std::unordered_map& albedo_map, + std::unordered_map& periodic_sense_map); + +void prepare_boundary_conditions(std::set>& periodic_pairs, + std::unordered_map& albedo_map, + std::unordered_map& periodic_sense_map); void free_memory_surfaces(); diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 8a145fb1f1e..a740740c1e6 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -55,8 +55,13 @@ void read_geometry_xml() void read_geometry_xml(pugi::xml_node root) { // Read surfaces, cells, lattice - read_surfaces(root); + std::set> periodic_pairs; + std::unordered_map albedo_map; + std::unordered_map periodic_sense_map; + + read_surfaces(root, periodic_pairs, albedo_map, periodic_sense_map); read_cells(root); + prepare_boundary_conditions(periodic_pairs, albedo_map, periodic_sense_map); read_lattices(root); // Check to make sure a boundary condition was applied to at least one diff --git a/src/surface.cpp b/src/surface.cpp index a75407251bb..0e276082d99 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1170,7 +1170,10 @@ Direction SurfaceZTorus::normal(Position r) const //============================================================================== -void read_surfaces(pugi::xml_node node) +void read_surfaces(pugi::xml_node node, + std::set>& periodic_pairs, + std::unordered_map& albedo_map, + std::unordered_map& periodic_sense_map) { // Count the number of surfaces int n_surfaces = 0; @@ -1181,9 +1184,6 @@ void read_surfaces(pugi::xml_node node) // Loop over XML surface elements and populate the array. Keep track of // periodic surfaces and their albedos. model::surfaces.reserve(n_surfaces); - std::set> periodic_pairs; - std::unordered_map albedo_map; - std::unordered_map periodic_sense_map; { pugi::xml_node surf_node; int i_surf; @@ -1278,25 +1278,22 @@ void read_surfaces(pugi::xml_node node) fmt::format("Two or more surfaces use the same unique ID: {}", id)); } } +} +void prepare_boundary_conditions(std::set>& periodic_pairs, + std::unordered_map& albedo_map, + std::unordered_map& periodic_sense_map) +{ // Fill the senses map for periodic surfaces - auto v = node.children("cell"); auto n_periodic = periodic_sense_map.size(); - for (auto it = v.begin(); (it != v.end()) && (n_periodic > 0); ++it) { - pugi::xml_node cell_node = *it; - // Read the region specification. - std::string region_spec; - if (check_for_node(cell_node, "region")) { - region_spec = get_node_value(cell_node, "region"); - // Get a tokenized representation of the region specification and apply - // De Morgans law - Region region(region_spec, 0); - for (auto s : region.surfaces()) { - auto id = model::surfaces[std::abs(s) - 1]->id_; - if (periodic_sense_map.find(id) != periodic_sense_map.end()) { - periodic_sense_map[id] = std::copysign(1, s); - --n_periodic; - } + for (auto it = model::cells.begin(); + (it != model::cells.end()) && (n_periodic > 0); ++it) { + auto cell = it*; + for (auto s : cell->region_) { + auto id = model::surfaces[std::abs(s) - 1]->id_; + if (periodic_sense_map.find(id) != periodic_sense_map.end()) { + periodic_sense_map[id] = std::copysign(1, s); + --n_periodic; } } } From f11391626af2a02d19f008b8070bd64fefc946b2 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 31 Dec 2025 01:08:06 +0200 Subject: [PATCH 17/21] fixed missing include --- include/openmc/surface.h | 1 + 1 file changed, 1 insertion(+) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 31c94bbe4b8..27a48387677 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -2,6 +2,7 @@ #define OPENMC_SURFACE_H #include // For numeric_limits +#include #include #include From 82880a61be571c69a7bb2a181874f79d264b158c Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 31 Dec 2025 01:14:08 +0200 Subject: [PATCH 18/21] fixed typo --- src/surface.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/surface.cpp b/src/surface.cpp index 0e276082d99..f8c4cc77c89 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1288,7 +1288,7 @@ void prepare_boundary_conditions(std::set>& periodic_pairs, auto n_periodic = periodic_sense_map.size(); for (auto it = model::cells.begin(); (it != model::cells.end()) && (n_periodic > 0); ++it) { - auto cell = it*; + auto cell = *it; for (auto s : cell->region_) { auto id = model::surfaces[std::abs(s) - 1]->id_; if (periodic_sense_map.find(id) != periodic_sense_map.end()) { From b33224ab8a16856d85fa3746954f131ecd713b03 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 31 Dec 2025 01:20:15 +0200 Subject: [PATCH 19/21] fixed another typo --- src/surface.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/surface.cpp b/src/surface.cpp index f8c4cc77c89..69c7818948c 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1288,8 +1288,8 @@ void prepare_boundary_conditions(std::set>& periodic_pairs, auto n_periodic = periodic_sense_map.size(); for (auto it = model::cells.begin(); (it != model::cells.end()) && (n_periodic > 0); ++it) { - auto cell = *it; - for (auto s : cell->region_) { + auto& cell = *it; + for (auto s : cell->surfaces()) { auto id = model::surfaces[std::abs(s) - 1]->id_; if (periodic_sense_map.find(id) != periodic_sense_map.end()) { periodic_sense_map[id] = std::copysign(1, s); From ebd5c8b76323a08f06e8a679233be0d75ebc1020 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 5 Jan 2026 23:49:36 -0600 Subject: [PATCH 20/21] Add comments in surface.h, simplify periodic_sense_map update --- include/openmc/surface.h | 25 +++++++++++++++++++++++++ src/surface.cpp | 13 ++++++++----- 2 files changed, 33 insertions(+), 5 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 27a48387677..2d8580345a4 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -379,11 +379,36 @@ class SurfaceZTorus : public Surface { // Non-member functions //============================================================================== +//! Read surface definitions from XML and populate the global surfaces vector. +//! +//! This function parses surface elements from the XML input, creates the +//! appropriate surface objects, and identifies periodic surfaces along with +//! their albedo values and sense information. +//! +//! \param node XML node containing surface definitions +//! \param[out] periodic_pairs Set of surface ID pairs representing periodic +//! boundary conditions +//! \param[out] albedo_map Map of surface IDs to albedo values for periodic +//! surfaces +//! \param[out] periodic_sense_map Map of surface IDs to their sense values +//! (used to determine orientation for periodic BCs) void read_surfaces(pugi::xml_node node, std::set>& periodic_pairs, std::unordered_map& albedo_map, std::unordered_map& periodic_sense_map); +//! Resolve periodic surface pairs and assign boundary conditions. +//! +//! This function completes the setup of periodic boundary conditions by +//! resolving unpaired periodic surfaces, determining whether each pair +//! represents translational or rotational periodicity based on surface +//! normals, and assigning the appropriate boundary condition objects. +//! +//! \param[inout] periodic_pairs Set of surface ID pairs representing periodic +//! boundary conditions; unpaired entries are resolved +//! \param albedo_map Map of surface IDs to albedo values for periodic surfaces +//! \param periodic_sense_map Map of surface IDs to their sense values (used to +//! determine orientation for periodic BCs) void prepare_boundary_conditions(std::set>& periodic_pairs, std::unordered_map& albedo_map, std::unordered_map& periodic_sense_map); diff --git a/src/surface.cpp b/src/surface.cpp index 69c7818948c..e1612fc00d3 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1286,12 +1286,15 @@ void prepare_boundary_conditions(std::set>& periodic_pairs, { // Fill the senses map for periodic surfaces auto n_periodic = periodic_sense_map.size(); - for (auto it = model::cells.begin(); - (it != model::cells.end()) && (n_periodic > 0); ++it) { - auto& cell = *it; + for (const auto& cell : model::cells) { + if (n_periodic == 0) + break; // Early exit once all periodic surfaces found + for (auto s : cell->surfaces()) { - auto id = model::surfaces[std::abs(s) - 1]->id_; - if (periodic_sense_map.find(id) != periodic_sense_map.end()) { + auto surf_idx = std::abs(s) - 1; + auto id = model::surfaces[surf_idx]->id_; + + if (periodic_sense_map.count(id)) { periodic_sense_map[id] = std::copysign(1, s); --n_periodic; } From e24a15d11d149794f31d1ab5d18c16e1ca4f90db Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 6 Jan 2026 00:20:23 -0600 Subject: [PATCH 21/21] Style fixes --- tests/regression_tests/periodic_6fold/test.py | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/tests/regression_tests/periodic_6fold/test.py b/tests/regression_tests/periodic_6fold/test.py index 7a72b3fad20..bded1c465df 100644 --- a/tests/regression_tests/periodic_6fold/test.py +++ b/tests/regression_tests/periodic_6fold/test.py @@ -10,7 +10,7 @@ @pytest.mark.parametrize("flip1", [False, True]) @pytest.mark.parametrize("flip2", [False, True]) def test_periodic(flip1, flip2): - model = openmc.model.Model() + model = openmc.Model() # Define materials water = openmc.Material() @@ -41,39 +41,39 @@ def test_periodic(flip1, flip2): x_max = openmc.XPlane(5., boundary_type='reflective') z_cyl = openmc.ZCylinder(x0=3*cos(pi/6), y0=3*sin(pi/6), r=2.0) - - match (flip1,flip2): - case (False,False): + + match (flip1, flip2): + case (False, False): outside_cyl = openmc.Cell(1, fill=water, region=( +plane1 & +plane2 & -x_max & +z_cyl)) inside_cyl = openmc.Cell(2, fill=fuel, region=( +plane1 & +plane2 & -z_cyl)) - case (False,True): + case (False, True): outside_cyl = openmc.Cell(1, fill=water, region=( +plane1 & -plane2 & -x_max & +z_cyl)) inside_cyl = openmc.Cell(2, fill=fuel, region=( - +plane1 & -plane2 & -z_cyl)) - case (True,False): + +plane1 & -plane2 & -z_cyl)) + case (True, False): outside_cyl = openmc.Cell(1, fill=water, region=( -plane1 & +plane2 & -x_max & +z_cyl)) inside_cyl = openmc.Cell(2, fill=fuel, region=( - -plane1 & +plane2 & -z_cyl)) - case (True,True): + -plane1 & +plane2 & -z_cyl)) + case (True, True): outside_cyl = openmc.Cell(1, fill=water, region=( -plane1 & -plane2 & -x_max & +z_cyl)) inside_cyl = openmc.Cell(2, fill=fuel, region=( - -plane1 & -plane2 & -z_cyl)) + -plane1 & -plane2 & -z_cyl)) root_universe = openmc.Universe(0, cells=(outside_cyl, inside_cyl)) model.geometry = openmc.Geometry(root_universe) # Define settings - model.settings = openmc.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), (5, 5, 0)) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box((0, 0, 0), (5, 5, 0)) ) + with change_directory(f'{flip1}-{flip2}'): harness = PyAPITestHarness('statepoint.4.h5', model) harness.main()