-
Notifications
You must be signed in to change notification settings - Fork 594
Fix a bug in rotational periodic boundary conditions #3692
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.
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
base: develop
Are you sure you want to change the base?
Conversation
paulromano
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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()|
In my testing, the normalization of the angle range to [-PI,PI] fix the warning problem. |
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 followed the style guidelines for Python source files (if applicable)I have made corresponding changes to the documentation (if applicable)