Moving surface

[ilhamv]

Description

This PR introduces a continuously moving surface capability for time-dependent simulations. The Surface class now has the method move that takes velocities and durations to support piecewise continuous linear movement. As an example:

# Surfaces
cylinder = openmc.ZCylinder(r=1.0)
top = openmc.ZPlane(z0=10.0)
bottom = openmc.ZPlane(z0=-10.0)

# Move surfaces
velocities = np.zeros((3,3))
velocities[:,2] = np.array([1.0, -2.0])
durations = np.array([4.0, 2.0])
top.move(velocities, durations)
bottom.move(velocities, durations)

# A moving region
moving_rod = -cylinder & -top & +bottom

The example above shows how we move a rod region by moving the top and bottom bounding surfaces. In the example, the rod region moves upward at the speed of 1 cm/s for 4 seconds, and then downward at the speed of 2 cm/s for 2 seconds, returning to the initial position.

The implementation in this PR is an extension of this conference paper.

Verification: Moving absorber

Let us consider a boron pellet, placed at a distance from a point source and moving in the air. Below is the resulting absorption reaction rate map over time, followed by the Python input script for the model.

import openmc
import numpy as np


# =============================================================================
# Materials
# =============================================================================

air = openmc.Material()
air.add_element('N', 0.784431)
air.add_element('O', 0.210748)
air.add_element('Ar',0.0046)
air.set_density('g/cc', 0.001205)

boron = openmc.Material()
boron.add_element('B', 1.0, enrichment=100.0, enrichment_target='B10')
boron.set_density('g/cm3', 2.34)

materials = openmc.Materials([air, boron])
materials.export_to_xml()

# =============================================================================
# Geometry
# =============================================================================

tank = openmc.model.RectangularParallelepiped(-10, 10, -10, 10, -10, 10, boundary_type='vacuum')

cylinder = openmc.ZCylinder(r=1.0, x0=-3.0, name='cylinder', surface_id=111)
top = openmc.ZPlane(z0=9.0, name='moving top', surface_id=333)
bottom = openmc.ZPlane(z0=6.0, name='moving bottom', surface_id=222)

# Z-direction movements
velocities = np.zeros((3,3))
velocities[:,2] = np.array([-1.0, 2.0, -10.0])
durations = np.array([5.0, 2.0, 1.0])
top.move(velocities, durations)
bottom.move(velocities, durations)

# X-direction movements
velocities = np.zeros((3,3))
velocities[:,0] = np.array([3.0, -4.0, 1.0])
durations = np.array([2.0, 2.0, 4.0])
cylinder.move(velocities, durations)

pellet_region = -cylinder & -top & +bottom
air_tank = openmc.Cell(region = -tank & ~pellet_region, fill=air, cell_id=1)
pellet = openmc.Cell(region = pellet_region, fill=boron, cell_id=2)

geometry = openmc.Geometry([air_tank, pellet])
geometry.export_to_xml()


# =============================================================================
# Settings
# =============================================================================

settings_file = openmc.Settings()
settings_file.run_mode = "fixed source"
settings_file.particles = 1000000000
settings_file.batches = 1
settings_file.output = {"tallies": False}
settings_file.cutoff = {"time_neutron": 9.0}

delta_dist = openmc.stats.Point((3.0, 0.0, 0.0))
isotropic = openmc.stats.Isotropic()
energy = openmc.stats.Discrete([1.4E7], [1.0])
time = openmc.stats.Uniform(0.0, 9.0)
settings_file.source = openmc.IndependentSource(space=delta_dist, angle=isotropic, energy=energy, time=time)
settings_file.export_to_xml()

# =============================================================================
# Tallies
# =============================================================================

mesh = openmc.RegularMesh()
mesh.dimension = (200, 1, 200)
mesh.lower_left = (-10, -10, -10)
mesh.upper_right = (10, 10, 10)
time_grid = np.linspace(0.0, 9.0, 181)

time_filter = openmc.TimeFilter(time_grid)
mesh_filter = openmc.MeshFilter(mesh)

tally = openmc.Tally(name='tally')
tally.estimator = "collision"
tally.filters = [time_filter, mesh_filter]
tally.scores = ["absorption"]

tallies = openmc.Tallies([tally])

tallies.export_to_xml()

Verification: Four-phase C5G7

Another form of verification is presented by the code-to-code comparison of the Four-phase C5G7 transient benchmark, which involves the movement of the control rod banks, exercising a range of transients.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[gonuke]

Thanks for this interesting addition @ilhamv. Can you clarify what time scales of movement this is meant to support? Specifically, are movements expected to be relevant in the time scale of an individual neutron history? Or is this intended for problems in which a snapshot of the geometry is used for a given simulation (set of histories)?

[ilhamv]

Thanks for the interest, @gonuke! It should support arbitrary time scales, even those comparable to neutron speed, as long as it is within OpenMC's precision tolerance. It is achieved by making both "surface evaluation" and "distance to surface" functions time-dependent, respecting the relative position and the relative velocity of the neutron and the surface.

[gonuke]

Do we worry about motion/source location in k-eigenvalue problems, esp. for delayed neutrons?

[ilhamv]

That's a great point. Delayed products should be drifted along with the containing, moving cell. This can be implemented by adding a drifting flag, as well as velocities and durations arrays as members of Cell, and then applying it whenever a delayed product is sampled by accordingly adjusting the position.

This PR, Surface.move(durations, velocities), does not include such a product drifting effect. A planned follow-up PR, Cell.drift(durations, velocities), would handle that. My original thought was that separating Surface.move and Cell.drift as two independent capabilities would give users more flexibility in building transient models (there may be cases where one only wants to use one of them). The best way to include these really is up for a discussion, though.

[cfichtlscherer]

This is very interesting—thank you so much for all your work!

I’m wondering: could I also use this feature to model a moving detector? For example, would it work if I define an absorption or pulse-height tally in a moving cell? Or is the mechanism that I tally a mesh?

[ilhamv]

This is very interesting—thank you so much for all your work!

I’m wondering: could I also use this feature to model a moving detector? For example, would it work if I define an absorption or pulse-height tally in a moving cell? Or is the mechanism that I tally a mesh?

Thanks, @cfichtlscherer!

It should be able to support moving detectors. To get the tally moving, we can filter by material or cell (MaterialFilter or CellFilter).