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