diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst
index 1210bc282..3a8a6d432 100644
--- a/docs/source/usersguide/tallies.rst
+++ b/docs/source/usersguide/tallies.rst
@@ -220,6 +220,16 @@ The following tables show all valid scores:
| |multiplicity from (n,2n), (n,3n), and (n,4n) |
| |reactions. |
+----------------------+---------------------------------------------------+
+ |H1-production |Total production of H1. |
+ +----------------------+---------------------------------------------------+
+ |H2-production |Total production of H2 (deuterium). |
+ +----------------------+---------------------------------------------------+
+ |H3-production |Total production of H1 (tritium). |
+ +----------------------+---------------------------------------------------+
+ |He3-production |Total production of He3. |
+ +----------------------+---------------------------------------------------+
+ |He4-production |Total production of He4 (alpha particles). |
+ +----------------------+---------------------------------------------------+
.. table:: **Miscellaneous scores: units are indicated for each.**
@@ -246,6 +256,10 @@ The following tables show all valid scores:
|inverse-velocity |The flux-weighted inverse velocity where the |
| |velocity is in units of centimeters per second. |
+----------------------+---------------------------------------------------+
+ |heating |Total neutron heating in units of eV per source |
+ | |particle. This corresponds to MT=301 produced by |
+ | |NJOY's HEATR module. |
+ +----------------------+---------------------------------------------------+
|kappa-fission |The recoverable energy production rate due to |
| |fission. The recoverable energy is defined as the |
| |fission product kinetic energy, prompt and delayed |
@@ -281,3 +295,7 @@ The following tables show all valid scores:
|decay-rate |The delayed-nu-fission-weighted decay rate where |
| |the decay rate is in units of inverse seconds. |
+----------------------+---------------------------------------------------+
+ |damage-energy |Damage energy production in units of eV per source |
+ | |particle. This corresponds to MT=444 produced by |
+ | |NJOY's HEATR module. |
+ +----------------------+---------------------------------------------------+
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index 067e6431a..87388c91a 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -225,6 +225,13 @@ constexpr int N_3P {197};
constexpr int N_N3P {198};
constexpr int N_3N2PA {199};
constexpr int N_5N2P {200};
+constexpr int N_XP {203};
+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 DAMAGE_ENERGY {444};
constexpr int COHERENT {502};
constexpr int INCOHERENT {504};
constexpr int PAIR_PROD_ELEC {515};
diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py
index eec083602..79761b0c9 100644
--- a/openmc/data/neutron.py
+++ b/openmc/data/neutron.py
@@ -448,11 +448,13 @@ class IncidentNeutron(EqualityMixin):
for rx in self.reactions.values():
# Skip writing redundant reaction if it doesn't have photon
# production or is a summed transmutation reaction. MT=4 is also
- # sometimes needed for probability tables.
+ # sometimes needed for probability tables. Also write gas
+ # production, heating, and damage energy production.
if rx.redundant:
photon_rx = any(p.particle == 'photon' for p in rx.products)
- transmutation_rx = (rx.mt in (16, 103, 104, 105, 106, 107))
- if not (photon_rx or transmutation_rx or rx.mt == 4):
+ keep_mts = (4, 16, 103, 104, 105, 106, 107,
+ 203, 204, 205, 206, 207, 301, 444)
+ if not (photon_rx or rx.mt in keep_mts):
continue
rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt))
@@ -615,13 +617,14 @@ class IncidentNeutron(EqualityMixin):
# Read energy grid
n_energy = ace.nxs[3]
- energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy]*EV_PER_MEV
+ i = ace.jxs[1]
+ energy = ace.xss[i : i + n_energy]*EV_PER_MEV
data.energy[strT] = energy
- total_xs = ace.xss[ace.jxs[1] + n_energy:ace.jxs[1] + 2 * n_energy]
- absorption_xs = ace.xss[ace.jxs[1] + 2 * n_energy:ace.jxs[1] +
- 3 * n_energy]
+ total_xs = ace.xss[i + n_energy : i + 2*n_energy]
+ absorption_xs = ace.xss[i + 2*n_energy : i + 3*n_energy]
+ heating_number = ace.xss[i + 3*n_energy : i + 4*n_energy]*EV_PER_MEV
- # Create redundant reactions (total and absorption)
+ # Create redundant reactions (total, absorption, and heating)
total = Reaction(1)
total.xs[strT] = Tabulated1D(energy, total_xs)
total.redundant = True
@@ -633,13 +636,15 @@ class IncidentNeutron(EqualityMixin):
absorption.redundant = True
data.reactions[101] = absorption
+ heating = Reaction(301)
+ heating.xs[strT] = Tabulated1D(energy, heating_number)
+ heating.redundant = True
+ data.reactions[301] = heating
+
# Read each reaction
n_reaction = ace.nxs[4] + 1
for i in range(n_reaction):
rx = Reaction.from_ace(ace, i)
- # Don't include gas production / damage cross sections
- if 200 < rx.mt < 219 or rx.mt == 444:
- continue
data.reactions[rx.mt] = rx
# Some photon production reactions may be assigned to MTs that don't
@@ -690,6 +695,8 @@ class IncidentNeutron(EqualityMixin):
mts = data.get_reaction_components(rx.mt)
if mts != [rx.mt]:
rx.redundant = True
+ if rx.mt in (203, 204, 205, 206, 207, 444):
+ rx.redundant = True
# Read unresolved resonance probability tables
urr = ProbabilityTables.from_ace(ace)
diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py
index 0c82f32f1..ab6ba54dc 100644
--- a/openmc/data/njoy.py
+++ b/openmc/data/njoy.py
@@ -72,8 +72,13 @@ broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
_TEMPLATE_HEATR = """
heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr} /
-{mat} 3 /
-302 318 402 /
+{mat} 4 /
+302 318 402 444 /
+"""
+
+_TEMPLATE_GASPR = """
+gaspr / %%%%%%%%%%%%%%%%%%%%%%%%% Add gas production %%%%%%%%%%%%%%%%%%%%%%%%%%%
+{nendf} {ngaspr_in} {ngaspr} /
"""
_TEMPLATE_PURR = """
@@ -211,8 +216,8 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
- error=0.001, broadr=True, heatr=True, purr=True, acer=True,
- **kwargs):
+ error=0.001, broadr=True, heatr=True, gaspr=True, purr=True,
+ acer=True, **kwargs):
"""Generate incident neutron ACE file from an ENDF file
Parameters
@@ -234,6 +239,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
Indicating whether to Doppler broaden XS when running NJOY
heatr : bool, optional
Indicating whether to add heating kerma when running NJOY
+ gaspr : bool, optional
+ Indicating whether to add gas production data when running NJOY
purr : bool, optional
Indicating whether to add probability table when running NJOY
acer : bool, optional
@@ -285,6 +292,13 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
commands += _TEMPLATE_HEATR
nlast = nheatr
+ # gaspr
+ if gaspr:
+ ngaspr_in = nlast
+ ngaspr = ngaspr_in + 1
+ commands += _TEMPLATE_GASPR
+ nlast = ngaspr
+
# purr
if purr:
npurr_in = nlast
diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py
index 11807cadc..a988f2d6d 100644
--- a/openmc/data/reaction.py
+++ b/openmc/data/reaction.py
@@ -51,7 +51,9 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
189: '(n,nta)', 190: '(n,2n2p)', 191: '(n,p3He)',
192: '(n,d3He)', 193: '(n,3Hea)', 194: '(n,4n2p)',
195: '(n,4n2a)', 196: '(n,4npa)', 197: '(n,3p)',
- 198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 444: '(n,damage)',
+ 198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 203: '(n,Xp)',
+ 204: '(n,Xd)', 205: '(n,Xt)', 206: '(n,X3He)', 207: '(n,Xa)',
+ 301: 'heating', 444: 'damage-energy',
649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)',
849: '(n,ac)', 891: '(n,2nc)'}
REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(50, 91)})
@@ -988,6 +990,10 @@ class Reaction(EqualityMixin):
# Read reaction cross section
xs = ace.xss[ace.jxs[7] + loc + 1:ace.jxs[7] + loc + 1 + n_energy]
+ # For damage energy production, convert to eV
+ if mt == 444:
+ xs *= EV_PER_MEV
+
# Fix negatives -- known issue for Y89 in JEFF 3.2
if np.any(xs < 0.0):
warn("Negative cross sections found for MT={} in {}. Setting "
diff --git a/src/reaction.cpp b/src/reaction.cpp
index 1f997f078..d2405d899 100644
--- a/src/reaction.cpp
+++ b/src/reaction.cpp
@@ -229,8 +229,10 @@ std::string reaction_name(int mt)
return "(n,X3He)";
} else if (mt == 207) {
return "(n,Xa)";
+ } else if (mt == 301) {
+ return "heating";
} else if (mt == 444) {
- return "(damage)";
+ return "damage-energy";
} else if (mt == COHERENT) {
return "coherent scatter";
} else if (mt == INCOHERENT) {
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 5bae26b75..682db3540 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -197,6 +197,20 @@ score_str_to_int(std::string score_str)
return N_PT;
if (score_str == "(n,da)")
return N_DA;
+ if (score_str == "(n,Xp)" || score_str == "H1-production")
+ return N_XP;
+ if (score_str == "(n,Xd)" || score_str == "H2-production")
+ return N_XD;
+ if (score_str == "(n,Xt)" || score_str == "H3-production")
+ return N_XT;
+ if (score_str == "(n,X3He)" || score_str == "He3-production")
+ 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;
// So far we have not identified this score string. Check to see if it is a
// deprecated score.
diff --git a/tests/regression_tests/tallies/inputs_true.dat b/tests/regression_tests/tallies/inputs_true.dat
index 80f55b034..3ec8d8d4c 100644
--- a/tests/regression_tests/tallies/inputs_true.dat
+++ b/tests/regression_tests/tallies/inputs_true.dat
@@ -511,4 +511,7 @@
total
tracklength
+
+ H1-production H2-production H3-production He3-production He4-production heating damage-energy
+
diff --git a/tests/regression_tests/tallies/results_true.dat b/tests/regression_tests/tallies/results_true.dat
index 08a3ba110..9830ff872 100644
--- a/tests/regression_tests/tallies/results_true.dat
+++ b/tests/regression_tests/tallies/results_true.dat
@@ -1 +1 @@
-b5edf87cb58db29aa1c38203141df2f055b88aca675aee9673879b579391412abd232cb3d404f60840ac86b936384050576cccb1e3a1fc5c8290b6e890dda74c
\ No newline at end of file
+de773a799f84348241ebd65bf7beae8114861d8674b652bfd443d578c146a7d37ca38f2efc5d05124fdbb1cf12829907768351e6b2d6b5f4bd2c121417757507
\ No newline at end of file
diff --git a/tests/regression_tests/tallies/test.py b/tests/regression_tests/tallies/test.py
index eeb74636f..543d0acd2 100644
--- a/tests/regression_tests/tallies/test.py
+++ b/tests/regression_tests/tallies/test.py
@@ -165,6 +165,10 @@ def test_tallies():
all_nuclide_tallies[3].filters = [mesh_filter]
all_nuclide_tallies[3].nuclides = ['U235']
+ fusion_tally = Tally()
+ fusion_tally.scores = ['H1-production', 'H2-production', 'H3-production',
+ 'He3-production', 'He4-production', 'heating', 'damage-energy']
+
model.tallies += [
azimuthal_tally1, azimuthal_tally2, azimuthal_tally3,
cellborn_tally, dg_tally, energy_tally, energyout_tally,
@@ -174,5 +178,6 @@ def test_tallies():
model.tallies += score_tallies
model.tallies += flux_tallies
model.tallies += all_nuclide_tallies
+ model.tallies.append(fusion_tally)
harness.main()