Keep the micro cross section caches addressed correctly

include/openmc/particle.h

@@ -125,6 +125,15 @@ class Particle : public ParticleData {
   //! \param[in] ncrystal_xs Thermal scattering xs from NCrystal
   void update_neutron_xs(int i_nuclide, int i_grid = C_NONE, int i_sab = C_NONE,
     double sab_frac = 0.0, double ncrystal_xs = -1.0);
+
+  //! Update the microscopic cross section cache of an element
+  //
+  //! Evaluates the element's cross sections at the particle's energy unless
+  //! the cache already holds them, which is what transport leaves behind for
+  //! the elements present where the particle last collided.
+  //!
+  //! \param[in] i_element Index in data::elements
+  void update_photon_xs(int i_element);
 };
 
 //============================================================================

include/openmc/particle_data.h

@@ -584,11 +584,14 @@ class ParticleData : public GeometryState {
 
   // Cross section caches
   // Microscopic neutron cross sections
-  NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; }
-  const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; }
+  NuclideMicroXS& neutron_xs(int i_nuclide) { return neutron_xs_[i_nuclide]; }
+  const NuclideMicroXS& neutron_xs(int i_nuclide) const
+  {
+    return neutron_xs_[i_nuclide];
+  }
 
   // Microscopic photon cross sections
-  ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; }
+  ElementMicroXS& photon_xs(int i_element) { return photon_xs_[i_element]; }
 
   // Macroscopic cross sections
   MacroXS& macro_xs() { return macro_xs_; }

include/openmc/photon.h

@@ -185,6 +185,14 @@ extern tensor::Tensor<double>
 extern std::unordered_map<std::string, int> element_map;
 extern vector<unique_ptr<PhotonInteraction>> elements;
 
+//! Index in \ref elements of the element each nuclide belongs to
+//
+//! Photon data is tabulated per element rather than per nuclide, so the photon
+//! micro cross section cache is indexed by element while the neutron one is
+//! indexed by nuclide. Anything holding a nuclide index needs this to reach the
+//! photon data. Entries are C_NONE when photon transport is off.
+extern vector<int> nuclide_to_element;
+
 } // namespace data
 
 } // namespace openmc

src/material.cpp

@@ -827,6 +827,8 @@ void Material::calculate_xs(Particle& p) const
 
 void Material::calculate_neutron_xs(Particle& p) const
 {
+  assert(p.type().is_neutron());
+
   // Find energy index on energy grid
   int neutron = ParticleType::neutron().transport_index();
   int i_grid =
@@ -901,6 +903,8 @@ void Material::calculate_neutron_xs(Particle& p) const
 
 void Material::calculate_photon_xs(Particle& p) const
 {
+  assert(p.type().is_photon());
+
   p.macro_xs().coherent = 0.0;
   p.macro_xs().incoherent = 0.0;
   p.macro_xs().photoelectric = 0.0;

src/nuclide.cpp

@@ -1168,6 +1168,14 @@ extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n)
         close_group(group);
         file_close(file_id);
       }
+
+      // Record which element this nuclide belongs to, so that a nuclide index
+      // can reach the photon data it shares with the other isotopes
+      int i_nuclide = data::nuclide_map.at(name);
+      if (data::nuclide_to_element.size() <= i_nuclide) {
+        data::nuclide_to_element.resize(i_nuclide + 1, C_NONE);
+      }
+      data::nuclide_to_element[i_nuclide] = data::element_map.at(element);
     }
   }
   return 0;

src/particle.cpp

@@ -1,5 +1,7 @@
 #include "openmc/particle.h"
 
+#include <cassert>
+
 #include <algorithm> // copy, min
 #include <cmath>     // log, abs
 
@@ -964,6 +966,8 @@ void Particle::write_restart() const
 void Particle::update_neutron_xs(
   int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
 {
+  assert(type().is_neutron());
+
   // Get microscopic cross section cache
   auto& micro = this->neutron_xs(i_nuclide);
 
@@ -982,6 +986,16 @@ void Particle::update_neutron_xs(
   }
 }
 
+void Particle::update_photon_xs(int i_element)
+{
+  assert(type().is_photon());
+
+  // If the cache doesn't match, recalculate micro xs
+  if (this->E() != this->photon_xs(i_element).last_E) {
+    data::elements[i_element]->calculate_xs(*this);
+  }
+}
+
 //==============================================================================
 // Non-method functions
 //==============================================================================

src/photon.cpp

@@ -1,5 +1,7 @@
 #include "openmc/photon.h"
 
+#include <cassert>
+
 #include "openmc/array.h"
 #include "openmc/bremsstrahlung.h"
 #include "openmc/constants.h"
@@ -36,6 +38,7 @@ tensor::Tensor<double> compton_profile_pz;
 
 std::unordered_map<std::string, int> element_map;
 vector<unique_ptr<PhotonInteraction>> elements;
+vector<int> nuclide_to_element;
 
 } // namespace data
 
