Adjust for secondary particle energy directly in heating scores (#3227)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

include/openmc/particle.h

@@ -53,6 +53,12 @@ class Particle : public ParticleData {
   //! \param type Particle type
   void create_secondary(double wgt, Direction u, double E, ParticleType type);
 
+  //! split a particle
+  //
+  //! creates a new particle with weight wgt
+  //! \param wgt Weight of the new particle
+  void split(double wgt);
+
   //! initialize from a source site
   //
   //! initializes a particle from data stored in a source site. The source

include/openmc/particle_data.h

@@ -430,7 +430,7 @@ class ParticleData : public GeometryState {
   int delayed_group_ {0};
 
   int n_bank_ {0};
-  int n_bank_second_ {0};
+  double bank_second_E_ {0.0};
   double wgt_bank_ {0.0};
   int n_delayed_bank_[MAX_DELAYED_GROUPS];
 
@@ -544,11 +544,13 @@ class ParticleData : public GeometryState {
   int& delayed_group() { return delayed_group_; } // delayed group
 
   // Post-collision data
-  int& n_bank() { return n_bank_; } // number of banked fission sites
-  int& n_bank_second()
+  double& bank_second_E()
   {
-    return n_bank_second_;
-  }                                        // number of secondaries banked
+    return bank_second_E_;
+  } // energy of last reaction secondaries
+  const double& bank_second_E() const { return bank_second_E_; }
+
+  int& n_bank() { return n_bank_; }        // number of banked fission sites
   double& wgt_bank() { return wgt_bank_; } // weight of banked fission sites
   int* n_delayed_bank()
   {

src/particle.cpp

@@ -91,17 +91,25 @@ void Particle::create_secondary(
     return;
   }
 
-  secondary_bank().emplace_back();
-
-  auto& bank {secondary_bank().back()};
+  auto& bank = secondary_bank().emplace_back();
   bank.particle = type;
   bank.wgt = wgt;
   bank.r = r();
   bank.u = u;
   bank.E = settings::run_CE ? E : g();
   bank.time = time();
+  bank_second_E() += bank.E;
+}
 
-  n_bank_second() += 1;
+void Particle::split(double wgt)
+{
+  auto& bank = secondary_bank().emplace_back();
+  bank.particle = type();
+  bank.wgt = wgt;
+  bank.r = r();
+  bank.u = u();
+  bank.E = settings::run_CE ? E() : g();
+  bank.time = time();
 }
 
 void Particle::from_source(const SourceSite* src)
@@ -356,7 +364,7 @@ void Particle::event_collide()
 
   // Reset banked weight during collision
   n_bank() = 0;
-  n_bank_second() = 0;
+  bank_second_E() = 0.0;
   wgt_bank() = 0.0;
   zero_delayed_bank();
 
@@ -417,6 +425,7 @@ void Particle::event_revive_from_secondary()
     from_source(&secondary_bank().back());
     secondary_bank().pop_back();
     n_event() = 0;
+    bank_second_E() = 0.0;
 
     // Subtract secondary particle energy from interim pulse-height results
     if (!model::active_pulse_height_tallies.empty() &&

src/physics.cpp

@@ -239,11 +239,10 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     }
 
     // Write fission particles to nuBank
-    p.nu_bank().emplace_back();
-    NuBank* nu_bank_entry = &p.nu_bank().back();
-    nu_bank_entry->wgt = site.wgt;
-    nu_bank_entry->E = site.E;
-    nu_bank_entry->delayed_group = site.delayed_group;
+    NuBank& nu_bank_entry = p.nu_bank().emplace_back();
+    nu_bank_entry.wgt = site.wgt;
+    nu_bank_entry.E = site.E;
+    nu_bank_entry.delayed_group = site.delayed_group;
   }
 
   // If shared fission bank was full, and no fissions could be added,

src/physics_mg.cpp

@@ -188,11 +188,10 @@ void create_fission_sites(Particle& p)
     }
 
     // Write fission particles to nuBank
-    p.nu_bank().emplace_back();
-    NuBank* nu_bank_entry = &p.nu_bank().back();
-    nu_bank_entry->wgt = site.wgt;
-    nu_bank_entry->E = site.E;
-    nu_bank_entry->delayed_group = site.delayed_group;
+    NuBank& nu_bank_entry = p.nu_bank().emplace_back();
+    nu_bank_entry.wgt = site.wgt;
+    nu_bank_entry.E = site.E;
+    nu_bank_entry.delayed_group = site.delayed_group;
   }
 
   // If shared fission bank was full, and no fissions could be added,

