26783873914

Committed 01 Jun 2026 09:43PM UTC coverage: 81.362% (+0.03%) from 81.333%

Build # 26783873914

Build Type

Pull #3948

github

Committed by

web-flow

Commit Message

Merge 6314ea576 into 111eb7706

Pull Request Pull Request #3948: Fix get_index_in_direction for regular meshes

Coverage Stats

18027 of 26121 branches covered (69.01%)

Branch coverage included in aggregate %.

43 of 45 new or added lines in 9 files covered. (95.56%)

443 existing lines in 12 files now uncovered.

59175 of 68766 relevant lines covered (86.05%)

48551677.44 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

91.08

/openmc/data/angle_distribution.py

from collections.abc import Iterable
from io import StringIO
from numbers import Real
from warnings import warn

import numpy as np

import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from openmc.stats import Univariate, Tabular, Uniform, Legendre
from .function import INTERPOLATION_SCHEME
from .data import EV_PER_MEV
from .endf import as_evaluation, get_head_record, get_cont_record, \
    get_tab1_record, get_list_record, get_tab2_record


class AngleDistribution(EqualityMixin):
    """Angle distribution as a function of incoming energy

    Parameters
    ----------
    energy : Iterable of float
        Incoming energies in eV at which distributions exist
    mu : Iterable of openmc.stats.Univariate
        Distribution of scattering cosines corresponding to each incoming energy

    Attributes
    ----------
    energy : Iterable of float
        Incoming energies in eV at which distributions exist
    mu : Iterable of openmc.stats.Univariate
        Distribution of scattering cosines corresponding to each incoming energy

    """

    def __init__(self, energy, mu):
        super().__init__()
        self.energy = energy
        self.mu = mu

    @property
    def energy(self):
        return self._energy

    @energy.setter
    def energy(self, energy):
        cv.check_type('angle distribution incoming energy', energy,
                      Iterable, Real)
        self._energy = energy

    @property
    def mu(self):
        return self._mu

    @mu.setter
    def mu(self, mu):
        cv.check_type('angle distribution scattering cosines', mu,
                      Iterable, Univariate)
        self._mu = mu

    def to_hdf5(self, group):
        """Write angle distribution to an HDF5 group

        Parameters
        ----------
        group : h5py.Group
            HDF5 group to write to

        """

        dset = group.create_dataset('energy', data=self.energy)

        # Make sure all data is tabular
        mu_tabular = [mu_i if isinstance(mu_i, Tabular) else
                      mu_i.to_tabular() for mu_i in self.mu]

        # Determine total number of (mu,p) pairs and create array
        n_pairs = sum([len(mu_i.x) for mu_i in mu_tabular])
        pairs = np.empty((3, n_pairs))

        # Create array for offsets
        offsets = np.empty(len(mu_tabular), dtype=int)
        interpolation = np.empty(len(mu_tabular), dtype=int)
        j = 0

        # Populate offsets and pairs array
        for i, mu_i in enumerate(mu_tabular):
            n = len(mu_i.x)
            offsets[i] = j
            interpolation[i] = 1 if mu_i.interpolation == 'histogram' else 2
            pairs[0, j:j+n] = mu_i.x
            pairs[1, j:j+n] = mu_i.p
            pairs[2, j:j+n] = mu_i.c
            j += n

        # Create dataset for distributions
        dset = group.create_dataset('mu', data=pairs)

        # Write interpolation as attribute
        dset.attrs['offsets'] = offsets
        dset.attrs['interpolation'] = interpolation

    @classmethod
    def from_hdf5(cls, group):
        """Generate angular distribution from HDF5 data

        Parameters
        ----------
        group : h5py.Group
            HDF5 group to read from

        Returns
        -------
        openmc.data.AngleDistribution
            Angular distribution

        """
        energy = group['energy'][()]
        data = group['mu']
        offsets = data.attrs['offsets']
        interpolation = data.attrs['interpolation']

        mu = []
        n_energy = len(energy)
        for i in range(n_energy):
            # Determine length of outgoing energy distribution and number of
            # discrete lines
            j = offsets[i]
            if i < n_energy - 1:
                n = offsets[i+1] - j
            else:
                n = data.shape[1] - j

            interp = INTERPOLATION_SCHEME[interpolation[i]]
            mu_i = Tabular(data[0, j:j+n], data[1, j:j+n], interp)
            mu_i.c = data[2, j:j+n]

            mu.append(mu_i)

        return cls(energy, mu)

    @classmethod
    def from_ace(cls, ace, location_dist, location_start):
        """Generate an angular distribution from ACE data

        Parameters
        ----------
        ace : openmc.data.ace.Table
            ACE table to read from
        location_dist : int
            Index in the XSS array corresponding to the start of a block,
            e.g. JXS(9).
        location_start : int
            Index in the XSS array corresponding to the start of an angle
            distribution array

        Returns
        -------
        openmc.data.AngleDistribution
            Angular distribution

        """
        # Set starting index for angle distribution
        idx = location_dist + location_start - 1

        # Number of energies at which angular distributions are tabulated
        n_energies = int(ace.xss[idx])
        idx += 1

        # Incoming energy grid
        energy = ace.xss[idx:idx + n_energies]*EV_PER_MEV
        idx += n_energies

        # Read locations for angular distributions
        lc = ace.xss[idx:idx + n_energies].astype(int)
        idx += n_energies

        mu = []
        for i in range(n_energies):
            if lc[i] > 0:
                # Equiprobable 32 bin distribution
                n_bins = 32
                idx = location_dist + abs(lc[i]) - 1
                cos = ace.xss[idx:idx + n_bins + 1]
                pdf = np.zeros(n_bins + 1)
                pdf[:n_bins] = 1.0/(n_bins*np.diff(cos))
                cdf = np.linspace(0.0, 1.0, n_bins + 1)

                mu_i = Tabular(cos, pdf, 'histogram', ignore_negative=True)
                mu_i.c = cdf
            elif lc[i] < 0:
                # Tabular angular distribution
                idx = location_dist + abs(lc[i]) - 1
                intt = int(ace.xss[idx])
                n_points = int(ace.xss[idx + 1])
                # Data is given as rows of (values, PDF, CDF)
                data = ace.xss[idx + 2:idx + 2 + 3*n_points]
                data.shape = (3, n_points)

                mu_i = Tabular(data[0], data[1], INTERPOLATION_SCHEME[intt])
                mu_i.c = data[2]
            else:
                # Isotropic angular distribution
                mu_i = Uniform(-1., 1.)

            mu.append(mu_i)

        return cls(energy, mu)

    @classmethod
    def from_endf(cls, ev, mt):
        """Generate an angular distribution from an ENDF evaluation

        Parameters
        ----------
        ev : openmc.data.endf.Evaluation or endf.Material
            ENDF evaluation
        mt : int
            The MT value of the reaction to get angular distributions for

        Returns
        -------
        openmc.data.AngleDistribution
            Angular distribution

        """
        ev = as_evaluation(ev)
        file_obj = StringIO(ev.section[4, mt])

        # Read HEAD record
        items = get_head_record(file_obj)
        lvt = items[2]
        ltt = items[3]

        # Read CONT record
        items = get_cont_record(file_obj)
        li = items[2]
        nk = items[4]

        # Check for obsolete energy transformation matrix. If present, just skip
        # it and keep reading
        if lvt > 0:
            warn('Obsolete energy transformation matrix in MF=4 angular '
                 'distribution.')
            for _ in range((nk + 5)//6):
                file_obj.readline()

        if ltt == 0 and li == 1:
            # Purely isotropic
            energy = np.array([0., ev.info['energy_max']])
            mu = [Uniform(-1., 1.), Uniform(-1., 1.)]

        elif ltt == 1 and li == 0:
            # Legendre polynomial coefficients
            params, tab2 = get_tab2_record(file_obj)
            n_energy = params[5]

            energy = np.zeros(n_energy)
            mu = []
            for i in range(n_energy):
                items, al = get_list_record(file_obj)
                energy[i] = items[1]
                coefficients = np.asarray([1.0] + al)
                mu.append(Legendre(coefficients))

        elif ltt == 2 and li == 0:
            # Tabulated probability distribution
            params, tab2 = get_tab2_record(file_obj)
            n_energy = params[5]

            energy = np.zeros(n_energy)
            mu = []
            for i in range(n_energy):
                params, f = get_tab1_record(file_obj)
                energy[i] = params[1]
                if f.n_regions > 1:
                    raise NotImplementedError('Angular distribution with multiple '
                                              'interpolation regions not supported.')
                mu.append(Tabular(f.x, f.y, INTERPOLATION_SCHEME[f.interpolation[0]]))

        elif ltt == 3 and li == 0:
            # Legendre for low energies / tabulated for high energies
            params, tab2 = get_tab2_record(file_obj)
            n_energy_legendre = params[5]

            energy_legendre = np.zeros(n_energy_legendre)
            mu = []
            for i in range(n_energy_legendre):
                items, al = get_list_record(file_obj)
                energy_legendre[i] = items[1]
                coefficients = np.asarray([1.0] + al)
                mu.append(Legendre(coefficients))

            params, tab2 = get_tab2_record(file_obj)
            n_energy_tabulated = params[5]

            energy_tabulated = np.zeros(n_energy_tabulated)
            for i in range(n_energy_tabulated):
                params, f = get_tab1_record(file_obj)
                energy_tabulated[i] = params[1]
                if f.n_regions > 1:
                    raise NotImplementedError('Angular distribution with multiple '
                                              'interpolation regions not supported.')
                mu.append(Tabular(f.x, f.y, INTERPOLATION_SCHEME[f.interpolation[0]]))

            energy = np.concatenate((energy_legendre, energy_tabulated))

        return AngleDistribution(energy, mu)

1	from collections.abc import Iterable	11✔
2	from io import StringIO	11✔
3	from numbers import Real	11✔
4	from warnings import warn	11✔
5
6	import numpy as np	11✔
7
8	import openmc.checkvalue as cv	11✔
9	from openmc.mixin import EqualityMixin	11✔
10	from openmc.stats import Univariate, Tabular, Uniform, Legendre	11✔
11	from .function import INTERPOLATION_SCHEME	11✔
12	from .data import EV_PER_MEV	11✔
13	from .endf import as_evaluation, get_head_record, get_cont_record, \	11✔
14	get_tab1_record, get_list_record, get_tab2_record
15
16
17	class AngleDistribution(EqualityMixin):	11✔
18	"""Angle distribution as a function of incoming energy
19
20	Parameters
21	----------
22	energy : Iterable of float
23	Incoming energies in eV at which distributions exist
24	mu : Iterable of openmc.stats.Univariate
25	Distribution of scattering cosines corresponding to each incoming energy
26
27	Attributes
28	----------
29	energy : Iterable of float
30	Incoming energies in eV at which distributions exist
31	mu : Iterable of openmc.stats.Univariate
32	Distribution of scattering cosines corresponding to each incoming energy
33
34	"""
35
36	def __init__(self, energy, mu):	11✔
37	super().__init__()	11✔
38	self.energy = energy	11✔
39	self.mu = mu	11✔
40
41	@property	11✔
42	def energy(self):	11✔
43	return self._energy	11✔
44
45	@energy.setter	11✔
46	def energy(self, energy):	11✔
47	cv.check_type('angle distribution incoming energy', energy,	11✔
48	Iterable, Real)
49	self._energy = energy	11✔
50
51	@property	11✔
52	def mu(self):	11✔
53	return self._mu	11✔
54
55	@mu.setter	11✔
56	def mu(self, mu):	11✔
57	cv.check_type('angle distribution scattering cosines', mu,	11✔
58	Iterable, Univariate)
59	self._mu = mu	11✔
60
61	def to_hdf5(self, group):	11✔
62	"""Write angle distribution to an HDF5 group
63
64	Parameters
65	----------
66	group : h5py.Group
67	HDF5 group to write to
68
69	"""
70
71	dset = group.create_dataset('energy', data=self.energy)	11✔
72
73	# Make sure all data is tabular
74	mu_tabular = [mu_i if isinstance(mu_i, Tabular) else	11✔
75	mu_i.to_tabular() for mu_i in self.mu]
76
77	# Determine total number of (mu,p) pairs and create array
78	n_pairs = sum([len(mu_i.x) for mu_i in mu_tabular])	11✔
79	pairs = np.empty((3, n_pairs))	11✔
80
81	# Create array for offsets
82	offsets = np.empty(len(mu_tabular), dtype=int)	11✔
83	interpolation = np.empty(len(mu_tabular), dtype=int)	11✔
84	j = 0	11✔
85
86	# Populate offsets and pairs array
87	for i, mu_i in enumerate(mu_tabular):	11✔
88	n = len(mu_i.x)	11✔
89	offsets[i] = j	11✔
90	interpolation[i] = 1 if mu_i.interpolation == 'histogram' else 2	11✔
91	pairs[0, j:j+n] = mu_i.x	11✔
92	pairs[1, j:j+n] = mu_i.p	11✔
93	pairs[2, j:j+n] = mu_i.c	11✔
94	j += n	11✔
95
96	# Create dataset for distributions
97	dset = group.create_dataset('mu', data=pairs)	11✔
98
99	# Write interpolation as attribute
100	dset.attrs['offsets'] = offsets	11✔
101	dset.attrs['interpolation'] = interpolation	11✔
102
103	@classmethod	11✔
104	def from_hdf5(cls, group):	11✔
105	"""Generate angular distribution from HDF5 data
106
107	Parameters
108	----------
109	group : h5py.Group
110	HDF5 group to read from
111
112	Returns
113	-------
114	openmc.data.AngleDistribution
115	Angular distribution
116
117	"""
118	energy = group['energy'][()]	11✔
119	data = group['mu']	11✔
120	offsets = data.attrs['offsets']	11✔
121	interpolation = data.attrs['interpolation']	11✔
122
123	mu = []	11✔
124	n_energy = len(energy)	11✔
125	for i in range(n_energy):	11✔
126	# Determine length of outgoing energy distribution and number of
127	# discrete lines
128	j = offsets[i]	11✔
129	if i < n_energy - 1:	11✔
130	n = offsets[i+1] - j	11✔
131	else:
132	n = data.shape[1] - j	11✔
133
134	interp = INTERPOLATION_SCHEME[interpolation[i]]	11✔
135	mu_i = Tabular(data[0, j:j+n], data[1, j:j+n], interp)	11✔
136	mu_i.c = data[2, j:j+n]	11✔
137
138	mu.append(mu_i)	11✔
139
140	return cls(energy, mu)	11✔
141
142	@classmethod	11✔
143	def from_ace(cls, ace, location_dist, location_start):	11✔
144	"""Generate an angular distribution from ACE data
145
146	Parameters
147	----------
148	ace : openmc.data.ace.Table
149	ACE table to read from
150	location_dist : int
151	Index in the XSS array corresponding to the start of a block,
152	e.g. JXS(9).
153	location_start : int
154	Index in the XSS array corresponding to the start of an angle
155	distribution array
156
157	Returns
158	-------
159	openmc.data.AngleDistribution
160	Angular distribution
161
162	"""
163	# Set starting index for angle distribution
164	idx = location_dist + location_start - 1	11✔
165
166	# Number of energies at which angular distributions are tabulated
167	n_energies = int(ace.xss[idx])	11✔
168	idx += 1	11✔
169
170	# Incoming energy grid
171	energy = ace.xss[idx:idx + n_energies]*EV_PER_MEV	11✔
172	idx += n_energies	11✔
173
174	# Read locations for angular distributions
175	lc = ace.xss[idx:idx + n_energies].astype(int)	11✔
176	idx += n_energies	11✔
177
178	mu = []	11✔
179	for i in range(n_energies):	11✔
180	if lc[i] > 0:	11✔
181	# Equiprobable 32 bin distribution
182	n_bins = 32	×
183	idx = location_dist + abs(lc[i]) - 1	×
184	cos = ace.xss[idx:idx + n_bins + 1]	×
185	pdf = np.zeros(n_bins + 1)	×
186	pdf[:n_bins] = 1.0/(n_bins*np.diff(cos))	×
187	cdf = np.linspace(0.0, 1.0, n_bins + 1)	×
188
189	mu_i = Tabular(cos, pdf, 'histogram', ignore_negative=True)	×
190	mu_i.c = cdf	×
191	elif lc[i] < 0:	11✔
192	# Tabular angular distribution
193	idx = location_dist + abs(lc[i]) - 1	11✔
194	intt = int(ace.xss[idx])	11✔
195	n_points = int(ace.xss[idx + 1])	11✔
196	# Data is given as rows of (values, PDF, CDF)
197	data = ace.xss[idx + 2:idx + 2 + 3*n_points]	11✔
198	data.shape = (3, n_points)	11✔
199
200	mu_i = Tabular(data[0], data[1], INTERPOLATION_SCHEME[intt])	11✔
201	mu_i.c = data[2]	11✔
202	else:
203	# Isotropic angular distribution
204	mu_i = Uniform(-1., 1.)	×
205
206	mu.append(mu_i)	11✔
207
208	return cls(energy, mu)	11✔
209
210	@classmethod	11✔
211	def from_endf(cls, ev, mt):	11✔
212	"""Generate an angular distribution from an ENDF evaluation
213
214	Parameters
215	----------
216	ev : openmc.data.endf.Evaluation or endf.Material
217	ENDF evaluation
218	mt : int
219	The MT value of the reaction to get angular distributions for
220
221	Returns
222	-------
223	openmc.data.AngleDistribution
224	Angular distribution
225
226	"""
227	ev = as_evaluation(ev)	11✔
228	file_obj = StringIO(ev.section[4, mt])	11✔
229
230	# Read HEAD record
231	items = get_head_record(file_obj)	11✔
232	lvt = items[2]	11✔
233	ltt = items[3]	11✔
234
235	# Read CONT record
236	items = get_cont_record(file_obj)	11✔
237	li = items[2]	11✔
238	nk = items[4]	11✔
239
240	# Check for obsolete energy transformation matrix. If present, just skip
241	# it and keep reading
242	if lvt > 0:	11✔
UNCOV 243	warn('Obsolete energy transformation matrix in MF=4 angular '	×
244	'distribution.')
245	for _ in range((nk + 5)//6):	×
UNCOV 246	file_obj.readline()	×
247
248	if ltt == 0 and li == 1:	11✔
249	# Purely isotropic
250	energy = np.array([0., ev.info['energy_max']])	11✔
251	mu = [Uniform(-1., 1.), Uniform(-1., 1.)]	11✔
252
253	elif ltt == 1 and li == 0:	11✔
254	# Legendre polynomial coefficients
255	params, tab2 = get_tab2_record(file_obj)	11✔
256	n_energy = params[5]	11✔
257
258	energy = np.zeros(n_energy)	11✔
259	mu = []	11✔
260	for i in range(n_energy):	11✔
261	items, al = get_list_record(file_obj)	11✔
262	energy[i] = items[1]	11✔
263	coefficients = np.asarray([1.0] + al)	11✔
264	mu.append(Legendre(coefficients))	11✔
265
266	elif ltt == 2 and li == 0:	11✔
267	# Tabulated probability distribution
268	params, tab2 = get_tab2_record(file_obj)	11✔
269	n_energy = params[5]	11✔
270
271	energy = np.zeros(n_energy)	11✔
272	mu = []	11✔
273	for i in range(n_energy):	11✔
274	params, f = get_tab1_record(file_obj)	11✔
275	energy[i] = params[1]	11✔
276	if f.n_regions > 1:	11✔
UNCOV 277	raise NotImplementedError('Angular distribution with multiple '	×
278	'interpolation regions not supported.')
279	mu.append(Tabular(f.x, f.y, INTERPOLATION_SCHEME[f.interpolation[0]]))	11✔
280
281	elif ltt == 3 and li == 0:	11✔
282	# Legendre for low energies / tabulated for high energies
283	params, tab2 = get_tab2_record(file_obj)	11✔
284	n_energy_legendre = params[5]	11✔
285
286	energy_legendre = np.zeros(n_energy_legendre)	11✔
287	mu = []	11✔
288	for i in range(n_energy_legendre):	11✔
289	items, al = get_list_record(file_obj)	11✔
290	energy_legendre[i] = items[1]	11✔
291	coefficients = np.asarray([1.0] + al)	11✔
292	mu.append(Legendre(coefficients))	11✔
293
294	params, tab2 = get_tab2_record(file_obj)	11✔
295	n_energy_tabulated = params[5]	11✔
296
297	energy_tabulated = np.zeros(n_energy_tabulated)	11✔
298	for i in range(n_energy_tabulated):	11✔
299	params, f = get_tab1_record(file_obj)	11✔
300	energy_tabulated[i] = params[1]	11✔
301	if f.n_regions > 1:	11✔
UNCOV 302	raise NotImplementedError('Angular distribution with multiple '	×
303	'interpolation regions not supported.')
304	mu.append(Tabular(f.x, f.y, INTERPOLATION_SCHEME[f.interpolation[0]]))	11✔
305
306	energy = np.concatenate((energy_legendre, energy_tabulated))	11✔
307
308	return AngleDistribution(energy, mu)	11✔

openmc-dev / openmc / 26783873914

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous