luca-fiorito-11 / sandy

Sampling nuclear data and uncertainty

SANDY

Sampling tool for nuclear data

SANDY is a python package that can read, write and perform a set of operations on nuclear data files in ENDF-6 format.

Stochastic sampling of nuclear data

The primary objective of the code, as it was originally conceived, is to produce perturbed files containing sampled parameters that represent the information stored in the evaluated nuclear data covariances. Such files can be ultimately used to propagate uncertainties through any given compatible system using a brute force technique.

Currently, SANDY can draw samples for:

• cross sections;
• angular distrbutions of outgoing particles;
• energy distrbutions of outgoing particles;
• fission neutron multiplicities;
• fission yields.

API for ENDF-6 files

The recent development on SANDY extended the original goal and focused on providing a simple interface for nuclear data files in ENDF-6 format. Nuclear data such as cross sections, fission yields, radioactive decay constants and so on can be imported into tabulated dataframes (making extensive use of pandas) for further post-processing, analysis, plotting, ...

Examples are available here.

---

🔧 Installation

SANDY can be installed both on Linux (recommended) or Windows (using Anaconda). The installation instructions are available here.

⌛ Development history and releases

The latest and older releases of SANDY are available here.

📔 Documentation

The official SANDY documentation can be found here.

🎮 Jupyter notebooks

Here you can find some cool Jupyter notebooks that kind of give an idea of what one can do with SANDY.

📞 Contacts

• Luca Fiorito - [email protected]

🔖 Acknowledgments

SANDY was conceived and developed as a part of the PhD thesis on Nuclear data uncertainty propagation and uncertainty quantification in nuclear codes in the framework of a collaboration between SCK CEN and ULB.

📋 Reference

Among the publications about SANDY, please use the following as references for citation.

• L. Fiorito, J. Dyrda and M. Fleming, JEFF-3.3 covariance application to ICSBEP using SANDY and NDAST, EPJ Web of Conferences 211, 07003 (2019)

• L. Fiorito, G. Žerovnik, A. Stankovskiy, G. Van den Eynde, P.E. Labeau, Nuclear data uncertainty propagation to integral responses using SANDY, Annals of Nuclear Energy, Volume 101, 2017, Pages 359-366, ISSN 0306-4549.

🌍 Publications

Here is a (incomplete) list of scientific studies citing SANDY.

If some info are not correct or missing, please let us know!

sandy's Issues

ACE file generation

ACE files are not correctly generated when using SANDY to sample ENDF files and afterwards automatically process the perturbed files with NJOY. Only the ACE and xsdir for the last perturbed file are generated, even though NJOY runs for each perturbed file, probably because each time, the perturbed file is written with the identifier of the number of samples instead of the identifier number of the perturbed file, so they get overwritten continuously.

ENDF-6 with MF32 and not MF33

ENDF-6 files with only MF32 and not MF33 are not processed with the default ERRORR input options.
This can solve the problem.

nuclides with such condition are 95-Am-241, 17-Cl- 35, 17-Cl- 37, 91-Pa-231 and 91-Pa-233 in JEFF-3.3

ENDF-6 file title

When concatenating ENDF-6 files, it would be nice to have a way to remove the file TITLE (first line) to all concatenated files except the first.
The same way the last line (FEND, or end-of-file line) should be removed from all files but the last.

Generating non-gaussian samples from TENDL

Hi,

I'm using sandy to generate some random samples of Fe56 from TENDL 2017, and it's giving me some non-gaussian samples.

Would you have any clues as to why this is happening?

Fortran integer issue in Windows

Integer definitions are inconsistent in rwfortran.f.

nubar MT selection with NJOY

When I select the nubar MT numbers and try to run get_errorr I get an error in NJOY.
The error is related to the selection of the MT numbers when NJOY runs groupr. An example is reported below.

mt=[456,455]
# ERRORR parameters
ek_errorr = sandy.energy_grids.ECCO33  # 33 group energy structure
err = 1  # RECONR recostruction tolerance, low for fast calculations
errorr = sandy.get_endf6_file("jeff_33", "xs", ZAMs[1]).get_errorr(err=err, ek_errorr=ek_errorr, mt=mt, nubar=True, xs=False, mubar=False, chi=False, verbose= True)

Read perturbations coefficients from file

Reading perturbation coefficients that were previously generated ensures reproducibility and could save disk space (no need to store all samples all the times). Also, only a fraction of samples could be generated from the perturbation coefficients, leaving a second fraction for later.

