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.

Sign up for GitHub

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

Jump to bottom

Lattice corner crossings leak particles #2445

Open
gridley opened this issue Mar 28, 2023 · 2 comments · May be fixed by #3355
Open

Lattice corner crossings leak particles #2445

gridley opened this issue Mar 28, 2023 · 2 comments · May be fixed by #3355

Comments

@gridley
Copy link
Contributor

gridley commented Mar 28, 2023 •
edited
Loading

Hey everyone,

The attached input (geometry shown) exhibits an incorrect handling of lattice crossing in OpenMC. If crossings happen sufficiently near the corner, the incorrect cell is calculated as containing the particle surrounded by a cylindrical boundary. This input is a 2x2 lattice of two pieces of iron and air. The particles are sourced in the lower left lattice cell, and cross the boundary at an angle phi.

(purple = iron, yellow = airlike substance; the arrow can be adjusted by changing phi in the input file)
plot_1

There should be practically zero attenuation if the cell is calculated correctly. By setting the variable "offset" (found below) to 1e-6, I get the follow tracklength-estimated flux tally:

no_atten

There is clearly no attenuation here, as expected. However, if the particle is in the iron region, some attenuation could be expected. By shifting the source offset to zero and crossing the lattice boundary within floating point precision, the following can be obtained:

atten

That's clearly incorrect! Of course, many "negative distance to lattice boundary" warnings come up here, because the tracking is completely thrown off at this point. The spatial coordinates and lattice coordinates simply do not agree.

Crossing a lattice at a corner is a rare event, true, but may occur in simulations employing large numbers of histories. It would be nice to make the geometry handling robust. I have encountered this scenario many times in the projection plots PR, because it's not uncommon to want to view some geometry at 45 degrees. If a centered lattice is present, it creates exactly this scenario.

It's hard for me to say what the best approach to fix this is. Should our lattice crossing code specifically handle corner crossings? These are kind of fundamentally numerically tricky. Alternatively, if a negative distance to lattice boundary is encountered, we could simply call exhaustive_find_cell on the particle, which should correctly reorient it and make its lattice and spatial coordinates agree.

Lastly, in the event that negative distances to lattice boundaries occur, the lattice itself could recalculate only the lattice coordinates based on the particle's current position, then return to a distance-to-boundary calculation.

Or, perhaps, we should simply give up on this literal corner case and accept that it happens with acceptably low probability in non-pathological simulations.

import openmc
import numpy as np

# length of lattice on each side
lat_size = 20.0

# angle that we're crossing the corner at
phi = np.pi/4.0

# number of lattice elements across
numlat = 2

air = openmc.Material(name='air')
air.set_density('g/cc', 0.001)
air.add_nuclide('N14', 1.)

metal = openmc.Material(name='metal')
metal.set_density('g/cc', 7)
metal.add_nuclide('Fe56', 1.)

mats = openmc.Materials([air, metal])
mats.export_to_xml()

left = openmc.XPlane(x0=-lat_size/2.0, name='left')
right = openmc.XPlane(x0=lat_size/2.0, name='right')
bottom = openmc.YPlane(y0=-lat_size/2.0, name='bottom')
top = openmc.YPlane(y0=lat_size/2.0, name='top')
cyl = openmc.ZCylinder(x0=0, y0=0, r=lat_size)
cyl.boundary_type = 'vacuum'

in_lattice_region = +left & -right & +bottom & -top

outside_lattice = openmc.Cell(name='outside_lattice')
outside_lattice.region = -cyl & ~in_lattice_region
outside_lattice.fill = air

inside_lattice = openmc.Cell(name='inside_lattice')
inside_lattice.region = in_lattice_region

root = openmc.Universe(name='root universe')

# can we straightforwardly define infinite universes?
metal_uni = openmc.Universe()
metal_cell = openmc.Cell()
metal_cell.fill = metal
metal_cell.region = -openmc.Sphere(r=1e90)
metal_uni.add_cell(metal_cell)

air_uni   = openmc.Universe()
air_cell = openmc.Cell()
air_cell.fill = air
air_cell.region = -openmc.Sphere(r=1e90)
air_uni.add_cell(air_cell)

# define a checkerboard lattice
lattice = openmc.RectLattice()
lattice.lower_left = [-lat_size/2.0, -lat_size/2.0]
lattice.pitch = [lat_size/numlat, lat_size/numlat]
lattice.universes = [[metal_uni if (i%2)==(j%2) else air_uni for i in range(numlat)]
                        for j in range(numlat)]
inside_lattice.fill = lattice
root.add_cells([outside_lattice, inside_lattice])

geometry = openmc.Geometry(root)
geometry.export_to_xml()

# Instantiate a Settings object, set all runtime parameters, and export to XML
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.batches = 10
settings.particles = int(1e3)
settings.max_lost_particles = int(1e8)
# settings.verbosity = 10

# define a source located outside the lattice and pointing straight into its
# corner at 45 degrees
angular_distr = openmc.stats.Monodirectional(reference_uvw=[np.cos(phi), np.sin(phi), 0.0])
offset = 0.0 #1e-6
spatial_distr = openmc.stats.Point(xyz=[-np.cos(phi)+offset, -np.sin(phi), 0.0])
settings.source = openmc.source.Source(space=spatial_distr, angle=angular_distr)
settings.export_to_xml()

plot = openmc.Plot(plot_id=1)
plot.origin = [0, 0, 0]

# encompass the surrounding cylinder and lattice
plot.width = [2 * lat_size, 2 * lat_size]
plot.pixels = [800, 800]
plot.color_by = 'material'

# Instantiate a Plots collection and export to XML
plot_file = openmc.Plots([plot])
plot_file.export_to_xml()

# Instantiate a tally mesh
mesh = openmc.RegularMesh()
mesh.dimension = [100, 100, 1]
mesh.lower_left = [-lat_size, -lat_size, -1e9]
mesh.upper_right = [lat_size, lat_size, 1e9]
mesh_filter = openmc.MeshFilter(mesh)

# Instantiate the Tally
tally = openmc.Tally(tally_id=1)
tally.filters = [mesh_filter]
tally.scores = ['flux']
tally.estimator = 'tracklength'
tallies_file = openmc.Tallies([tally])
tallies_file.export_to_xml()
@pshriwise
Copy link
Contributor

If there were an official geometry debugging season, it seems this would be it! Thanks for reporting this @gridley. And thanks for providing an MWE as well 💯 I'll have a look as soon as I can carve out some time.

@gridley
Copy link
Contributor Author

gridley commented Mar 28, 2023 •
edited
Loading

Cool cool. I'm thinking that it would be best to just explicitly check if the particle is crossing through a corner as defined by some appropriate tolerance (and this value may require some tweaking). This should be straightforward to implement, in practice. I think it would be something like checking if pairs of lattice edges both have a sufficiently close distance to intersection.

This should be easy to handle using some additional logic in RectLattice::distance where we include a possibility of values of lattice_trans being something like {1, 1, 0} and so on and so forth.

@gridley gridley mentioned this issue Mar 15, 2025
Open
@nuclearkevin nuclearkevin linked a pull request Mar 18, 2025 that will close this issue
Draft
5 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Add corner checks for rectangular lattice crossings nuclearkevin/openmc
2 participants
@pshriwise @gridley