Merge pull request #1449 from paulromano/photon-energydep-fixes

Fix/improve energy deposition for coupled neutron-photon calculations

docs/source/io_formats/settings.rst

@@ -96,6 +96,17 @@ When the DAGMC mode is enabled, the OpenMC geometry will be read from the file
 ``dagmc.h5m``. If a :ref:`geometry.xml <io_geometry>` file is present with
 ``dagmc`` set to ``true``, it will be ignored.
 
+----------------------------
+``<delayed_photon_scaling>``
+----------------------------
+
+Determines whether to scale the fission photon yield to account for delayed
+photon energy. The photon yields are scaled as (EGP + EGD)/EGP where EGP and EGD
+are the prompt and delayed photon components of energy release, respectively,
+from MF=1, MT=458 on an ENDF evaluation.
+
+  *Default*: true
+
 --------------------------------
 ``<electron_treatment>`` Element
 --------------------------------

docs/source/publications.rst

@@ -71,6 +71,10 @@ Coupling and Multi-physics
   Coupling of OpenMC and Nek5000 within the MOOSE Framework," *Proc. PHYSOR*,
   Cancun, Mexico, Apr. 22-26 (2018).
 
+- Sterling Harper, Kord Smith, and Benoit Forget, "Faster Monte Carlo
+  multiphysics using temperature derivatives," *Proc. PHYSOR*, Cancun, Mexico,
+  Apr. 22-26 (2018).
+
 - Jun Chen, Liangzhi Cao, Chuanqi Zhao, and Zhouyu Liu, "`Development of
   Subchannel Code SUBSC for high-fidelity multi-physics coupling application
   <https://doi.org/10.1016/j.egypro.2017.08.121>`_", *Energy Procedia*, **127**,

include/openmc/nuclide.h

@@ -97,7 +97,7 @@ public:
   std::vector<int> index_inelastic_scatter_;
 
 private:
-  void create_derived();
+  void create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons);
 
   static int XS_TOTAL;
   static int XS_ABSORPTION;

include/openmc/settings.h

@@ -29,6 +29,7 @@ extern bool confidence_intervals;     //!< use confidence intervals for results?
 extern bool create_fission_neutrons;  //!< create fission neutrons (fixed source)?
 extern "C" bool cmfd_run;             //!< is a CMFD run?
 extern "C" bool dagmc;                //!< indicator of DAGMC geometry
+extern bool delayed_photon_scaling;   //!< Scale fission photon yield to include delayed
 extern "C" bool entropy_on;           //!< calculate Shannon entropy?
 extern bool legendre_to_tabular;      //!< convert Legendre distributions to tabular?
 extern bool material_cell_offsets;    //!< create material cells offsets?

openmc/settings.py

@@ -40,6 +40,10 @@ class Settings(object):
         below which particle type will be killed.
     dagmc : bool
         Indicate that a CAD-based DAGMC geometry will be used.
+    delayed_photon_scaling : bool
+        Indicate whether to scale the fission photon yield by (EGP + EGD)/EGP
+        where EGP is the energy release of prompt photons and EGD is the energy
+        release of delayed photons.
     electron_treatment : {'led', 'ttb'}
         Whether to deposit all energy from electrons locally ('led') or create
         secondary bremsstrahlung photons ('ttb').
@@ -219,6 +223,7 @@ class Settings(object):
             VolumeCalculation, 'volume calculations')
 
         self._create_fission_neutrons = None
+        self._delayed_photon_scaling = None
         self._material_cell_offsets = None
         self._log_grid_bins = None
 
@@ -356,6 +361,10 @@ class Settings(object):
     def create_fission_neutrons(self):
         return self._create_fission_neutrons
 
+    @property
+    def delayed_photon_scaling(self):
+        return self._delayed_photon_scaling
+
     @property
     def material_cell_offsets(self):
         return self._material_cell_offsets
@@ -691,6 +700,11 @@ class Settings(object):
                       create_fission_neutrons, bool)
         self._create_fission_neutrons = create_fission_neutrons
 