warning for negative eigenvalues

A warning should be issued when large negative eigenvalues are found in a covariance matrix.

Print covariance eigenvalues and eigenvectors in sampling procedure

The sampling procedure shall give the user the possibility to access the sampled coefficient, but also the eigenvalues and eigenvectors of the original covariance matrix.

reading TEMP and ERR in PENDF

PENDF files produced by NJOY add two parameters in the file header MF1/MT451:

• TEMP: temperature of the evaluation (position C1)
• ERR: reconstruction tolerance used by RECONR (position C2)

SANDY should be able to read and write these parameters and it currently does not do it.

Need for additional energy grids

More tabulated multi-group energy grids should be provided as a part of the code.

Cholesky decomposition

Discussed in #133

^{Originally posted by EnricaBelfiore February 18, 2022}
@luca-fiorito-11 @AitorBengoechea

I and @GrimFe noticed that the last modification introduced in the sampling procedure does not work for all nuclear data covariance matrices. More in the details, we tested it for covariances of Pu239, U235, U233 and H1.
First of all, was this method already tested with real test cases? In case it was, which nuclides did you use for testing @AitorBengoechea?

We think that this problem is related to the fact that the covarinace for nuclear data are not always positive defined. As far as we understand, to_positive method handles this by transforming such matrices into SEMI-positive defined matrices. LU decomposition implemented inside sparse_table_cholesky method does not properly work for singular matrix (det=0).

We propose to back up to the previous decomposition scheme (cov = UDU.T) and to return L=U*sqrt(D), so that cov=L*L.T, where cov is the semipositive approximation of the nuclear data covariance matrix.
Unfortunately, scipy.linalg.ldl (implementing UDU.T decomposition) does not handle sparse matrix. Therefore our suggestion is to directly implement scipy.linalg.ldl decomposition in get_L. Further possible memory optimization can and should be considered.

Add an extra (optional) point at the end of the perturbation interval

The perturbation procedure from Sandy sometimes produces interpolation artifacts when the energy grid of the initial cross section is too sparse. This can be seen on the attached images, where the lines labeled "Broadr + perturbed" and "Broadr + perturbed + proc" feature a smooth tail from the maxiumum value to zero at the edge of the perturbation. From the code sample from the custom perturbation function:

            enew = np.union1d(self.data.index.values, pert.right.index.values)
            u_xs = self.reshape(enew)
            u_pert = pert.reshape(enew)
            u_xs.data[(mat, mt)] = u_xs.data[(mat, mt)] * u_pert.right.values

it is clear that the point marking the edge of the perturbation interval is added to the total grid, but another point right next to it could be added to define the edge of the perturbation. This additional point can be optional and selected by the user. I am currently working on this feature and will submit a pull request soon.

