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

src/weight_windows.cpp

@@ -52,7 +52,6 @@ WeightWindows::WeightWindows(int32_t id)
 {
   index_ = variance_reduction::weight_windows.size();
   set_id(id);
-  set_defaults();
 }
 
 WeightWindows::WeightWindows(pugi::xml_node node)
@@ -70,9 +69,9 @@ WeightWindows::WeightWindows(pugi::xml_node node)
   int32_t id = std::stoi(get_node_value(node, "id"));
   this->set_id(id);
 
-  // get the particle type
+  // Get the particle type
   auto particle_type_str = std::string(get_node_value(node, "particle_type"));
-  particle_type_ = ParticleType {particle_type_str};
+  set_particle_type(ParticleType {particle_type_str});
 
   // Determine associated mesh
   int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
@@ -117,8 +116,6 @@ WeightWindows::WeightWindows(pugi::xml_node node)
   // read the lower/upper weight bounds
   this->set_bounds(get_node_array<double>(node, "lower_ww_bounds"),
     get_node_array<double>(node, "upper_ww_bounds"));
-
-  set_defaults();
 }
 
 WeightWindows::~WeightWindows()
@@ -251,6 +248,10 @@ void WeightWindows::set_particle_type(ParticleType p_type)
     fatal_error(fmt::format(
       "Particle type '{}' cannot be applied to weight windows.", p_type.str()));
   particle_type_ = p_type;
+
+  // The default energy grid is particle dependent, so derive it now that the
+  // particle type is known
+  set_defaults();
 }
 
 void WeightWindows::set_mesh(int32_t mesh_idx)
@@ -847,7 +848,6 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
   if (e_bounds.size() > 0)
     wws->set_energy_bounds(e_bounds);
   wws->set_particle_type(particle_type);
-  wws->set_defaults();
 }
 
 void WeightWindowsGenerator::create_tally()
@@ -1341,6 +1341,9 @@ extern "C" int openmc_weight_windows_export(const char* filename)
   std::vector<int32_t> mesh_ids;
   std::vector<int32_t> ww_ids;
   for (const auto& ww : variance_reduction::weight_windows) {
+    // Backstop for objects built through the C API whose particle type was
+    // never set explicitly, so an empty energy grid is never written out
+    ww->set_defaults();
 
     ww->to_hdf5(weight_windows_group);
     ww_ids.push_back(ww->id());

tests/unit_tests/weightwindows/test_ww_defaults.py

@@ -0,0 +1,93 @@
+"""Default weight window energy bounds must follow the particle type."""
+
+import numpy as np
+import pytest
+import openmc
+import openmc.lib
+
+
+@pytest.fixture
+def lib_model(run_in_tmpdir):
+    """Minimal model with both neutron and photon data available."""
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    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')
+    cell = openmc.Cell(fill=water, region=-sphere)
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 100
+    model.settings.batches = 1
+    model.settings.photon_transport = True
+
+    model.export_to_model_xml()
+    openmc.lib.init()
+    # Required: initialize_data() runs here, not in openmc_init(), and it is
+    # what narrows data::energy_min/max to the loaded data for each particle
+    openmc.lib.simulation_init()
+    yield model
+    openmc.lib.simulation_finalize()
+    openmc.lib.finalize()
+
+
+def _lib_mesh():
+    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))
+    return mesh
+
+
+@pytest.mark.parametrize('particle', ('neutron', 'photon'))
+def test_default_energy_bounds_follow_particle(lib_model, particle):
+    """Defaults are derived after the particle type is known, not before."""
+    ww = openmc.lib.WeightWindows(300 if particle == 'neutron' else 301)
+    ww.mesh = _lib_mesh()
+    ww.particle = particle
+
+    bounds = np.asarray(ww.energy_bounds)
+    assert bounds.size == 2
+
+    # Compare against the range the library reports for this particle. Using
+    # the other particle's range would be the symptom of deriving defaults in
+    # the constructor.
+    other = 'photon' if particle == 'neutron' else 'neutron'
+    other_ww = openmc.lib.WeightWindows(400 if particle == 'neutron' else 401)
+    other_ww.mesh = _lib_mesh()
+    other_ww.particle = other
+    other_bounds = np.asarray(other_ww.energy_bounds)
+
+    assert not np.allclose(bounds, other_bounds), (
+        f'{particle} and {other} weight windows have identical default energy '
+        'bounds, which suggests the defaults were not derived from the '
+        'particle type'
+    )
+
+
+def test_explicit_energy_bounds_survive_particle_type(lib_model):
+    """Setting the particle type must not overwrite an explicit grid."""
+    ww = openmc.lib.WeightWindows(302)
+    ww.mesh = _lib_mesh()
+    ww.energy_bounds = (1.0e3, 1.0e5, 1.0e7)
+    ww.particle = 'photon'
+
+    np.testing.assert_allclose(ww.energy_bounds, (1.0e3, 1.0e5, 1.0e7))
+
+
+def test_bounds_survive_particle_type(lib_model):
+    """Setting the particle type must not discard weight window bounds."""
+    ww = openmc.lib.WeightWindows(303)
+    ww.mesh = _lib_mesh()
+    ww.energy_bounds = (0.0, 1.0e7)
+    lower = np.arange(1.0, 9.0)
+    ww.bounds = lower, 5.0 * lower
+    ww.particle = 'photon'
+
+    np.testing.assert_allclose(ww.bounds[0], lower)
+    np.testing.assert_allclose(ww.bounds[1], 5.0 * lower)