+    @delayed_photon_scaling.setter
+    def delayed_photon_scaling(self, value):
+        cv.check_type('delayed photon scaling', value, bool)
+        self._delayed_photon_scaling = value
+
     @material_cell_offsets.setter
     def material_cell_offsets(self, value):
         cv.check_type('material cell offsets', value, bool)
@@ -930,6 +944,11 @@ class Settings(object):
             elem = ET.SubElement(root, "create_fission_neutrons")
             elem.text = str(self._create_fission_neutrons).lower()
 
+    def _create_delayed_photon_scaling_subelement(self, root):
+        if self._delayed_photon_scaling is not None:
+            elem = ET.SubElement(root, "delayed_photon_scaling")
+            elem.text = str(self._delayed_photon_scaling).lower()
+
     def _create_material_cell_offsets_subelement(self, root):
         if self._material_cell_offsets is not None:
             elem = ET.SubElement(root, "material_cell_offsets")
@@ -1166,6 +1185,11 @@ class Settings(object):
         if text is not None:
             self.create_fission_neutrons = text in ('true', '1')
 
+    def _delayed_photon_scaling_from_xml_element(self, root):
+        text = get_text(root, 'delayed_photon_scaling')
+        if text is not None:
+            self.delayed_photon_scaling = text in ('true', '1')
+
     def _material_cell_offsets_from_xml_element(self, root):
         text = get_text(root, 'material_cell_offsets')
         if text is not None:
@@ -1225,6 +1249,7 @@ class Settings(object):
         self._create_resonance_scattering_subelement(root_element)
         self._create_volume_calcs_subelement(root_element)
         self._create_create_fission_neutrons_subelement(root_element)
+        self._create_delayed_photon_scaling_subelement(root_element)
         self._create_material_cell_offsets_subelement(root_element)
         self._create_log_grid_bins_subelement(root_element)
         self._create_dagmc_subelement(root_element)
@@ -1291,6 +1316,7 @@ class Settings(object):
         settings._ufs_mesh_from_xml_element(root)
         settings._resonance_scattering_from_xml_element(root)
         settings._create_fission_neutrons_from_xml_element(root)
+        settings._delayed_photon_scaling_from_xml_element(root)
         settings._material_cell_offsets_from_xml_element(root)
         settings._log_grid_bins_from_xml_element(root)
         settings._dagmc_from_xml_element(root)

src/error.cpp

@@ -57,7 +57,7 @@ void output(const std::string& message, std::ostream& out, int indent)
   while (i_start < length) {
     if (length - i_start < line_len) {
       // Remainder of message will fit on line
-      out << message.substr(i_start) << '\n';
+      out << message.substr(i_start) << std::endl;
       break;
 
     } else {

src/finalize.cpp

@@ -69,6 +69,7 @@ int openmc_finalize()
   settings::check_overlaps = false;
   settings::confidence_intervals = false;
   settings::create_fission_neutrons = true;
+  settings::delayed_photon_scaling = true;
   settings::electron_treatment = ELECTRON_LED;
   settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0};
   settings::entropy_on = false;

src/nuclide.cpp

@@ -259,17 +259,24 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature, int i_nucl
   }
 
   // Read fission energy release data if present
+  std::unique_ptr<Function1D> prompt_photons;
+  std::unique_ptr<Function1D> delayed_photons;
   if (object_exists(group, "fission_energy_release")) {
     hid_t fer_group = open_group(group, "fission_energy_release");
     fission_q_prompt_ = read_function(fer_group, "q_prompt");
     fission_q_recov_ = read_function(fer_group, "q_recoverable");
+
+    // We need prompt/delayed photon energy release for scaling fission photon
+    // production
+    prompt_photons = read_function(fer_group, "prompt_photons");
+    delayed_photons = read_function(fer_group, "delayed_photons");
     close_group(fer_group);
   }
 
-  this->create_derived();
+  this->create_derived(prompt_photons.get(), delayed_photons.get());
 }
 
