Merge branch 'openmc-dev:develop' into develop

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} \

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

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 <https://cython.org/>`_
-      Cython is used for resonance reconstruction for ENDF data converted to
-      :class:`openmc.data.IncidentNeutron`.
-
    `vtk <https://vtk.org/>`_
       The Python VTK bindings are needed to convert voxel and track files to VTK
       format.

openmc/data/endf.py

@@ -14,11 +14,7 @@
 
 from .data import gnds_name
 from .function import Tabulated1D
-try:
-    from ._endf import float_endf
-    _CYTHON = True
-except ImportError:
-    _CYTHON = False
+from endf.records import float_endf
 
 
 _LIBRARY = {0: 'ENDF/B', 1: 'ENDF/A', 2: 'JEFF', 3: 'EFF',
@@ -91,10 +87,6 @@ def py_float_endf(s):
     return float(ENDF_FLOAT_RE.sub(r'\1e\2\3', s))
 
 
-if not _CYTHON:
-    float_endf = py_float_endf
-
-
 def int_endf(s):
     """Convert string of integer number in ENDF to int.
 

openmc/data/function.py

@@ -708,28 +708,6 @@ def __init__(self, resonances, background, mt):
         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

openmc/data/neutron.py

@@ -16,8 +16,7 @@
     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 @@ def add_temperature_from_ace(self, ace_or_filename, metastable_scheme='nndc'):
         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 @@ def add_elastic_0K_from_endf(self, filename, overwrite=False):
             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 @@ def add_elastic_0K_from_endf(self, filename, overwrite=False):
         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
-
-                # 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)
+            # Run NJOY to create a pointwise ENDF file
+            make_pendf(filename, **kwargs)
 
-        # 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.

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,20 +230,19 @@ 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
     ------
     subprocess.CalledProcessError
         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,