Merge pull request #1776 from paulromano/urr-capture-fix

Ensure (n,gamma) reaction rate tally uses sampled URR capture cross section

src/nuclide.cpp

@@ -896,18 +896,21 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
   if (fission < 0.) {fission = 0.;}
   if (capture < 0.) {capture = 0.;}
 
-  // Set elastic, absorption, fission, and total x/s. Note that the total x/s
-  // is calculated as a sum of partials instead of the table-provided value
+  // Set elastic, absorption, fission, total, and capture x/s. Note that the
+  // total x/s is calculated as a sum of partials instead of the table-provided
+  // value
   micro.elastic = elastic;
   micro.absorption = capture + fission;
   micro.fission = fission;
   micro.total = elastic + inelastic + capture + fission;
+  if (simulation::need_depletion_rx) {
+    micro.reaction[0] = capture;
+  }
 
   // Determine nu-fission cross-section
   if (fissionable_) {
     micro.nu_fission = nu(p.E_, EmissionMode::total) * micro.fission;
   }
-
 }
 
 std::pair<gsl::index, double> Nuclide::find_temperature(double T) const

src/tallies/tally_scoring.cpp

@@ -479,6 +479,13 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) {
   const auto& rxn {*nuc.reactions_[m]};
   const auto& micro {p.neutron_xs_[i_nuclide]};
 
+  // In the URR, the (n,gamma) cross section is sampled randomly from
+  // probability tables. Make sure we use the sampled value (which is equal to
+  // absorption - fission) rather than the dilute average value
+  if (micro.use_ptable && score_bin == N_GAMMA) {
+    return micro.absorption - micro.fission;
+  }
+
   auto i_temp = micro.index_temp;
   if (i_temp >= 0) { // Can be false due to multipole
     // Get index on energy grid and interpolation factor

tests/regression_tests/filter_energyfun/results_true.dat

@@ -1,2 +1,2 @@
    energyfunction nuclide                    score     mean  std. dev.
-0  d2effa26cb3cf2   Am241  ((n,gamma) / (n,gamma)) 1.74e-01   3.55e-03
+0  d2effa26cb3cf2   Am241  ((n,gamma) / (n,gamma)) 1.74e-01   3.52e-03

tests/regression_tests/mgxs_library_no_nuclides/results_true.dat

@@ -184,7 +184,7 @@ nu-diffusion-coefficient
 0         1         2   total  0.000000   0.000000
 (n,gamma)
    material  group in nuclide      mean  std. dev.
-1         1         1   total  0.019727   0.002262
+1         1         1   total  0.019720   0.002260
 0         1         2   total  0.071719   0.006262
 (n,a)
    material  group in nuclide      mean  std. dev.
@@ -490,7 +490,7 @@ nu-diffusion-coefficient
 0         2         2   total  0.000000   0.000000
 (n,gamma)
    material  group in nuclide      mean  std. dev.
-1         2         1   total  0.001569   0.000322
+1         2         1   total  0.001571   0.000323
 0         2         2   total  0.005400   0.000618
 (n,a)
    material  group in nuclide          mean     std. dev.

tests/regression_tests/mgxs_library_nuclides/results_true.dat

@@ -1 +1 @@
-a9999488e2aa2ad0f1d694afedb71fbb56f1d0c8e8abdb6616698505a9d85302bf90586866fdb58eac96f76712a2dfbc884cc03087df7b5d148fb29d14dd2bbb
+54a6351f3444c5ea85f972ed05f48c81da2de27a6295a660c321e7c595ec3865e5882fb7fa1f2e4a3bfd616e4b637482e2a16930dd59c3968d04c59d4b1c08b1

tests/regression_tests/multipole/results_true.dat

@@ -5,8 +5,8 @@ tally 1:
 2.951721E+00
 2.785273E+00
 1.554128E+00
-5.349716E-01
-5.732174E-02
+5.352493E-01
+5.738169E-02
 4.499834E-01
 4.055011E-02
 0.000000E+00
@@ -25,8 +25,8 @@ tally 1:
 2.628739E+04
 2.785273E+00
 1.554128E+00
-2.176478E+00
-9.480405E-01
+2.176755E+00
+9.482820E-01
 3.574110E+02
 2.555993E+04
 1.146136E-04

tests/regression_tests/tallies/results_true.dat

@@ -1 +1 @@
-ba4369130ac65812939cb417c7088ae6c8cdec118fd7dce5645b40a20a3ae47072cd72c24aa091ccc09f833a5bfc0bd72407a0dcde7d3120e39379942614a75b
+0f036e2c34f3e9cf350d4f6bac6b0c5933d77d0665a5d461a3911c29d8f1d49b74100e680920da286c63538fabb86cbcd17532c6aad6048ee78e6cef75dd76ae

tests/unit_tests/test_deplete_resultslist.py

@@ -22,7 +22,7 @@ def test_get_atoms(res):
 
     t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0])
     n_ref = np.array(
-        [6.67473282e+08, 3.76986925e+14, 3.68587383e+14, 3.91338675e+14])
+        [6.67473282e+08, 3.76987065e+14, 3.68586723e+14, 3.91338392e+14])
 
     np.testing.assert_allclose(t, t_ref)
     np.testing.assert_allclose(n, n_ref)
