Add test for flux bias with weight windows in multigroup mode (#3202)

Co-authored-by: Paul Romano <[email protected]>

src/physics_mg.cpp

@@ -28,12 +28,12 @@ void collision_mg(Particle& p)
   // Add to the collision counter for the particle
   p.n_collision()++;
 
-  if (settings::weight_window_checkpoint_collision)
-    apply_weight_windows(p);
-
   // Sample the reaction type
   sample_reaction(p);
 
+  if (settings::weight_window_checkpoint_collision)
+    apply_weight_windows(p);
+
   // Display information about collision
   if ((settings::verbosity >= 10) || p.trace()) {
     write_message(fmt::format("    Energy Group = {}", p.g()), 1);

tests/unit_tests/weightwindows/test_ww_mg.py

@@ -0,0 +1,51 @@
+import numpy as np
+import openmc
+
+
+def test_weight_windows_mg(request, run_in_tmpdir):
+    # import basic random ray model
+    model = openmc.examples.random_ray_three_region_cube()
+
+    # create a mesh tally
+    mesh = openmc.RegularMesh.from_domain(model.geometry, (3, 3, 3))
+    mesh_tally = openmc.Tally()
+    mesh_tally.filters = [openmc.MeshFilter(mesh)]
+    mesh_tally.scores = ['flux']
+    model.tallies = [mesh_tally]
+
+    # replace random ray settings with fixed source settings
+    settings = openmc.Settings()
+    settings.particles = 5000
+    settings.batches = 10
+    settings.energy_mode = 'multi-group'
+    settings.run_mode = 'fixed source'
+    space = openmc.stats.Point((1, 1, 1))
+    energy = openmc.stats.delta_function(1e6)
+    settings.source = openmc.IndependentSource(space=space, energy=energy)
+    model.settings = settings
+
+    # perform analog simulation
+    statepoint = model.run()
+
+    # extract flux from analog simulation
+    with openmc.StatePoint(statepoint) as sp:
+        tally_out = sp.get_tally(id=mesh_tally.id)
+        flux_analog = tally_out.mean
+
+    # load the weight windows for this problem and apply them
+    ww_lower_bnds = np.loadtxt(request.path.parent / 'ww_mg.txt')
+    weight_windows = openmc.WeightWindows(mesh, lower_ww_bounds=ww_lower_bnds, upper_bound_ratio=5.0)
+    model.settings.weight_windows = weight_windows
+    model.settings.weight_windows_on = True
+
+    # re-run with weight windows
+    statepoint = model.run()
+    with openmc.StatePoint(statepoint) as sp:
+        tally_out = sp.get_tally(id=mesh_tally.id)
+        flux_ww = tally_out.mean
+
+    # the sum of the fluxes should approach the same value (no bias introduced)
+    analog_sum = flux_analog.sum()
+    ww_sum = flux_ww.sum()
+    assert np.allclose(analog_sum, ww_sum, rtol=1e-2)
+

tests/unit_tests/weightwindows/ww_mg.txt

@@ -0,0 +1,27 @@
+5.000000000000000278e-02
+1.023184121346435924e-02
+3.006624096325660397e-03
+1.030178532774538719e-02
+6.058877789444589400e-03
+2.216914234856166583e-03
+3.061186967456217597e-03
+2.148267952185671601e-03
+1.026171712186230980e-03
+1.040022203443005666e-02
+6.040633813799485378e-03
+2.364143118752318716e-03
+6.119726639841410569e-03
+4.329097093078606955e-03
+1.873104469085542763e-03
+2.246957229279350661e-03
+1.851165248260521617e-03
+7.825824911598703530e-04
+3.021300894848398706e-03
+2.286420236345795311e-03
+9.318583473482396160e-04
+2.234702678114806364e-03
+1.813664566119888152e-03
+7.969287848384389462e-04
+1.017895970086981662e-03
+7.707144532950136471e-04
+3.386087166633241791e-04