Fix a bug in rotational periodic boundary conditions

[GuySten]

Description

This PR simplify and fix logic of periodic boundary conditions.
It prevents particles getting lost when the periodic surfaces have normals in any direction.
I also added tests that check that the code run correctly regardless of the periodic plane normal direction.

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)

[paulromano]

Thanks for this fix and simplification! The one thing that still doesn't seem quite right is the warning about the angle not evenly subdividing 360. Depending on the orientation of the surfaces and whether positive/negative half-spaces are used to define the region, sometimes that warning is shown even when the angle is fine. To illustrate that, I've put together this example with two planes that are created using Plane.from_points, where each set of three points is randomly permuted. The region being plotted is the same in each case but some permutations of the points end up with the warning and some don't.

import openmc
from math import cos, sin, pi, radians
import numpy as np
import random

start = 0.
degrees = 45.
ang1 = radians(start)
ang2 = radians(start + degrees)

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)

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

(r1 & r2 & -zcyl).plot()

[GuySten]

In my testing, the normalization of the angle range to [-PI,PI] fix the warning problem.

[paulromano]

Thanks for the update. The update seems to have fixed the warning, but now cases that are 90 or 120 degrees rotationally periodic result in lost particles. I added a test to check for this. Can you look into this further?

[GuySten]

I had to compute the periodic surfaces senses to compute the angle using the canonical normal direction.

[lewisgross1296]

@GuySten thanks for this! Could you explain exactly what the problem was and what the solution was?

I think I might have experienced an error related to this. After using the new periodic BC, everything seemed to be working, but during one of my transport steps (as part of a Cardinal multiphysics simulation where OpenMC runs many times), the below happened

 Running OpenMC with 50000 particles per batch...
Executing filter editors...
Executing tally editors...
 Maximum neutron transport energy: 20000000 eV for He3
 Initializing source particles...

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

  Bat./Gen.      k            Average k
  =========   ========   ====================
        1/1    0.72088
...
       28/1    1.00152
       29/1    1.01043
       30/1   -nan    

This happened after four successful OpenMC runs, so it must be caused by a rare event. The error I got was

 WARNING: Particle 20017 underwent maximum number of events.
 ERROR: No fission sites banked on MPI rank 31
 ERROR: No fission sites banked on MPI rank 23
 ERROR: No fission sites banked on MPI rank 27
 ERROR: No fission sites banked on MPI rank 34
 ERROR: No fission sites banked on MPI rank 10
 ERROR: No fission sites banked on MPI rank 41
 ERROR: No fission sites banked on MPI rank 1
 ERROR: No fission sites banked on MPI rank 35
 ERROR: No fission sites banked on MPI rank 9
 ERROR: No fission sites banked on MPI rank 15
 ERROR: No fission sites banked on MPI rank 19
 ERROR: No fission sites banked on MPI rank 44
 ERROR: No fission sites banked on MPI rank 38
 ERROR: No fission sites banked on MPI rank 13
 ERROR: No fission sites banked on MPI rank 28
 ERROR: No fission sites banked on MPI rank 11
 ERROR: No fission sites banked on MPI rank 14
 ERROR: No fission sites banked on MPI rank 8
 ERROR: No fission sites banked on MPI rank 26
 ERROR: No fission sites banked on MPI rank 12
 ERROR: No fission sites banked on MPI rank 42
 ERROR: No fission sites banked on MPI rank 5
 ERROR: No fission sites banked on MPI rank 21
 ERROR: No fission sites banked on MPI rank 20
 ERROR: No fission sites banked on MPI rank 29
 ERROR: No fission sites banked on MPI rank 22
 ERROR: No fission sites banked on MPI rank 18
 ERROR: No fission sites banked on MPI rank 7
 ERROR: No fission sites banked on MPI rank 25
 ERROR: No fission sites banked on MPI rank 46
 ERROR: No fission sites banked on MPI rank 45
 ERROR: No fission sites banked on MPI rank 3
 ERROR: No fission sites banked on MPI rank 6
 ERROR: No fission sites banked on MPI rank 24
 ERROR: No fission sites banked on MPI rank 47
 ERROR: No fission sites banked on MPI rank 30
 ERROR: No fission sites banked on MPI rank 37
 ERROR: No fission sites banked on MPI rank 33
 ERROR: No fission sites banked on MPI rank 40
 ERROR: No fission sites banked on MPI rank 39
 ERROR: No fission sites banked on MPI rank 16
 ERROR: No fission sites banked on MPI rank 36
 ERROR: No fission sites banked on MPI rank 4
 ERROR: No fission sites banked on MPI rank 0
 ERROR: No fission sites banked on MPI rank 32
 ERROR: No fission sites banked on MPI rank 43
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 26 in communicator MPI_COMM_WORLD
with errorcode -1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
 ERROR: No fission sites banked on MPI rank 2
 ERROR: No fission sites banked on MPI rank 17

Does this seem similar to the problems you noticed in your testing?

[GuySten]

The only problems I encountered are missing particles.

The problem was that If you have two intersecting planes the angle between them is not well defined (there are two possible angles).
The angle between them is well defined only if they are halfspaces.
And you have to know how the halfspace is used in the geometry to get the correct orientation of the halfspace to compute the correct angle.
The rest is just math.

[lewisgross1296]

I see, yea in this case I do get the warning about 299.999 degrees not evenly dividing into 360, so in my case it's measuring the other angle.

Do you think the -nan I encountered would be fixed by the changes in this PR? I thought from the title of this PR and reply on discource that this PR might to address the -nan as well as the angle warning.

Though, I wonder if

 WARNING: Particle 20017 underwent maximum number of events.

might be related to an issue with the current implementation of PeriodicBCs. Perhaps it's a separate problem? And maybe it's causing the no fission sites banked because all the other successful runs had enough fission sites per MPI process

[GuySten]

I don't know but you can install this branch and check if that solves the problem.

[lewisgross1296]

Since my issue was happening within a Cardinal simulation, and after multiple OpenMC solves on a computationally expensive run, it's not super straightforward to just build Cardinal with your OpenMC change and test to see if it fixes it. (It might take order of days to find out if this worked). However, I'm going to do some testing to see if I can prove my issue is unrelated / a coincidence. From what I can tell from your comments/source code changes, this PR is only aimed to fix the warning, and is unrelated to my -nan issue.

Back to the original discourse forum post I made (re-sharing the image here)

Should we add a test that using flip_normal method on a surface involved with a periodic boundary condition does not change the eigenvalue? I think there may be something else strange with my model causing the -nan but at the very least, I don't think using the flip_normal method on a surface should change the actual eigenvalue as I observed.

[lewisgross1296]

Using openmc.Geometry.merge_surfaces = True ended up removing lines from my XML. I've now had one case get past the -nan where it originally failed, so some duplicated surfaces could have caused my problem. This suggests my -nan issue is specific to my model and shouldn't block anything here. Sorry, I originally couldn't tell if this PR was meant to fix both the angle and the nan from your response on discourse. Thanks for this contribution!

[paulromano]

Thanks for working through this bug @GuySten!