Merge pull request #2433 from paulromano/photon-from-ace-fix

Handle zero photon cross sections in `IncidentPhoton.from_ace`

openmc/data/photon.py

@@ -514,13 +514,13 @@ def from_ace(cls, ace_or_filename):
 
         # Read each reaction
         data = cls(Z)
-        for mt in (502, 504, 515, 522, 525):
+        for mt in (502, 504, 517, 522, 525):
             data.reactions[mt] = PhotonReaction.from_ace(ace, mt)
 
         # Get heating cross sections [eV-barn] from factors [eV per collision]
         # by multiplying with total xs
         data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in
-                                         (502, 504, 515, 522)])
+                                         (502, 504, 517, 522)])
 
         # Compton profiles
         n_shell = ace.nxs[5]
@@ -1000,7 +1000,7 @@ def from_ace(cls, ace, mt):
         elif mt == 504:
             # Incoherent scattering
             idx = ace.jxs[1] + n
-        elif mt == 515:
+        elif mt == 517:
             # Pair production
             idx = ace.jxs[1] + 4*n
         elif mt == 522:
@@ -1021,6 +1021,9 @@ def from_ace(cls, ace, mt):
         else:
             nonzero = (xs != 0.0)
             xs[nonzero] = np.exp(xs[nonzero])
+
+            # Replace zero elements to small non-zero to enable log-log
+            xs[~nonzero] = np.exp(-500.0)
         rx.xs = Tabulated1D(energy, xs, [n], [5])
 
         # Get form factors for incoherent/coherent scattering

src/photon.cpp

@@ -91,9 +91,13 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   close_group(rgroup);
 
   // Read pair production
-  rgroup = open_group(group, "pair_production_electron");
-  read_dataset(rgroup, "xs", pair_production_electron_);
-  close_group(rgroup);
+  if (object_exists(group, "pair_production_electron")) {
+    rgroup = open_group(group, "pair_production_electron");
+    read_dataset(rgroup, "xs", pair_production_electron_);
+    close_group(rgroup);
+  } else {
+    pair_production_electron_ = xt::zeros_like(energy_);
+  }
 
   // Read pair production
   if (object_exists(group, "pair_production_nuclear")) {