@@ -710,6 +713,8 @@ void PhotonInteraction::compton_doppler(
 
 void PhotonInteraction::calculate_xs(Particle& p) const
 {
+  assert(p.type().is_photon());
+
   // Perform binary search on the element energy grid in order to determine
   // which points to interpolate between
   int n_grid = energy_.size();
@@ -1116,6 +1121,7 @@ std::pair<double, double> klein_nishina(double alpha, uint64_t* seed)
 void free_memory_photon()
 {
   data::elements.clear();
+  data::nuclide_to_element.clear();
   data::compton_profile_pz.resize({0});
   data::ttb_e_grid.resize({0});
   data::ttb_k_grid.resize({0});

src/tallies/tally_scoring.cpp

@@ -177,6 +177,33 @@ void score_fission_delayed_dg(int i_tally, int d_bin, double score,
   dg_match.bins_[i_bin] = original_bin;
 }
 
+//! Helper function to refresh the micro xs cache of an absent nuclide
+//
+//! A tally binned by nuclide with multiply_density off scores a nuclide the
+//! scoring region need not contain, and transport only evaluates the cross
+//! sections of what is actually there. Each cache is read back solely by
+//! scoring of the particle type it holds cross sections for, so there is
+//! nothing to refresh for any other kind of particle.
+//
+//! \param[in,out] p Particle being scored
+//! \param[in] i_nuclide Index in data::nuclides
+//! \param[in,out] i_log_union Index on the log union grid, or C_NONE if it has
+//!   not been determined yet
+void update_absent_nuclide_xs(Particle& p, int i_nuclide, int& i_log_union)
+{
+  if (p.type().is_neutron()) {
+    // Determine log union grid index
+    if (i_log_union == C_NONE) {
+      int neutron = ParticleType::neutron().transport_index();
+      i_log_union =
+        std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing;
+    }
+    p.update_neutron_xs(i_nuclide, i_log_union);
+  } else if (p.type().is_photon()) {
+    p.update_photon_xs(data::nuclide_to_element[i_nuclide]);
+  }
+}
+
 //! Helper function to retrieve fission q value from a nuclide
 
 double get_nuc_fission_q(const Nuclide& nuc, const Particle& p, int score_bin)
@@ -586,6 +613,13 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
   // Get the pre-collision energy of the particle.
   auto E = p.E_last();
 
+  // Photon data is tabulated per element, so the cache holding photon micro
+  // cross sections is indexed by element while i_nuclide is not. Resolve it
+  // once here; C_NONE whenever no score below can ask for one.
+  int i_element = (i_nuclide >= 0 && p.type().is_photon())
+                    ? data::nuclide_to_element[i_nuclide]
+                    : C_NONE;
+
   for (auto i = 0; i < tally.scores_.size(); ++i) {
     auto score_bin = tally.scores_[i];
     auto score_index = start_index + i;
@@ -601,7 +635,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
         if (p.type().is_neutron()) {
           score = p.neutron_xs(i_nuclide).total * atom_density * flux;
         } else if (p.type().is_photon()) {
-          score = p.photon_xs(i_nuclide).total * atom_density * flux;
+          score = p.photon_xs(i_element).total * atom_density * flux;
         }
       } else {
         score = p.macro_xs().total * flux;
@@ -625,7 +659,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
           const auto& micro = p.neutron_xs(i_nuclide);
           score = (micro.total - micro.absorption) * atom_density * flux;
         } else {
-          const auto& micro = p.photon_xs(i_nuclide);
+          const auto& micro = p.photon_xs(i_element);
           score = (micro.coherent + micro.incoherent) * atom_density * flux;
         }
       } else {
@@ -645,7 +679,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
         if (p.type().is_neutron()) {
           score = p.neutron_xs(i_nuclide).absorption * atom_density * flux;
         } else {
-          const auto& xs = p.photon_xs(i_nuclide);
+          const auto& xs = p.photon_xs(i_element);
           score =
             (xs.total - xs.coherent - xs.incoherent) * atom_density * flux;
         }
@@ -1052,7 +1086,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
         continue;
 
       if (i_nuclide >= 0) {
-        const auto& micro = p.photon_xs(i_nuclide);
+        const auto& micro = p.photon_xs(i_element);
         double xs = (score_bin == COHERENT)        ? micro.coherent
                     : (score_bin == INCOHERENT)    ? micro.incoherent
                     : (score_bin == PHOTOELECTRIC) ? micro.photoelectric
@@ -2447,14 +2481,8 @@ void score_tracklength_tally_general(
             if (!tally.multiply_density())
               atom_density = 1.0;
           } else if (!tally.multiply_density()) {
-            // Determine log union grid index
-            if (i_log_union == C_NONE) {
-              int neutron = ParticleType::neutron().transport_index();
-              i_log_union = std::log(p.E() / data::energy_min[neutron]) /
-                            simulation::log_spacing;
-            }
             // Update micro xs cache
-            p.update_neutron_xs(i_nuclide, i_log_union);
+            update_absent_nuclide_xs(p, i_nuclide, i_log_union);
             atom_density = 1.0;
           }
         }
@@ -2577,14 +2605,8 @@ void score_collision_tally(Particle& p)
             if (!tally.multiply_density())
               atom_density = 1.0;
           } else if (!tally.multiply_density()) {
-            // Determine log union grid index
-            if (i_log_union == C_NONE) {
-              int neutron = ParticleType::neutron().transport_index();
-              i_log_union = std::log(p.E() / data::energy_min[neutron]) /
-                            simulation::log_spacing;
-            }
             // Update micro xs cache
-            p.update_neutron_xs(i_nuclide, i_log_union);
+            update_absent_nuclide_xs(p, i_nuclide, i_log_union);
             atom_density = 1.0;
           }
         }

tests/unit_tests/test_photon_micro_xs.py

@@ -0,0 +1,88 @@
+"""Photon micro cross section tallies are looked up per element.
+
+Photon data does not distinguish isotopes, so the cache holding photon micro
+cross sections is indexed by element while a tally's nuclide bins are indexed by
+nuclide. Two isotopes of the same element present at the same atom density must
+therefore score the same microscopic cross section, whatever that cross section
+happens to be -- which makes this check independent of the data library.
+"""
+
+import openmc
+import pytest
+
+
+def test_isotopes_of_one_element_score_alike(run_in_tmpdir):
+    openmc.reset_auto_ids()
+
+    mat = openmc.Material()
+    mat.add_nuclide('U235', 1.0)
+    mat.add_nuclide('U238', 1.0)
+    mat.set_density('g/cm3', 10.0)
+
+    sphere = openmc.Sphere(r=5.0, boundary_type='vacuum')
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([openmc.Cell(fill=mat, region=-sphere)])
+    model.settings.run_mode = 'fixed source'
+    model.settings.photon_transport = True
+    model.settings.particles = 200
+    model.settings.batches = 2
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point(),
+        energy=openmc.stats.Discrete([1.0e6], [1.0]),
+        particle='photon')
+
+    tally = openmc.Tally()
+    tally.nuclides = ['U235', 'U238']
+    tally.scores = ['total']
+    model.tallies = openmc.Tallies([tally])
+
+    model.run(apply_tally_results=True)
+
+    u235, u238 = tally.mean.ravel()
+    assert u235 > 0.0
+    assert u238 == pytest.approx(u235, rel=1e-12)
+
+
+def test_nuclide_absent_from_the_material(run_in_tmpdir):
+    """A tally may name a nuclide no material contains.
+
+    Doing so adds it to the global nuclide list without adding an element, so
+    its index reaches past the end of the element-indexed photon cache unless
+    the lookup goes through the nuclide-to-element mapping.
+
+    multiply_density has to be off for the score to say anything: an absent
+    nuclide has an atom density of zero, which multiplies away whatever came
+    out of the cache. With it off the score is the microscopic cross section
+    itself, and it is also the branch that refreshes the cache for a nuclide
+    the scoring region does not contain.
+    """
+    openmc.reset_auto_ids()
+
+    mat = openmc.Material()
+    mat.add_nuclide('U235', 1.0)
+    mat.set_density('g/cm3', 10.0)
+
+    sphere = openmc.Sphere(r=5.0, boundary_type='vacuum')
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([openmc.Cell(fill=mat, region=-sphere)])
+    model.settings.run_mode = 'fixed source'
+    model.settings.photon_transport = True
+    model.settings.particles = 200
+    model.settings.batches = 2
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point(),
+        energy=openmc.stats.Discrete([1.0e6], [1.0]),
+        particle='photon')
+
+    tally = openmc.Tally()
+    tally.nuclides = ['U235', 'U238']  # U238 is not in the material
+    tally.scores = ['total']
+    tally.multiply_density = False
+    model.tallies = openmc.Tallies([tally])
+
+    model.run(apply_tally_results=True)
+
+    # Same element as U235, so the same microscopic cross section
+    u235, u238 = tally.mean.ravel()
+    assert u235 > 0.0
+    assert u238 == pytest.approx(u235, rel=1e-12)