-void Nuclide::create_derived()
+void Nuclide::create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons)
 {
   for (const auto& grid : grid_) {
     // Allocate and initialize cross section
@@ -294,7 +301,21 @@ void Nuclide::create_derived()
           auto pprod = xt::view(xs_[t], xt::range(j, j+n), XS_PHOTON_PROD);
           for (int k = 0; k < n; ++k) {
             double E = grid_[t].energy[k+j];
-            pprod[k] += xs[k] * (*p.yield_)(E);
+
+            // For fission, artificially increase the photon yield to account
+            // for delayed photons
+            double f = 1.0;
+            if (settings::delayed_photon_scaling) {
+              if (is_fission(rx->mt_)) {
+                if (prompt_photons && delayed_photons) {
+                  double energy_prompt = (*prompt_photons)(E);
+                  double energy_delayed = (*delayed_photons)(E);
+                  f = (energy_prompt + energy_delayed)/(energy_prompt);
+                }
+              }
+            }
+
+            pprod[k] += f * xs[k] * (*p.yield_)(E);
           }
         }
       }

src/output.cpp

@@ -92,7 +92,7 @@ void title()
   // Write number of OpenMP threads
   std::cout << "    OpenMP Threads | " << omp_get_max_threads() << '\n';
 #endif
-  std::cout << '\n';
+  std::cout << std::endl;
 }
 
 //==============================================================================
@@ -126,7 +126,7 @@ header(const char* msg, int level) {
 
   // Print header based on verbosity level.
   if (settings::verbosity >= level)
-    std::cout << '\n' << out << "\n\n";
+    std::cout << '\n' << out << "\n" << std::endl;
 }
 
 //==============================================================================

src/relaxng/settings.rnc

@@ -14,6 +14,8 @@ element settings {
     (element energy_positron { xsd:double } | attribute energy_positron { xsd:double })?
   }? &
 
+  element delayed_photon_scaling { xsd:boolean }? &
+
   element electron_treatment { ( "led" | "ttb" ) }? &
 
   element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" | "material-union" | "union" ) }? &

src/relaxng/settings.rng

@@ -82,6 +82,11 @@
         </interleave>
       </element>
     </optional>
+    <optional>
+      <element name="delayed_photon_scaling">
+        <data type="boolean"/>
+      </element>
+    </optional>
     <optional>
       <element name="electron_treatment">
         <choice>

src/settings.cpp

@@ -44,6 +44,7 @@ bool cmfd_run                {false};
 bool confidence_intervals    {false};
 bool create_fission_neutrons {true};
 bool dagmc                   {false};
+bool delayed_photon_scaling  {true};
 bool entropy_on              {false};
 bool legendre_to_tabular     {true};
 bool material_cell_offsets   {true};
@@ -778,6 +779,11 @@ void read_settings_xml()
     }
   }
 
