Skip to content

Fix weight window energy defaults being derived before the particle type is known - #4091

Merged
paulromano merged 3 commits into
openmc-dev:developfrom
GuySten:ww-particle-type-fix
Sep 1, 2026
Merged

Fix weight window energy defaults being derived before the particle type is known#4091
paulromano merged 3 commits into
openmc-dev:developfrom
GuySten:ww-particle-type-fix

Conversation

@GuySten

@GuySten GuySten commented Aug 31, 2026

Copy link
Copy Markdown
Contributor

Fix weight window energy defaults being derived before the particle type is known

Scope: this affects weight windows created through the C API during a
simulation
without an explicit energy_bounds. Photon ones previously received
neutron default energy bounds; they now receive photon ones. Weight windows
defined in XML are unaffected, as are any with an explicit energy_bounds, as
are any created before openmc_simulation_init().

Problem

WeightWindows::set_defaults() derives the default energy grid from
data::energy_min/max for the object's particle type, and only does so when the
grid is empty:

void WeightWindows::set_defaults()
{
  if (energy_bounds_.size() == 0) {
    int p_type = particle_type_.transport_index();
    ...
    energy_bounds_.push_back(data::energy_min[p_type]);
    energy_bounds_.push_back(data::energy_max[p_type]);
  }
}

The WeightWindows(int32_t id) constructor, which is the one the C API and
openmc.lib use, called it immediately:

WeightWindows::WeightWindows(int32_t id)
{
  index_ = variance_reduction::weight_windows.size();
  set_id(id);
  set_defaults();   // particle type, mesh and energy grid all unknown here
}

ParticleType default-constructs to PDG_NEUTRON, so this does not error. It
silently derives neutron bounds. set_particle_type() does not re-derive them,
and because the grid is no longer empty every subsequent set_defaults() call
is a no-op. So the early call does not merely compute the wrong answer, it
prevents the right one from ever being computed.

The XML constructor is not affected: it assigns particle_type_ before calling
set_defaults() at the end.

The window in which this is observable is narrower than it first appears.
data::energy_min/max are populated by initialize_data(), which runs from
openmc_simulation_init(), not openmc_init(). Before simulation
initialization both arrays hold their static defaults of 0.0 and INFTY for
every particle, so the neutron and photon defaults coincide and nothing is
visibly wrong. The bug bites weight windows constructed through the C API once a
simulation is running -- which is the regime weight window generation operates
in, and is why WeightWindowsGenerator needs its workaround.

Two things in the existing code point at this:

  • WeightWindowsGenerator::WeightWindowsGenerator calls set_defaults() again
    after set_mesh(), set_energy_bounds() and set_particle_type(). That is a
    workaround, and it only works because the generator constructs through a path
    where the grid is still empty.
  • The energy_bounds_.size() == 0 guard is what turns an early call into a
    permanent one.

MWE

import numpy as np
import openmc
import openmc.lib

water = openmc.Material()
water.set_density('g/cm3', 1.0)
water.add_nuclide('H1', 2.0)
water.add_nuclide('O16', 1.0)

sphere = openmc.Sphere(r=10.0, boundary_type='vacuum')
model = openmc.Model()
model.geometry = openmc.Geometry([openmc.Cell(fill=water, region=-sphere)])
model.settings.run_mode = 'fixed source'
model.settings.particles = 100
model.settings.batches = 1
model.settings.photon_transport = True     # so photon data is loaded too
model.export_to_model_xml()


def make_ww(uid, particle):
    """Create a weight window through the C API, as openmc.lib does."""
    mesh = openmc.lib.RegularMesh()
    mesh.dimension = (2, 2, 2)
    mesh.set_parameters(lower_left=(-1.0, -1.0, -1.0),
                        upper_right=(1.0, 1.0, 1.0))
    ww = openmc.lib.WeightWindows(uid)
    ww.mesh = mesh
    ww.particle = particle
    return np.asarray(ww.energy_bounds)


openmc.lib.init()
openmc.lib.simulation_init()   # populates data::energy_min/max
try:
    neutron = make_ww(900, 'neutron')
    photon = make_ww(901, 'photon')
finally:
    openmc.lib.simulation_finalize()
    openmc.lib.finalize()

print(f'\n  neutron default energy bounds  {neutron}')
print(f'  photon  default energy bounds  {photon}')

if np.allclose(neutron, photon):
    print('\n  FAIL: the photon weight window has neutron energy bounds.')
    print('  Photon data extends far beyond the neutron limit, so an upper')
    print('  bound equal to the neutron one cannot be right.')
    raise SystemExit(1)

print('\n  PASS: each particle got its own default energy range.')

Changes

  • WeightWindows(int32_t id) no longer calls set_defaults(). A comment
    explains why, so it is not reinstated.

  • set_particle_type() calls set_defaults(), which is the point at which the
    defaults are actually determined. It is a no-op when an explicit grid has
    already been supplied.

  • openmc_weight_windows_export() calls set_defaults() once per weight window
    before writing, so objects built through the C API whose particle type was
    never set do not export an empty energy grid.

  • WeightWindowsGenerator drops its trailing set_defaults() call, which is
    now redundant. This is the workaround the bug forced, so its removal is the
    clearest demonstration that the fix works.

  • The XML constructor assigns the particle type through set_particle_type()
    instead of writing particle_type_ directly, and drops its own trailing
    set_defaults(). Two side benefits: the XML path now gets the same
    neutron/photon validation as the C API path, and every existing XML-based
    weight window test exercises the new call site.

After this, set_defaults() has exactly two call sites -- set_particle_type()
and the export backstop -- so every path that assigns a particle type derives
its defaults the same way.

bounds_size() already treats an empty energy grid as one bin:

int num_energy_bins = energy_bounds_.size() > 0 ? energy_bounds_.size() - 1 : 1;

so set_mesh() continues to allocate a correctly shaped 1 x n_spatial bounds
tensor in the window between construction and the particle type being set.

Testing

New tests/unit_tests/weightwindows/test_ww_defaults.py, which initializes a
model with both neutron and photon data and calls simulation_init() so that
data::energy_min/max are actually populated:

  • Neutron and photon weight windows must get different default energy bounds.
    Identical bounds are the signature of the defaults being derived from the
    default particle type at construction.
  • An explicit energy_bounds grid survives a later set_particle_type().
  • Populated weight window bounds survive a later set_particle_type().

Notes

Found while reviewing #4057, but it is independent of that PR and reproduces on
develop. It may allow #4057 to drop its reset_energy_bounds() addition,
since a freshly constructed openmc.lib.WeightWindows will now pick up correct
defaults on its own once the particle type is assigned.

Checklist

  • I have performed a self-review of my own code
  • I have run clang-format (version 18) 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)

@GuySten
GuySten force-pushed the ww-particle-type-fix branch from aaf8210 to 78adf76 Compare August 31, 2026 19:05
@GuySten GuySten added the Bugs label Aug 31, 2026
@GuySten
GuySten marked this pull request as ready for review August 31, 2026 19:54
@GuySten
GuySten requested a review from pshriwise as a code owner August 31, 2026 19:54
@GuySten
GuySten requested a review from paulromano August 31, 2026 19:54
@paulromano
paulromano enabled auto-merge (squash) September 1, 2026 13:26
@paulromano
paulromano merged commit c6f9187 into openmc-dev:develop Sep 1, 2026
16 checks passed
@GuySten
GuySten deleted the ww-particle-type-fix branch September 1, 2026 14:52
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants