From 9483cce0bc334467748cd842cb122fefc6f3305d Mon Sep 17 00:00:00 2001 From: Jon Shimwell Date: Wed, 14 Aug 2024 15:34:17 +0100 Subject: [PATCH] Remove resonance reconstruction and Cython dependency (#3111) Co-authored-by: Paul Romano --- Dockerfile | 2 +- MANIFEST.in | 2 - docs/source/usersguide/install.rst | 4 - openmc/data/function.py | 22 -- openmc/data/neutron.py | 88 +---- openmc/data/njoy.py | 15 +- openmc/data/reconstruct.pyx | 522 -------------------------- pyproject.toml | 3 +- setup.py | 14 - tests/unit_tests/test_data_neutron.py | 27 +- tools/ci/gha-install.sh | 3 +- 11 files changed, 33 insertions(+), 669 deletions(-) delete mode 100644 openmc/data/reconstruct.pyx delete mode 100755 setup.py diff --git a/Dockerfile b/Dockerfile index 35b9cf578d..4a94e3c0b1 100644 --- a/Dockerfile +++ b/Dockerfile @@ -95,7 +95,7 @@ RUN cd $HOME \ RUN if [ "$build_dagmc" = "on" ]; then \ # Install addition packages required for DAGMC apt-get -y install libeigen3-dev libnetcdf-dev libtbb-dev libglfw3-dev \ - && pip install --upgrade numpy "cython<3.0" \ + && pip install --upgrade numpy \ # Clone and install EMBREE && mkdir -p $HOME/EMBREE && cd $HOME/EMBREE \ && git clone --single-branch -b ${EMBREE_TAG} --depth 1 ${EMBREE_REPO} \ diff --git a/MANIFEST.in b/MANIFEST.in index afd016cb02..b73218af0d 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -26,8 +26,6 @@ recursive-include include *.h recursive-include include *.h.in recursive-include include *.hh recursive-include man *.1 -recursive-include openmc *.pyx -recursive-include openmc *.c recursive-include src *.cc recursive-include src *.cpp recursive-include src *.rnc diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 130e96c0ae..1c0b7fa5b3 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -584,10 +584,6 @@ distributions. parallel runs. This package is needed if you plan on running depletion simulations in parallel using MPI. - `Cython `_ - Cython is used for resonance reconstruction for ENDF data converted to - :class:`openmc.data.IncidentNeutron`. - `vtk `_ The Python VTK bindings are needed to convert voxel and track files to VTK format. diff --git a/openmc/data/function.py b/openmc/data/function.py index 23fd5e9d4f..c5914f513d 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -708,28 +708,6 @@ class ResonancesWithBackground(EqualityMixin): self.background = background self.mt = mt - def __call__(self, x): - # Get background cross section - xs = self.background(x) - - for r in self.resonances: - if not isinstance(r, openmc.data.resonance._RESOLVED): - continue - - if isinstance(x, Iterable): - # Determine which energies are within resolved resonance range - within = (r.energy_min <= x) & (x <= r.energy_max) - - # Get resonance cross sections and add to background - resonant_xs = r.reconstruct(x[within]) - xs[within] += resonant_xs[self.mt] - else: - if r.energy_min <= x <= r.energy_max: - resonant_xs = r.reconstruct(x) - xs += resonant_xs[self.mt] - - return xs - @property def background(self): return self._background diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 894be18717..95a3424ea4 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -16,8 +16,7 @@ from .endf import ( Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_evaluations) from .fission_energy import FissionEnergyRelease from .function import Tabulated1D, Sum, ResonancesWithBackground -from .grid import linearize, thin -from .njoy import make_ace +from .njoy import make_ace, make_pendf from .product import Product from .reaction import Reaction, _get_photon_products_ace, FISSION_MTS from . import resonance as res @@ -286,7 +285,7 @@ class IncidentNeutron(EqualityMixin): if strT in data.urr: self.urr[strT] = data.urr[strT] - def add_elastic_0K_from_endf(self, filename, overwrite=False): + def add_elastic_0K_from_endf(self, filename, overwrite=False, **kwargs): """Append 0K elastic scattering cross section from an ENDF file. Parameters @@ -297,6 +296,8 @@ class IncidentNeutron(EqualityMixin): If existing 0 K data is present, this flag can be used to indicate that it should be overwritten. Otherwise, an exception will be thrown. + **kwargs + Keyword arguments passed to :func:`openmc.data.njoy.make_pendf` Raises ------ @@ -309,75 +310,22 @@ class IncidentNeutron(EqualityMixin): if '0K' in self.energy and not overwrite: raise ValueError('0 K data already exists for this nuclide.') - data = type(self).from_endf(filename) - if data.resonances is not None: - x = [] - y = [] - for rr in data.resonances: - if isinstance(rr, res.RMatrixLimited): - raise TypeError('R-Matrix Limited not supported.') - elif isinstance(rr, res.Unresolved): - continue + with tempfile.TemporaryDirectory() as tmpdir: + # Set arguments for make_pendf + pendf_path = os.path.join(tmpdir, 'pendf') + kwargs.setdefault('output_dir', tmpdir) + kwargs.setdefault('pendf', pendf_path) - # Get energies/widths for resonances - e_peak = rr.parameters['energy'].values - if isinstance(rr, res.MultiLevelBreitWigner): - gamma = rr.parameters['totalWidth'].values - elif isinstance(rr, res.ReichMoore): - df = rr.parameters - gamma = (df['neutronWidth'] + - df['captureWidth'] + - abs(df['fissionWidthA']) + - abs(df['fissionWidthB'])).values + # Run NJOY to create a pointwise ENDF file + make_pendf(filename, **kwargs) - # Determine peak energies and widths - e_min, e_max = rr.energy_min, rr.energy_max - in_range = (e_peak > e_min) & (e_peak < e_max) - e_peak = e_peak[in_range] - gamma = gamma[in_range] - - # Get midpoints between resonances (use min/max energy of - # resolved region as absolute lower/upper bound) - e_mid = np.concatenate( - ([e_min], (e_peak[1:] + e_peak[:-1])/2, [e_max])) - - # Add grid around each resonance that includes the peak +/- the - # width times each value in _RESONANCE_ENERGY_GRID. Values are - # constrained so that points around one resonance don't overlap - # with points around another. This algorithm is from Fudge - # (https://doi.org/10.1063/1.1945057). - energies = [] - for e, g, e_lower, e_upper in zip(e_peak, gamma, e_mid[:-1], - e_mid[1:]): - e_left = e - g*_RESONANCE_ENERGY_GRID - energies.append(e_left[e_left > e_lower][::-1]) - e_right = e + g*_RESONANCE_ENERGY_GRID[1:] - energies.append(e_right[e_right < e_upper]) - - # Concatenate all points - energies = np.concatenate(energies) - - # Create 1000 equal log-spaced energies over RRR, combine with - # resonance peaks and half-height energies - e_log = np.logspace(log10(e_min), log10(e_max), 1000) - energies = np.union1d(e_log, energies) - - # Linearize and thin cross section - xi, yi = linearize(energies, data[2].xs['0K']) - xi, yi = thin(xi, yi) - - # If there are multiple resolved resonance ranges (e.g. Pu239 in - # ENDF/B-VII.1), combine them - x = np.concatenate((x, xi)) - y = np.concatenate((y, yi)) - else: - energies = data[2].xs['0K'].x - x, y = linearize(energies, data[2].xs['0K']) - x, y = thin(x, y) - - # Set 0K energy grid and elastic scattering cross section - self.energy['0K'] = x - self[2].xs['0K'] = Tabulated1D(x, y) + # Add 0K elastic scattering cross section + pendf = Evaluation(pendf_path) + file_obj = StringIO(pendf.section[3, 2]) + get_head_record(file_obj) + params, xs = get_tab1_record(file_obj) + self.energy['0K'] = xs.x + self[2].xs['0K'] = xs def get_reaction_components(self, mt): """Determine what reactions make up redundant reaction. diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index ac1b5e345e..1bf44891ef 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -221,7 +221,7 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False, shutil.move(tmpfilename, str(filename)) -def make_pendf(filename, pendf='pendf', error=0.001, stdout=False): +def make_pendf(filename, pendf='pendf', **kwargs): """Generate pointwise ENDF file from an ENDF file Parameters @@ -230,10 +230,9 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False): Path to ENDF file pendf : str, optional Path of pointwise ENDF file to write - error : float, optional - Fractional error tolerance for NJOY processing - stdout : bool - Whether to display NJOY standard output + **kwargs + Keyword arguments passed to :func:`openmc.data.njoy.make_ace`. All NJOY + module arguments other than pendf default to False. Raises ------ @@ -241,9 +240,9 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False): If the NJOY process returns with a non-zero status """ - - make_ace(filename, pendf=pendf, error=error, broadr=False, - heatr=False, purr=False, acer=False, stdout=stdout) + for key in ('broadr', 'heatr', 'gaspr', 'purr', 'acer'): + kwargs.setdefault(key, False) + make_ace(filename, pendf=pendf, **kwargs) def make_ace(filename, temperatures=None, acer=True, xsdir=None, diff --git a/openmc/data/reconstruct.pyx b/openmc/data/reconstruct.pyx deleted file mode 100644 index cd0bbc38b9..0000000000 --- a/openmc/data/reconstruct.pyx +++ /dev/null @@ -1,522 +0,0 @@ -from libc.stdlib cimport malloc, calloc, free -from libc.math cimport cos, sin, sqrt, atan, M_PI - -cimport numpy as np -import numpy as np -from numpy.linalg import inv -cimport cython - - -cdef extern from "complex.h": - double cabs(double complex) - double complex conj(double complex) - double creal(complex double) - double cimag(complex double) - double complex cexp(double complex) - -# Physical constants are from CODATA 2014 -cdef double NEUTRON_MASS_ENERGY = 939.5654133e6 # eV/c^2 -cdef double HBAR_C = 197.3269788e5 # eV-b^0.5 - - -@cython.cdivision(True) -def wave_number(double A, double E): - r"""Neutron wave number in center-of-mass system. - - ENDF-102 defines the neutron wave number in the center-of-mass system in - Equation D.10 as - - .. math:: - k = \frac{2m_n}{\hbar} \frac{A}{A + 1} \sqrt{|E|} - - Parameters - ---------- - A : double - Ratio of target mass to neutron mass - E : double - Energy in eV - - Returns - ------- - double - Neutron wave number in b^-0.5 - - """ - return A/(A + 1)*sqrt(2*NEUTRON_MASS_ENERGY*abs(E))/HBAR_C - -@cython.cdivision(True) -cdef double _wave_number(double A, double E): - return A/(A + 1)*sqrt(2*NEUTRON_MASS_ENERGY*abs(E))/HBAR_C - - -@cython.cdivision(True) -cdef double phaseshift(int l, double rho): - """Calculate hardsphere phase shift as given in ENDF-102, Equation D.13 - - Parameters - ---------- - l : int - Angular momentum quantum number - rho : float - Product of the wave number and the channel radius - - Returns - ------- - double - Hardsphere phase shift - - """ - if l == 0: - return rho - elif l == 1: - return rho - atan(rho) - elif l == 2: - return rho - atan(3*rho/(3 - rho**2)) - elif l == 3: - return rho - atan((15*rho - rho**3)/(15 - 6*rho**2)) - elif l == 4: - return rho - atan((105*rho - 10*rho**3)/(105 - 45*rho**2 + rho**4)) - - -@cython.cdivision(True) -def penetration_shift(int l, double rho): - r"""Calculate shift and penetration factors as given in ENDF-102, Equations D.11 - and D.12. - - Parameters - ---------- - l : int - Angular momentum quantum number - rho : float - Product of the wave number and the channel radius - - Returns - ------- - double - Penetration factor for given :math:`l` - double - Shift factor for given :math:`l` - - """ - cdef double den - - if l == 0: - return rho, 0. - elif l == 1: - den = 1 + rho**2 - return rho**3/den, -1/den - elif l == 2: - den = 9 + 3*rho**2 + rho**4 - return rho**5/den, -(18 + 3*rho**2)/den - elif l == 3: - den = 225 + 45*rho**2 + 6*rho**4 + rho**6 - return rho**7/den, -(675 + 90*rho**2 + 6*rho**4)/den - elif l == 4: - den = 11025 + 1575*rho**2 + 135*rho**4 + 10*rho**6 + rho**8 - return rho**9/den, -(44100 + 4725*rho**2 + 270*rho**4 + 10*rho**6)/den - - -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -def reconstruct_mlbw(mlbw, double E): - """Evaluate cross section using MLBW data. - - Parameters - ---------- - mlbw : openmc.data.MultiLevelBreitWigner - Multi-level Breit-Wigner resonance parameters - E : double - Energy in eV at which to evaluate the cross section - - Returns - ------- - elastic : double - Elastic scattering cross section in barns - capture : double - Radiative capture cross section in barns - fission : double - Fission cross section in barns - - """ - cdef int i, nJ, ij, l, n_res, i_res - cdef double elastic, capture, fission - cdef double A, k, rho, rhohat, I - cdef double P, S, phi, cos2phi, sin2phi - cdef double Ex, Q, rhoc, rhochat, P_c, S_c - cdef double jmin, jmax, j, Dl - cdef double E_r, gt, gn, gg, gf, gx, P_r, S_r, P_rx - cdef double gnE, gtE, Eprime, x, f - cdef double *g - cdef double (*s)[2] - cdef double [:,:] params - - I = mlbw.target_spin - A = mlbw.atomic_weight_ratio - k = _wave_number(A, E) - - elastic = 0. - capture = 0. - fission = 0. - - for i, l in enumerate(mlbw._l_values): - params = mlbw._parameter_matrix[l] - - rho = k*mlbw.channel_radius[l](E) - rhohat = k*mlbw.scattering_radius[l](E) - P, S = penetration_shift(l, rho) - phi = phaseshift(l, rhohat) - cos2phi = cos(2*phi) - sin2phi = sin(2*phi) - - # Determine shift and penetration at modified energy - if mlbw._competitive[i]: - Ex = E + mlbw.q_value[l]*(A + 1)/A - rhoc = mlbw.channel_radius[l](Ex) - rhochat = mlbw.scattering_radius[l](Ex) - P_c, S_c = penetration_shift(l, rhoc) - if Ex < 0: - P_c = 0 - - # Determine range of total angular momentum values based on equation - # 41 in LA-UR-12-27079 - jmin = abs(abs(I - l) - 0.5) - jmax = I + l + 0.5 - nJ = int(jmax - jmin + 1) - - # Determine Dl factor using Equation 43 in LA-UR-12-27079 - Dl = 2*l + 1 - g = malloc(nJ*sizeof(double)) - for ij in range(nJ): - j = jmin + ij - g[ij] = (2*j + 1)/(4*I + 2) - Dl -= g[ij] - - s = calloc(2*nJ, sizeof(double)) - for i_res in range(params.shape[0]): - # Copy resonance parameters - E_r = params[i_res, 0] - j = params[i_res, 2] - ij = int(j - jmin) - gt = params[i_res, 3] - gn = params[i_res, 4] - gg = params[i_res, 5] - gf = params[i_res, 6] - gx = params[i_res, 7] - P_r = params[i_res, 8] - S_r = params[i_res, 9] - P_rx = params[i_res, 10] - - # Calculate neutron and total width at energy E - gnE = P*gn/P_r # ENDF-102, Equation D.7 - gtE = gnE + gg + gf - if gx > 0: - gtE += gx*P_c/P_rx - - Eprime = E_r + (S_r - S)/(2*P_r)*gn # ENDF-102, Equation D.9 - x = 2*(E - Eprime)/gtE # LA-UR-12-27079, Equation 26 - f = 2*gnE/(gtE*(1 + x*x)) # Common factor in Equation 40 - s[ij][0] += f # First sum in Equation 40 - s[ij][1] += f*x # Second sum in Equation 40 - capture += f*g[ij]*gg/gtE - if gf > 0: - fission += f*g[ij]*gf/gtE - - for ij in range(nJ): - # Add all but last term of LA-UR-12-27079, Equation 40 - elastic += g[ij]*((1 - cos2phi - s[ij][0])**2 + - (sin2phi + s[ij][1])**2) - - # Add final term with Dl from Equation 40 - elastic += 2*Dl*(1 - cos2phi) - - # Free memory - free(g) - free(s) - - capture *= 2*M_PI/(k*k) - fission *= 2*M_PI/(k*k) - elastic *= M_PI/(k*k) - - return (elastic, capture, fission) - - -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -def reconstruct_slbw(slbw, double E): - """Evaluate cross section using SLBW data. - - Parameters - ---------- - slbw : openmc.data.SingleLevelBreitWigner - Single-level Breit-Wigner resonance parameters - E : double - Energy in eV at which to evaluate the cross section - - Returns - ------- - elastic : double - Elastic scattering cross section in barns - capture : double - Radiative capture cross section in barns - fission : double - Fission cross section in barns - - """ - cdef int i, l, i_res - cdef double elastic, capture, fission - cdef double A, k, rho, rhohat, I - cdef double P, S, phi, cos2phi, sin2phi, sinphi2 - cdef double Ex, rhoc, rhochat, P_c, S_c - cdef double E_r, J, gt, gn, gg, gf, gx, P_r, S_r, P_rx - cdef double gnE, gtE, Eprime, f - cdef double x, theta, psi, chi - cdef double [:,:] params - - I = slbw.target_spin - A = slbw.atomic_weight_ratio - k = _wave_number(A, E) - - elastic = 0. - capture = 0. - fission = 0. - - for i, l in enumerate(slbw._l_values): - params = slbw._parameter_matrix[l] - - rho = k*slbw.channel_radius[l](E) - rhohat = k*slbw.scattering_radius[l](E) - P, S = penetration_shift(l, rho) - phi = phaseshift(l, rhohat) - cos2phi = cos(2*phi) - sin2phi = sin(2*phi) - sinphi2 = sin(phi)**2 - - # Add potential scattering -- first term in ENDF-102, Equation D.2 - elastic += 4*M_PI/(k*k)*(2*l + 1)*sinphi2 - - # Determine shift and penetration at modified energy - if slbw._competitive[i]: - Ex = E + slbw.q_value[l]*(A + 1)/A - rhoc = k*slbw.channel_radius[l](Ex) - rhochat = k*slbw.scattering_radius[l](Ex) - P_c, S_c = penetration_shift(l, rhoc) - if Ex < 0: - P_c = 0 - - for i_res in range(params.shape[0]): - # Copy resonance parameters - E_r = params[i_res, 0] - J = params[i_res, 2] - gt = params[i_res, 3] - gn = params[i_res, 4] - gg = params[i_res, 5] - gf = params[i_res, 6] - gx = params[i_res, 7] - P_r = params[i_res, 8] - S_r = params[i_res, 9] - P_rx = params[i_res, 10] - - # Calculate neutron and total width at energy E - gnE = P*gn/P_r # Equation D.7 - gtE = gnE + gg + gf - if gx > 0: - gtE += gx*P_c/P_rx - - Eprime = E_r + (S_r - S)/(2*P_r)*gn # Equation D.9 - gJ = (2*J + 1)/(4*I + 2) # Mentioned in section D.1.1.4 - - # Calculate common factor for elastic, capture, and fission - # cross sections - f = M_PI/(k*k)*gJ*gnE/((E - Eprime)**2 + gtE**2/4) - - # Add contribution to elastic per Equation D.2 - elastic += f*(gnE*cos2phi - 2*(gg + gf)*sinphi2 - + 2*(E - Eprime)*sin2phi) - - # Add contribution to capture per Equation D.3 - capture += f*gg - - # Add contribution to fission per Equation D.6 - if gf > 0: - fission += f*gf - - return (elastic, capture, fission) - - -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -def reconstruct_rm(rm, double E): - """Evaluate cross section using Reich-Moore data. - - Parameters - ---------- - rm : openmc.data.ReichMoore - Reich-Moore resonance parameters - E : double - Energy in eV at which to evaluate the cross section - - Returns - ------- - elastic : double - Elastic scattering cross section in barns - capture : double - Radiative capture cross section in barns - fission : double - Fission cross section in barns - - """ - cdef int i, l, m, n, i_res - cdef int i_s, num_s, i_J, num_J - cdef double elastic, capture, fission, total - cdef double A, k, rho, rhohat, I - cdef double P, S, phi - cdef double smin, smax, s, Jmin, Jmax, J, j - cdef double E_r, gn, gg, gfa, gfb, P_r - cdef double E_diff, abs_value, gJ - cdef double Kr, Ki, x - cdef double complex Ubar, U_, factor - cdef bint hasfission - cdef np.ndarray[double, ndim=2] one - cdef np.ndarray[double complex, ndim=2] K, Imat, U - cdef double [:,:] params - - # Get nuclear spin - I = rm.target_spin - - elastic = 0. - fission = 0. - total = 0. - A = rm.atomic_weight_ratio - k = _wave_number(A, E) - one = np.eye(3) - K = np.zeros((3,3), dtype=complex) - - for i, l in enumerate(rm._l_values): - # Check for l-dependent scattering radius - rho = k*rm.channel_radius[l](E) - rhohat = k*rm.scattering_radius[l](E) - - # Calculate shift and penetrability - P, S = penetration_shift(l, rho) - - # Calculate phase shift - phi = phaseshift(l, rhohat) - - # Calculate common factor on collision matrix terms (term outside curly - # braces in ENDF-102, Eq. D.27) - Ubar = cexp(-2j*phi) - - # The channel spin is the vector sum of the target spin, I, and the - # neutron spin, 1/2, so can take on values of |I - 1/2| < s < I + 1/2 - smin = abs(I - 0.5) - smax = I + 0.5 - num_s = int(smax - smin + 1) - - for i_s in range(num_s): - s = i_s + smin - - # Total angular momentum is the vector sum of l and s and can assume - # values between |l - s| < J < l + s - Jmin = abs(l - s) - Jmax = l + s - num_J = int(Jmax - Jmin + 1) - - for i_J in range(num_J): - J = i_J + Jmin - - # Initialize K matrix - for m in range(3): - for n in range(3): - K[m,n] = 0.0 - - hasfission = False - if (l, J) in rm._parameter_matrix: - params = rm._parameter_matrix[l, J] - - for i_res in range(params.shape[0]): - # Sometimes, the same (l, J) quantum numbers can occur - # for different values of the channel spin, s. In this - # case, the sign of the channel spin indicates which - # spin is to be used. If the spin is negative assume - # this resonance comes from the I - 1/2 channel and vice - # versa. - j = params[i_res, 2] - if l > 0: - if (j < 0 and s != smin) or (j > 0 and s != smax): - continue - - # Copy resonance parameters - E_r = params[i_res, 0] - gn = params[i_res, 3] - gg = params[i_res, 4] - gfa = params[i_res, 5] - gfb = params[i_res, 6] - P_r = params[i_res, 7] - - # Calculate neutron width at energy E - gn = sqrt(P*gn/P_r) - - # Calculate j/2 * inverse of denominator of K matrix terms - factor = 0.5j/(E_r - E - 0.5j*gg) - - # Upper triangular portion of K matrix -- see ENDF-102, - # Equation D.28 - K[0,0] = K[0,0] + gn*gn*factor - if gfa != 0.0 or gfb != 0.0: - # Negate fission widths if necessary - gfa = (-1 if gfa < 0 else 1)*sqrt(abs(gfa)) - gfb = (-1 if gfb < 0 else 1)*sqrt(abs(gfb)) - - K[0,1] = K[0,1] + gn*gfa*factor - K[0,2] = K[0,2] + gn*gfb*factor - K[1,1] = K[1,1] + gfa*gfa*factor - K[1,2] = K[1,2] + gfa*gfb*factor - K[2,2] = K[2,2] + gfb*gfb*factor - hasfission = True - - # Get collision matrix - gJ = (2*J + 1)/(4*I + 2) - if hasfission: - # Copy upper triangular portion of K to lower triangular - K[1,0] = K[0,1] - K[2,0] = K[0,2] - K[2,1] = K[1,2] - - Imat = inv(one - K) - U = Ubar*(2*Imat - one) # ENDF-102, Eq. D.27 - elastic += gJ*cabs(1 - U[0,0])**2 # ENDF-102, Eq. D.24 - total += 2*gJ*(1 - creal(U[0,0])) # ENDF-102, Eq. D.23 - - # Calculate fission from ENDF-102, Eq. D.26 - fission += 4*gJ*(cabs(Imat[1,0])**2 + cabs(Imat[2,0])**2) - else: - U_ = Ubar*(2/(1 - K[0,0]) - 1) - if abs(creal(K[0,0])) < 3e-4 and abs(phi) < 3e-4: - # If K and phi are both very small, the calculated cross - # sections can lose precision because the real part of U - # ends up very close to unity. To get around this, we - # use Euler's formula to express Ubar by real and - # imaginary parts, expand cos(2phi) = 1 - 2phi^2 + - # O(phi^4), and then simplify - Kr = creal(K[0,0]) - Ki = cimag(K[0,0]) - x = 2*(-Kr + (Kr*Kr + Ki*Ki)*(1 - phi*phi) + phi*phi - - sin(2*phi)*Ki)/((1 - Kr)*(1 - Kr) + Ki*Ki) - total += 2*gJ*x - elastic += gJ*(x*x + cimag(U_)**2) - else: - total += 2*gJ*(1 - creal(U_)) # ENDF-102, Eq. D.23 - elastic += gJ*cabs(1 - U_)**2 # ENDF-102, Eq. D.24 - - # Calculate capture as difference of other cross sections as per ENDF-102, - # Equation D.25 - capture = total - elastic - fission - - elastic *= M_PI/(k*k) - capture *= M_PI/(k*k) - fission *= M_PI/(k*k) - - return (elastic, capture, fission) diff --git a/pyproject.toml b/pyproject.toml index 98b1d152fb..39aa261c1b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,6 @@ [build-system] -requires = ["setuptools", "wheel", "numpy", "cython"] +requires = ["setuptools", "wheel"] +build-backend = "setuptools.build_meta" [project] name = "openmc" diff --git a/setup.py b/setup.py deleted file mode 100755 index 88a45a3609..0000000000 --- a/setup.py +++ /dev/null @@ -1,14 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -from setuptools import setup -from Cython.Build import cythonize - - -kwargs = { - # Cython is used to add resonance reconstruction - 'ext_modules': cythonize('openmc/data/*.pyx'), - 'include_dirs': [np.get_include()] -} - -setup(**kwargs) diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index c0d6a1f154..db6ae1eb85 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -282,10 +282,6 @@ def test_slbw(xe135): s = resolved.parameters.iloc[0] assert s['energy'] == pytest.approx(0.084) - xs = resolved.reconstruct([10., 30., 100.]) - assert sorted(xs.keys()) == [2, 18, 102] - assert np.all(xs[18] == 0.0) - def test_mlbw(sm150): resolved = sm150.resonances.resolved @@ -294,10 +290,6 @@ def test_mlbw(sm150): assert resolved.energy_max == pytest.approx(1570.) assert resolved.target_spin == 0.0 - xs = resolved.reconstruct([10., 100., 1000.]) - assert sorted(xs.keys()) == [2, 18, 102] - assert np.all(xs[18] == 0.0) - def test_reichmoore(gd154): res = gd154.resonances @@ -319,7 +311,6 @@ def test_reichmoore(gd154): elastic = gd154.reactions[2].xs['0K'] assert isinstance(elastic, openmc.data.ResonancesWithBackground) - assert elastic(0.0253) == pytest.approx(5.7228949796394524) def test_rml(cl35): @@ -347,8 +338,6 @@ def test_mlbw_cov_lcomp0(cf252): assert not subset.parameters.empty assert (subset.file2res.parameters['energy'] < 100).all() samples = cov.sample(1) - xs = samples[0].reconstruct([10., 100., 1000.]) - assert sorted(xs.keys()) == [2, 18, 102] def test_mlbw_cov_lcomp1(ti50): @@ -365,9 +354,7 @@ def test_mlbw_cov_lcomp1(ti50): subset = cov.subset('L', [1, 1]) assert not subset.parameters.empty assert (subset.file2res.parameters['L'] == 1).all() - samples = cov.sample(1) - xs = samples[0].reconstruct([10., 100., 1000.]) - assert sorted(xs.keys()) == [2, 18, 102] + cov.sample(1) def test_mlbw_cov_lcomp2(na23): @@ -384,9 +371,7 @@ def test_mlbw_cov_lcomp2(na23): subset = cov.subset('L', [1, 1]) assert not subset.parameters.empty assert (subset.file2res.parameters['L'] == 1).all() - samples = cov.sample(1) - xs = samples[0].reconstruct([10., 100., 1000.]) - assert sorted(xs.keys()) == [2, 18, 102] + cov.sample(1) def test_rmcov_lcomp1(gd154): @@ -403,9 +388,7 @@ def test_rmcov_lcomp1(gd154): subset = cov.subset('energy', [0, 100]) assert not subset.parameters.empty assert (subset.file2res.parameters['energy'] < 100).all() - samples = cov.sample(1) - xs = samples[0].reconstruct([10., 100., 1000.]) - assert sorted(xs.keys()) == [2, 18, 102] + cov.sample(1) def test_rmcov_lcomp2(th232): @@ -422,9 +405,7 @@ def test_rmcov_lcomp2(th232): subset = cov.subset('energy', [0, 100]) assert not subset.parameters.empty assert (subset.file2res.parameters['energy'] < 100).all() - samples = cov.sample(1) - xs = samples[0].reconstruct([10., 100., 1000.]) - assert sorted(xs.keys()) == [2, 18, 102] + cov.sample(1) def test_madland_nix(am241): diff --git a/tools/ci/gha-install.sh b/tools/ci/gha-install.sh index 87952fda9c..cff7dc834f 100755 --- a/tools/ci/gha-install.sh +++ b/tools/ci/gha-install.sh @@ -40,8 +40,7 @@ if [[ $MPI == 'y' ]]; then export CC=mpicc export HDF5_MPI=ON export HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/mpich - pip install wheel "cython<3.0" - pip install --no-binary=h5py --no-build-isolation h5py + pip install --no-binary=h5py h5py fi # Build and install OpenMC executable