Commit 8907be42 authored by William A. Wieselquist's avatar William A. Wieselquist
Browse files

Add example ii.json processing script

parent 274550ca
Loading
Loading
Loading
Loading
+54 −0
Original line number Diff line number Diff line
@@ -55,3 +55,57 @@ presented results and support with modeling challenges.
## Inventory Interface File

This work has developed a new inventory interface (`ii`) file format for exchanging nuclide inventory data with the MELCOR team. It is a heirarchical description of the inventory in one or more zones, at one or more time points. Currently the `ii` (pronounced aye-aye) data is emitted in standard [JSON](https://www.json.org/json-en.html) format so that one can read the files outside of SCALE. We currently use `ii.json` as the extension for these files. In the future we may provide an [HDF5](https://www.hdfgroup.org/solutions/hdf5/) version as well to save disk space--these would have the extension `.ii.h5`. We made the decision to package relevant nuclear data with the inventory in the `ii` format because typically downstream operations would like to calculate decay heat, mass, activity, and other derived quantities from the number of atoms/moles of each nuclide. Instead of providing number of moles, grams, and activity for each nuclide, we provide the molar mass (grams/mol) and decay constant (1/s) for each nuclide. Including nuclear data also aids us in nuclear data uncertainty studies where we would generate hundreds of `ii.json` files for each realization of the nuclear data. Each one of these realizations is the initial condition for a MELCOR severe accident analysis realization. Although the variation of the inventory in moles has the implicit effect of variations in nuclear data, the direct effect on decay heat and activity requires the perturbed recoverable energy from decay values and the decay constants. Additional details on the format are provided in this copy of the SCALE Quality Assurance (QA) case for this feature: [SCL-2020-003](supplemental/scl-2020-003_inventory-interface.md).

## Example of Processing the Inventory Interface File

An example script `supplemental/example_out.py` shows how to process the `ii.json`
file. It takes the file name and the label for the response. The `case(0)` 
response corresponds to a core total inventory.

```
python3 supplemental/example_out.py INL-A/full_core/depletion/INL-A_core_t6-depl_ce_v7.1.ii.json 'case(0)'
```

The screen output should be:
```
available responses for file INL-A/full_core/depletion/INL-A_core_t6-depl_ce_v7.1.ii.json : case(-1) case(-2) case(0) case(100) case(101) case(102) case(103) case(104) case(105) case(106) case(107) case(108) case(109) case(110) case(111) case(112) case(113) case(114) case(115) case(116) case(117) case(118) case(119) case(120) case(121) case(122) case(123) case(124) case(125) case(126) case(127) case(128) case(129) case(130) case(131) case(132) case(133) case(134) case(135) case(136) case(137) case(138) case(139) case(140) case(141) case(142) case(143) case(144) case(145) case(146) case(147) case(148) case(149) case(150) case(151) case(152) case(153) case(154) case(155) case(156) case(157) case(158) case(159) case(160) case(161) case(162) case(163) case(164) case(165) case(166) case(167) case(168) case(169) case(170) case(171) case(172) case(173) case(174) case(175) case(176) case(177) case(178) case(179) case(180) case(61) case(62) case(71) case(72) case(73) case(74) case(75) case(76) case(81) case(82) case(83) case(84) case(85) case(86) case(87) case(88) case(89) case(90) case(91) case(92) case(93) case(94) case(95) case(96) case(97) case(98) case(99)
   found response case(0) !
   found time list [0.0, 259200.0, 13365490.0, 26471780.0, 39578068.0, 52684360.0, 65790648.0, 78896936.0, 92046424.0, 105195920.0, 118345408.0, 131494896.0, 144644384.0, 157793872.0] in SECONDS !

SYSTEM MASS
     time(s)      mass(g)
 0.00000e+00  2.99799e+06
 2.59200e+05  2.99799e+06
 1.33655e+07  2.99799e+06
 2.64718e+07  2.99799e+06
 3.95781e+07  2.99799e+06
 5.26844e+07  2.99799e+06
 6.57906e+07  2.99799e+06
 7.88969e+07  2.99799e+06
 9.20464e+07  2.99799e+06
 1.05196e+08  2.99799e+06
 1.18345e+08  2.99799e+06
 1.31495e+08  2.99799e+06
 1.44644e+08  2.99799e+06
 1.57794e+08  2.99799e+06

TARGET NUCLIDE Cs137
     time(s)      mass(g) activity(Ci)
 0.00000e+00  1.00450e-14  1.37523e-12
 2.59200e+05  9.02606e-04  1.23573e-01
 1.33655e+07  4.63803e-02  6.34979e+00
 2.64718e+07  9.14244e-02  1.25166e+01
 3.95781e+07  1.36039e-01  1.86247e+01
 5.26844e+07  1.80228e-01  2.46746e+01
 6.57906e+07  2.23996e-01  3.06667e+01
 7.88969e+07  2.67347e-01  3.66017e+01
 9.20464e+07  3.10425e-01  4.24993e+01
 1.05196e+08  3.53090e-01  4.83405e+01
 1.18345e+08  3.95348e-01  5.41259e+01
 1.31495e+08  4.37201e-01  5.98559e+01
 1.44644e+08  4.78653e-01  6.55310e+01
 1.57794e+08  5.19709e-01  7.11519e+01
 ```

The script can be modified to iterate through nuclides and output any desired
quantity.
+94 −0
Original line number Diff line number Diff line
#import matplotlib.pyplot as plt
import numpy as np
import json
import sys    
import pprint as pp
import os

def canonical_name(izzzaaa,elem):
    a = izzzaaa % 1000
    name = elem+str(a)
    i = izzzaaa // 1000000
    if( i>0 ):
        name+='m'+str(i)
    return name

def read_data_set(file, name):

    # Load JSON file and data set name
    with open(file) as json_file:
        x = json.load(json_file)

    # Pull out data and response definitions first
    nucdb = x['data'].pop('nuclides')
    defines = x['responses'].pop('definitions')

    # Then print available data sets
    print('available responses for file',file,':',*x['responses'])
    
    # Then extract this response data set
    resp = x['responses'][name]
    vecdef_name = resp['nuclideVector']
    vecdef = defines['nuclideVectors'][vecdef_name]

    print('   found response',name,'!')

    # Initialize dictionary for processed data
    data=dict()
    data['name']=os.path.basename(file)+':'+name
    
    times = resp['timeList']['list']
    units = resp['timeList']['units']  
    data['seconds']=times
    data['elements']=[]
    data['nuclides']=[]
    print('   found time list',times,'in',units,'!')
 
    # Store the mass here
    data['grams']=[0] * len(times)
    
    for inuc in range(0,len(vecdef)):
        id=vecdef[inuc]
        nucdata=nucdb[id]
        # Add unit conversions to nucdata
        na = 6.022e23; #atoms/mole
        lam = nucdata['decayConstant'] #1/seconds
        q = nucdata['decayEnergy'] #Joules/decay
        ci_per_bq = 3.7e10; #Curies/Becquerel
        nucdata['moles_to_watts']=na*lam*q;
        nucdata['moles_to_curies']=lam*na/ci_per_bq;
        nucdata['moles_to_grams']=nucdata['mass'];
        izzzaaa=nucdata['IZZZAAA']
        elem=nucdata['element']  
        name=canonical_name(izzzaaa,elem)
        data['elements'].append(elem)
        data['nuclides'].append(name)
        data[name]=dict()
        data[name]['data']=nucdata
        data[name]['moles']=[0] * len(times)
        for itime in range(0,len(times)):
            amount=resp['amountList'][itime][inuc]
            data[name]['moles'][itime]=amount
            data['grams'][itime] += nucdata['moles_to_grams']*amount

    # Uniquify
    data['elements']=list(set(data['elements']))
    data['nuclides']=list(set(data['nuclides']))
    
    return data

d = read_data_set(sys.argv[1],sys.argv[2])

# print the whole data structure
# pp.pprint(d)
print('\nSYSTEM MASS\n{:>12s} {:>12s}'.format('time(s)','mass(g)'))
for t,m in zip(d['seconds'],d['grams']):
    print('{:12.5e} {:12.5e}'.format(t,m))

target='Cs137'
print('\nTARGET NUCLIDE {:12s}\n{:>12s} {:>12s} {:>12s}'.format(target,'time(s)','mass(g)','activity(Ci)'))
times=d['seconds']
nucdata = d[target]['data']
for itime in range(0,len(times)):
    moles = d[target]['moles'][itime]
    print('{:12.5e} {:12.5e} {:12.5e}'.format(times[itime],moles,moles*nucdata['moles_to_grams'],moles*nucdata['moles_to_curies']))