Skip to content

Conversation

@lewisgross1296
Copy link
Contributor

Description

Currently only periodicity around the Z-axis is supported. This PR extends this capability for X and Y periodicity. Fixes #3559.

Checklist

  • I have performed a self-review of my own code
  • I have run clang-format (version 15) on any C++ source files (if applicable)
  • I have followed the style guidelines for Python source files (if applicable)
  • I have made corresponding changes to the documentation (if applicable)
  • I have added tests that prove my fix is effective or that my feature works (if applicable)

@lewisgross1296 lewisgross1296 marked this pull request as draft October 2, 2025 03:27
@lewisgross1296
Copy link
Contributor Author

I think this is getting close. I'm currently running into a problem when trying to generate the test files

 WARNING: Couldn't find particle after hitting periodic boundary on surface 4.
          The normal vector of one periodic surface may need to be reversed.

I think I understand conceptually what's going on, i.e. the particle is moved to the other periodic surface and continues in the same direction (to where there is no geometry) when it should continue in the opposite direction (inside the geometry). I'm trying to find the correct line to fix this, but I didn't think I changed this behavior from how it worked when there was only rotational periodicity around the z-axis.

@lewisgross1296 lewisgross1296 marked this pull request as ready for review November 6, 2025 16:33
Comment on lines 264 to 267
Position new_r = {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[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]};
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just realized that my "generalization" was based on the z-case that existed here before. The x and z cases should be similar but the y rotation matrix has the sign flipped, which I think explains why I'm seeing one test failure in this way. If we condition on zero_axis_idx I think I can get it right
rotation_matrices

@lewisgross1296
Copy link
Contributor Author

lewisgross1296 commented Nov 6, 2025

So I just pushed my attempt at accounting for the signs in the transformation of r->new_r and u->new_u, but I'm still failing the tests. I did some scaffolding and it seems like my transformations are not actually transforming anything at all and I'm confused why. For example, the ycyl_model test

                                %%%%%%%%%%%%%%%
                           %%%%%%%%%%%%%%%%%%%%%%%%
                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %%%%%%%%%%%%%%%%%%%%%%%%
                                     %%%%%%%%%%%%%%%%%%%%%%%%
                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%
                ##################     %%%%%%%%%%%%%%%%%%%%%%%
                ###################     %%%%%%%%%%%%%%%%%%%%%%%
                ####################     %%%%%%%%%%%%%%%%%%%%%%
                #####################     %%%%%%%%%%%%%%%%%%%%%
                ######################     %%%%%%%%%%%%%%%%%%%%
                #######################     %%%%%%%%%%%%%%%%%%
                 #######################     %%%%%%%%%%%%%%%%%
                 ######################     %%%%%%%%%%%%%%%%%
                  ####################     %%%%%%%%%%%%%%%%%
                    #################     %%%%%%%%%%%%%%%%%
                     ###############     %%%%%%%%%%%%%%%%
                       ############     %%%%%%%%%%%%%%%
                          ########     %%%%%%%%%%%%%%
                                      %%%%%%%%%%%

                 | The OpenMC Monte Carlo Code
       Copyright | 2011-2025 MIT, UChicago Argonne LLC, and contributors
         License | https://docs.openmc.org/en/latest/license.html
         Version | 0.13.1-dev2826
     Commit Hash | dbd1a421d063fa78ff93e9c5e6995f84f7177a06
       Date/Time | 2025-11-06 19:45:53
   MPI Processes | 1
  OpenMP Threads | 1

 Reading model XML file 'model.xml' ...
 Reading cross sections XML file...
 Reading U235 from /root/nndc_hdf5/U235.h5
 Reading U238 from /root/nndc_hdf5/U238.h5
 Minimum neutron data temperature: 294 K
 Maximum neutron data temperature: 294 K
 Preparing distributed cell instances...
 Writing summary.h5 file...
 Maximum neutron transport energy: 20000000 eV for U235
 Initializing source particles...

 ====================>     K EIGENVALUE SIMULATION     <====================

  Bat./Gen.      k            Average k
  =========   ========   ====================
r = (2.58686, 2.82336, 1.49352)
new r =(2.58686, 2.82336, 1.49352)
u = (-0.454864, 0.605916, -0.65266)
new u= (-0.454863, 0.605916, -0.65266)
 WARNING: Couldn't find particle after hitting periodic boundary on surface 5.
          The normal vector of one periodic surface may need to be reversed.
 ERROR: Maximum number of lost particles has been reached.
Abort(-1) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0

Am I missing some C++ thing about scope or something like that?

  // rotations around the y-axis have sign flipped for the sin_theta terms
  int flip;
  if (zero_axis_idx_ == 1) {
    flip = -1;
  } else {
    flip = 1;
  }

  Position new_r;
  new_r[zero_axis_idx_] = r[zero_axis_idx_];
  new_r[axis_1_idx_] =
    cos_theta * r[axis_1_idx_] - flip * sin_theta * r[axis_2_idx_];
  new_r[axis_2_idx_] =
    flip * sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_];

  Direction new_u;
  new_u[zero_axis_idx_] = u[zero_axis_idx_];
  new_u[axis_1_idx_] =
    cos_theta * u[axis_1_idx_] - flip * sin_theta * u[axis_2_idx_];
  new_u[axis_2_idx_] =
    flip * sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];

