Making preprocessing more efficient by reducing parsing instances

tools/ALARAJOYWrapper/README.md

@@ -18,7 +18,11 @@ This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Dat
   * [Warnings](https://docs.python.org/3/library/warnings.html)
 - Generic Python packages
   * [Matplotlib.pyplot](https://matplotlib.org/3.5.3/api/_as_gen/matplotlib.pyplot.html)
-  * [NumPy](https://numpy.org/install/)
+  * [NumPy](https://numpy.org/install/) ([numpy.f2py](https://numpy.org/doc/stable/f2py/) is used for Fortran interfacing, which depends on the following):
+    * C-compiler (i.e. [gcc](https://hprc.tamu.edu/kb/Software/GNU-Compiler-Collection/#gcc-versions))
+    * Fortran-compiler (i.e. [gfortran](https://fortran-lang.org/learn/os_setup/install_gfortran/))
+    * [Meson](https://mesonbuild.com/)
+    * [Ninja](https://github.com/ninja-build/ninja.git)
   * [Pandas](https://pandas.pydata.org/docs/getting_started/install.html)
   * [PyYAML](https://pyyaml.org/wiki/PyYAMLDocumentation) (only needed if running `xs_plotting.py` as an independent script.)
 - Domain-specific packages

tools/ALARAJOYWrapper/njoy_endf_wrapper.f90

@@ -0,0 +1,59 @@
+! ===========================
+! This lightweight NJOY wrapper accesses the NJOY ENDF module's utility
+! function endf.terpa() to make use of NJOY's built-in interpolation 
+! functionality to apply the appropriate interpolation schemes. The
+! designation for these schemes are encoded within the TAB1 record header 
+! variable INT, representing one of the following interpolation schemes, as
+! laid out in the ENDF-6 manual
+! (https://www.nndc.bnl.gov/endfdocs/ENDF-102-2023.pdf), under section 0.5.2
+! ("Interpolation Laws"):
+!
+!     INT   | Interpolation Scheme
+!     1     | y is constant in x (constant, histogram)
+!     2     | y is linear in x (linear-linear)
+!     3     | y is linear in ln(x) (linear-log)
+!     4     | ln(y) is linear in x (log-linear)
+!     5     | ln(y) is linear in ln(x) (log-log)
+!     6     | special one-dimensional interpolation law, used for charged-
+!           | particle cross sections only
+!     11-25 | method of corresponding points (follow interpolation laws 1-5)
+!     21-25 | unit base interpolation (follow interpolation laws of 1-5)
+!
+! This module's incorporation and usage within ALARAJOYWrapper is managed by
+! njoy_tools.import_njoy_endf_wrapper(), which conditionally compiles this
+! Fortran file to an executable using numpy.f2py (if such an executable has
+! not already been created) and importing njoy_endf_wrapper as a Python
+! package. This allows the subroutine interpolate_tab1() to be callable within
+! ALARAJOYWrapper, which is necessary for the construction of pathway-specific
+! reaction cross-section from MF9 multiplicities multiplied by MF3 cumulative
+! cross-sections (see xs_plotting.extract_continuous_data() for specific use-
+! case implementation).
+! ===========================
+
+module njoy_endf_wrapper
+
+   use endf
+   implicit none
+
+contains
+
+   subroutine interpolate_tab1(tab1, x, y)
+
+      real(kind=8), intent(in) :: tab1(:)
+      real(kind=8), intent(in) :: x(:)
+      real(kind=8), intent(out) :: y(size(x))
+
+      integer :: i
+      integer :: ip, ir, idis
+      real(kind=8) :: xnext
+
+      ip = 2
+      ir = 1
+
+      do i = 1, size(x)
+         call terpa(y(i), x(i), xnext, idis, tab1, ip, ir)
+      end do
+
+   end subroutine interpolate_tab1
+
+end module njoy_endf_wrapper

tools/ALARAJOYWrapper/njoy_tools.py

@@ -4,6 +4,8 @@
 from pathlib import Path
 import re
 import numpy as np
+from sys import executable
+from shutil import which
 
 def set_directory():
     '''
@@ -699,4 +701,49 @@ def cleanup_njoy_files(element, A):
 
     output_dir = dir / 'njoy_outputs'
     output_dir.mkdir(exist_ok=True)
-    Path('output').rename(output_dir / f'njoy_output_{element}{A}.out')
+    Path('output').rename(output_dir / f'njoy_output_{element}{A}.out')
+
+def import_njoy_endf_wrapper():
+    """
+    Import the njoy_endf_wrapper module defined by the targeted NJOY-wrapper
+        njoy_endf_wrapper.f90. This wrapper utilizes the NJOY function
+        `endf.terpa()`, which interpolates TAB1 data according to the encoded
+        interpolation scheme(s). If the module cannot be accessed, compile
+        njoy_endf_wrapper.f90 to a CPython executable using NumPy.f2py and
+        subsequently import the module.
+
+    Arguments:
+        None
+
+    Returns:
+        njoy_endf_wrapper.njoy_endf_wrapper (fortran object): Python module of
+            the compiled njoy_endf_wrapper.f90 NJOY wrapper containing the
+            subroutine `interpolate_tab1()`, which can be used to apply
+            `endf.terpa()` to interpolate TAB1 according to the encoded
+            interpolation scheme(s).
+    """
+
+    try:
+        import njoy_endf_wrapper
+
+    except (ModuleNotFoundError, ImportError):
+        njoy_dir = Path(which('njoy')).parent
+        subprocess.run(
+            [
+                executable,
+                '-m',
+                'numpy.f2py',
+                '-c',
+                '-m',
+                'njoy_endf_wrapper',
+                str(Path(__file__).parent / 'njoy_endf_wrapper.f90'),
+                f'-I{njoy_dir / "fortran_modules"}',
+                f'-L{njoy_dir}',
+                '-lnjoy'
+            ],
+            check=True,
+            cwd=Path.cwd()
+        )
+        import njoy_endf_wrapper
+
+    return njoy_endf_wrapper.njoy_endf_wrapper

tools/ALARAJOYWrapper/preprocess_fendl3.py

@@ -115,7 +115,7 @@ def configure_logging(redirect_warnings=False):
 
 def process_pendf(
     material_id, MTs, pKZA, mt_dict, temperature,
-    tendl_path, tendl_dir, unresr_err_cases
+    endf_dict, tendl_dir, unresr_err_cases
 ):
     """
     Prepare and run initial NJOY run with MODER, RECONR, BROADR, UNRESR, and
@@ -133,8 +133,8 @@ def process_pendf(
         MTs (set): Set of all MTs from the original TENDL file.
         mt_dict (dict): Dictionary formatted data structure for mt_table.csv.
         temperature (float): Temperature at which to run NJOY modules.
-        tendl_path (pathlib._local.PosixPath): Path to the original,
-            unmodified TENDL file.
+        endf_dict (dict): Nested EndfParserPy-formatted dictionary containing
+            all parsed nuclear data from a TENDL file.
         tendl_dir (pathlib._local.PosixPath): Path to the directory in which
             the original TENDL nuclide files are contained.
         unresr_err_cases (list of str): List of all nuclides that required an
@@ -216,11 +216,9 @@ def process_pendf(
     else:
         njoy_error += prep_error
 
-    pendf_MTs = set(tp.parse_endf_file_level_data(
-        pendf_path, endf_format='pendf'
-    )[0])
+    pendf_MTs = set(tp.parse_endf_data(pendf_path, endf_format='pendf')[1])
     MTs |= pendf_MTs.intersection(set(rxd.GAS_DF['total_mt']))
-    isomer_dict = tp.determine_all_excitations(tendl_path, MTs)
+    isomer_dict = tp.determine_all_excitations(endf_dict, MTs)
 
     return MTs, isomer_dict, njoy_error, unresr_err_cases
 
@@ -429,7 +427,7 @@ def rxn_to_str(parent, daughter, MT, rxn):
     return dsv_row + ' '.join(str(xs) for xs in rxn['xsections'])
 
 def store_results(
-    dsv_path, all_rxns, nGroups, tendl_dir, group_name, plotting
+    dsv_path, all_rxns, nGroups, tendl_dir, group_name, plotting, endf_dict
 ):
     """
     Save groupwise-converted cross-section data to a space-delimited DSV file
@@ -473,7 +471,9 @@ def store_results(
             cross-sections.
         plotting (bool): Boolean to set whether to produce cross-section
             plots.
-
+        endf_dict (dict): Nested EndfParserPy-formatted dictionary containing
+            all parsed nuclear data from a TENDL file.
+
     Returns:
         None 
     """
@@ -501,8 +501,7 @@ def store_results(
                                 )
 
                                 continuous_dict = xp.extract_continuous_data(
-                                    tendl_dir / f'{element}{A}.tendl',
-                                    xp.flagged_num_to_int(MT)
+                                    endf_dict, MT
                                 )
 
                                 energies = njt.load_external_group_struct(
@@ -521,7 +520,7 @@ def store_results(
                                 plot_path = xp.set_plot_save_path(
                                     element, A, emitted, tendl_dir, group_name
                                 )
-                                
+
                                 plt.savefig(plot_path)
 
         if plotting:
@@ -581,8 +580,8 @@ def main():
     for file_properties in tp.search_for_files(search_dir):
         element, A, pKZA, endf_path = tuple(file_properties.values())
         TAPE20.write_bytes(endf_path.read_bytes())
-        endf_file_dict, material_id = tp.parse_endf_file_level_data(TAPE20)
-        MTs = set(endf_file_dict)
+        endf_dict, mf3_file_dict, material_id = tp.parse_endf_data(TAPE20)
+        MTs = set(mf3_file_dict)
 
         if len((MTs - rxd.SPEC_MTS) - endf6_MTs) > 0:
             invalid_MTs = sorted((MTs - rxd.SPEC_MTS) - endf6_MTs)
@@ -594,7 +593,7 @@ def main():
 
         MTs, isomer_dict, njoy_prep_error, unresr_err_cases = process_pendf(
             material_id, MTs, pKZA, mt_dict, temperature,
-            TAPE20, search_dir, unresr_err_cases
+            endf_dict, search_dir, unresr_err_cases
         )
 
         if not njoy_prep_error:
@@ -628,7 +627,7 @@ def main():
     dsv_path = dir / 'cumulative_gendf_data.dsv'
     store_results(
         dsv_path, gas_filtered, nGroups,
-        search_dir, group_name, args.xs_plotting 
+        search_dir, group_name, args.xs_plotting, endf_dict
     )
     print(
         f'Neutron activation cross-sections converted to {nGroups} groups ' \

tools/ALARAJOYWrapper/tendl_processing.py

@@ -23,8 +23,9 @@
     for val in arr
 }
 ISOMERIC_STATES = 'mnopqrstuvwxyz'
+PATH_SPECIFIC_MFS = (9,10)
 
-def parse_endf_file_level_data(endf_path, MF=3, endf_format='endf6-ext'):
+def parse_endf_data(endf_path, MF=3, endf_format='endf6-ext'):
     """
     For an ENDF-formatted TENDL file containing neutron activation data for a
         single nuclide, parse and store the file's data into a nested
@@ -52,7 +53,29 @@ def parse_endf_file_level_data(endf_path, MF=3, endf_format='endf6-ext'):
 
     endf_dict = EndfParserPy(endf_format=endf_format).parsefile(endf_path)
 
-    return endf_dict.get(MF, {}), endf_dict[1][451]['MAT']
+    return endf_dict, endf_dict.get(MF, {}), endf_dict[1][451]['MAT']
+
+def get_section_dict(endf_dict, MF, MT):
+    """
+    Produce a reaction (MT)-specific subdictionary from a nested EndfParserPy-
+        formatted nested dictionary containing a whole TENDL file's parsed
+        nuclear data. Will return an empty dictionary if either the MF or MT
+        are not present in the provided dictionary.
+
+    Arguments:
+        endf_dict (dict): Nested EndfParserPy-formatted dictionary containing
+            all parsed nuclear data from a TENDL file.
+        MF (int): ENDF file number.
+        MT (int): Unique reaction identifier.
+
+    Returns:
+        section (dict): Sub-dictionary containing nuclear data for a
+            given MF/MT combination from a parsed TENDL file. Will return an
+            empty dictionary if either the MF or MT is not present in
+            `endf_dict`.
+    """
+
+    return endf_dict.get(MF, {}).get(MT, {})
 
 def calculate_KZA_from_ENDF(filepath, MF=1, MT=451):
     """
@@ -76,7 +99,7 @@ def calculate_KZA_from_ENDF(filepath, MF=1, MT=451):
         KZA (int): Unique ZZZAAAM identifier for a given nuclide.
     """
 
-    nuc_data = parse_endf_file_level_data(filepath, MF)[0][MT]
+    nuc_data = parse_endf_data(filepath, MF)[1][MT]
 
     return int(nuc_data['ZA'] * 10 + nuc_data['LISO'])
 
@@ -161,7 +184,7 @@ def search_for_files(dir = Path.cwd()):
 
     return file_info
 
-def determine_all_excitations(endf_path, MTs):
+def determine_all_excitations(endf_dict, MTs):
     """
     Reference an ENDF file's MF9 and MF10 file data and explicitly defined
         excitation reactions to construct a nested dictionary keyed by
@@ -172,8 +195,8 @@ def determine_all_excitations(endf_path, MTs):
         cross-section data.
 
     Arguments:
-        endf_path (pathlib._local.PosixPath): Path to the ENDF (TENDL) file to
-            be processed.
+        endf_dict (dict): Nested EndfParserPy-formatted dictionary containing
+            all parsed nuclear data from a TENDL file.
         MTs (set): Set of all MT reaction numbers contained in the TENDL file.
 
     Returns:
@@ -187,20 +210,16 @@ def determine_all_excitations(endf_path, MTs):
 
     isomer_dict = defaultdict(lambda: defaultdict(list))
 
-    path_specific_MFs = (9,10)
-    mf_dict = {
-        MF: parse_endf_file_level_data(endf_path, MF)[0]
-        for MF in path_specific_MFs
-    }
-
     for MT in MTs:
         cumulative_MT = REVERSE_EXCITATION_DICT.get(MT)
         if MT not in EXCITATION_REACTIONS:
             # Isomer pathways contained either in MF 9 ("Multiplicities for
             # Production of Radioactive Nuclides") and MF 10 ("Cross Sections
             # for Production of Radioactive Nuclides").
-            for MF in path_specific_MFs:
-                pathways = mf_dict[MF].get(MT, {}).get('subsection', {})
+            for MF in PATH_SPECIFIC_MFS:
+                pathways = get_section_dict(
+                    endf_dict, MF, MT
+                ).get('subsection', {})
                 for pathway_data in pathways.values():
                     isomer_dict[MT][MF].append(pathway_data['LFS'])