diff --git a/tools/ALARAJOYWrapper/README.md b/tools/ALARAJOYWrapper/README.md index ca2c080c..e81e042e 100644 --- a/tools/ALARAJOYWrapper/README.md +++ b/tools/ALARAJOYWrapper/README.md @@ -24,7 +24,7 @@ This preprocessor uses [NJOY 2016](https://github.com/njoy/NJOY2016) Nuclear Dat - Domain-specific packages * [Endf-parserpy](https://github.com/IAEA-NDS/endf-parserpy) * [NJOY 2016](https://github.com/njoy/NJOY2016) - * [OpenMC](https://docs.openmc.org/en/stable/quickinstall.html) (only needed if specifying a multigroup energy structure by name from the dictionary `openmc.mgxs.GROUP_STRUCTURES`) + * [OpenMC](https://docs.openmc.org/en/stable/quickinstall.html) (needed if specifying a multigroup energy structure by name from the dictionary `openmc.mgxs.GROUP_STRUCTURES` or utilizing the `xs_plotting` module) diff --git a/tools/ALARAJOYWrapper/preprocess_fendl3.py b/tools/ALARAJOYWrapper/preprocess_fendl3.py index f71fe285..4146be88 100644 --- a/tools/ALARAJOYWrapper/preprocess_fendl3.py +++ b/tools/ALARAJOYWrapper/preprocess_fendl3.py @@ -11,6 +11,7 @@ from pathlib import Path from collections import defaultdict from subprocess import TimeoutExpired +from openmc.data import endf def make_argparser(): parser = argparse.ArgumentParser() @@ -482,6 +483,7 @@ def store_results( dsv.write(f'{nGroups} {group_name}\n') for parent in sorted(all_rxns): element, A = tp.interpret_KZA(parent) + endf_obj = endf.Evaluation(tendl_dir / f'{element}{A}.tendl') for daughter in all_rxns[parent]: if parent != daughter: for MT, rxn in all_rxns[parent][daughter].items(): @@ -501,8 +503,7 @@ def store_results( ) continuous_dict = xp.extract_continuous_data( - tendl_dir / f'{element}{A}.tendl', - xp.flagged_num_to_int(MT) + endf_obj, MT ) energies = njt.load_external_group_struct( @@ -521,7 +522,7 @@ def store_results( plot_path = xp.set_plot_save_path( element, A, emitted, tendl_dir, group_name ) - + plt.savefig(plot_path) if plotting: diff --git a/tools/ALARAJOYWrapper/xs_plotting.py b/tools/ALARAJOYWrapper/xs_plotting.py index b9c007ae..4f0cd01b 100644 --- a/tools/ALARAJOYWrapper/xs_plotting.py +++ b/tools/ALARAJOYWrapper/xs_plotting.py @@ -5,7 +5,9 @@ import njoy_tools as njt import reaction_data as rxd import matplotlib.pyplot as plt +import pandas as pd from pathlib import Path +from openmc.data import Reaction, endf def flagged_num_to_int(num): """ @@ -26,14 +28,15 @@ def flagged_num_to_int(num): instance of '*' contained in the original value. """ - re_match = re.match(r'^-?\d+', str(num)) + num = str(num) + re_match = re.match(r'^-?\d+', num) if not re_match: raise ValueError( f'Invalid flagged number {num}. Must be formatted with numeric ' \ 'characters before non-numeric characters.' ) - - return int(re_match.group()) + + return int(re_match.group()), num.count('*') def ensure_emission_specificity(emitted, dKZA): """ @@ -61,14 +64,13 @@ def ensure_emission_specificity(emitted, dKZA): return emitted -def extract_continuous_data(tendl_path, MT): +def extract_continuous_data(endf_obj, MT): """ For a given nuclide and reaction, extract its continuous-energy cross- sections and corresponding energies from its original TENDL file. Arguments: - tendl_path (pathlib._local.PosixPath): Path to the nuclide's original - TENDL file. + endf_obj (openmc.data.endf.Evaluation): OpenMC parsed-ENDF object. MT (int): Reaction identifying number. Returns: @@ -83,16 +85,39 @@ def extract_continuous_data(tendl_path, MT): lists. """ - xs_table = ( - tp.parse_endf_file_level_data(tendl_path)[0] - .get(MT, {}) - .get('xstable', {'E' : [], 'xs' : []}) - ) + continuous_dict = {'energies' : [], 'xs' : []} + MT, isomeric_state = flagged_num_to_int(MT) + rxn = Reaction.from_endf(endf_obj, MT) - return { - 'xs' : xs_table['xs'], - 'energies' : xs_table['E'] - } + # For excitation reactions, calculate specific pathway reactions by + # multiplying reaction multiplicities by MF3 cumulative cross-sections + # interpolated by the multiplicities' energy array + if isomeric_state > 0: + + pathways = [] + for product in rxn.products: + if product.particle not in {'neutron', 'photon', 'electron'}: + iso_flag = re.compile(r'_e(\d+)$').search(product.particle) + excited_state = int(iso_flag.group(1)) if iso_flag else 0 + pathways.append((excited_state, product)) + + pathways.sort(key=lambda pathway: pathway[0]) + + if pathways and isomeric_state < len(pathways): + product = pathways[isomeric_state][1] + energies = product.yield_.x + continuous_dict['energies'].extend(energies) + continuous_dict['xs'].extend( + product.yield_.y * rxn.xs['0K'](energies) + ) + + else: + mf3_xs_table = rxn.xs.get('0K') + if mf3_xs_table: + continuous_dict['energies'].extend(mf3_xs_table.x) + continuous_dict['xs'].extend(mf3_xs_table.y) + + return continuous_dict def extract_groupwise_data_from_DSV(dsv_list, KZA, MT): """ @@ -141,14 +166,14 @@ def extract_groupwise_data_from_DSV(dsv_list, KZA, MT): dsv_pKZA, dsv_dKZA, dsv_MT, emitted = rxn[:4] emitted = ensure_emission_specificity(emitted, dsv_dKZA) - if KZA == dsv_pKZA and MT == flagged_num_to_int(dsv_MT)[0]: + if KZA == dsv_pKZA and str(MT) == dsv_MT: groupwise_dict[group_name] = { 'xs' : np.array(rxn[4:]).astype(float), 'energies' : energy_bounds } break - return groupwise_dict, ensure_emission_specificity(emitted, dsv_dKZA) + return groupwise_dict, emitted def set_plot_save_path( element, A, emitted, tendl_dir, group_names, img_ext='png' @@ -348,7 +373,30 @@ def find_all_mass_nums(tendl_dir, element): mass_nums.add(nuc_match.group(1)) return mass_nums - + +def find_all_MTs(dsv_list, pKZA): + """ + Given a list of preprocessed groupwise DSV files and a parent nuclide + identified by its KZA, compile all reaction identifiers (MTs) that + exist for that nuclide in any of the DSVs. + + Arguments: + dsv_list (list): List of filepaths to DSV files containing + ALARAJOYWrapper-processed groupwise TENDL data. + pKZA (int): ZZZAAAM identifier of the parent nuclide. + + Returns: + MTs (set): Set of all reaction types for the parent nuclide present in + any of the DSV files provided. + """ + + MTs = set() + for dsv in dsv_list: + df = pd.read_csv(dsv, sep=r'\s+', skiprows=1, header=None) + MTs.update(df.loc[df[0] == pKZA, 2]) + + return MTs + def main(): # Only load in yaml module when executing xs_plotting.py as a script, @@ -357,6 +405,7 @@ def main(): from yaml import safe_load plt.rcParams.update({'figure.max_open_warning': 0}) + plot_path = None parser = argparse.ArgumentParser() parser.add_argument('--yaml', '-y') @@ -394,7 +443,7 @@ def main(): for A in mass_nums: KZA = str(( - njt.elements[element] * 1000 + flagged_num_to_int(A) + njt.elements[element] * 1000 + flagged_num_to_int(A)[0] ) * 10 + tp.ISOMERIC_STATES.find(str(A)[-1]) + 1) MTs = adjust_dict_for_all_tag(element_dict, A) @@ -403,17 +452,14 @@ def main(): MTs = [MTs] if check_all_tag(MTs): - MTs = rxd.process_mt_data(rxd.load_mt_table( - njt.set_directory() / 'mt_table.csv' - )).keys() + MTs = find_all_MTs(dsv_list, KZA) - for MT in [flagged_num_to_int(MT) for MT in MTs]: + for MT in MTs: fig, ax = plt.subplots(figsize=(10,6)) continuous_dict = extract_continuous_data( - tendl_dir / f'{element}{A}.tendl', flagged_num_to_int(MT) + endf.Evaluation(tendl_dir / f'{element}{A}.tendl'), MT ) - groupwise_dict, emitted = extract_groupwise_data_from_DSV( dsv_list, KZA, MT ) @@ -428,10 +474,11 @@ def main(): ) plt.savefig(plot_path) - print( - f'Cross-section plots saved to {plot_path.parents[2]}/, ' \ - 'organized by element, nuclide, reaction.' - ) + if plot_path: + print( + f'Cross-section plots saved to {plot_path.parents[2]}/, ' \ + 'organized by element, nuclide, reaction.' + ) if __name__ == '__main__':