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). |
|
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 Though, I wonder if 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 |
|
I don't know but you can install this branch and check if that solves the problem. |
|
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 Back to the original discourse forum post I made (re-sharing the image here) |
| 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) |
There was a problem hiding this comment.
Can you explain where the -1 comes from? This seems like a significant change. Tagging @gonuke who was also curious.
| axis_1_idx_ = 0; // x component independent | ||
| axis_2_idx_ = 2; // z component dependent | ||
| axis_1_idx_ = 2; // z component independent | ||
| axis_2_idx_ = 0; // x component dependent |
There was a problem hiding this comment.
@GuySten I think lines 176 - 183 were originally your suggestion in my PR (i.e. switching the independent and dependent axis for the y-case to make them confirm to the other cases). Can you explain why this is now being changed / how making this change does not break the assumptions of the rotation matrix code when computing new_r and new_u?
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)