From 822879cc4838e52bf87202d11429756dbe234474 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 12 Dec 2021 16:27:12 -0600
Subject: [PATCH] Ensure secondary particles below energy cutoff are not
 created

---
 src/particle.cpp                       |  6 +++
 tests/unit_tests/test_energy_cutoff.py | 63 ++++++++++++++++++++++++++
 2 files changed, 69 insertions(+)
 create mode 100644 tests/unit_tests/test_energy_cutoff.py

diff --git a/src/particle.cpp b/src/particle.cpp
index 7579f5f5a8..fa35f634af 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -39,6 +39,12 @@ namespace openmc {
 void Particle::create_secondary(
   double wgt, Direction u, double E, ParticleType type)
 {
+  // If energy is below cutoff for this particle, don't create secondary
+  // particle
+  if (E < settings::energy_cutoff[static_cast<int>(type)]) {
+    return;
+  }
+
   secondary_bank().emplace_back();
 
   auto& bank {secondary_bank().back()};
diff --git a/tests/unit_tests/test_energy_cutoff.py b/tests/unit_tests/test_energy_cutoff.py
new file mode 100644
index 0000000000..fb82c579a1
--- /dev/null
+++ b/tests/unit_tests/test_energy_cutoff.py
@@ -0,0 +1,63 @@
+from random import uniform
+
+import pytest
+import openmc
+
+
+def inf_medium_model(cutoff_energy, source_energy):
+    """Infinite medium problem with a monoenergetic photon source"""
+    model = openmc.Model()
+
+    m = openmc.Material()
+    m.add_nuclide('Zr90', 1.0)
+    m.set_density('g/cm3', 1.0)
+
+    sph = openmc.Sphere(r=100.0, boundary_type='reflective')
+    cell = openmc.Cell(fill=m, region=-sph)
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.source = openmc.Source(
+        particle='photon',
+        energy=openmc.stats.Discrete([source_energy], [1.0]),
+    )
+    model.settings.particles = 100
+    model.settings.batches = 10
+    model.settings.cutoff = {'energy_photon': cutoff_energy}
+
+    tally_flux = openmc.Tally(name='flux')
+    tally_flux.filters = [
+        openmc.EnergyFilter([0.0, cutoff_energy, source_energy]),
+        openmc.ParticleFilter(['photon'])
+    ]
+    tally_flux.scores = ['flux']
+    tally_heating = openmc.Tally(name='heating')
+    tally_heating.scores = ['heating']
+    model.tallies = openmc.Tallies([tally_flux, tally_heating])
+
+    return model
+
+
+def test_energy_cutoff(run_in_tmpdir):
+    # Pick a random cutoff energy between 1 and 5 keV
+    cutoff_energy = uniform(1e3, 5e3)
+
+    # Pick a random source energy some factor higher than cutoff energy
+    source_energy = uniform(10, 20) * cutoff_energy
+
+    # Create model and run simulation
+    model = inf_medium_model(cutoff_energy, source_energy)
+    statepoint_path = model.run()
+
+    # Get resulting flux and heating values
+    with openmc.StatePoint(statepoint_path) as sp:
+        flux = sp.get_tally(name='flux').mean.ravel()
+        heating = sp.get_tally(name='heating').mean.ravel()
+
+    # There should be no flux below the cutoff energy (first bin in the tally)
+    assert flux[0] == 0.0
+    assert flux[1] > 0.0
+
+    # Despite killing particles below the cutoff, the total heating should be
+    # equal to the source energy
+    assert heating[0] == pytest.approx(source_energy)