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