Skip to content
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

Cylindrical prism with optional filleted edges #2309

Merged
merged 26 commits into from
Jul 18, 2023

Conversation

yardasol
Copy link
Contributor

@yardasol yardasol commented Nov 23, 2022

This PR adds a new function, cylindrical_prism, that creates a finite cylindrical prism region with the option to fillet the cylinder's edges

@yardasol yardasol changed the title Rounded cyl corners Cylindrical prism with rounded edges option Nov 23, 2022
Copy link
Contributor

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like a useful addition @yardasol. I have a few suggestions related mostly to documentation and code reuse/maintainability.

openmc/model/funcs.py Outdated Show resolved Hide resolved
boundary_type='transmission',
upper_edge_radius=0.,
lower_edge_radius=0.):
"""Get a finite cylindrical prism.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps more information here about the rounded edges. From reviewing the code, it seems that this option adds what I would call a fillet at each end.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1 for fillet terminology

openmc/model/funcs.py Outdated Show resolved Hide resolved
Comment on lines 213 to 235
args['a'] = r - upper_edge_radius
args['b'] = upper_edge_radius
args['c'] = upper_edge_radius
args[axis + '0'] = origin[axcoord] + height/2 - upper_edge_radius
tor_upper_max = tor(name='{} max'.format(axis), **args)

axis_upper_min = plane(axis, 'upper min', height/2 + origin[axcoord] - upper_edge_radius)
cyl_upper_min = cylinder(axis, 'upper min', x1, origin[axcoord1], x2, origin[axcoord2], r - upper_edge_radius)

corners_upper = (+cyl_upper_min & +tor_upper_max & + axis_upper_min)

if lower_edge_radius > 0.:
args['a'] = r - lower_edge_radius
args['b'] = lower_edge_radius
args['c'] = lower_edge_radius
args[axis + '0'] = origin[axcoord] - height/2 + lower_edge_radius
tor_lower_min = tor(name='{} min'.format(axis), **args)

axis_lower_max = plane(axis, 'lower max', origin[axcoord] - height/2 + lower_edge_radius)
cyl_lower_min = cylinder(axis, 'lower min', x1, origin[axcoord1], x2, origin[axcoord2], r - lower_edge_radius)

corners_lower = (+cyl_lower_min & +tor_lower_min & - axis_lower_max)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems as though there is an opportunity to reuse at least some of this code in a modular way (e.g. a function or loop) rather than cut/paste for the upper/lower. The benefit of this kind of reuse is for consistency/maintainability and not for performance.

Copy link
Contributor

@paulromano paulromano left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this contribution @yardasol! In addition to the line comments, can you also add this function to the documentation in docs/source/pythonapi/model.rst?

openmc/model/funcs.py Outdated Show resolved Hide resolved
openmc/model/funcs.py Outdated Show resolved Hide resolved
Comment on lines 160 to 163
cls = globals()['{}Plane'.format(axis.upper())]
return cls(name='{} {}'.format(name, axis),
boundary_type=boundary_type,
**{axis + '0': value})
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The use of globals() feels a little hacky. I'd rather we import openmc and then get the object from that. Also, if you put the value of the axis first when creating the plane, it will simplify it a bit:

Suggested change
cls = globals()['{}Plane'.format(axis.upper())]
return cls(name='{} {}'.format(name, axis),
boundary_type=boundary_type,
**{axis + '0': value})
cls = getattr(openmc, f'{axis.upper()}Plane')
return cls(value, name=f'{name} {axis}',
boundary_type=boundary_type)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The globals() approach is used in rectangular_prism and hexagonal_prism so these should probably get changed as well.

openmc/model/funcs.py Outdated Show resolved Hide resolved
openmc/model/funcs.py Outdated Show resolved Hide resolved
openmc/model/funcs.py Outdated Show resolved Hide resolved
openmc/model/funcs.py Outdated Show resolved Hide resolved
openmc/model/funcs.py Outdated Show resolved Hide resolved
openmc/model/funcs.py Outdated Show resolved Hide resolved
openmc/model/funcs.py Show resolved Hide resolved
Copy link
Contributor

@pshriwise pshriwise left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a few smalll comments from me for now, @yardasol. Sorry it's taken me a bit to get to.

One thought that came to mind is the case where the fillet radius is larger than the cylinder radius. This would result in a degenerate torus (a minor radius larger than the major radius). I think we can do transport on those (@paulromano would know for sure), but I'm curious about your thoughts on what the behavior of this function should be in that case. My initial feeling is that if the corner radius is greater than the cylinder radius, we should limit the fillet radius to the radius of the cylinder and display a warning.

openmc/model/funcs.py Outdated Show resolved Hide resolved
5,
origin=(-1,1,3))
# clear checks
assert (-1, 1, 3) in cyl_prism # center
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be beneficial to use the openmc.lib.find_cell method to check that the C++ code is also getting the correct result too since it's what will affect transport.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great idea!

boundary_type='transmission',
upper_edge_radius=0.,
lower_edge_radius=0.):
"""Get a finite cylindrical prism.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1 for fillet terminology