+  // Check whether to scale fission photon yields
+  if (check_for_node(root, "delayed_photon_scaling")) {
+    delayed_photon_scaling = get_node_value_bool(root, "delayed_photon_scaling");
+  }
+
   // Check whether material cell offsets should be generated
   if (check_for_node(root, "material_cell_offsets")) {
     material_cell_offsets = get_node_value_bool(root, "material_cell_offsets");

src/tallies/tally_scoring.cpp

@@ -1267,6 +1267,10 @@ score_general_ce(Particle* p, int i_tally, int start_index,
       break;
 
     default:
+      // The default block is really only meant for redundant neutron reactions
+      // (e.g. 444, 901)
+      if (p->type_ != Particle::Type::neutron) continue;
+
       if (tally.estimator_ == ESTIMATOR_ANALOG) {
         // Any other score is assumed to be a MT number. Thus, we just need
         // to check if it matches the MT number of the event

tests/regression_tests/photon_production/inputs_true.dat

@@ -1,16 +1,16 @@
 <?xml version='1.0' encoding='utf-8'?>
 <geometry>
-  <cell id="13" material="void" region="-9 10 -11" universe="9" />
-  <cell id="14" material="13" region="-9 11 -12" universe="9" />
-  <cell id="15" material="void" region="~(-9 10 -12)" universe="9" />
-  <surface boundary="vacuum" coeffs="0.0 0.0 1.0" id="9" type="x-cylinder" />
-  <surface boundary="vacuum" coeffs="-1.0" id="10" type="x-plane" />
-  <surface coeffs="1.0" id="11" type="x-plane" />
-  <surface boundary="vacuum" coeffs="1000000000.0" id="12" type="x-plane" />
+  <cell id="1" material="void" region="-1 2 -3" universe="1" />
+  <cell id="2" material="1" region="-1 3 -4" universe="1" />
+  <cell id="3" material="void" region="~(-1 2 -4)" universe="1" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 1.0" id="1" type="x-cylinder" />
+  <surface boundary="vacuum" coeffs="-1.0" id="2" type="x-plane" />
+  <surface coeffs="1.0" id="3" type="x-plane" />
+  <surface boundary="vacuum" coeffs="1000000000.0" id="4" type="x-plane" />
 </geometry>
 <?xml version='1.0' encoding='utf-8'?>
 <materials>
-  <material id="13">
+  <material id="1">
     <density units="g/cm3" value="2.6989" />
     <nuclide ao="1.0" name="Al27" />
   </material>
@@ -38,7 +38,7 @@
 <?xml version='1.0' encoding='utf-8'?>
 <tallies>
   <filter id="1" type="surface">
-    <bins>9</bins>
+    <bins>1</bins>
   </filter>
   <filter id="2" type="particle">
     <bins>photon</bins>

tests/regression_tests/photon_production/test.py

@@ -1,90 +1,72 @@
-from math import pi
-
-import numpy as np
 import openmc
+import pytest
 
 from tests.testing_harness import PyAPITestHarness
 
 
-class SourceTestHarness(PyAPITestHarness):
-    def _build_inputs(self):
-        mat = openmc.Material()
-        mat.set_density('g/cm3', 2.6989)
-        mat.add_nuclide('Al27', 1.0)
-        materials = openmc.Materials([mat])
-        materials.export_to_xml()
+@pytest.fixture
+def model():
+    model = openmc.model.Model()
+    mat = openmc.Material()
+    mat.set_density('g/cm3', 2.6989)
+    mat.add_nuclide('Al27', 1.0)
+    model.materials.append(mat)
 
-        cyl = openmc.XCylinder(r=1.0, boundary_type='vacuum')
-        x_plane_left = openmc.XPlane(-1.0, 'vacuum')
-        x_plane_center = openmc.XPlane(1.0)
-        x_plane_right = openmc.XPlane(1.0e9, 'vacuum')
+    cyl = openmc.XCylinder(r=1.0, boundary_type='vacuum')
+    x_plane_left = openmc.XPlane(-1.0, 'vacuum')
+    x_plane_center = openmc.XPlane(1.0)
+    x_plane_right = openmc.XPlane(1.0e9, 'vacuum')
 
-        inner_cyl_left = openmc.Cell()
-        inner_cyl_right = openmc.Cell()
-        outer_cyl = openmc.Cell()
+    inner_cyl_left = openmc.Cell()
+    inner_cyl_right = openmc.Cell()
+    outer_cyl = openmc.Cell()
 
-        inner_cyl_left.region = -cyl & +x_plane_left & -x_plane_center
-        inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right
-        outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right)
-        inner_cyl_right.fill = mat
-        geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl])
-        geometry.export_to_xml()
+    inner_cyl_left.region = -cyl & +x_plane_left & -x_plane_center
+    inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right
+    outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right)
+    inner_cyl_right.fill = mat
+    model.geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl])
 
-        source = openmc.Source()
-        source.space = openmc.stats.Point((0, 0, 0))
-        source.angle = openmc.stats.Monodirectional()
-        source.energy = openmc.stats.Discrete([14.0e6], [1.0])
-        source.particle = 'neutron'
+    source = openmc.Source()
+    source.space = openmc.stats.Point((0, 0, 0))
+    source.angle = openmc.stats.Monodirectional()
+    source.energy = openmc.stats.Discrete([14.0e6], [1.0])
+    source.particle = 'neutron'
 
