Do not transport source particles below the energy cutoff

src/event.cpp

@@ -69,7 +69,9 @@ void process_init_events(int64_t n_particles, int64_t source_offset)
   for (int64_t i = 0; i < n_particles; i++) {
     initialize_particle_track(
       simulation::particles[i], source_offset + i + 1, false);
-    dispatch_xs_event(i);
+    if (simulation::particles[i].alive()) {
+      dispatch_xs_event(i);
+    }
   }
   simulation::time_event_init.stop();
 }

src/particle.cpp

@@ -279,9 +279,23 @@ void Particle::event_advance()
   // Find the distance to the nearest boundary
   boundary() = distance_to_boundary(*this);
 
+  // Photons below the energy cutoff (e.g., sampled from a source) are not
+  // transported. They collide at their current position, where the energy
+  // cutoff check in sample_photon_reaction deposits their energy locally. In a
+  // void there is nowhere to deposit the energy, so they are simply killed.
+  bool photon_below_cutoff =
+    type().is_photon() &&
+    E() < settings::energy_cutoff[type().transport_index()];
+  if (photon_below_cutoff && material() == MATERIAL_VOID) {
+    wgt() = 0.0;
+    return;
+  }
+
   // Sample a distance to collision
-  if (type() == ParticleType::electron() ||
-      type() == ParticleType::positron()) {
+  if (photon_below_cutoff) {
+    collision_distance() = 0.0;
+  } else if (type() == ParticleType::electron() ||
+             type() == ParticleType::positron()) {
     collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
   } else if (macro_xs().total == 0.0) {
     collision_distance() = INFINITY;

src/simulation.cpp

@@ -754,8 +754,18 @@ void initialize_particle_track(
   // Set particle track.
   p.write_track() = check_track_criteria(p);
 
+  // Determine whether a source neutron is below the energy cutoff, in which
+  // case it is killed before being transported. Photons, electrons, and
+  // positrons below the cutoff instead deposit their energy locally at their
+  // first collision, which occurs at the birth site.
+  bool below_cutoff = false;
+  if (!is_secondary && p.type().is_neutron()) {
+    int idx = p.type().transport_index();
+    below_cutoff = p.E() < settings::energy_cutoff[idx];
+  }
+
   // Set the particle's initial weight window value.
-  if (!is_secondary) {
+  if (!is_secondary && !below_cutoff) {
     p.wgt_ww_born() = -1.0;
     apply_weight_windows(p);
   }
@@ -771,6 +781,11 @@ void initialize_particle_track(
     simulation::total_weight += p.wgt();
   }
 
+  // Kill source neutrons below the energy cutoff
+  if (below_cutoff) {
+    p.wgt() = 0.0;
+  }
+
   // Force calculation of cross-sections by setting last energy to zero
   if (settings::run_CE) {
     p.invalidate_neutron_xs();

tests/unit_tests/test_energy_cutoff.py

@@ -61,3 +61,25 @@ def test_energy_cutoff(run_in_tmpdir):
     # Despite killing particles below the cutoff, the total heating should be
     # equal to the source energy
     assert heating[0] == pytest.approx(source_energy)
+
+
+def test_source_below_energy_cutoff(run_in_tmpdir):
+    # Source photons born below the energy cutoff should not be transported
+    # but should deposit their energy locally, as photons that fall below the
+    # cutoff during transport do
+    source_energy = 1e3
+    cutoff_energy = 2e3
+    model = inf_medium_model(cutoff_energy, 2*cutoff_energy)
+    model.settings.source = openmc.IndependentSource(
+        particle='photon',
+        energy=openmc.stats.Discrete([source_energy], [1.0]),
+    )
+    model.tallies[0].filters = [openmc.ParticleFilter(['photon'])]
+    statepoint_path = model.run()
+
+    with openmc.StatePoint(statepoint_path) as sp:
+        flux = sp.get_tally(name='flux').mean.ravel()
+        heating = sp.get_tally(name='heating').mean.ravel()
+
+    assert flux[0] == 0.0
+    assert heating[0] == pytest.approx(source_energy)