The same images show another potential issue in the sampling procedure, where the cross section can change because broadening happens after the perturbation, but I will open a separate issue for that and submit a separate pull request (edit: reference: #254 ) .

Do Doppler broadening before perturbation

The current default sampling procedure for cross sections read from ENDF files in Sandy can be summerized (according to my understanding, please correct if I misunderstood) as:

1. Read the ENDF file using the read_formatted_file(file, listmat=None, listmf=None, listmt=None) function
2. A PENDF with linearized cross sections is constructed and merged with the endf6 file
3. Samples are extracted from the covariance matrix in the same grid as the covariance matrix
4. An union energy grid is constructed based on the energy grids of the linearized cross sections and the samples is constructed
5. The cross section is interpolated to the new grid
6. The perturbations are "interpolated" to the new grid by filling empty values with ones (P.reindex(P.index.union(frame[mat,mt].index)).ffill().fillna(1).reindex(frame[mat,mt].index))
7. Perturbations are capped to (0,2)
8. Perturbations are multiplied with the cross section
9. The tape is updated and dumped to a file
10. (optional) a NJOY sequence is run to produce ace files

I have written out this procedure because I would like you to correct me if I misunderstood something.
The issue is that to use the cross sections, doppler broadening happens after the perturbations. This happens when using the default Sandy procedure or a custom processing procedure. The perturbations are seen as "small resonances" and get flattened by doppler broadening. I have attached two graphs of the results of sampling procedures where the MT18 cross section for U235 was perturbed by 30% in the WIMS group structure. The graphs show perturbations of the first two groups starting from the lowest energies by 30% (to better see the effects). According to my results, the perturbation of the thermal group flattens the perturbation, leading to possible underestimation of the perturbation effect. The same graph shows what happens when the sampling procedure is modified so the doppler broadening happens before the perturbation and other processing steps (reconstruction of resonances, unresolved resonances, ..) take place after. Using the modified procedure, the perturbation remains intact.

I am currently working on a patch that would implement this modified procedure in Sandy.

RDD data not found

with the dafault sandy installation python setup.py install, the local decay data libraries in sandy/appendix/decay_data are not included in the package.
The results is that an error is raised when running, say, sandy.get_endf6_file("jeff_33", "decay_data", "all").

convert fortran parser to python

Not in control format

With the newer endfb (and tendl) libraries I keep getting this error:

> sandy n-026_Fe_056.endf 

ERROR:  line ' $Rev:: 1518     $  $Date:: 2018-02-02#$                             1 0  0
' is not in CONTROL format

Issue in `CategoryCov.invert`

CategoryCov.invert fails for the following example:

sandy.get_endf6_file("jeff_33", "xs", 10010).get_errorr(err=1, mt=[102]).get_cov(multigroup=False)

with error:

Scipy.sparse RuntimeError: Factor is exactly singular

A similar issue can be found on stackoverflow at Scipy.sparse RuntimeError: Factor is exactly singular

Extracting MF=31 (nubar) covariance with SANDY API

Dear Sandy developers,
I am currently using SANDY as an API to generate the ERRORR files for different applications (e.g., singular value decomposition, sandwhich rule calculations using the sensitivity coeffs. computed by other codes, and so on). So far, I was able to generate the ERRORR files for the MF=32,33 and the related MT using some lines of code as follows,

endf = sandy.get_endf6_file(the_library_i_want, "xs", za)
pendf = endf.get_pendf(temperature=temperature)
myerrorr = endf.get_errorr(njoy=os.environ['NJOY'], pendf=pendf, temperature=temperature, ign=1,ek=group_structure)

Now I am interested in doing the same thing but for the nubar covariance (MF=31). How can I do that? Looking at the sandy.njoy.process method and sub-methods I do not actually see how to specify the mfcov option of the card 7 in the NJOY deck (which allows to choose among MF 31,32,33,...).

Could you please help me?
Thank you so much
Nicolo'

add njoy wrapper

Add njoy module at root level to produce predefined NJOY2016 inputs and to run the code.
The following function must be implemented:

• produce ace files for fast evaluations
• produce pendf files for fast evaluations
• produce ace files for thermal evaluations
• produce pendf files for thermal evaluations
• produce ace files for proton evaluations
• produce ace files for photon evaluations

Parsing KEFF in ALEPH output in NPS calculations

ALEPH-2 outputs for NPS calculations contain estimates of KEFF preceded by the following;

Keff estimate NPS problem

This is different than for KCODE calculations, where the following is reported

Keff  eff. mult. factor

These variations should be reflected also in the parsing algorithm in sandy.aleph2.OutputFile.

ENDF6 Cross Sections with many data points does not extract all data points from endf file

Hi, I am trying to extract all the datapoints from cross section files for things like U-235 and I am running into an issue where it seems like some part of the reader/extraction method in your code stops processing information after a certain amount of data points OR its like its missing a part of the file when reading and only getting a chunk of it.

The following is what I am running with:

import sandy
import pandas as pd
import plotly.graph_objects as go
import os
from scripts import functions as fn

Read ENDF-6 data for a specific nuclide (e.g., U-235)

#nuclide = fn.format_input(input("Input nuclide of interest\n"))
fy = sandy.get_endf6_file("endfb_80", "xs", 922350)

fyr = fy.get_records()

fydict = fy.read_section(mat=9228, mf=3, mt=18)

here you can see that it gets most of the upper energy datapoints in detail but then jumps to 1e-5 all of a sudden.

print(fydict['E'], fydict['XS'])

Create a scatter plot with log-log scale

fig = go.Figure()

fig.add_trace(go.Scatter(x=fydict['E'], y=fydict['XS'], mode="lines", name=f"{nuclide}", text="test"))

fig.show()

---

irespr=1

when running ERRORR, option IRESPR=1 should be preferred.
Although faster, IRESPR=0 cannot process the covariance matrix of the resonance parameters of many fission products (using JEFF-3.3) for whatever reason.

Fix docstrings in decay.py

few typos and inconsistencies in the docstrings of decay.py that were not tackles when PR #229 was merged.

Not urgent!

Parsing MF32 LRF=3 - Reich Moore

this feature is not implemented

multi-threading pandas / numpy /scipy

Set maximum number of threads created by some of the pandas/numpy/scipy functions.

export OMP_NUM_THREADS=$N
export NUMEXPR_NUM_THREADS=$N
export MKL_NUM_THREADS=$N
export MKL_DOMAIN_NUM_THREADS=$N

where N is the number of threads.

Behaviour for `sandy.sampling.run()` for libraries without covaraince data

Please note that I'm not too knowledgeable with the code, so feel free to comment angry messages and close the issue if this is in error :)

For example, for ENDF-B-VII.1-neutrons n-008-O_017.endf,

sandy.sampling.run() errors because the get_perturbations() will return an empty dictionary (i.e., smps={}). This will cause the apply_perturbations() function to error with the following error message:

Traceback (most recent call last):
  File "go.py", line 10, in <module>
    obj.run()
  File "/home/4ib/anaconda3/envs/env_neutronics/lib/python3.7/site-packages/fermi-0.1-py3.7.egg/fermi/neutronics/nuclear_data_perturbation/openmc_data_pert.py", line 40, in run
    self.generate_perturbed_data()
  File "/home/4ib/anaconda3/envs/env_neutronics/lib/python3.7/site-packages/fermi-0.1-py3.7.egg/fermi/neutronics/nuclear_data_perturbation/openmc_data_pert.py", line 95, in generate_perturbed_data
    sandy.sampling.run(' '.join(cli))
  File "/home/4ib/anaconda3/envs/env_neutronics/lib/python3.7/site-packages/sandy-1.0.37-py3.7.egg/sandy/sampling.py", line 185, in inner
    foo(iargs)
  File "/home/4ib/anaconda3/envs/env_neutronics/lib/python3.7/site-packages/sandy-1.0.37-py3.7.egg/sandy/sampling.py", line 305, in run
    verbose=iargs.debug,
  File "/home/4ib/anaconda3/envs/env_neutronics/lib/python3.7/site-packages/sandy-1.0.37-py3.7.egg/sandy/core/endf6.py", line 2204, in apply_perturbations
    if not item:
UnboundLocalError: local variable 'item' referenced before assignment

Is there a way to add an argument so that if the file does not have any covariance data it will just output something helpful? Thanks!

Integer initialization in rwfortran.f

It seems that one of the "integer*4" declarations did actually need to be specifically set to length 4; 8 causes it to crash. (line 156 inside of "rilist")

Apply samples to doppler-broadened cross sections

This issue was already raised in #254 .
Here more detailes are given on the implementation, not on the physics.

Applying perturbations to doppler-boradedn cross sections rquires:

• to use the method sandy.Endf6.apply_perturbations on a pregenerated PENDF file (already doppler-broadened)
• possibly to concatenate ENDF-6 objects at multiple temperatures in a single file

Describe sampling procedure for fission yields

No module named 'rwf' and module 'pandas' has no attribute 'Panel'

I tried to install SANDY via cloning the repository and running all the conda/pip lines quite fine up to the "on windows" lines, which I didn't do, of course. I am on a conda virtual environment (sandy-devel) as suggested. I met two problems: after having set the pythonpath I use ipython to try "import sandy" and I get the following error:

No module named 'rwf'

While, if I try the import sandy command on a VScode notebook (still in my sandy-devel environment), I get the following error:

AttributeError: module 'pandas' has no attribute 'Panel'

In any case I am stuck at the "import sandy" command.
I am on WSL2 (Ubuntu on Windows).

It takes forever to load decay data!

Loading decay data with sandy.get_endf6_file("jeff_33", "decay", "all") takes forever. The reason is that each decay data file is retrieved, unzipped and read from the web. This operation is repeated in series about 3000-4000 times to get all the nuclide inventory

sorting LPC when writing

Legendre polynomial coefficients in file MF4 are not sorted by energy when writing to file.

>5 digit numbers in SEND section

SEND pointers larger than 99999 should be handled properly

Remove dependence from pytables

functions explicitly using pytables are obsolete and shall be removed.

wrong implementation of `CategoryCov.sandwich`

CategoryCov.sandwich is not doing what it is supposed to, that is, an inner product.

'Edistr' object has no attribute 'add_points'

Hi,

I'm trying to generate random samples from an ENDF file of U-235, running SANDY from the Linux command line. There seems to be an issue when sampling the outgoing-particle energy distributions, as I get the following error:

edistr = sandy.Edistr.from_endf6(endfnew).add_points(extra_points)
AttributeError: 'Edistr' object has no attribute 'add_points'

It seems that there was an 'add_points' function in an earlier version of the Edistr class definition in pfns.py, but that this was removed at some point (perhaps replaced by add_energy_points, but the reference in sampling.py is not updated?).

/Markus

Add sampling procedure for Pu239 FYs

Use covariance data provided by OECD/NEA

reading xs from JENDL-4.0u

Handle e+0 Case in write_float Function

The write_float function handles the e-0 case but does not address the e+0 case, leading to inconsistent formatting in scientific notation. This inconsistency results in improperly formatted ENDF files, causing processing to fail. The function should be updated to replace both e-0 with - and e+0 with + to ensure proper formatting and successful processing.

Example where the old behavior fails:
sandy.write_float(1-1e-8)

Possible issue in getting covariance matrix eigenvalues

Discussed in #135

^{Originally posted by GrimFe February 18, 2022}
@luca-fiorito-11 @AitorBengoechea

When trying to identify the relevant eigenvalues of the Pu239 covariance matrix, @EnricaBelfiore and I noticed a possible issue with the latest implementation.
Only six eigenvalues (the larger ones) are now returned by eig method of CategoryCov.
This is of course not coherent with the matrix shape and should be clarified.
Even if in the reported case further eigenvalues are low (despite them not being zero), this truncation to the sixth eigenvalue happens also when considering U233 covariance matrix. There the seventh eigenvalue is of the same order of magnitude as the sixth, suggesting the truncation is not order of magnitude-based. Moreover, in our opinion all the eigenvalues should be returned to get the best possible matrix representation.

It is worth noticing that np.linalg.eig(cov.data) reports all the eigenvalues, meaning there is no lack of information, but rather an error in eigenvalue calculation. Browsing for this issue, we found this scipy/scipy#11198 reporting a possible error in scipy.sparse.linalg.eigs(), which might be a good starting point for debugging.

Example of what reported in the following - considering U233:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
tape = sandy.get_endf6_file("jeff_33", "xs", 922330)
Cov = sandy.XsCov.from_endf6(tape)
cov = sandy.CategoryCov(Cov)
len(cov.eig()[0])
len(pd.Series(np.linalg.eig(cov.data)[0]).sort_values(ascending=False))
len(pd.Series(scipy.linalg.eig(cov.data)[0]).sort_values(ascending=False))

EE = pd.Series(np.linalg.eig(cov.data)[0]).sort_values(ascending=False)
plt.plot(np.array(range(len(EE))), EE, 'o-')
plt.grid()

functions for q-matrix

Hi @m-fleming

please, review the list of items that you need.

Luca

HTTP Error 403: Forbidden

see #325

incompatibility with pandas>=1

Installation in Windows

Due to MSVC compiler not being supported by setuptools from version 65.0.0, sandy installation fails when doing python setup.py install with the following message:

from distutils.msvccompiler import get_build_version as get_build_msvc_version ModuleNotFoundError: No module named 'distutils.msvccompiler'

To fix this error, setuptools should be downgraded to the latest version compatible with MSVC:

conda install setuptools=64.0.3

Furthermore, to avoid future problems, it is recommended that the file requirements_conda.txt reflects this.

ALEPH parse issue for "Major contributors"

The ALEPH-2 parser cannot process correctly Table 3 and Table 4 (section "major contributors") when there is only 1 major contributor

Large number of points for E_in and E_out in MF5 can lead to allocation problems in NJOY (acer)

To better reproduce the MF35 covariance block structure for incoming neutron energy, SANDY adds extra incident energy points for which the outgoing particle distribution is interpolated.
To simplify the interpolation, all the distributions are restructured on a union grid containing all the points provided for each of them.
Consequently, it is not rare to see hundreds of incident incoming energies that have outgoing particle distributions defined over thousand of energy points.
This situation is not accepted by NJOY2016 (acer) that will raise a memory allocation error.

API update

Correct, document and test the API.

negative samples

When sampling from a Normal distribution with large standard deviations it is likely to draw negative perturbations.
Physically many quantities such as cross sections or energy distributions do not take negative values.

SANDY must handle the negative samples in a consistent manner.