-        settings = openmc.Settings()
-        settings.particles = 10000
-        settings.run_mode = 'fixed source'
-        settings.batches = 1
-        settings.photon_transport = True
-        settings.electron_treatment = 'ttb'
-        settings.cutoff = {'energy_photon' : 1000.0}
-        settings.source = source
-        settings.export_to_xml()
+    model.settings.particles = 10000
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.photon_transport = True
+    model.settings.electron_treatment = 'ttb'
+    model.settings.cutoff = {'energy_photon' : 1000.0}
+    model.settings.source = source
 
-        surface_filter = openmc.SurfaceFilter(cyl)
-        particle_filter = openmc.ParticleFilter('photon')
-        current_tally = openmc.Tally()
-        current_tally.filters = [surface_filter, particle_filter]
-        current_tally.scores = ['current']
-        tally_tracklength = openmc.Tally()
-        tally_tracklength.filters = [particle_filter]
-        tally_tracklength.scores = ['total', 'heating']
-        tally_tracklength.nuclides = ['Al27', 'total']
-        tally_tracklength.estimator = 'tracklength'
-        tally_collision = openmc.Tally()
-        tally_collision.filters = [particle_filter]
-        tally_collision.scores = ['total', 'heating']
-        tally_collision.nuclides = ['Al27', 'total']
-        tally_collision.estimator = 'collision'
-        tally_analog = openmc.Tally()
-        tally_analog.filters = [particle_filter]
-        tally_analog.scores = ['total', 'heating']
-        tally_analog.nuclides = ['Al27', 'total']
-        tally_analog.estimator = 'analog'
-        tallies = openmc.Tallies([current_tally, tally_tracklength,
-                                  tally_collision, tally_analog])
-        tallies.export_to_xml()
+    surface_filter = openmc.SurfaceFilter(cyl)
+    particle_filter = openmc.ParticleFilter('photon')
+    current_tally = openmc.Tally()
+    current_tally.filters = [surface_filter, particle_filter]
+    current_tally.scores = ['current']
+    tally_tracklength = openmc.Tally()
+    tally_tracklength.filters = [particle_filter]
+    tally_tracklength.scores = ['total', 'heating']
+    tally_tracklength.nuclides = ['Al27', 'total']
+    tally_tracklength.estimator = 'tracklength'
+    tally_collision = openmc.Tally()
+    tally_collision.filters = [particle_filter]
+    tally_collision.scores = ['total', 'heating']
+    tally_collision.nuclides = ['Al27', 'total']
+    tally_collision.estimator = 'collision'
+    tally_analog = openmc.Tally()
+    tally_analog.filters = [particle_filter]
+    tally_analog.scores = ['total', 'heating']
+    tally_analog.nuclides = ['Al27', 'total']
+    tally_analog.estimator = 'analog'
+    model.tallies.extend([current_tally, tally_tracklength,
+                          tally_collision, tally_analog])
 
-    def _get_results(self):
-        with openmc.StatePoint(self._sp_name) as sp:
-            outstr = ''
-            for i, tally_ind in enumerate(sp.tallies):
-                tally = sp.tallies[tally_ind]
-                results = np.zeros((tally.sum.size * 2, ))
-                results[0::2] = tally.sum.ravel()
-                results[1::2] = tally.sum_sq.ravel()
-                results = ['{0:12.6E}'.format(x) for x in results]
-
-                outstr += 'tally {}:\n'.format(i + 1)
-                outstr += '\n'.join(results) + '\n'
-            return outstr
+    return model
 
 
-def test_photon_production():
-    harness = SourceTestHarness('statepoint.1.h5')
+def test_photon_production(model):
+    harness = PyAPITestHarness('statepoint.1.h5', model)
     harness.main()

tests/regression_tests/photon_production_fission/inputs_true.dat

