diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index f005b50606..1676417c66 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -263,6 +263,11 @@ The following tables show all valid scores: | |deposition (analog estimator) or pre-generated | | |photon heating number. | +----------------------+---------------------------------------------------+ + |heating-local |Total nuclear heating in units of eV per source | + | |particle assuming energy from secondary photons is | + | |deposited locally. Note that this score should only| + | |be used for incident neutrons. | + +----------------------+---------------------------------------------------+ |kappa-fission |The recoverable energy production rate due to | | |fission. The recoverable energy is defined as the | | |fission product kinetic energy, prompt and delayed | diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 394f90740f..51df59accf 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -232,7 +232,7 @@ constexpr int N_XD {204}; constexpr int N_XT {205}; constexpr int N_X3HE {206}; constexpr int N_XA {207}; -constexpr int NEUTRON_HEATING {301}; +constexpr int HEATING {301}; constexpr int DAMAGE_ENERGY {444}; constexpr int COHERENT {502}; constexpr int INCOHERENT {504}; @@ -252,6 +252,7 @@ constexpr int N_A0 {800}; constexpr int N_AC {849}; constexpr int N_2N0 {875}; constexpr int N_2NC {891}; +constexpr int HEATING_LOCAL {901}; constexpr std::array DEPLETION_RX {N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N}; @@ -369,7 +370,6 @@ 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}; diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index a936ac0102..a1571ea58e 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -93,7 +93,7 @@ public: std::vector urr_data_; std::vector> reactions_; //!< Reactions - std::array reaction_index_; //!< Index of each reaction + std::array reaction_index_; //!< Index of each reaction std::vector index_inelastic_scatter_; private: diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index dfd44a4d07..022cf8283d 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -452,7 +452,7 @@ class IncidentNeutron(EqualityMixin): if rx.redundant: photon_rx = any(p.particle == 'photon' for p in rx.products) keep_mts = (4, 16, 103, 104, 105, 106, 107, - 203, 204, 205, 206, 207, 301, 318, 444) + 203, 204, 205, 206, 207, 301, 444, 901) if not (photon_rx or rx.mt in keep_mts): continue @@ -553,20 +553,6 @@ class IncidentNeutron(EqualityMixin): fer_group = group['fission_energy_release'] data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group) - # Rebuild non-fission heating - total_heating = data.reactions.get(301) - fission_heating = data.reactions.get(318) - if total_heating is not None and fission_heating is not None: - non_fission_heating = Reaction(999) - non_fission_heating.redundant = True - for strT, total in total_heating.xs.items(): - fission = fission_heating.xs.get(strT) - if fission is None: - continue - non_fission_heating.xs[strT] = Tabulated1D( - total.x, total.y - fission(total.x)) - data.reactions[999] = non_fission_heating - return data @classmethod @@ -744,7 +730,7 @@ class IncidentNeutron(EqualityMixin): # Instantiate incident neutron data data = cls(name, atomic_number, mass_number, metastable, - atomic_weight_ratio, temperature) + atomic_weight_ratio, [temperature]) if (2, 151) in ev.section: data.resonances = res.Resonances.from_endf(ev) @@ -777,9 +763,10 @@ class IncidentNeutron(EqualityMixin): for mt, rx in data.reactions.items(): if mt in (19, 20, 21, 38): if (5, mt) not in ev.section: - neutron = data.reactions[18].products[0] - rx.products[0].applicability = neutron.applicability - rx.products[0].distribution = neutron.distribution + if rx.products: + neutron = data.reactions[18].products[0] + rx.products[0].applicability = neutron.applicability + rx.products[0].distribution = neutron.distribution # Read fission energy release (requires that we already know nu for # fission) @@ -826,37 +813,6 @@ class IncidentNeutron(EqualityMixin): for table in lib.tables[1:]: data.add_temperature_from_ace(table) - # Add fission energy release data - ev = evaluation if evaluation is not None else Evaluation(filename) - if (1, 458) in ev.section: - data.fission_energy = FissionEnergyRelease.from_endf(ev, data) - # Add 318 fission heating data from heatr - non_fission_heating = Reaction(999) - non_fission_heating.redundant = True - fission_heating = Reaction(318) - fission_heating.redundant = True - - heatr_evals = get_evaluations(kwargs["heatr"]) - for heatr in heatr_evals: - temp = "{}K".format(round(heatr.target["temperature"])) - f318 = StringIO(heatr.section[3, 318]) - get_head_record(f318) - _params, fission_kerma = get_tab1_record(f318) - total_heating_xs = data.reactions[301].xs.get(temp) - if total_heating_xs is None: - fission_heating.xs[temp] = fission_kerma - continue - # Cast fission heating to same grid as total heating - new_fission_heat_xs = fission_kerma(total_heating_xs.x) - fission_heating.xs[temp] = Tabulated1D( - total_heating_xs.x, new_fission_heat_xs) - non_fission_heating.xs[temp] = Tabulated1D( - total_heating_xs.x, - total_heating_xs.y - new_fission_heat_xs) - - data.reactions[318] = fission_heating - data.reactions[999] = non_fission_heating - # Add 0K elastic scattering cross section if '0K' not in data.energy: pendf = Evaluation(kwargs['pendf']) @@ -866,6 +822,74 @@ class IncidentNeutron(EqualityMixin): data.energy['0K'] = xs.x data[2].xs['0K'] = xs + # Add fission energy release data + ev = evaluation if evaluation is not None else Evaluation(filename) + if (1, 458) in ev.section: + data.fission_energy = f = FissionEnergyRelease.from_endf(ev, data) + else: + f = None + + # For energy deposition, we want to store two different KERMAs: + # one calculated assuming outgoing photons deposit their energy + # locally, and one calculated assuming they carry their energy + # away. This requires two HEATR runs (which make_ace does by + # default). Here, we just need to correct for the fact that NJOY + # uses a fission heating number of h = EFR, whereas we want: + # + # 1) h = EFR + EGP + EGD + EB (for local case) + # 2) h = EFR + EB (for non-local case) + # + # The best way to handle this is to subtract off the fission + # KERMA that NJOY calculates and add back exactly what we want. + + # If NJOY is not run with HEATR at all, skip everything below + if not kwargs["heatr"]: + return data + + # Helper function to get a cross section from an ENDF file on a + # given energy grid + def get_file3_xs(ev, mt, E): + file_obj = StringIO(ev.section[3, mt]) + get_head_record(file_obj) + _, xs = get_tab1_record(file_obj) + return xs(E) + + heating_local = Reaction(901) + heating_local.redundant = True + + heatr_evals = get_evaluations(kwargs["heatr"]) + heatr_local_evals = get_evaluations(kwargs["heatr"] + "_local") + for ev, ev_local in zip(heatr_evals, heatr_local_evals): + temp = "{}K".format(round(ev.target["temperature"])) + + # Get total KERMA (originally from ACE file) and energy grid + kerma = data.reactions[301].xs[temp] + E = kerma.x + + if f is not None: + # Replace fission KERMA with (EFR + EB)*sigma_f + fission = data.reactions[18].xs[temp] + kerma_fission = get_file3_xs(ev, 318, E) + kerma.y = kerma.y - kerma_fission + ( + f.fragments(E) + f.betas(E)) * fission.y + + # For local KERMA, we first need to get the values from the + # HEATR run with photon energy deposited locally and put + # them on the same energy grid + kerma_local = get_file3_xs(ev_local, 301, E) + + if f is not None: + # When photons deposit their energy locally, we replace the + # fission KERMA with (EFR + EGP + EGD + EB)*sigma_f + kerma_fission_local = get_file3_xs(ev_local, 318, E) + kerma_local = kerma_local - kerma_fission_local + ( + f.fragments(E) + f.prompt_photons(E) + + f.delayed_photons(E) + f.betas(E))*fission.y + + heating_local.xs[temp] = Tabulated1D(E, kerma_local) + + data.reactions[901] = heating_local + return data def _get_redundant_reaction(self, mt, mts): diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index 7993ec4004..4983098f73 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -71,7 +71,14 @@ broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%% _TEMPLATE_HEATR = """ heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%% {nendf} {nheatr_in} {nheatr} / -{mat} 4 / +{mat} 4 0 0 0 / +302 318 402 444 / +""" + +_TEMPLATE_HEATR_LOCAL = """ +heatr / %%%%%%%%%%%%%%%%% Add heating kerma (local photons) %%%%%%%%%%%%%%%%%%%% +{nendf} {nheatr_in} {nheatr_local} / +{mat} 4 0 0 1 / 302 318 402 444 / """ @@ -216,8 +223,7 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False): def make_ace(filename, temperatures=None, acer=True, xsdir=None, output_dir=None, pendf=False, error=0.001, broadr=True, - heatr=True, gaspr=True, purr=True, evaluation=None, - **kwargs): + heatr=True, gaspr=True, purr=True, evaluation=None, **kwargs): """Generate incident neutron ACE file from an ENDF file File names can be passed to @@ -320,7 +326,11 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, # heatr if heatr: nheatr_in = nlast - nheatr = nheatr_in + 1 + nheatr_local = nheatr_in + 1 + tapeout[nheatr_local] = (output_dir / "heatr_local") if heatr is True \ + else heatr + '_local' + commands += _TEMPLATE_HEATR_LOCAL + nheatr = nheatr_local + 1 tapeout[nheatr] = (output_dir / "heatr") if heatr is True else heatr commands += _TEMPLATE_HEATR nlast = nheatr diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 696310ccaf..83e773f204 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -54,9 +54,8 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)', 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', - 318: "fission-heating", 999: "non-fission-heating", 649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)', - 849: '(n,ac)', 891: '(n,2nc)'} + 849: '(n,ac)', 891: '(n,2nc)', 901: 'heating-local'} REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(50, 91)}) REACTION_NAME.update({i: '(n,p{})'.format(i - 600) for i in range(600, 649)}) REACTION_NAME.update({i: '(n,d{})'.format(i - 650) for i in range(650, 699)}) diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index 0001b2c4b6..97b7437ea3 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/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', -17: 'heating' + -15: 'fission-q-recoverable', -16: 'decay-rate' } _ESTIMATORS = { 1: 'analog', 2: 'tracklength', 3: 'collision' diff --git a/src/output.cpp b/src/output.cpp index 6516b48b61..40cd1dbfe9 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -620,7 +620,6 @@ const std::unordered_map 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. diff --git a/src/reaction.cpp b/src/reaction.cpp index 31bf8efd33..f7f1e87acb 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -88,193 +88,189 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) std::string reaction_name(int mt) { if (mt == SCORE_FLUX) { - return "flux"; + return "flux"; } else if (mt == SCORE_TOTAL) { - return "total"; + return "total"; } else if (mt == SCORE_SCATTER) { - return "scatter"; + return "scatter"; } else if (mt == SCORE_NU_SCATTER) { - return "nu-scatter"; + return "nu-scatter"; } else if (mt == SCORE_ABSORPTION) { - return "absorption"; + return "absorption"; } else if (mt == SCORE_FISSION) { - return "fission"; + return "fission"; } else if (mt == SCORE_NU_FISSION) { - return "nu-fission"; + return "nu-fission"; } else if (mt == SCORE_DECAY_RATE) { - return "decay-rate"; + return "decay-rate"; } else if (mt == SCORE_DELAYED_NU_FISSION) { - return "delayed-nu-fission"; + return "delayed-nu-fission"; } else if (mt == SCORE_PROMPT_NU_FISSION) { - return "prompt-nu-fission"; + return "prompt-nu-fission"; } else if (mt == SCORE_KAPPA_FISSION) { - return "kappa-fission"; + return "kappa-fission"; } else if (mt == SCORE_CURRENT) { - return "current"; + return "current"; } else if (mt == SCORE_EVENTS) { - return "events"; + return "events"; } else if (mt == SCORE_INVERSE_VELOCITY) { - return "inverse-velocity"; + return "inverse-velocity"; } else if (mt == SCORE_FISS_Q_PROMPT) { - return "fission-q-prompt"; + return "fission-q-prompt"; } else if (mt == SCORE_FISS_Q_RECOV) { - return "fission-q-recoverable"; - } else if (mt == SCORE_HEATING) { - return "heating"; + return "fission-q-recoverable"; // Normal ENDF-based reactions } else if (mt == TOTAL_XS) { - return "(n,total)"; + return "(n,total)"; } else if (mt == ELASTIC) { - return "(n,elastic)"; + return "(n,elastic)"; } else if (mt == N_LEVEL) { - return "(n,level)"; + return "(n,level)"; } else if (mt == N_2ND) { - return "(n,2nd)"; + return "(n,2nd)"; } else if (mt == N_2N) { - return "(n,2n)"; + return "(n,2n)"; } else if (mt == N_3N) { - return "(n,3n)"; + return "(n,3n)"; } else if (mt == N_FISSION) { - return "(n,fission)"; + return "(n,fission)"; } else if (mt == N_F) { - return "(n,f)"; + return "(n,f)"; } else if (mt == N_NF) { - return "(n,nf)"; + return "(n,nf)"; } else if (mt == N_2NF) { - return "(n,2nf)"; + return "(n,2nf)"; } else if (mt == N_NA) { - return "(n,na)"; + return "(n,na)"; } else if (mt == N_N3A) { - return "(n,n3a)"; + return "(n,n3a)"; } else if (mt == N_2NA) { - return "(n,2na)"; + return "(n,2na)"; } else if (mt == N_3NA) { - return "(n,3na)"; + return "(n,3na)"; } else if (mt == N_NP) { - return "(n,np)"; + return "(n,np)"; } else if (mt == N_N2A) { - return "(n,n2a)"; + return "(n,n2a)"; } else if (mt == N_2N2A) { - return "(n,2n2a)"; + return "(n,2n2a)"; } else if (mt == N_ND) { - return "(n,nd)"; + return "(n,nd)"; } else if (mt == N_NT) { - return "(n,nt)"; + return "(n,nt)"; } else if (mt == N_N3HE) { - return "(n,nHe-3)"; + return "(n,nHe-3)"; } else if (mt == N_ND2A) { - return "(n,nd2a)"; + return "(n,nd2a)"; } else if (mt == N_NT2A) { - return "(n,nt2a)"; + return "(n,nt2a)"; } else if (mt == N_4N) { - return "(n,4n)"; + return "(n,4n)"; } else if (mt == N_3NF) { - return "(n,3nf)"; + return "(n,3nf)"; } else if (mt == N_2NP) { - return "(n,2np)"; + return "(n,2np)"; } else if (mt == N_3NP) { - return "(n,3np)"; + return "(n,3np)"; } else if (mt == N_N2P) { - return "(n,n2p)"; + return "(n,n2p)"; } else if (mt == N_NPA) { - return "(n,npa)"; + return "(n,npa)"; } else if (N_N1 <= mt && mt <= N_N40) { - return "(n,n" + std::to_string(mt-50) + ")"; + return "(n,n" + std::to_string(mt-50) + ")"; } else if (mt == N_NC) { - return "(n,nc)"; + return "(n,nc)"; } else if (mt == N_DISAPPEAR) { - return "(n,disappear)"; + return "(n,disappear)"; } else if (mt == N_GAMMA) { - return "(n,gamma)"; + return "(n,gamma)"; } else if (mt == N_P) { - return "(n,p)"; + return "(n,p)"; } else if (mt == N_D) { - return "(n,d)"; + return "(n,d)"; } else if (mt == N_T) { - return "(n,t)"; + return "(n,t)"; } else if (mt == N_3HE) { - return "(n,3He)"; + return "(n,3He)"; } else if (mt == N_A) { - return "(n,a)"; + return "(n,a)"; } else if (mt == N_2A) { - return "(n,2a)"; + return "(n,2a)"; } else if (mt == N_3A) { - return "(n,3a)"; + return "(n,3a)"; } else if (mt == N_2P) { - return "(n,2p)"; + return "(n,2p)"; } else if (mt == N_PA) { - return "(n,pa)"; + return "(n,pa)"; } else if (mt == N_T2A) { - return "(n,t2a)"; + return "(n,t2a)"; } else if (mt == N_D2A) { - return "(n,d2a)"; + return "(n,d2a)"; } else if (mt == N_PD) { - return "(n,pd)"; + return "(n,pd)"; } else if (mt == N_PT) { - return "(n,pt)"; + return "(n,pt)"; } else if (mt == N_DA) { - return "(n,da)"; + return "(n,da)"; } else if (mt == 201) { - return "(n,Xn)"; + return "(n,Xn)"; } else if (mt == 202) { - return "(n,Xgamma)"; - } else if (mt == 203) { - return "(n,Xp)"; - } else if (mt == 204) { - return "(n,Xd)"; - } else if (mt == 205) { - return "(n,Xt)"; - } else if (mt == 206) { - return "(n,X3He)"; - } else if (mt == 207) { - return "(n,Xa)"; - } else if (mt == 301) { + return "(n,Xgamma)"; + } else if (mt == N_XP) { + return "(n,Xp)"; + } else if (mt == N_XD) { + return "(n,Xd)"; + } else if (mt == N_XT) { + return "(n,Xt)"; + } else if (mt == N_X3HE) { + return "(n,X3He)"; + } else if (mt == N_XA) { + return "(n,Xa)"; + } else if (mt == HEATING) { return "heating"; - } else if (mt == 318) { - return "fission-heating"; - } else if (mt == 999) { - return "non-fission-heating"; - } else if (mt == 444) { - return "damage-energy"; + } else if (mt == DAMAGE_ENERGY) { + return "damage-energy"; } else if (mt == COHERENT) { - return "coherent scatter"; + return "coherent scatter"; } else if (mt == INCOHERENT) { - return "incoherent scatter"; + return "incoherent scatter"; } else if (mt == PAIR_PROD_ELEC) { - return "pair production, electron"; + return "pair production, electron"; } else if (mt == PAIR_PROD) { - return "pair production"; + return "pair production"; } else if (mt == PAIR_PROD_NUC) { - return "pair production, nuclear"; + return "pair production, nuclear"; } else if (mt == PHOTOELECTRIC) { - return "photoelectric"; + return "photoelectric"; } else if (534 <= mt && mt <= 572) { std::stringstream name; name << "photoelectric, " << SUBSHELLS[mt - 534] << " subshell"; return name.str(); } else if (600 <= mt && mt <= 648) { - return "(n,p" + std::to_string(mt-600) + ")"; + return "(n,p" + std::to_string(mt-600) + ")"; } else if (mt == 649) { - return "(n,pc)"; + return "(n,pc)"; } else if (650 <= mt && mt <= 698) { - return "(n,d" + std::to_string(mt-650) + ")"; + return "(n,d" + std::to_string(mt-650) + ")"; } else if (mt == 699) { - return "(n,dc)"; + return "(n,dc)"; } else if (700 <= mt && mt <= 748) { - return "(n,t" + std::to_string(mt-700) + ")"; + return "(n,t" + std::to_string(mt-700) + ")"; } else if (mt == 749) { - return "(n,tc)"; + return "(n,tc)"; } else if (750 <= mt && mt <= 798) { - return "(n,3He" + std::to_string(mt-750) + ")"; + return "(n,3He" + std::to_string(mt-750) + ")"; } else if (mt == 799) { - return "(n,3Hec)"; + return "(n,3Hec)"; } else if (800 <= mt && mt <= 848) { - return "(n,a" + std::to_string(mt-800) + ")"; + return "(n,a" + std::to_string(mt-800) + ")"; } else if (mt == 849) { - return "(n,ac)"; + return "(n,ac)"; + } else if (mt == HEATING_LOCAL) { + return "heating-local"; } else { - return "MT=" + std::to_string(mt); + return "MT=" + std::to_string(mt); } } diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 27d6cc58c0..563beee766 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -111,7 +111,10 @@ score_str_to_int(std::string score_str) return SCORE_FISS_Q_RECOV; if (score_str == "heating") - return SCORE_HEATING; + return HEATING; + + if (score_str == "heating-local") + return HEATING_LOCAL; if (score_str == "current") return SCORE_CURRENT; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 1f6a8010a5..37310070fe 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -192,7 +192,7 @@ double get_nuc_fission_q(const Nuclide& nuc, const Particle* p, int score_bin) //! Pulled out to support both the fission_q scores and energy deposition //! score -double score_fission_q(const Particle* p, int score_bin, const Tally& tally, +double score_fission_q(const Particle* p, int score_bin, const Tally& tally, double flux, int i_nuclide, double atom_density) { if (tally.estimator_ == ESTIMATOR_ANALOG) { @@ -231,7 +231,7 @@ double score_fission_q(const Particle* p, int score_bin, const Tally& tally, auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); const Nuclide& nuc {*data::nuclides[j_nuclide]}; - score += get_nuc_fission_q(nuc, p, score_bin) * atom_density + score += get_nuc_fission_q(nuc, p, score_bin) * atom_density * p->neutron_xs_[j_nuclide].fission; } return score * flux; @@ -268,7 +268,7 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux, // Get heating macroscopic "cross section" double heating_xs; if (i_nuclide >= 0) { - const Nuclide& nuc {*data::nuclides[i_nuclide]}; + const Nuclide& nuc {*data::nuclides[i_nuclide]}; heating_xs = get_nuclide_neutron_heating(p, nuc, rxn_bin, i_nuclide); if (tally.estimator_ == ESTIMATOR_ANALOG) { heating_xs /= p->neutron_xs_[i_nuclide].total; @@ -1210,10 +1210,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, break; - case SCORE_HEATING: + case HEATING: score = 0.; if (p->type_ == Particle::Type::neutron) { - score = score_neutron_heating(p, tally, flux, NEUTRON_HEATING, + score = score_neutron_heating(p, tally, flux, HEATING, i_nuclide, atom_density); } else if (p->type_ == Particle::Type::photon) { if (tally.estimator_ == ESTIMATOR_ANALOG) { diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index ed83956031..d8639cc5a1 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -205,20 +205,22 @@ def test_derived_products(am244): assert total_neutron.yield_(6e6) == pytest.approx(4.2558) -def test_heating(run_in_tmpdir, am244): - assert 318 in am244.reactions - assert 999 in am244.reactions - # compare values in 999 reaction - total_heating = am244.reactions[301].xs["294K"] - fission_heating = am244.reactions[318].xs["294K"] - expected_y = total_heating.y - fission_heating(total_heating.x) - assert am244.reactions[999].xs["294K"].y == pytest.approx(expected_y) +def test_kerma(run_in_tmpdir, am244, h2): + # Make sure kerma w/ local photon is >= regular kerma + for nuc in (am244, h2): + assert 301 in nuc + assert 901 in nuc + for T in nuc.temperatures: + k, k_local = nuc[301].xs[T], nuc[901].xs[T] + assert np.all(k.x == k_local.x) + assert np.all(k_local.y >= k.y) - am244.export_to_hdf5("am244.h5") - read_in = openmc.data.IncidentNeutron.from_hdf5("am244.h5") - assert 318 in read_in.reactions - assert 999 in read_in.reactions - assert read_in.reactions[999].xs["294K"].y == pytest.approx(expected_y) + # Make sure 301/901 get exported/imported correctly + h2.export_to_hdf5("H2.h5") + read_in = openmc.data.IncidentNeutron.from_hdf5("H2.h5") + assert 301 in read_in + assert 901 in read_in + assert np.all(read_in[901].xs['300K'].y == h2[901].xs['300K'].y) def test_urr(pu239):