openmc-dev
BoundingBox inconsistency between C++ and Python
Bug Description
I'm getting slightly different bounding boxes out of openmc.Cell.bounding_box and openmc.lib.cells[cell_id].bounding_box. It would be nice if these two methods returned the same thing.
Steps to Reproduce
I have attached a MWE below which is a small portion of an OpenMC model that was generated by the openmc_mcnp_adapter
import openmc
import openmc.lib
s19 = openmc.YCylinder(0, 0, 17.7)
s48 = openmc.Plane(0.7071067811865476, 6.123233995736766e-17, 0.7071067811865476, 11.45)
s58 = openmc.Plane(0.7071067811865476, 6.123233995736766e-17, 0.7071067811865476, 14.35)
s62 = openmc.Plane(6.123233995736766e-17, 1.0, 6.123233995736766e-17, 1.5999999999999999)
s64 = openmc.Plane(6.123233995736766e-17, 1.0, 6.123233995736766e-17, -1.3)
s65 = openmc.Plane(-0.7071067811865475, 6.123233995736766e-17, 0.7071067811865476, 1.45)
s68 = openmc.Plane(-0.7071067811865475, 6.123233995736766e-17, 0.7071067811865476, -1.45)
s101 = openmc.YPlane(5.6, boundary_type='vacuum')
region = +s64 & -s101 & -s19 & +s48 & (+s58 | +s62 | -s64 | +s65 | -s68)
cell = openmc.Cell(region=region)
univ = openmc.Universe(cells=[cell])
geometry = openmc.Geometry(root=univ)
settings = openmc.Settings(particles=100, batches=100)
model = openmc.model.Model(geometry=geometry, settings=settings)
model.export_to_model_xml()
openmc.lib.init()
print(cell.bounding_box)
print(openmc.lib.cells[1].bounding_box)
openmc.lib.finalize()
And here is the associated output:
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
############### %%%%%%%%%%%%%%%%%%%%%%%%
################## %%%%%%%%%%%%%%%%%%%%%%%
################### %%%%%%%%%%%%%%%%%%%%%%%
#################### %%%%%%%%%%%%%%%%%%%%%%
##################### %%%%%%%%%%%%%%%%%%%%%
###################### %%%%%%%%%%%%%%%%%%%%
####################### %%%%%%%%%%%%%%%%%%
####################### %%%%%%%%%%%%%%%%%
###################### %%%%%%%%%%%%%%%%%
#################### %%%%%%%%%%%%%%%%%
################# %%%%%%%%%%%%%%%%%
############### %%%%%%%%%%%%%%%%
############ %%%%%%%%%%%%%%%
######## %%%%%%%%%%%%%%
%%%%%%%%%%%
| The OpenMC Monte Carlo Code
Copyright | 2011-2023 MIT, UChicago Argonne LLC, and contributors
License | https://docs.openmc.org/en/latest/license.html
Version | 0.13.4-dev
Git SHA1 | cd13d3471da30213bf68b91739994b73bb07cb90
Date/Time | 2023-08-07 15:13:03
OpenMP Threads | 32
Reading model XML file './model.xml' ...
Reading cross sections XML file...
Minimum neutron data temperature: 1.7976931348623157e+308 K
Maximum neutron data temperature: 0 K
Preparing distributed cell instances...
Writing summary.h5 file...
BoundingBox(lower_left=(-17.7, -1.3, -17.7), upper_right=(17.7, 5.6, 17.7))
BoundingBox(lower_left=(-17.7, -inf, -17.7), upper_right=(17.7, 5.6, 17.7))
Environment
OpenMC version 0.13.4-dev on Ubuntu 20.04 on WSL
Just had a look at this, and it turns out it's due to some very clever Python coding by one @eepeterson: https://github.com/openmc-dev/openmc/blob/8e77a1ffb60de4f1f63dce319f2bc7fb818462a7/openmc/surface.py#L527-L542
Right now we don't have such an equivalent on the C++ side, where we only provide bounding boxes for axis-aligned planes (and otherwise default to infinite BBs). So, the Python side is smart enough to figure out that s64 above is basically a YPlane, but the C++ side can't.
@paulromano @eepeterson if no one else has been working on this I'd be happy to take it on
@caderache2014 I don't think either of us is going to tackle this in the near term so feel free! if you have any questions let me know.
@eepeterson @paulromano hi all. As far as either of you know, is this issue still active? (I know I said I'd take it on but shortly after my last exchange with @eepeterson I had a bicycle accident in Oxford and had to finish my internship with multiple ribs and my wrist broken)
Oh no!! So sorry to hear that @caderache2014! Hope you are recovering well at this point. As far as I'm aware, no one has been working on this so it's still up for grabs.
@paulromano thanks! I'll see what I can do with it. (And also, a belated thanks to @shimwell. I was suffered the bike accident during my morning commute and walked into the First Light office with a severe laceration on my arm; John personally helped me dress the wound using the office first aide kit ;-)
@paulromano and @eepeterson I have a question. for this bug, we're trying to get the C++ code to recognize a bounding box on cell object even if said boundaries are not axis aligned. It seems to me that we could add this logic to C++ in openmc/src/surface.cpp
https://github.com/openmc-dev/openmc/blob/24e1c951614d71fb7014a914615e55a697eef04c/src/surface.cpp#L287-L295
But that would require us to implement the logic in multiples types (e.g. SurfaceZPlane, SurfaceYPlane SurfaceXCylinder, etc)
Alternatively, we could implement this in the Cell class itself in openmc/src/cell.cpp
https://github.com/openmc-dev/openmc/blob/24e1c951614d71fb7014a914615e55a697eef04c/src/cell.cpp#L909C1-L953C1
BoundingBox Region::bounding_box(int32_t cell_id) const { if (simple_) { return bounding_box_simple(); } else { auto postfix = generate_postfix(cell_id); return bounding_box_complex(postfix); } }
//==============================================================================
BoundingBox Region::bounding_box_simple() const { BoundingBox bbox; for (int32_t token : expression_) { bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); } return bbox; }
//==============================================================================
BoundingBox Region::bounding_box_complex(vector<int32_t> postfix) const { vector<BoundingBox> stack(postfix.size()); int i_stack = -1;
for (auto& token : postfix) { if (token == OP_UNION) { stack[i_stack - 1] = stack[i_stack - 1] | stack[i_stack]; i_stack--; } else if (token == OP_INTERSECTION) { stack[i_stack - 1] = stack[i_stack - 1] & stack[i_stack]; i_stack--; } else { i_stack++; stack[i_stack] = model::surfaces[abs(token) - 1]->bounding_box(token > 0); } }
Ensures(i_stack == 0); return stack.front(); }
If either of you have any thoughts what makes the most sense, please let me know. Chris
Hi @paulromano (and @eepeterson ) I have been looking through the C++ code base for this bug and I had some questions about the class hierarchy.
I understand that bounding boxes play a significant role in these simulations. I also understand that determining a bounding box around region (i.e. a boolean composition of planes, half-spaces, spheres, etc) is challenging.
In the example shown above, @eepeterson implemented his clever logic in Python that can surmise that a region is actually YPlane and derive a bounding box appropriately. He implemented this functionality using a Mixin class for surfaces.
I'm definitely able to reproduce the inconsistency implementing the same code above. It is my understanding that calling bounding_box on openmc.lib.cell[cell_idx] triggers the C++ functionality in src/cell.cpp. However, looking through the path of execution on the C++ code, I'm not sure I see an appropriate point to introduce additional logic. Would we add functionality to the C++ Region class in src/surface.cpp? Or to the Surface class?
If I'm totally off base please let me know. I've tried modifying the bounding_box functions (in src/surface.cpp) for both Region and YSurface on my local branch and I don't think it's affecting the output of openmc.lib.cell[1].bounding_box. If the relevant C++ functionality is somewhere else, I'm still looking for it.
@caderache2014 What you would need to do is create a new SurfacePlane::bounding_box method that overrides the default Surface::bounding_box method (which just returns an infinite bounding box). The SurfacePlane::bounding_box method should have the same logic that PlaneMixin.bounding_box has on the Python side.