@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="1" region="-1" universe="1" />
+  <surface boundary="reflective" coeffs="0.0 0.0 0.0 100.0" id="1" type="sphere" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material depletable="true" id="1">
+    <density units="g/cm3" value="10.0" />
+    <nuclide ao="1.0" name="U235" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>eigenvalue</run_mode>
+  <particles>1000</particles>
+  <batches>5</batches>
+  <inactive>2</inactive>
+  <source strength="1.0">
+    <space type="point">
+      <parameters>0 0 0</parameters>
+    </space>
+  </source>
+  <photon_transport>true</photon_transport>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <filter id="1" type="particle">
+    <bins>neutron photon</bins>
+  </filter>
+  <tally id="1">
+    <filters>1</filters>
+    <nuclides>U235 total</nuclides>
+    <scores>fission heating heating-local</scores>
+    <estimator>tracklength</estimator>
+  </tally>
+  <tally id="2">
+    <filters>1</filters>
+    <nuclides>U235 total</nuclides>
+    <scores>fission heating heating-local</scores>
+    <estimator>collision</estimator>
+  </tally>
+  <tally id="3">
+    <filters>1</filters>
+    <nuclides>U235 total</nuclides>
+    <scores>fission heating heating-local</scores>
+    <estimator>analog</estimator>
+  </tally>
+</tallies>

tests/regression_tests/photon_production_fission/results_true.dat

@@ -0,0 +1,77 @@
+k-combined:
+2.223956E+00 2.598313E-03
+tally 1:
+2.600054E+00
+2.253433E+00
+4.390171E+08
+6.424554E+16
+0.000000E+00
+0.000000E+00
+2.600054E+00
+2.253433E+00
+4.390171E+08
+6.424554E+16
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.764667E+07
+7.568912E+14
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.764667E+07
+7.568912E+14
+0.000000E+00
+0.000000E+00
+tally 2:
+2.602512E+00
+2.258176E+00
+4.394176E+08
+6.437709E+16
+0.000000E+00
+0.000000E+00
+2.602512E+00
+2.258176E+00
+4.394176E+08
+6.437709E+16
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.749088E+07
+7.519602E+14
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.749088E+07
+7.519602E+14
+0.000000E+00
+0.000000E+00
+tally 3:
+2.663535E+00
+2.364845E+00
+4.394176E+08
+6.437709E+16
+0.000000E+00
+0.000000E+00
+2.663535E+00
+2.364845E+00
+4.394176E+08
+6.437709E+16
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.729796E+07
+7.458273E+14
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.743827E+07
+7.502589E+14
+0.000000E+00
+0.000000E+00

tests/regression_tests/photon_production_fission/test.py

@@ -0,0 +1,48 @@
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def model():
+    model = openmc.model.Model()
+    mat = openmc.Material()
+    mat.set_density('g/cm3', 10.0)
+    mat.add_nuclide('U235', 1.0)
+    model.materials.append(mat)
+
+    sph = openmc.Sphere(r=100.0, boundary_type='reflective')
+    cell = openmc.Cell(fill=mat, region=-sph)
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings.particles = 1000
+    model.settings.batches = 5
+    model.settings.inactive = 2
+    model.settings.photon_transport = True
+    model.settings.source = openmc.Source(space=openmc.stats.Point((0, 0, 0)))
+
+    particle_filter = openmc.ParticleFilter(['neutron', 'photon'])
+    tally_tracklength = openmc.Tally()
+    tally_tracklength.filters = [particle_filter]
+    tally_tracklength.scores = ['fission', 'heating', 'heating-local']
+    tally_tracklength.nuclides = ['U235', 'total']
+    tally_tracklength.estimator = 'tracklength'
+    tally_collision = openmc.Tally()
+    tally_collision.filters = [particle_filter]
+    tally_collision.scores = ['fission', 'heating', 'heating-local']
+    tally_collision.nuclides = ['U235', 'total']
+    tally_collision.estimator = 'collision'
+    tally_analog = openmc.Tally()
+    tally_analog.filters = [particle_filter]
+    tally_analog.scores = ['fission', 'heating', 'heating-local']
+    tally_analog.nuclides = ['U235', 'total']
+    tally_analog.estimator = 'analog'
+    model.tallies.extend([tally_tracklength, tally_collision, tally_analog])
+
+    return model
+
+
+def test_photon_production(model):
+    harness = PyAPITestHarness('statepoint.5.h5', model)
+    harness.main()