EDIT: The above problem derived from the fact that we didn't actually apply the x or y cases in the constructor and so we need to add detection from the two periodic planes to determine which case we are (x,y, or z periodic). Essentially, it was computing a periodic angle of 2pi in the above case

@lewisgross1296 lewisgross1296 force-pushed the generalize-periodic-bc branch 2 times, most recently from 9febb99 to 17ec3eb Compare November 19, 2025 20:33
@lewisgross1296
Copy link
Contributor Author

I'm now passing the xcyl test but not the ycyl.

                                %%%%%%%%%%%%%%%
                           %%%%%%%%%%%%%%%%%%%%%%%%
                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %%%%%%%%%%%%%%%%%%%%%%%%
                                     %%%%%%%%%%%%%%%%%%%%%%%%
                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%
                ##################     %%%%%%%%%%%%%%%%%%%%%%%
                ###################     %%%%%%%%%%%%%%%%%%%%%%%
                ####################     %%%%%%%%%%%%%%%%%%%%%%
                #####################     %%%%%%%%%%%%%%%%%%%%%
                ######################     %%%%%%%%%%%%%%%%%%%%
                #######################     %%%%%%%%%%%%%%%%%%
                 #######################     %%%%%%%%%%%%%%%%%
                 ######################     %%%%%%%%%%%%%%%%%
                  ####################     %%%%%%%%%%%%%%%%%
                    #################     %%%%%%%%%%%%%%%%%
                     ###############     %%%%%%%%%%%%%%%%
                       ############     %%%%%%%%%%%%%%%
                          ########     %%%%%%%%%%%%%%
                                      %%%%%%%%%%%

                 | The OpenMC Monte Carlo Code
       Copyright | 2011-2025 MIT, UChicago Argonne LLC, and contributors
         License | https://docs.openmc.org/en/latest/license.html
         Version | 0.13.1-dev2842
     Commit Hash | 17ec3ebdcb532b32e67010edc794e35c30376a66
       Date/Time | 2025-11-19 21:03:08
   MPI Processes | 1
  OpenMP Threads | 1

 Reading model XML file 'model.xml' ...
 Reading cross sections XML file...
angle after call =1.0472
angle after call =1.0472
 Reading U235 from /root/nndc_hdf5/U235.h5
 Reading U238 from /root/nndc_hdf5/U238.h5
 Minimum neutron data temperature: 294 K
 Maximum neutron data temperature: 294 K
 Preparing distributed cell instances...
 Writing summary.h5 file...
 Maximum neutron transport energy: 20000000 eV for U235
 Initializing source particles...

 ====================>     K EIGENVALUE SIMULATION     <====================

  Bat./Gen.      k            Average k
  =========   ========   ====================
r = (2.58686, 2.82336, 1.49352)
new r =(2.58686, 2.82336, -1.49352)
u = (-0.454864, 0.605916, -0.65266)
new u= (-0.792652, 0.605916, 0.0675933)
 WARNING: Couldn't find particle after hitting periodic boundary on surface 10.
          The normal vector of one periodic surface may need to be reversed.

I did some printing of details, and it looks like the rotation angle currently being detected is 60 degrees (1.0472 radians). When I do the rotation matrix by hand, new_r is what I would expect, but it is outside my geometry. If the angle is -60 instead, it would move the particle to the other surface location like I'd want, see diagram below
geom_diagram
I tried the openmc.Plane.flip_normal(), but what happens instead is that it detects 240 degrees and still translates particles to a LOST PARTICLE ZONE.

I'm thinking this is another consequence of the y-rotation matrix / the orientation of the axes and how it corresponds with our diagrams. See below
differing_coordinate_axes
In the top two cases, the third axis goes out of the page, but in the bottom case, the y-axis goes into the page. Since all three cases compute angles the same way, but the axes are flipped for the y-case, there needs to be another negative to get the right rotation angle between two planes for a RotationalPeriodicBC::PeriodicAxis axis = y. I'm going to test this out and then report back

break;
case y:
zero_axis_idx_ = 1; // y component of plane must be zero
axis_1_idx_ = 2; // z component independent
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you can simplify this PR significantly by switching axis_1_idx , axis_2_idx here.
This switch should make every special case because of y superfluous.

Copy link
Contributor Author

@lewisgross1296 lewisgross1296 Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I get what you mean, i.e. if we switch x and z, then it will behave like the others as long as we're careful about everything. I think it's a little less readable that way and would need a lengthy comment to explain. I kind of like it as is, though. It makes the most physical sense to me, and the math works that way (the y being a special case in the rotation matrix and the compute angle). I'm curious if others have feelings about leaving it or using your clever trick

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Generalize RotationalPeriodicBC to allow periodicity around x and y axes

2 participants