diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 8174108e1ae..09eb7a21270 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -415,6 +415,7 @@ Kalbach-Mann :Object type: Group :Attributes: - **type** (*char[]*) -- 'kalbach-mann' + - **particle** (*char[]*) -- Type of particle :Datasets: - **energy** (*double[]*) -- Incoming energies at which distributions exist :Attributes: diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index 83806d35248..a1a27dc85b9 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -44,6 +44,7 @@ class KalbachMann : public AngleEnergy { xt::xtensor a; //!< Parameterized function }; + bool is_photon_; //!< Whether the projectile is a photon int n_region_; //!< Number of interpolation regions vector breakpoints_; //!< Breakpoints between regions vector interpolation_; //!< Interpolation laws diff --git a/openmc/data/data.py b/openmc/data/data.py index 5ecadd37be0..e075ebbe145 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -274,6 +274,7 @@ # Unit conversions EV_PER_MEV = 1.0e6 JOULE_PER_EV = 1.602176634e-19 +EV_PER_AMU = 931.49410372e6 # eV/c^2 per amu # Avogadro's constant AVOGADRO = 6.02214076e23 @@ -281,6 +282,9 @@ # Neutron mass in units of amu NEUTRON_MASS = 1.00866491595 +# Neutron mass in units of eV/c^2 +NEUTRON_MASS_EV = NEUTRON_MASS*EV_PER_AMU + # Used in atomic_mass function as a cache _ATOMIC_MASS: dict[str, float] = {} diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index d92bf9c213a..b612110a410 100644 --- a/openmc/data/kalbach_mann.py +++ b/openmc/data/kalbach_mann.py @@ -9,7 +9,7 @@ from openmc.stats import Tabular, Univariate, Discrete, Mixture from .function import Tabulated1D, INTERPOLATION_SCHEME from .angle_energy import AngleEnergy -from .data import EV_PER_MEV +from .data import EV_PER_MEV, NEUTRON_MASS_EV from .endf import get_list_record, get_tab2_record @@ -204,13 +204,14 @@ def kalbach_slope(energy_projectile, energy_emitted, za_projectile, Kalbach-Mann slope given with the same format as ACE file. """ - # TODO: develop for photons as projectile - # TODO: test for other particles than neutron - if za_projectile != 1: - raise NotImplementedError( - "Developed and tested for neutron projectile only." - ) + if za_projectile == 0: + # Calculate slope for photons using Eq. 3 in doi:10.1080/18811248.1995.9731830 + # or ENDF-6 Formats Manual section 6.2.3.2 + slope_n = kalbach_slope(energy_projectile, energy_emitted, 1, + za_emitted, za_target) + return slope_n * np.sqrt(0.5*energy_projectile/NEUTRON_MASS_EV)*np.minimum(4,np.maximum(1,9.3/np.sqrt(energy_emitted/EV_PER_MEV))) + # Special handling of elemental carbon if za_emitted == 6000: za_emitted = 6012 @@ -268,6 +269,8 @@ class KalbachMann(AngleEnergy): slope : Iterable of openmc.data.Tabulated1D Kalbach-Chadwick angular distribution slope value 'a' as a function of outgoing energy for each incoming energy + particle : {'neutron', 'photon'} + Incident particle type, defaults to neutron Attributes ---------- @@ -285,11 +288,13 @@ class KalbachMann(AngleEnergy): slope : Iterable of openmc.data.Tabulated1D Kalbach-Chadwick angular distribution slope value 'a' as a function of outgoing energy for each incoming energy + particle : {'neutron', 'photon'} + incident particle type, defaults to neutron """ def __init__(self, breakpoints, interpolation, energy, energy_out, - precompound, slope): + precompound, slope, particle = 'neutron'): super().__init__() self.breakpoints = breakpoints self.interpolation = interpolation @@ -297,6 +302,7 @@ def __init__(self, breakpoints, interpolation, energy, energy_out, self.energy_out = energy_out self.precompound = precompound self.slope = slope + self.particle = particle @property def breakpoints(self): @@ -317,6 +323,15 @@ def interpolation(self, interpolation): cv.check_type('Kalbach-Mann interpolation', interpolation, Iterable, Integral) self._interpolation = interpolation + + @property + def particle(self): + return self._particle + + @particle.setter + def particle(self, particle): + cv.check_value('Kalbach-Mann incident particle', particle, ['neutron', 'photon']) + self._particle = particle @property def energy(self): @@ -367,6 +382,7 @@ def to_hdf5(self, group): """ group.attrs['type'] = np.bytes_('kalbach-mann') + group.attrs['particle'] = np.bytes_(self.particle) dset = group.create_dataset('energy', data=self.energy) dset.attrs['interpolation'] = np.vstack((self.breakpoints, @@ -436,6 +452,7 @@ def from_hdf5(cls, group): Kalbach-Mann energy distribution """ + particle = group.attrs.get("particle", b"neutron").decode() interp_data = group['energy'].attrs['interpolation'] energy_breakpoints = interp_data[0, :] energy_interpolation = interp_data[1, :] @@ -491,7 +508,7 @@ def from_hdf5(cls, group): slope.append(km_a) return cls(energy_breakpoints, energy_interpolation, - energy, energy_out, precompound, slope) + energy, energy_out, precompound, slope, particle = particle) @classmethod def from_ace(cls, ace, idx, ldis): @@ -514,6 +531,7 @@ def from_ace(cls, ace, idx, ldis): Kalbach-Mann energy-angle distribution """ + particle = {'u':'photon', 'c':'neutron'}[ace.data_type.value] # Read number of interpolation regions and incoming energies n_regions = int(ace.xss[idx]) n_energy_in = int(ace.xss[idx + 1 + 2*n_regions]) @@ -586,10 +604,10 @@ def from_ace(cls, ace, idx, ldis): km_r.append(Tabulated1D(data[0], data[3])) km_a.append(Tabulated1D(data[0], data[4])) - return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a) + return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a, particle = particle) @classmethod - def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): + def from_endf(cls, file_obj, za_emitted, za_target, za_projectile): """Generate Kalbach-Mann distribution from an ENDF evaluation. If the projectile is a neutron, the slope is calculated when it is @@ -606,14 +624,8 @@ def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): ZA identifier of the emitted particle za_target : int ZA identifier of the target - projectile_mass : float - Mass of the projectile - - Warns - ----- - UserWarning - If the mass of the projectile is not equal to 1 (other than - a neutron), the slope is not calculated and set to 0 if missing. + za_projectile : int + ZA identifier of the projectile Returns ------- @@ -621,6 +633,7 @@ def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): Kalbach-Mann energy-angle distribution """ + particle = {0: 'photon', 1: 'neutron'}[za_projectile] params, tab2 = get_tab2_record(file_obj) lep = params[3] ne = params[5] @@ -654,31 +667,19 @@ def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): a_i = values[:, 3] calculated_slope.append(False) else: - # Check if the projectile is not a neutron - if not np.isclose(projectile_mass, 1.0, atol=1.0e-12, rtol=0.): - warn( - "Kalbach-Mann slope calculation is only available with " - "neutrons as projectile. Slope coefficients are set to 0." - ) - a_i = np.zeros_like(r_i) - calculated_slope.append(False) - - else: - # TODO: retrieve ZA of the projectile - za_projectile = 1 - a_i = [kalbach_slope(energy_projectile=energy[i], - energy_emitted=e, - za_projectile=za_projectile, - za_emitted=za_emitted, - za_target=za_target) - for e in eout_i] - calculated_slope.append(True) + a_i = [kalbach_slope(energy_projectile=energy[i], + energy_emitted=e, + za_projectile=za_projectile, + za_emitted=za_emitted, + za_target=za_target) + for e in eout_i] + calculated_slope.append(True) precompound.append(Tabulated1D(eout_i, r_i)) slope.append(Tabulated1D(eout_i, a_i)) km_distribution = cls(tab2.breakpoints, tab2.interpolation, energy, - energy_out, precompound, slope) + energy_out, precompound, slope, particle = particle) # List of bool to indicate slope calculation by OpenMC km_distribution._calculated_slope = calculated_slope diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 65b59582cf8..218a5eac8e3 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -87,6 +87,8 @@ def _get_products(ev, mt): is not defined in the 'center-of-mass' system. The breakup logic is not implemented which can lead to this error being raised while the definition of the product is correct. + NotImplementedError + When the projectile is not a neutron and not a photon. Returns ------- @@ -163,11 +165,16 @@ def _get_products(ev, mt): ) zat = ev.target["atomic_number"] * 1000 + ev.target["mass_number"] - projectile_mass = ev.projectile["mass"] + if ev.projectile['mass'] == 0.0: + za_projectile = 0 + elif np.isclose(ev.projectile['mass'], 1.0, atol=1.0e-12, rtol=0.): + za_projectile = 1 + else: + raise NotImplementedError('Unknown projectile') p.distribution = [KalbachMann.from_endf(file_obj, za, zat, - projectile_mass)] + za_projectile)] elif law == 2: # Discrete two-body scattering diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 4a7878343c6..15a3bab2045 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -22,6 +22,14 @@ namespace openmc { KalbachMann::KalbachMann(hid_t group) { + is_photon_ = false; + // Check if projectile is a photon + if (attribute_exists(group, "particle")) { + std::string temp; + read_attribute(group, "particle", temp); + if (temp == "photon") + is_photon_ = true; + } // Open incoming energy dataset hid_t dset = open_dataset(group, "energy"); @@ -229,12 +237,19 @@ void KalbachMann::sample( } // Sampled correlated angle from Kalbach-Mann parameters - if (prn(seed) > km_r) { - double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a); - mu = std::log(T + std::sqrt(T * T + 1.0)) / km_a; + if (is_photon_) { + double T = uniform_distribution(0., 1., seed); + double sinha = std::sinh(km_a); + mu = sinha * ((1 + T) - 2 * km_r) / + (sinha * T + std::cosh(km_a) - std::exp(km_a * T)); } else { - double r1 = prn(seed); - mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a; + if (prn(seed) > km_r) { + double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a); + mu = std::log(T + std::sqrt(T * T + 1.0)) / km_a; + } else { + double r1 = prn(seed); + mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a; + } } } diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py index 5d06669f7b9..fca2d34e2b3 100644 --- a/tests/unit_tests/test_data_kalbach_mann.py +++ b/tests/unit_tests/test_data_kalbach_mann.py @@ -103,16 +103,13 @@ def test_kalbach_slope(): energy_projectile = 10.2 # [eV] energy_emitted = 5.4 # [eV] - # Check that NotImplementedError is raised if the projectile is not - # a neutron - with pytest.raises(NotImplementedError): - kalbach_slope( - energy_projectile=energy_projectile, - energy_emitted=energy_emitted, - za_projectile=1000, - za_emitted=1, - za_target=6012 - ) + kalbach_slope( + energy_projectile=energy_projectile, + energy_emitted=energy_emitted, + za_projectile=0, + za_emitted=1, + za_target=6012 + ) assert kalbach_slope( energy_projectile=energy_projectile,