1from collections import namedtuple 2from io import StringIO 3import os 4import shutil 5from subprocess import Popen, PIPE, STDOUT, CalledProcessError 6import tempfile 7from pathlib import Path 8 9from . import endf 10 11 12# For a given MAT number, give a name for the ACE table and a list of ZAID 13# identifiers. This is based on Appendix C in the ENDF manual. 14ThermalTuple = namedtuple('ThermalTuple', ['name', 'zaids', 'nmix']) 15_THERMAL_DATA = { 16 1: ThermalTuple('hh2o', [1001], 1), 17 2: ThermalTuple('parah', [1001], 1), 18 3: ThermalTuple('orthoh', [1001], 1), 19 5: ThermalTuple('hyh2', [1001], 1), 20 7: ThermalTuple('hzrh', [1001], 1), 21 8: ThermalTuple('hcah2', [1001], 1), 22 10: ThermalTuple('hice', [1001], 1), 23 11: ThermalTuple('dd2o', [1002], 1), 24 12: ThermalTuple('parad', [1002], 1), 25 13: ThermalTuple('orthod', [1002], 1), 26 14: ThermalTuple('dice', [1002], 1), 27 26: ThermalTuple('be', [4009], 1), 28 27: ThermalTuple('bebeo', [4009], 1), 29 28: ThermalTuple('bebe2c', [4009], 1), 30 30: ThermalTuple('graph', [6000, 6012, 6013], 1), 31 31: ThermalTuple('grph10', [6000, 6012, 6013], 1), 32 32: ThermalTuple('grph30', [6000, 6012, 6013], 1), 33 33: ThermalTuple('lch4', [1001], 1), 34 34: ThermalTuple('sch4', [1001], 1), 35 35: ThermalTuple('sch4p2', [1001], 1), 36 37: ThermalTuple('hch2', [1001], 1), 37 38: ThermalTuple('mesi00', [1001], 1), 38 39: ThermalTuple('lucite', [1001], 1), 39 40: ThermalTuple('benz', [1001, 6000, 6012], 2), 40 42: ThermalTuple('tol00', [1001], 1), 41 43: ThermalTuple('sisic', [14028, 14029, 14030], 1), 42 44: ThermalTuple('csic', [6000, 6012, 6013], 1), 43 45: ThermalTuple('ouo2', [8016, 8017, 8018], 1), 44 46: ThermalTuple('obeo', [8016, 8017, 8018], 1), 45 47: ThermalTuple('sio2-a', [8016, 8017, 8018, 14028, 14029, 14030], 3), 46 48: ThermalTuple('osap00', [92238], 1), 47 49: ThermalTuple('sio2-b', [8016, 8017, 8018, 14028, 14029, 14030], 3), 48 50: ThermalTuple('oice', [8016, 8017, 8018], 1), 49 51: ThermalTuple('od2o', [8016, 8017, 8018], 1), 50 52: ThermalTuple('mg24', [12024], 1), 51 53: ThermalTuple('al27', [13027], 1), 52 55: ThermalTuple('yyh2', [39089], 1), 53 56: ThermalTuple('fe56', [26056], 1), 54 58: ThermalTuple('zrzrh', [40000, 40090, 40091, 40092, 40094, 40096], 1), 55 59: ThermalTuple('si00', [14028], 1), 56 60: ThermalTuple('asap00', [13027], 1), 57 71: ThermalTuple('n-un', [7014, 7015], 1), 58 72: ThermalTuple('u-un', [92238], 1), 59 75: ThermalTuple('uuo2', [8016, 8017, 8018], 1), 60} 61 62 63def _get_thermal_data(ev, mat): 64 """Return appropriate ThermalTuple, accounting for bugs.""" 65 66 # JEFF assigns MAT=59 to Ca in CaH2 (which is supposed to be silicon). 67 if ev.info['library'][0] == 'JEFF': 68 if ev.material == 59: 69 if 'CaH2' in ''.join(ev.info['description']): 70 zaids = [20040, 20042, 20043, 20044, 20046, 20048] 71 return ThermalTuple('cacah2', zaids, 1) 72 73 # Before ENDF/B-VIII.0, crystalline graphite was MAT=31 74 if ev.info['library'] != ('ENDF/B', 8, 0): 75 if ev.material == 31: 76 return _THERMAL_DATA[30] 77 78 # ENDF/B incorrectly assigns MAT numbers for UO2 79 # 80 # Material | ENDF Manual | VII.0 | VII.1 | VIII.0 81 # ---------|-------------|-------|-------|------- 82 # O in UO2 | 45 | 75 | 75 | 75 83 # U in UO2 | 75 | 76 | 48 | 48 84 if ev.info['library'][0] == 'ENDF/B': 85 if ev.material == 75: 86 return _THERMAL_DATA[45] 87 version = ev.info['library'][1:] 88 if version in ((7, 1), (8, 0)) and ev.material == 48: 89 return _THERMAL_DATA[75] 90 if version == (7, 0) and ev.material == 76: 91 return _THERMAL_DATA[75] 92 93 # If not a problematic material, use the dictionary as is 94 return _THERMAL_DATA[mat] 95 96 97_TEMPLATE_RECONR = """ 98reconr / %%%%%%%%%%%%%%%%%%% Reconstruct XS for neutrons %%%%%%%%%%%%%%%%%%%%%%% 99{nendf} {npendf} 100'{library} PENDF for {zsymam}'/ 101{mat} 2/ 102{error}/ err 103'{library}: {zsymam}'/ 104'Processed by NJOY'/ 1050/ 106""" 107 108_TEMPLATE_BROADR = """ 109broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%% 110{nendf} {npendf} {nbroadr} 111{mat} {num_temp} 0 0 0. / 112{error}/ errthn 113{temps} 1140/ 115""" 116 117_TEMPLATE_HEATR = """ 118heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%% 119{nendf} {nheatr_in} {nheatr} / 120{mat} 4 0 0 0 / 121302 318 402 444 / 122""" 123 124_TEMPLATE_HEATR_LOCAL = """ 125heatr / %%%%%%%%%%%%%%%%% Add heating kerma (local photons) %%%%%%%%%%%%%%%%%%%% 126{nendf} {nheatr_in} {nheatr_local} / 127{mat} 4 0 0 1 / 128302 318 402 444 / 129""" 130 131_TEMPLATE_GASPR = """ 132gaspr / %%%%%%%%%%%%%%%%%%%%%%%%% Add gas production %%%%%%%%%%%%%%%%%%%%%%%%%%% 133{nendf} {ngaspr_in} {ngaspr} / 134""" 135 136_TEMPLATE_PURR = """ 137purr / %%%%%%%%%%%%%%%%%%%%%%%% Add probability tables %%%%%%%%%%%%%%%%%%%%%%%%% 138{nendf} {npurr_in} {npurr} / 139{mat} {num_temp} 1 20 64 / 140{temps} 1411.e10 1420/ 143""" 144 145_TEMPLATE_ACER = """ 146acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%% 147{nendf} {nacer_in} 0 {nace} {ndir} 1481 0 1 .{ext} / 149'{library}: {zsymam} at {temperature}'/ 150{mat} {temperature} 1511 1/ 152/ 153""" 154 155_THERMAL_TEMPLATE_THERMR = """ 156thermr / %%%%%%%%%%%%%%%% Add thermal scattering data (free gas) %%%%%%%%%%%%%%% 1570 {nthermr1_in} {nthermr1} 1580 {mat} 12 {num_temp} 1 0 {iform} 1 221 1/ 159{temps} 160{error} {energy_max} 161thermr / %%%%%%%%%%%%%%%% Add thermal scattering data (bound) %%%%%%%%%%%%%%%%%% 162{nthermal_endf} {nthermr2_in} {nthermr2} 163{mat_thermal} {mat} 16 {num_temp} {inelastic} {elastic} {iform} {natom} 222 1/ 164{temps} 165{error} {energy_max} 166""" 167 168_THERMAL_TEMPLATE_ACER = """ 169acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%% 170{nendf} {nthermal_acer_in} 0 {nace} {ndir} 1712 0 1 .{ext}/ 172'{library}: {zsymam_thermal} processed by NJOY'/ 173{mat} {temperature} '{data.name}' / 174{zaids} / 175222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} {iwt}/ 176""" 177 178 179def run(commands, tapein, tapeout, input_filename=None, stdout=False, 180 njoy_exec='njoy'): 181 """Run NJOY with given commands 182 183 Parameters 184 ---------- 185 commands : str 186 Input commands for NJOY 187 tapein : dict 188 Dictionary mapping tape numbers to paths for any input files 189 tapeout : dict 190 Dictionary mapping tape numbers to paths for any output files 191 input_filename : str, optional 192 File name to write out NJOY input commands 193 stdout : bool, optional 194 Whether to display output when running NJOY 195 njoy_exec : str, optional 196 Path to NJOY executable 197 198 Raises 199 ------ 200 subprocess.CalledProcessError 201 If the NJOY process returns with a non-zero status 202 203 """ 204 205 if input_filename is not None: 206 with open(str(input_filename), 'w') as f: 207 f.write(commands) 208 209 with tempfile.TemporaryDirectory() as tmpdir: 210 # Copy evaluations to appropriates 'tapes' 211 for tape_num, filename in tapein.items(): 212 tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num)) 213 shutil.copy(str(filename), tmpfilename) 214 215 # Start up NJOY process 216 njoy = Popen([njoy_exec], cwd=tmpdir, stdin=PIPE, stdout=PIPE, 217 stderr=STDOUT, universal_newlines=True) 218 219 njoy.stdin.write(commands) 220 njoy.stdin.flush() 221 lines = [] 222 while True: 223 # If process is finished, break loop 224 line = njoy.stdout.readline() 225 if not line and njoy.poll() is not None: 226 break 227 228 lines.append(line) 229 if stdout: 230 # If user requested output, print to screen 231 print(line, end='') 232 233 # Check for error 234 if njoy.returncode != 0: 235 raise CalledProcessError(njoy.returncode, njoy_exec, 236 ''.join(lines)) 237 238 # Copy output files back to original directory 239 for tape_num, filename in tapeout.items(): 240 tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num)) 241 if os.path.isfile(tmpfilename): 242 shutil.move(tmpfilename, str(filename)) 243 244 245def make_pendf(filename, pendf='pendf', error=0.001, stdout=False): 246 """Generate pointwise ENDF file from an ENDF file 247 248 Parameters 249 ---------- 250 filename : str 251 Path to ENDF file 252 pendf : str, optional 253 Path of pointwise ENDF file to write 254 error : float, optional 255 Fractional error tolerance for NJOY processing 256 stdout : bool 257 Whether to display NJOY standard output 258 259 Raises 260 ------ 261 subprocess.CalledProcessError 262 If the NJOY process returns with a non-zero status 263 264 """ 265 266 make_ace(filename, pendf=pendf, error=error, broadr=False, 267 heatr=False, purr=False, acer=False, stdout=stdout) 268 269 270def make_ace(filename, temperatures=None, acer=True, xsdir=None, 271 output_dir=None, pendf=False, error=0.001, broadr=True, 272 heatr=True, gaspr=True, purr=True, evaluation=None, **kwargs): 273 """Generate incident neutron ACE file from an ENDF file 274 275 File names can be passed to 276 ``[acer, xsdir, pendf, broadr, heatr, gaspr, purr]`` 277 to specify the exact output for the given module. 278 Otherwise, the files will be writen to the current directory 279 or directory specified by ``output_dir``. Default file 280 names mirror the variable names, e.g. ``heatr`` output 281 will be written to a file named ``heatr`` unless otherwise 282 specified. 283 284 Parameters 285 ---------- 286 filename : str 287 Path to ENDF file 288 temperatures : iterable of float, optional 289 Temperatures in Kelvin to produce ACE files at. If omitted, data is 290 produced at room temperature (293.6 K). 291 acer : bool or str, optional 292 Flag indicating if acer should be run. If a string is give, write the 293 resulting ``ace`` file to this location. Path of ACE file to write. 294 Defaults to ``"ace"`` 295 xsdir : str, optional 296 Path of xsdir file to write. Defaults to ``"xsdir"`` in the same 297 directory as ``acer`` 298 output_dir : str, optional 299 Directory to write output for requested modules. If not provided 300 and at least one of ``[pendf, broadr, heatr, gaspr, purr, acer]`` 301 is ``True``, then write output files to current directory. If given, 302 must be a path to a directory. 303 pendf : str, optional 304 Path of pendf file to write. If omitted, the pendf file is not saved. 305 error : float, optional 306 Fractional error tolerance for NJOY processing 307 broadr : bool or str, optional 308 Indicating whether to Doppler broaden XS when running NJOY. If string, 309 write the output tape to this file. 310 heatr : bool or str, optional 311 Indicating whether to add heating kerma when running NJOY. If string, 312 write the output tape to this file. 313 gaspr : bool or str, optional 314 Indicating whether to add gas production data when running NJOY. 315 If string, write the output tape to this file. 316 purr : bool or str, optional 317 Indicating whether to add probability table when running NJOY. 318 If string, write the output tape to this file. 319 evaluation : openmc.data.endf.Evaluation, optional 320 If the ENDF file contains multiple material evaluations, this argument 321 indicates which evaluation should be used. 322 **kwargs 323 Keyword arguments passed to :func:`openmc.data.njoy.run` 324 325 Raises 326 ------ 327 subprocess.CalledProcessError 328 If the NJOY process returns with a non-zero status 329 IOError 330 If ``output_dir`` does not point to a directory 331 332 """ 333 if output_dir is None: 334 output_dir = Path() 335 else: 336 output_dir = Path(output_dir) 337 if not output_dir.is_dir(): 338 raise IOError("{} is not a directory".format(output_dir)) 339 340 ev = evaluation if evaluation is not None else endf.Evaluation(filename) 341 mat = ev.material 342 zsymam = ev.target['zsymam'] 343 344 # Determine name of library 345 library = '{}-{}.{}'.format(*ev.info['library']) 346 347 if temperatures is None: 348 temperatures = [293.6] 349 num_temp = len(temperatures) 350 temps = ' '.join(str(i) for i in temperatures) 351 352 # Create njoy commands by modules 353 commands = "" 354 355 nendf, npendf = 20, 21 356 tapein = {nendf: filename} 357 tapeout = {} 358 if pendf: 359 tapeout[npendf] = (output_dir / "pendf") if pendf is True else pendf 360 361 # reconr 362 commands += _TEMPLATE_RECONR 363 nlast = npendf 364 365 # broadr 366 if broadr: 367 nbroadr = nlast + 1 368 tapeout[nbroadr] = (output_dir / "broadr") if broadr is True else broadr 369 commands += _TEMPLATE_BROADR 370 nlast = nbroadr 371 372 # heatr 373 if heatr: 374 nheatr_in = nlast 375 nheatr_local = nheatr_in + 1 376 tapeout[nheatr_local] = (output_dir / "heatr_local") if heatr is True \ 377 else heatr + '_local' 378 commands += _TEMPLATE_HEATR_LOCAL 379 nheatr = nheatr_local + 1 380 tapeout[nheatr] = (output_dir / "heatr") if heatr is True else heatr 381 commands += _TEMPLATE_HEATR 382 nlast = nheatr 383 384 # gaspr 385 if gaspr: 386 ngaspr_in = nlast 387 ngaspr = ngaspr_in + 1 388 tapeout[ngaspr] = (output_dir / "gaspr") if gaspr is True else gaspr 389 commands += _TEMPLATE_GASPR 390 nlast = ngaspr 391 392 # purr 393 if purr: 394 npurr_in = nlast 395 npurr = npurr_in + 1 396 tapeout[npurr] = (output_dir / "purr") if purr is True else purr 397 commands += _TEMPLATE_PURR 398 nlast = npurr 399 400 commands = commands.format(**locals()) 401 402 # acer 403 if acer: 404 nacer_in = nlast 405 for i, temperature in enumerate(temperatures): 406 # Extend input with an ACER run for each temperature 407 nace = nacer_in + 1 + 2*i 408 ndir = nace + 1 409 ext = '{:02}'.format(i + 1) 410 commands += _TEMPLATE_ACER.format(**locals()) 411 412 # Indicate tapes to save for each ACER run 413 tapeout[nace] = output_dir / "ace_{:.1f}".format(temperature) 414 tapeout[ndir] = output_dir / "xsdir_{:.1f}".format(temperature) 415 commands += 'stop\n' 416 run(commands, tapein, tapeout, **kwargs) 417 418 if acer: 419 ace = (output_dir / "ace") if acer is True else Path(acer) 420 xsdir = (ace.parent / "xsdir") if xsdir is None else xsdir 421 with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file: 422 for temperature in temperatures: 423 # Get contents of ACE file 424 text = (output_dir / "ace_{:.1f}".format(temperature)).read_text() 425 426 # If the target is metastable, make sure that ZAID in the ACE 427 # file reflects this by adding 400 428 if ev.target['isomeric_state'] > 0: 429 mass_first_digit = int(text[3]) 430 if mass_first_digit <= 2: 431 text = text[:3] + str(mass_first_digit + 4) + text[4:] 432 433 # Concatenate into destination ACE file 434 ace_file.write(text) 435 436 # Concatenate into destination xsdir file 437 xsdir_in = output_dir / "xsdir_{:.1f}".format(temperature) 438 xsdir_file.write(xsdir_in.read_text()) 439 440 # Remove ACE/xsdir files for each temperature 441 for temperature in temperatures: 442 (output_dir / "ace_{:.1f}".format(temperature)).unlink() 443 (output_dir / "xsdir_{:.1f}".format(temperature)).unlink() 444 445 446def make_ace_thermal(filename, filename_thermal, temperatures=None, 447 ace='ace', xsdir=None, output_dir=None, error=0.001, 448 iwt=2, evaluation=None, evaluation_thermal=None, **kwargs): 449 """Generate thermal scattering ACE file from ENDF files 450 451 Parameters 452 ---------- 453 filename : str 454 Path to ENDF neutron sublibrary file 455 filename_thermal : str 456 Path to ENDF thermal scattering sublibrary file 457 temperatures : iterable of float, optional 458 Temperatures in Kelvin to produce data at. If omitted, data is produced 459 at all temperatures given in the ENDF thermal scattering sublibrary. 460 ace : str, optional 461 Path of ACE file to write 462 xsdir : str, optional 463 Path of xsdir file to write. Defaults to ``"xsdir"`` in the same 464 directory as ``ace`` 465 output_dir : str, optional 466 Directory to write ace and xsdir files. If not provided, then write 467 output files to current directory. If given, must be a path to a 468 directory. 469 error : float, optional 470 Fractional error tolerance for NJOY processing 471 iwt : int 472 `iwt` parameter used in NJOR/ACER card 9 473 evaluation : openmc.data.endf.Evaluation, optional 474 If the ENDF neutron sublibrary file contains multiple material 475 evaluations, this argument indicates which evaluation to use. 476 evaluation_thermal : openmc.data.endf.Evaluation, optional 477 If the ENDF thermal scattering sublibrary file contains multiple 478 material evaluations, this argument indicates which evaluation to use. 479 **kwargs 480 Keyword arguments passed to :func:`openmc.data.njoy.run` 481 482 Raises 483 ------ 484 subprocess.CalledProcessError 485 If the NJOY process returns with a non-zero status 486 487 """ 488 if output_dir is None: 489 output_dir = Path() 490 else: 491 output_dir = Path(output_dir) 492 if not output_dir.is_dir(): 493 raise IOError("{} is not a directory".format(output_dir)) 494 495 ev = evaluation if evaluation is not None else endf.Evaluation(filename) 496 mat = ev.material 497 zsymam = ev.target['zsymam'] 498 499 ev_thermal = (evaluation_thermal if evaluation_thermal is not None 500 else endf.Evaluation(filename_thermal)) 501 mat_thermal = ev_thermal.material 502 zsymam_thermal = ev_thermal.target['zsymam'] 503 504 # Determine name, isotopes based on MAT number 505 data = _get_thermal_data(ev_thermal, mat_thermal) 506 zaids = ' '.join(str(zaid) for zaid in data.zaids[:3]) 507 508 # Determine name of library 509 library = '{}-{}.{}'.format(*ev_thermal.info['library']) 510 511 # Determine if thermal elastic is present 512 if (7, 2) in ev_thermal.section: 513 elastic = 1 514 mt_elastic = 223 515 516 # Determine whether elastic is incoherent (0) or coherent (1) 517 file_obj = StringIO(ev_thermal.section[7, 2]) 518 elastic_type = endf.get_head_record(file_obj)[2] - 1 519 else: 520 elastic = 0 521 mt_elastic = 0 522 elastic_type = 0 523 524 # Determine number of principal atoms 525 file_obj = StringIO(ev_thermal.section[7, 4]) 526 items = endf.get_head_record(file_obj) 527 items, values = endf.get_list_record(file_obj) 528 energy_max = values[3] 529 natom = int(values[5]) 530 531 # Note that the 'iform' parameter is omitted in NJOY 99. We assume that the 532 # user is using NJOY 2012 or later. 533 iform = 0 534 inelastic = 2 535 536 # Determine temperatures from MF=7, MT=4 if none were specified 537 if temperatures is None: 538 file_obj = StringIO(ev_thermal.section[7, 4]) 539 endf.get_head_record(file_obj) 540 endf.get_list_record(file_obj) 541 endf.get_tab2_record(file_obj) 542 params = endf.get_tab1_record(file_obj)[0] 543 temperatures = [params[0]] 544 for i in range(params[2]): 545 temperatures.append(endf.get_list_record(file_obj)[0][0]) 546 547 num_temp = len(temperatures) 548 temps = ' '.join(str(i) for i in temperatures) 549 550 # Create njoy commands by modules 551 commands = "" 552 553 nendf, nthermal_endf, npendf = 20, 21, 22 554 tapein = {nendf: filename, nthermal_endf: filename_thermal} 555 tapeout = {} 556 557 # reconr 558 commands += _TEMPLATE_RECONR 559 nlast = npendf 560 561 # broadr 562 nbroadr = nlast + 1 563 commands += _TEMPLATE_BROADR 564 nlast = nbroadr 565 566 # thermr 567 nthermr1_in = nlast 568 nthermr1 = nthermr1_in + 1 569 nthermr2_in = nthermr1 570 nthermr2 = nthermr2_in + 1 571 commands += _THERMAL_TEMPLATE_THERMR 572 nlast = nthermr2 573 574 commands = commands.format(**locals()) 575 576 # acer 577 nthermal_acer_in = nlast 578 for i, temperature in enumerate(temperatures): 579 # Extend input with an ACER run for each temperature 580 nace = nthermal_acer_in + 1 + 2*i 581 ndir = nace + 1 582 ext = '{:02}'.format(i + 1) 583 commands += _THERMAL_TEMPLATE_ACER.format(**locals()) 584 585 # Indicate tapes to save for each ACER run 586 tapeout[nace] = output_dir / "ace_{:.1f}".format(temperature) 587 tapeout[ndir] = output_dir / "xsdir_{:.1f}".format(temperature) 588 commands += 'stop\n' 589 run(commands, tapein, tapeout, **kwargs) 590 591 ace = output_dir / ace 592 xsdir = (ace.parent / "xsdir") if xsdir is None else Path(xsdir) 593 with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file: 594 # Concatenate ACE and xsdir files together 595 for temperature in temperatures: 596 ace_in = output_dir / "ace_{:.1f}".format(temperature) 597 ace_file.write(ace_in.read_text()) 598 599 xsdir_in = output_dir / "xsdir_{:.1f}".format(temperature) 600 xsdir_file.write(xsdir_in.read_text()) 601 602 # Remove ACE/xsdir files for each temperature 603 for temperature in temperatures: 604 (output_dir / "ace_{:.1f}".format(temperature)).unlink() 605 (output_dir / "xsdir_{:.1f}".format(temperature)).unlink() 606