src/tallies/tally_scoring.cpp

@@ -964,14 +964,9 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
           // The energy deposited is the difference between the pre-collision
           // and post-collision energy...
           score = E - p.E();
-
           // ...less the energy of any secondary particles since they will be
           // transported individually later
-          const auto& bank = p.secondary_bank();
-          for (auto it = bank.end() - p.n_bank_second(); it < bank.end();
-               ++it) {
-            score -= it->E;
-          }
+          score -= p.bank_second_E();
 
           score *= p.wgt_last();
         } else {
@@ -1500,13 +1495,9 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
         // The energy deposited is the difference between the pre-collision and
         // post-collision energy...
         score = E - p.E();
-
         // ...less the energy of any secondary particles since they will be
         // transported individually later
-        const auto& bank = p.secondary_bank();
-        for (auto it = bank.end() - p.n_bank_second(); it < bank.end(); ++it) {
-          score -= it->E;
-        }
+        score -= p.bank_second_E();
 
         score *= p.wgt_last();
       }

src/track_output.cpp

@@ -31,8 +31,8 @@ int n_tracks_written; //! Number of tracks written
 
 void add_particle_track(Particle& p)
 {
-  p.tracks().emplace_back();
-  p.tracks().back().particle = p.type();
+  auto& track = p.tracks().emplace_back();
+  track.particle = p.type();
 }
 
 void write_particle_track(Particle& p)

src/weight_windows.cpp

@@ -114,7 +114,7 @@ void apply_weight_windows(Particle& p)
     // Create secondaries and divide weight among all particles
     int i_split = std::round(n_split);
     for (int l = 0; l < i_split - 1; l++) {
-      p.create_secondary(weight / n_split, p.u(), p.E(), p.type());
+      p.split(weight / n_split);
     }
     // remaining weight is applied to current particle
     p.wgt() = weight / n_split;

tests/unit_tests/weightwindows/test.py

@@ -4,7 +4,6 @@
 import numpy as np
 import pytest
 from uncertainties import ufloat
-
 import openmc
 import openmc.lib
 from openmc.stats import Discrete, Point
@@ -224,6 +223,47 @@ def test_lower_ww_bounds_shape():
     assert ww.lower_ww_bounds.shape == (2, 3, 4, 1)
 
 
+def test_photon_heating(run_in_tmpdir):
+    water = openmc.Material()
+    water.add_nuclide('H1', 1.0)
+    water.add_nuclide('O16', 2.0)
+    water.set_density('g/cm3', 1.0)
+
+    box = openmc.model.RectangularParallelepiped(
+        -300, 300, -300, 300, -300, 300, boundary_type='reflective')
+    cell = openmc.Cell(region=-box, fill=water)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+
+    mesh = openmc.RegularMesh.from_domain(model.geometry, dimension=(5, 5, 5))
+    wwg = openmc.WeightWindowGenerator(mesh, particle_type='photon')
+    model.settings.weight_window_generators = [wwg]
+
+    space = openmc.stats.Point((0, 0, 0))
+    energy = openmc.stats.delta_function(5e6)
+    model.settings.source = openmc.IndependentSource(
+        space=space, energy=energy, particle='photon')
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 5
+    model.settings.particles = 100
+
+    tally = openmc.Tally()
+    tally.scores = ['heating']
+    tally.filters = [
+        openmc.ParticleFilter(['photon']),
+        openmc.MeshFilter(mesh)
+    ]
+    model.tallies = [tally]
+
+    sp_file = model.run()
+    with openmc.StatePoint(sp_file) as sp:
+        tally_mean = sp.tallies[tally.id].mean
+
+    # these values should be nearly identical
+    assert np.all(tally_mean >= 0)
+
+
 def test_roundtrip(run_in_tmpdir, model, wws):
     model.settings.weight_windows = wws
 
@@ -249,11 +289,11 @@ def test_ww_attrs_python(model):
     # is successful
     wws = openmc.WeightWindows(mesh, lower_bounds, upper_bound_ratio=10.0)
 
-    assert wws.energy_bounds == None
+    assert wws.energy_bounds is None
 
     wwg = openmc.WeightWindowGenerator(mesh)
 
-    assert wwg.energy_bounds == None
+    assert wwg.energy_bounds is None
 
 def test_ww_attrs_capi(run_in_tmpdir, model):
     model.export_to_xml()