Fix a bug in rotational periodic boundary conditions #3692
Conversation
GuySten
changed the title
paulromano
left a comment
paulromano
left a comment
There was a problem hiding this comment.
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. |
paulromano
left a comment
paulromano
left a comment
There was a problem hiding this comment.
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?
|
I had to compute the periodic surfaces senses to compute the angle using the canonical normal direction. |
src/surface.cpp
Outdated
| 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; |
There was a problem hiding this comment.
I understand the need to get the sense information from cell region definitions, but the code you've added here will result in us iterating over all cells (and reading their regions) twice, which is not ideal. Instead, I would suggest deferring the creation of the periodic BC objects until after cells have been read so that we can get the sense information without having to redundantly parse the XML.
|
@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 This happened after four successful OpenMC runs, so it must be caused by a rare event. The error I got was Does this seem similar to the problems you noticed in your testing? |
|
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). |
3 participants
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)