@@ -48,8 +48,8 @@ def test_get_reaction_rate(res):
     t, r = res.get_reaction_rate("1", "Xe135", "(n,gamma)")
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
-    n_ref = [6.67473282e+08, 3.76986925e+14, 3.68587383e+14, 3.91338675e+14]
-    xs_ref = [3.32282266e-05, 2.76207120e-05, 4.10986677e-05, 3.72453665e-05]
+    n_ref = [6.67473282e+08, 3.76987065e+14, 3.68586723e+14, 3.91338392e+14]
+    xs_ref = [3.32282064e-05, 2.76208092e-05, 4.10987995e-05, 3.72454755e-05]
 
     np.testing.assert_allclose(t, t_ref)
     np.testing.assert_allclose(r, np.array(n_ref) * xs_ref)
@@ -60,8 +60,8 @@ def test_get_eigenvalue(res):
     t, k = res.get_eigenvalue()
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
-    k_ref = [1.16984322, 1.19097427, 1.03012572, 1.20045627]
-    u_ref = [0.0375587, 0.0347639, 0.07216021, 0.02839642]
+    k_ref = [1.16984322, 1.19097429, 1.03012517, 1.20045563]
+    u_ref = [0.0375587, 0.03476389, 0.07215969, 0.02839639]
 
     np.testing.assert_allclose(t, t_ref)
     np.testing.assert_allclose(k[:, 0], k_ref)

tests/unit_tests/test_urr_capture.py

@@ -0,0 +1,42 @@
+import openmc
+import pytest
+
+
+@pytest.fixture
+def th232_model():
+    # URR boundaries for Th232
+    e_min, e_max = 4000.0, 100000.0
+
+    model = openmc.model.Model()
+    th232 = openmc.Material()
+    th232.add_nuclide('Th232', 1.0)
+
+    surf = openmc.Sphere(r=100.0, boundary_type='reflective')
+    cell = openmc.Cell(fill=th232, region=-surf)
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings.particles = 100
+    model.settings.batches = 10
+    model.settings.run_mode = 'fixed source'
+    energies = openmc.stats.Uniform(e_min, e_max)
+    model.settings.source = openmc.Source(energy=energies)
+
+    tally = openmc.Tally(name='rates')
+    tally.filters = [openmc.EnergyFilter([e_min, e_max])]
+    tally.scores = ['(n,gamma)', 'absorption', 'fission']
+    model.tallies.append(tally)
+    return model
+
+
+def test_urr_capture(run_in_tmpdir, th232_model):
+    # Export and run model
+    th232_model.export_to_xml()
+    openmc.run()
+
+    # Get reaction rates from tally
+    with openmc.StatePoint('statepoint.10.h5') as sp:
+        t = sp.get_tally(name='rates')
+        ngamma, absorption, fission = t.mean.flatten()
+
+    # In URR, the (n,gamma) rate should be equal to absorption - fission
+    assert ngamma == pytest.approx(absorption - fission)