|
| 1 | +#!/usr/bin/env python3 |
| 2 | + |
| 3 | +import glob |
| 4 | +import os |
| 5 | +from zipfile import ZipFile |
| 6 | +from tarfile import TarFile |
| 7 | +from collections import OrderedDict, defaultdict |
| 8 | +from io import StringIO |
| 9 | +from itertools import chain |
| 10 | + |
| 11 | +try: |
| 12 | + import lxml.etree as ET |
| 13 | + _have_lxml = True |
| 14 | +except ImportError: |
| 15 | + import xml.etree.ElementTree as ET |
| 16 | + _have_lxml = False |
| 17 | + |
| 18 | +import openmc.data |
| 19 | +import openmc.deplete |
| 20 | +from openmc._xml import clean_indentation |
| 21 | +from openmc.deplete.chain import REACTIONS, replace_missing_fpy |
| 22 | +from openmc.deplete.nuclide import Nuclide, DecayTuple, ReactionTuple, \ |
| 23 | + FissionYieldDistribution |
| 24 | + |
| 25 | +from casl_chain import CASL_CHAIN, UNMODIFIED_DECAY_BR |
| 26 | +from utils import download |
| 27 | + |
| 28 | + |
| 29 | +URLS = [ |
| 30 | + 'https://www.oecd-nea.org/dbdata/jeff/jeff33/downloads/JEFF33-n.tgz', |
| 31 | + 'https://www.oecd-nea.org/dbdata/jeff/jeff33/downloads/JEFF33-rdd.zip', |
| 32 | + 'https://www.oecd-nea.org/dbdata/jeff/jeff33/downloads/JEFF33-nfy.asc' |
| 33 | +] |
| 34 | + |
| 35 | + |
| 36 | +def replace_missing_decay_product(product, decay_data, all_decay_data): |
| 37 | + # Determine atomic number, mass number, and metastable state |
| 38 | + Z, A, state = openmc.data.zam(product) |
| 39 | + symbol = openmc.data.ATOMIC_SYMBOL[Z] |
| 40 | + |
| 41 | + # Iterate until we find an existing nuclide in the chain |
| 42 | + while product not in decay_data: |
| 43 | + # If product has no decay data in the library, nothing further can be done |
| 44 | + if product not in all_decay_data: |
| 45 | + product = None |
| 46 | + break |
| 47 | + |
| 48 | + # If the current product is not in the chain but is stable, there's |
| 49 | + # nothing further we can do. Also, we only want to continue down the |
| 50 | + # decay chain if the half-life is short, so we also make a cutoff here |
| 51 | + # to terminate if the half-life is more than 1 day. |
| 52 | + decay_obj = all_decay_data[product] |
| 53 | + if decay_obj.nuclide['stable'] or decay_obj.half_life.n > 24*60*60: |
| 54 | + product = None |
| 55 | + break |
| 56 | + |
| 57 | + dominant_mode = max(decay_obj.modes, key=lambda x: x.branching_ratio) |
| 58 | + product = dominant_mode.daughter |
| 59 | + |
| 60 | + return product |
| 61 | + |
| 62 | + |
| 63 | +def main(): |
| 64 | + if os.path.isdir('./decay') and os.path.isdir('./nfy') and os.path.isdir('./neutrons'): |
| 65 | + endf_dir = '.' |
| 66 | + elif 'OPENMC_ENDF_DATA' in os.environ: |
| 67 | + endf_dir = os.environ['OPENMC_ENDF_DATA'] |
| 68 | + else: |
| 69 | + for url in URLS: |
| 70 | + basename = download(url) |
| 71 | + if url.endswith('.tgz'): |
| 72 | + tar = TarFile.open(basename, 'r') |
| 73 | + print('Extracting {}...'.format(basename)) |
| 74 | + tar.extractall() |
| 75 | + os.system('mv endf6 neutrons') |
| 76 | + elif url.endswith('.zip'): |
| 77 | + with ZipFile(basename, 'r') as zf: |
| 78 | + print('Extracting {}...'.format(basename)) |
| 79 | + os.mkdir('decay') |
| 80 | + zf.extractall('decay') |
| 81 | + else: |
| 82 | + os.system('./nfy_to_endf.sh ' + str(basename)) |
| 83 | + endf_dir = '.' |
| 84 | + |
| 85 | + decay_files = glob.glob(os.path.join(endf_dir, 'decay', '*.ASC')) |
| 86 | + fpy_files = glob.glob(os.path.join(endf_dir, 'nfy', '*.endf')) |
| 87 | + neutron_files = glob.glob(os.path.join(endf_dir, 'neutrons', '*.jeff33')) |
| 88 | + |
| 89 | + # Create a Chain |
| 90 | + chain = openmc.deplete.Chain() |
| 91 | + |
| 92 | + print('Reading ENDF nuclear data from "{}"...'.format(os.path.abspath(endf_dir))) |
| 93 | + |
| 94 | + # Create dictionary mapping target to filename |
| 95 | + print('Processing neutron sub-library files...') |
| 96 | + reactions = {} |
| 97 | + for f in neutron_files: |
| 98 | + evaluation = openmc.data.endf.Evaluation(f) |
| 99 | + nuc_name = evaluation.gnds_name |
| 100 | + if nuc_name in CASL_CHAIN: |
| 101 | + reactions[nuc_name] = {} |
| 102 | + for mf, mt, nc, mod in evaluation.reaction_list: |
| 103 | + # Q value for each reaction is given in MF=3 |
| 104 | + if mf == 3: |
| 105 | + file_obj = StringIO(evaluation.section[3, mt]) |
| 106 | + openmc.data.endf.get_head_record(file_obj) |
| 107 | + q_value = openmc.data.endf.get_cont_record(file_obj)[1] |
| 108 | + reactions[nuc_name][mt] = q_value |
| 109 | + |
| 110 | + # Determine what decay and FPY nuclides are available |
| 111 | + print('Processing decay sub-library files...') |
| 112 | + decay_data = {} |
| 113 | + all_decay_data = {} |
| 114 | + for f in decay_files: |
| 115 | + decay_obj = openmc.data.Decay(f) |
| 116 | + nuc_name = decay_obj.nuclide['name'] |
| 117 | + all_decay_data[nuc_name] = decay_obj |
| 118 | + if nuc_name in CASL_CHAIN: |
| 119 | + decay_data[nuc_name] = decay_obj |
| 120 | + |
| 121 | + for nuc_name in CASL_CHAIN: |
| 122 | + if nuc_name not in decay_data: |
| 123 | + print('WARNING: {} has no decay data!'.format(nuc_name)) |
| 124 | + |
| 125 | + print('Processing fission product yield sub-library files...') |
| 126 | + fpy_data = {} |
| 127 | + for f in fpy_files: |
| 128 | + fpy_obj = openmc.data.FissionProductYields(f) |
| 129 | + name = fpy_obj.nuclide['name'] |
| 130 | + if name in CASL_CHAIN: |
| 131 | + fpy_data[name] = fpy_obj |
| 132 | + |
| 133 | + print('Creating depletion_chain...') |
| 134 | + missing_daughter = [] |
| 135 | + missing_rx_product = [] |
| 136 | + missing_fpy = [] |
| 137 | + |
| 138 | + for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)): |
| 139 | + data = decay_data[parent] |
| 140 | + |
| 141 | + nuclide = Nuclide(parent) |
| 142 | + |
| 143 | + chain.nuclides.append(nuclide) |
| 144 | + chain.nuclide_dict[parent] = idx |
| 145 | + |
| 146 | + if not CASL_CHAIN[parent][0] and \ |
| 147 | + not data.nuclide['stable'] and data.half_life.nominal_value != 0.0: |
| 148 | + nuclide.half_life = data.half_life.nominal_value |
| 149 | + nuclide.decay_energy = data.decay_energy.nominal_value |
| 150 | + nuclide.sources = data.sources |
| 151 | + sum_br = 0.0 |
| 152 | + for mode in data.modes: |
| 153 | + decay_type = ','.join(mode.modes) |
| 154 | + if mode.daughter in decay_data: |
| 155 | + target = mode.daughter |
| 156 | + else: |
| 157 | + missing_daughter.append((parent, mode)) |
| 158 | + continue |
| 159 | + |
| 160 | + # Append decay mode |
| 161 | + br = mode.branching_ratio.nominal_value |
| 162 | + nuclide.add_decay_mode(decay_type, target, br) |
| 163 | + |
| 164 | + # Ensure sum of branching ratios is unity by slightly modifying last |
| 165 | + # value if necessary |
| 166 | + sum_br = sum(m.branching_ratio for m in nuclide.decay_modes) |
| 167 | + if sum_br != 1.0 and nuclide.decay_modes and parent not in UNMODIFIED_DECAY_BR: |
| 168 | + decay_type, target, br = nuclide.decay_modes.pop() |
| 169 | + br = 1.0 - sum(m.branching_ratio for m in nuclide.decay_modes) |
| 170 | + nuclide.add_decay_mode(decay_type, target, br) |
| 171 | + |
| 172 | + # If nuclide has incident neutron data, we need to list what |
| 173 | + # transmutation reactions are possible |
| 174 | + fissionable = False |
| 175 | + transmutation_reactions = ('(n,2n)', '(n,3n)', '(n,4n)', '(n,gamma)', |
| 176 | + '(n,p)', '(n,a)') |
| 177 | + if parent in reactions: |
| 178 | + reactions_available = reactions[parent].keys() |
| 179 | + for name in transmutation_reactions: |
| 180 | + mts, changes, _ = REACTIONS[name] |
| 181 | + if mts & reactions_available: |
| 182 | + delta_A, delta_Z = changes |
| 183 | + A = data.nuclide['mass_number'] + delta_A |
| 184 | + Z = data.nuclide['atomic_number'] + delta_Z |
| 185 | + daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A) |
| 186 | + |
| 187 | + if daughter not in decay_data: |
| 188 | + daughter = replace_missing_decay_product( |
| 189 | + daughter, decay_data, all_decay_data) |
| 190 | + if daughter is None: |
| 191 | + missing_rx_product.append((parent, name, daughter)) |
| 192 | + |
| 193 | + # Store Q value -- use sorted order so we get summation |
| 194 | + # reactions (e.g., MT=103) first |
| 195 | + for mt in sorted(mts): |
| 196 | + if mt in reactions[parent]: |
| 197 | + q_value = reactions[parent][mt] |
| 198 | + break |
| 199 | + else: |
| 200 | + q_value = 0.0 |
| 201 | + |
| 202 | + nuclide.add_reaction(name, daughter, q_value, 1.0) |
| 203 | + |
| 204 | + # Check for fission reactions |
| 205 | + if any(mt in reactions_available for mt in [18, 19, 20, 21, 38]): |
| 206 | + q_value = reactions[parent][18] |
| 207 | + nuclide.add_reaction('fission', None, q_value, 1.0) |
| 208 | + fissionable = True |
| 209 | + |
| 210 | + if fissionable: |
| 211 | + if parent in fpy_data: |
| 212 | + fpy = fpy_data[parent] |
| 213 | + else: |
| 214 | + nuclide._fpy = replace_missing_fpy(parent, fpy_data, decay_data) |
| 215 | + missing_fpy.append((parent, nuclide._fpy)) |
| 216 | + continue |
| 217 | + |
| 218 | + if fpy.energies is not None: |
| 219 | + yield_energies = fpy.energies |
| 220 | + else: |
| 221 | + yield_energies = [0.0] |
| 222 | + |
| 223 | + yield_data = {} |
| 224 | + for E, table_yd, table_yc in zip(yield_energies, fpy.independent, fpy.cumulative): |
| 225 | + yields = defaultdict(float) |
| 226 | + for product in table_yd: |
| 227 | + if product in decay_data: |
| 228 | + # identifier |
| 229 | + ifpy = CASL_CHAIN[product][2] |
| 230 | + # 1 for independent |
| 231 | + if ifpy == 1: |
| 232 | + if product not in table_yd: |
| 233 | + print('No independent fission yields found for {} in {}'.format(product, parent)) |
| 234 | + else: |
| 235 | + yields[product] += table_yd[product].nominal_value |
| 236 | + # 2 for cumulative |
| 237 | + elif ifpy == 2: |
| 238 | + if product not in table_yc: |
| 239 | + print('No cumulative fission yields found for {} in {}'.format(product, parent)) |
| 240 | + else: |
| 241 | + yields[product] += table_yc[product].nominal_value |
| 242 | + # 3 for special treatment with weight fractions |
| 243 | + elif ifpy == 3: |
| 244 | + for name_i, weight_i, ifpy_i in CASL_CHAIN[product][3]: |
| 245 | + if name_i not in table_yd: |
| 246 | + print('No fission yields found for {} in {}'.format(name_i, parent)) |
| 247 | + else: |
| 248 | + if ifpy_i == 1: |
| 249 | + yields[product] += weight_i * table_yd[name_i].nominal_value |
| 250 | + elif ifpy_i == 2: |
| 251 | + yields[product] += weight_i * table_yc[name_i].nominal_value |
| 252 | + |
| 253 | + yield_data[E] = yields |
| 254 | + |
| 255 | + nuclide.yield_data = FissionYieldDistribution(yield_data) |
| 256 | + |
| 257 | + # Replace missing FPY data |
| 258 | + for nuclide in chain.nuclides: |
| 259 | + if hasattr(nuclide, '_fpy'): |
| 260 | + nuclide.yield_data = chain[nuclide._fpy].yield_data |
| 261 | + |
| 262 | + # Display warnings |
| 263 | + if missing_daughter: |
| 264 | + print('The following decay modes have daughters with no decay data:') |
| 265 | + for parent, mode in missing_daughter: |
| 266 | + print(' {} -> {} ({})'.format(parent, mode.daughter, ','.join(mode.modes))) |
| 267 | + print('') |
| 268 | + |
| 269 | + if missing_rx_product: |
| 270 | + print('The following reaction products have no decay data:') |
| 271 | + for vals in missing_rx_product: |
| 272 | + print('{} {} -> {}'.format(*vals)) |
| 273 | + print('') |
| 274 | + |
| 275 | + if missing_fpy: |
| 276 | + print('The following fissionable nuclides have no fission product yields:') |
| 277 | + for parent, replacement in missing_fpy: |
| 278 | + print(' {}, replaced with {}'.format(parent, replacement)) |
| 279 | + print('') |
| 280 | + |
| 281 | + chain.export_to_xml('chain_casl.xml') |
| 282 | + |
| 283 | + |
| 284 | +if __name__ == '__main__': |
| 285 | + main() |
0 commit comments