do heating tallies for both neutron and photon

include/openmc/constants.h

@@ -230,7 +230,7 @@ constexpr int N_XD    {204};
 constexpr int N_XT    {205};
 constexpr int N_X3HE  {206};
 constexpr int N_XA    {207};
-constexpr int HEATING {301};
+constexpr int NEUTRON_HEATING {301};
 constexpr int DAMAGE_ENERGY {444};
 constexpr int COHERENT {502};
 constexpr int INCOHERENT {504};
@@ -366,6 +366,7 @@ constexpr int SCORE_INVERSE_VELOCITY   {-13}; // flux-weighted inverse velocity
 constexpr int SCORE_FISS_Q_PROMPT      {-14}; // prompt fission Q-value
 constexpr int SCORE_FISS_Q_RECOV       {-15}; // recoverable fission Q-value
 constexpr int SCORE_DECAY_RATE         {-16}; // delayed neutron precursor decay rate
+constexpr int SCORE_HEATING            {-17}; // nuclear heating (neutron or photon)
 
 // Tally map bin finding
 constexpr int NO_BIN_FOUND {-1};

openmc/capi/tally.py

@@ -89,7 +89,7 @@ _SCORES = {
     -5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission',
     -9: 'current', -10: 'events', -11: 'delayed-nu-fission',
     -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt',
-    -15: 'fission-q-recoverable', -16: 'decay-rate'
+    -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'heating'
 }
 _ESTIMATORS = {
     1: 'analog', 2: 'tracklength', 3: 'collision'

src/output.cpp

@@ -618,6 +618,7 @@ const std::unordered_map<int, const char*> score_names = {
   {SCORE_FISS_Q_PROMPT,      "Prompt fission power"},
   {SCORE_FISS_Q_RECOV,       "Recoverable fission power"},
   {SCORE_CURRENT,            "Current"},
+  {SCORE_HEATING,            "Heating"},
 };
 
 //! Create an ASCII output file showing all tally results.

src/photon.cpp

@@ -97,9 +97,13 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
   close_group(rgroup);
 
   // Read heating
-  rgroup = open_group(group, "heating");
-  read_dataset(rgroup, "xs", heating_);
-  close_group(rgroup);
+  if (object_exists(group, "heating")) {
+    rgroup = open_group(group, "heating");
+    read_dataset(rgroup, "xs", heating_);
+    close_group(rgroup);
+  } else {
+    heating_ = xt::zeros_like(energy_);
+  }
 
   // Read subshell photoionization cross section and atomic relaxation data
   rgroup = open_group(group, "subshells");

src/reaction.cpp

@@ -121,6 +121,8 @@ std::string reaction_name(int mt)
     return  "fission-q-prompt";
   } else if (mt == SCORE_FISS_Q_RECOV) {
     return  "fission-q-recoverable";
+  } else if (mt == SCORE_HEATING) {
+    return  "heating";
 
   // Normal ENDF-based reactions
   } else if (mt == TOTAL_XS) {

src/tallies/tally.cpp

@@ -109,6 +109,9 @@ score_str_to_int(std::string score_str)
   if (score_str == "fission-q-recoverable")
     return SCORE_FISS_Q_RECOV;
 
+  if (score_str == "heating")
+    return SCORE_HEATING;
+
   if (score_str == "current")
     return SCORE_CURRENT;
 
@@ -207,8 +210,6 @@ score_str_to_int(std::string score_str)
     return N_X3HE;
   if (score_str == "(n,Xa)" || score_str == "He4-production")
     return N_XA;
-  if (score_str == "heating")
-    return HEATING;
   if (score_str == "damage-energy")
     return DAMAGE_ENERGY;
 

src/tallies/tally_scoring.cpp

@@ -7,6 +7,7 @@
 #include "openmc/material.h"
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
+#include "openmc/photon.h"
 #include "openmc/reaction_product.h"
 #include "openmc/search.h"
 #include "openmc/settings.h"
@@ -1148,6 +1149,46 @@ score_general_ce(Particle* p, int i_tally, int start_index,
       break;
 
 
+    case SCORE_HEATING:
+      if (p->type_ == Particle::Type::neutron) {
+        score_bin = NEUTRON_HEATING;
+        // No break here, continue to run the default neutron MT scoring
+      } else if (p->type_ == Particle::Type::photon) {
+        if (tally.estimator_ != ESTIMATOR_ANALOG) {
+          // Calculate photon heating cross section on-the-fly
+          score = 0.;
+          if (i_nuclide >= 0) {
+            // Find the element corresponding to the nuclide
+            auto name = data::nuclides[i_nuclide]->name_;
+            int pos = name.find_first_of("0123456789");
+            std::string element = name.substr(0, pos);
+            int i_element = data::element_map[element];
+            auto& heating {data::elements[i_element].heating_};
+            auto i_grid = simulation::micro_photon_xs[i_element].index_grid;
+            auto f = simulation::micro_photon_xs[i_element].interp_factor;
+            score = std::exp(heating(i_grid) + f * (heating(i_grid+1) -
+              heating(i_grid))) * atom_density * flux;
+          } else {
+            if (p->material_ != MATERIAL_VOID) {
+              const Material& material {*model::materials[p->material_]};
+              for (auto i = 0; i < material.nuclide_.size(); ++i) {
+                auto i_element = material.element_[i];
+                auto atom_density = material.atom_density_(i);
+                auto& heating {data::elements[i_element].heating_};
+                auto i_grid = simulation::micro_photon_xs[i_element].index_grid;
+                auto f = simulation::micro_photon_xs[i_element].interp_factor;
+                score += std::exp(heating(i_grid) + f * (heating(i_grid+1) -
+                  heating(i_grid))) * atom_density * flux;
+              }
+            }
+          }
+        }
+        break;
+      } else {
+        // Nuclear heating of any other particles not implmented
+        break;
+      }
+
     default:
       if (tally.estimator_ == ESTIMATOR_ANALOG) {
         // Any other score is assumed to be a MT number. Thus, we just need