Ensure secondary particles below energy cutoff are not created

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()};

tests/regression_tests/photon_production/results_true.dat

@@ -1,8 +1,8 @@
 tally 1:
 8.610000E-01
 7.413210E-01
-9.503000E-01
-9.030701E-01
+9.493000E-01
+9.011705E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -16,12 +16,12 @@ tally 2:
 1.573004E+00
 4.296434E-04
 1.845934E-07
-8.124327E-01
-6.600469E-01
+2.337049E-01
+5.461796E-02
 0.000000E+00
 0.000000E+00
-8.124327E-01
-6.600469E-01
+2.337049E-01
+5.461796E-02
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -53,40 +53,40 @@ tally 3:
 4.104374E+12
 4.296582E-04
 1.846062E-07
-8.034000E-01
-6.454516E-01
-6.951747E+03
-4.832678E+07
+2.286000E-01
+5.225796E-02
+7.054033E+03
+4.975938E+07
 0.000000E+00
 0.000000E+00
-8.034000E-01
-6.454516E-01
-6.951747E+03
-4.832678E+07
+2.286000E-01
+5.225796E-02
+7.054033E+03
+4.975938E+07
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.762393E+05
-3.106028E+10
+1.764573E+05
+3.113718E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.762393E+05
-3.106028E+10
+1.764573E+05
+3.113718E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.631073E+03
-5.823328E+07
+7.691658E+03
+5.916160E+07
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.631073E+03
-5.823328E+07
+7.691658E+03
+5.916160E+07
 0.000000E+00
 0.000000E+00
 tally 4:
@@ -102,16 +102,16 @@ tally 4:
 4.104374E+12
 0.000000E+00
 0.000000E+00
-2.294000E-01
-5.262436E-02
-6.947074E+03
-4.826184E+07
+2.286000E-01
+5.225796E-02
+7.054033E+03
+4.975938E+07
 0.000000E+00
 0.000000E+00
-8.034000E-01
-6.454516E-01
-6.951747E+03
-4.832678E+07
+2.286000E-01
+5.225796E-02
+7.054033E+03
+4.975938E+07
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -122,8 +122,8 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.762393E+05
-3.106028E+10
+1.764573E+05
+3.113718E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -134,7 +134,7 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.631073E+03
-5.823328E+07
+7.691658E+03
+5.916160E+07
 0.000000E+00
 0.000000E+00

tests/regression_tests/photon_production_fission/results_true.dat

@@ -1,12 +1,12 @@
 k-combined:
-2.268512E+00 8.628052E-03
+2.294874E+00 2.454640E-02
 tally 1:
-2.657963E+00
-2.355039E+00
+2.691351E+00
+2.415076E+00
 0.000000E+00
 0.000000E+00
-2.657963E+00
-2.355039E+00
+2.691351E+00
+2.415076E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -18,52 +18,52 @@ tally 1:
 0.000000E+00
 0.000000E+00
 tally 2:
-2.640348E+00
-2.323851E+00
-4.246717E+08
-6.011634E+16
+2.654218E+00
+2.348806E+00
+4.260252E+08
+6.051306E+16
 0.000000E+00
 0.000000E+00
-2.640348E+00
-2.323851E+00
-4.246717E+08
-6.011634E+16
+2.654218E+00
+2.348806E+00
+4.260252E+08
+6.051306E+16
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.355438E+06
-1.849734E+12
+2.507232E+06
+2.096890E+12
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.355438E+06
-1.849734E+12
+2.507232E+06
+2.096890E+12
 0.000000E+00
 0.000000E+00
 tally 3:
-2.654437E+00
-2.348704E+00
-4.246717E+08
-6.011634E+16
+2.656221E+00
+2.351840E+00
+4.260252E+08
+6.051306E+16
 0.000000E+00
 0.000000E+00
-2.654437E+00
-2.348704E+00
-4.246717E+08
-6.011634E+16
+2.656221E+00
+2.351840E+00
+4.260252E+08
+6.051306E+16
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.173264E+06
-1.574708E+12
+2.507232E+06
+2.096890E+12
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.355438E+06
-1.849734E+12
+2.507232E+06
+2.096890E+12
 0.000000E+00
 0.000000E+00

tests/regression_tests/photon_source/results_true.dat

@@ -1,5 +1,5 @@
 tally 1:
-2.256359E+02
-5.091156E+04
+2.263938E+02
+5.125417E+04
 0.000000E+00
 0.000000E+00

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)