@yardasol
Copy link
Contributor Author

yardasol commented Dec 8, 2022

Thanks for the reviews everyone! I believe I have addressed everyone's comments.

@pshriwise I added the cell check via openmc.lib like you asked, but I don't think it's behaving correctly. The rest of the test_geometry.py test seems to get skipped after test_cylindrical_prism(). I'll try to take a look at it with a profiler but I'd also like to know your thoughts on my implementation and also on this behavior.

One thought that came to mind is the case where the fillet radius is larger than the cylinder radius. This would result in a degenerate torus (a minor radius larger than the major radius). I think we can do transport on those (@paulromano would know for sure), but I'm curious about your thoughts on what the behavior of this function should be in that case. My initial feeling is that if the corner radius is greater than the cylinder radius, we should limit the fillet radius to the radius of the cylinder and display a warning.

I believe that I enforce the fillet radii must be less than the cylinder radius using the check_less_than() function near the beginning of the function

@yardasol yardasol changed the title Cylindrical prism with rounded edges option Cylindrical prism with optional filled edges Dec 8, 2022
@pshriwise
Copy link
Contributor

I believe that I enforce the fillet radii must be less than the cylinder radius using the check_less_than() function near the beginning of the function

Ah, yeah so you do. Missed it. Sorry about that.

(-2.4, 1, 0.51)]]

openmc.lib.init()
for prism, inside_points, outside_points in zip(prisms, inside_points, outside_points):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are a few things to address here:

  • When openmc.lib.init is called, it will read the XML files and initialize OpenMC in memory, so any subsequent XML files that are generated will be ignored. The init and finalize calls should be inside of the loop so OpenMC is started and shutdown for each prism being tested here. I'd even look into the openmc.lib.run_in_memory() context manager too to simplify this further.
  • The name cell is going to be undefined in the loop here. Might be good to return the cell created from the make_geometry function.
  • Going to need a boundary condition on at least one surface in the geometry.
  • Need to wrap the string passed to match in pytest.raises in an re.escape call to ask the regex matching not to treat the parentheses characters as a pattern. I've been hit by that one before.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, got it working!

Copy link
Contributor

@gonuke gonuke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few more ideas for streamlining the code. I confess some are based on speculation of how things work and incomplete expertise on my part, but looked like opportunities for simplifying things.

# Handle rounded corners if given
if upper_fillet_radius > 0. or lower_fillet_radius > 0.:
# Get torus class corresponding to given axis
tor = globals()['{}Torus'.format(axis.upper())]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Note: this is outside my expertise, but based on design patterns....) don't you want to access this the same way you access planes and cylinders through the OpenMC object?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, can this be moved into your create_fillet function?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, and yes!

cyl_name = 'upper_min'
axial_bound = +plane(axis, 'upper min', coord)
else:
coord = origin[axcoord] - height / 2 + lower_fillet_radius
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
coord = origin[axcoord] - height / 2 + lower_fillet_radius
coord = origin[axcoord] - height / 2 + fillet_radius

openmc/model/funcs.py Outdated Show resolved Hide resolved
Comment on lines 248 to 255
if fillet_lower is not None and fillet_upper is not None:
fillet = fillet_lower | fillet_upper
elif fillet_lower is None:
fillet = fillet_upper
else:
fillet = fillet_lower

prism = prism & ~fillet
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Speculation?) If neither fillet is defined then fillet will be None in this last operation. Is that geometry operation well defined if fillet is None?

If it is.... can you just use the fillet_upper and fillet_lower without the above logic to see if they are defined?

