Add TimedMeshFilter

[ilhamv]

Description

There has been a known tallying issue in using a combination of TimeFilter and MeshFilter with the track-length estimator (for example, as mentioned in [1], [2], and [3]). This is because neutron tracks span 4 dimensions (x, y, z, and t), so filtering both in time and space needs to be done simultaneously, not separately, as in the case of using TimeFilter + MeshFilter.

This PR adds TimedMeshFilter that accurately filters events with respect to a given time grid and a given mesh, even with the track-length estimator.

Verification

The 1D 1G AZURV1 problem (a supercritical version, input script attached below) is used for verification. Note that there are known bugs in the OpenMC's time-dependent MG mode (see [3]), so to make it works correctly, some of the fixes in [3] are implemented in this PR as well.

Below are OpenMC results (10 batches, 100k particles/batch) for
(1) TimeFilter + MeshFilter with collision estimator,
(2) TimeFilter + MeshFilter with track-length estimator, and
(3) TimedMeshFilter with track-length estimator.

TimeFilter + MeshFilter results in an accurate solution if used with collition estimator, but it gives a wrong solution if used with track-length estimator. In particular, note that TimeFilter + MeshFilter with tracklength estimator produces nonphysical solution (non-zero flux beyond the physical wavefront of the neutrons)! TimedMeshFilter with track-lengh estimator, however, produces accurate results.

To drive home the verification, here are the error convergences of the three cases as a function of the number of particles per batch $N$:

Here is the input script for the verification test:

import openmc
import numpy as np
import h5py

# ===========================================================================
# Set Library
# ===========================================================================

SigmaC = 1.0/3.0
SigmaF = 1.0/3.0
nu = 2.3
SigmaA = SigmaC + SigmaF
SigmaS = 1.0/3.0
SigmaT = SigmaA + SigmaS
v = 1.0

groups = openmc.mgxs.EnergyGroups([0.0, 2e7])

xsdata = openmc.XSdata("mat", groups)
xsdata.order = 0

xsdata.set_inverse_velocity([1.0 / v], temperature=294.0)

xsdata.set_total([SigmaT], temperature=294.0)
xsdata.set_absorption([SigmaA], temperature=294.0)
xsdata.set_scatter_matrix(np.ones((1, 1, 1)) * SigmaS, temperature=294.0)

xsdata.set_nu_fission([nu * SigmaF], temperature=294.0)
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.add_xsdata(xsdata)
mg_cross_sections_file.export_to_hdf5("mgxs.h5")

# ===========================================================================
# Exporting to OpenMC materials.xml file
# ===========================================================================

materials = {}
materials["mat"] = openmc.Material(name="mat")
materials["mat"].set_density("macro", 1.0)
materials["mat"].add_macroscopic("mat")
materials_file = openmc.Materials(materials.values())
materials_file.cross_sections = "mgxs.h5"
materials_file.export_to_xml()

# ===========================================================================
# Exporting to OpenMC geometry.xml file
# ===========================================================================

# Instantiate ZCylinder surfaces
surf_Z1 = openmc.XPlane(surface_id=1, x0=-1e10, boundary_type="reflective")
surf_Z2 = openmc.XPlane(surface_id=2, x0=1e10, boundary_type="reflective")

# Instantiate Cells
cell_F = openmc.Cell(cell_id=1, name="F")

# Use surface half-spaces to define regions
cell_F.region = +surf_Z1 & -surf_Z2

# Register Materials with Cells
cell_F.fill = materials["mat"]

# Instantiate Universes
root = openmc.Universe(universe_id=0, name="root universe", cells=[cell_F])

# Instantiate a Geometry, register the root Universe, and export to XML
geometry = openmc.Geometry(root)
geometry.export_to_xml()

# ===========================================================================
# Exporting to OpenMC settings.xml file
# ===========================================================================

# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.run_mode = "fixed source"
settings_file.particles = 100000
settings_file.batches = 10
settings_file.output = {"tallies": False}
settings_file.cutoff = {"time_neutron": 20}
settings_file.energy_mode = "multi-group"

# Create an initial uniform spatial source distribution over fissionable zones
delta_dist = openmc.stats.Point()
isotropic = openmc.stats.Isotropic()
settings_file.source = openmc.IndependentSource(space=delta_dist, angle=isotropic)
settings_file.export_to_xml()


# ===========================================================================
# Exporting to OpenMC tallies.xml file
# ===========================================================================

# Create a mesh filter that can be used in a tally
mesh = openmc.RegularMesh()
mesh.dimension = (201, 1, 1)
mesh.lower_left = (-20.5, -1e10, -1e10)
mesh.upper_right = (20.5, 1e10, 1e10)
time_grid = np.linspace(0.0, 20.0, 41)

mesh_filter = openmc.MeshFilter(mesh)
time_filter = openmc.TimeFilter(time_grid)
timed_mesh_filter = openmc.TimedMeshFilter(mesh, time_grid)

# Now use the mesh filter in a tally and indicate what scores are desired
tally1 = openmc.Tally(name="old-collision")
tally1.estimator = "collision"
tally1.filters = [time_filter, mesh_filter]
tally1.scores = ["flux"]

tally2 = openmc.Tally(name="old-tracklength")
tally2.estimator = "tracklength"
tally2.filters = [time_filter, mesh_filter]
tally2.scores = ["flux"]

tally3 = openmc.Tally(name="new-tracklength")
tally3.estimator = "tracklength"
tally3.filters = [timed_mesh_filter]
tally3.scores = ["flux"]

# Instantiate a Tallies collection and export to XML
tallies = openmc.Tallies([tally1, tally2, tally3])
tallies.export_to_xml()

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)