Suggested change
if fillet_lower is not None and fillet_upper is not None:
fillet = fillet_lower | fillet_upper
elif fillet_lower is None:
fillet = fillet_upper
else:
fillet = fillet_lower
prism = prism & ~fillet
prism = prism & ~(fillet_upper | fillet_lower)

Are all those operation well defined if one or both are None?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a great question. My guess is that it doesn't work that way based on the documentation for the | operator:

class Union(Region, MutableSequence):
     r"""Union of two or more regions.
 
     Instances of Union are generally created via the | operator applied to two
     instances of :class:`openmc.Region`. This is illustrated in the following
     example:
 
     >>> s1 = openmc.ZPlane(z0=0.0)
     >>> s2 = openmc.Sphere(r=637.1e6)
     >>> type(-s2 | +s1)
     <class 'openmc.region.Union'>
 
     Instances of this class behave like a mutable sequence, e.g., they can be
     indexed and have an append() method.
 
     Parameters
     ----------
     nodes : iterable of openmc.Region

but maybe @paulromano or @pshriwise can verify?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If neither fillet is defined then fillet will be None in this last operation. Is that geometry operation well defined if fillet is None?

At least one fillet is guaranteed to exist because this code block is contained within an if statement that requires at least one of uppet_fillet_radius and lower_fillet_radius to be nonzero.

@yardasol yardasol requested a review from gonuke December 9, 2022 21:06
@paulromano
Copy link
Contributor

@yardasol Now I'm realizing that we already have openmc.model.RightCircularCylinder, which provides the same thing but without fillets. What do you think about just extending that class to support fillets rather than having two different ways to create cylindrical prisms?

@yardasol yardasol changed the title Cylindrical prism with optional filled edges Cylindrical prism with optional filleted edges Jan 9, 2023
@yardasol
Copy link
Contributor Author

yardasol commented Jan 9, 2023

@yardasol Now I'm realizing that we already have openmc.model.RightCircularCylinder, which provides the same thing but without fillets. What do you think about just extending that class to support fillets rather than having two different ways to create cylindrical prisms?

This is a good idea. I've been busy doing other things, which is why I added the "on hold" tag, but I'll be able to get back to this soon (sometime this week most likely) after my thesis work is done

@yardasol
Copy link
Contributor Author

yardasol commented Jun 11, 2023

Finally got around to dealing with this one. @paulromano @gonuke @pshriwise

docs/source/pythonapi/model.rst Outdated Show resolved Hide resolved
openmc/model/surface_composite.py Outdated Show resolved Hide resolved
openmc/model/surface_composite.py Outdated Show resolved Hide resolved
openmc/model/surface_composite.py Outdated Show resolved Hide resolved
@yardasol
Copy link
Contributor Author

@paulromano bump

Copy link
Contributor

@paulromano paulromano left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for the long delay @yardasol! Thanks for refactoring this into the RightCircularCylinder class -- it's looking really good

openmc/model/funcs.py Outdated Show resolved Hide resolved
openmc/model/surface_composite.py Outdated Show resolved Hide resolved
openmc/model/surface_composite.py Outdated Show resolved Hide resolved
openmc/model/surface_composite.py Outdated Show resolved Hide resolved
openmc/model/surface_composite.py Outdated Show resolved Hide resolved
tests/unit_tests/test_geometry.py Outdated Show resolved Hide resolved
Copy link
Contributor

@paulromano paulromano left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the quick updates @yardasol. Good to go!

@paulromano paulromano enabled auto-merge (squash) July 14, 2023 17:40
@paulromano
Copy link
Contributor

@pshriwise do you have any further requests on this PR?

@pshriwise
Copy link
Contributor

@pshriwise do you have any further requests on this PR?

No further requests from me 🚀

@paulromano
Copy link
Contributor

@pshriwise Thanks! need you to hit the approve button for this one to land

Copy link
Contributor

@pshriwise pshriwise left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @yardasol!

@paulromano paulromano merged commit 922b3f9 into openmc-dev:develop Jul 18, 2023
stchaker pushed a commit to stchaker/openmc that referenced this pull request Oct 25, 2023
Co-authored-by: Paul Romano <[email protected]>
Co-authored-by: Patrick Shriwise <[email protected]>
Co-authored-by: Paul Wilson <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants