import os
import numpy as np
import pandas as pd
import sandy
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
In this notebook we apply a constant relative perturbation, say 5%, to an energy-dependent cross section. The perturbation is applied over an energy bin of choice, e.g. between 10 eV and 100 eV.
As an example we use the JEFF-3.3 evaluation for Pu-239 downloaded from the OECD/NEA website and processed with NJOY-2016 into PENDF format using the following input file:
tape = sandy.get_endf6_file("jeff_33", "xs", 942390).get_pendf(
temperature=900,
purr=False,
heatr=False,
gaspr=False,
verbose=True,
)
moder 20 -21 / reconr -21 -22 / 'sandy runs njoy'/ 9437 0 0 / 0.001 0. / 0/ broadr -21 -22 -23 / 9437 1 0 0 0. / 0.001 / 900.0 / 0 / thermr 0 -23 -24 / 0 9437 20 1 1 0 0 1 221 0 / 900.0 / 0.001 10 / moder -24 30 / stop njoy 2016.70 04Apr23 07/05/23 15:16:06 ***************************************************************************** moder... 0.0s reconr... 0.1s ---message from rdf2bw---calculation of angular distribution not installed. broadr... 52.3s thermr... 75.3s wrote thermal data for temp = 9.0000E+02 79.7s moder... 79.7s 83.8s *****************************************************************************
xs = sandy.Xs.from_endf6(tape)
Now, we apply perturbations according to the WIMS-69 energy grid, available on the Serpent Wiki. Still, any other energy grids would work just fine.
egrid = sandy.energy_grids.WIMS69
egrid
array([1.0000e-05, 5.0000e-03, 1.0000e-02, 1.5000e-02, 2.0000e-02, 2.5000e-02, 3.0000e-02, 3.5000e-02, 4.2000e-02, 5.0000e-02, 5.8000e-02, 6.7000e-02, 8.0000e-02, 1.0000e-01, 1.4000e-01, 1.8000e-01, 2.2000e-01, 2.5000e-01, 2.8000e-01, 3.0000e-01, 3.2000e-01, 3.5000e-01, 4.0000e-01, 5.0000e-01, 6.2500e-01, 7.8000e-01, 8.5000e-01, 9.1000e-01, 9.5000e-01, 9.7200e-01, 9.9600e-01, 1.0200e+00, 1.0450e+00, 1.0710e+00, 1.0970e+00, 1.1230e+00, 1.1500e+00, 1.3000e+00, 1.5000e+00, 2.1000e+00, 2.6000e+00, 3.3000e+00, 4.0000e+00, 9.8770e+00, 1.5968e+01, 2.7700e+01, 4.8052e+01, 7.5501e+01, 1.4873e+02, 3.6726e+02, 9.0690e+02, 1.4251e+03, 2.2395e+03, 3.5191e+03, 5.5300e+03, 9.1180e+03, 1.5030e+04, 2.4780e+04, 4.0850e+04, 6.7340e+04, 1.1100e+05, 1.8300e+05, 3.0250e+05, 5.0000e+05, 8.2100e+05, 1.3530e+06, 2.2310e+06, 3.6790e+06, 6.0655e+06, 1.0000e+07])
The energy grid was converted in eV to be consistent with the cross section data.
sandy.Pert
object¶SANDY uses a dedicated class to store (relative) perturbation coefficients.
sandy.Pert([1, 1.05], index=[10, 100])
ENERGY (0.0, 10.0] 1.00000e+00 (10.0, 100.0] 1.05000e+00 dtype: float64
In this example the Pert
instance defines a multiplicative perturbation factor 1 (0%) for all xs values between 0 and 10 eV, and 1.05 (5%) for all xs between 10 and 100 eV.
Below we apply perturbations to the fission cross section (MT=18) of Pu-239 (MAT=9437). For each energy bin in the WIMS grid, we increase all xs points in that energy range by 30%. The process is repreated iteratively and 69 perturbed xs objects are created.
mat = 9437
mt = 18
pert_coeff = 30 / 100 # large perturbationfor better visualization
perturbed_xs = []
for i in range(1, egrid.size):
e_start = egrid[i-1]
e_stop = egrid[i]
index = egrid[i-1: i+1]
pert = sandy.Pert([1, 1 + pert_coeff], index=index)
print(f"perturbed xs in energy bin #{i} [{e_start:.5e}, {e_stop:.5e}]")
xspert = xs.custom_perturbation(mat, mt, pert)
perturbed_xs.append(xspert)
perturbed xs in energy bin #1 [1.00000e-05, 5.00000e-03] perturbed xs in energy bin #2 [5.00000e-03, 1.00000e-02] perturbed xs in energy bin #3 [1.00000e-02, 1.50000e-02] perturbed xs in energy bin #4 [1.50000e-02, 2.00000e-02] perturbed xs in energy bin #5 [2.00000e-02, 2.50000e-02] perturbed xs in energy bin #6 [2.50000e-02, 3.00000e-02] perturbed xs in energy bin #7 [3.00000e-02, 3.50000e-02] perturbed xs in energy bin #8 [3.50000e-02, 4.20000e-02] perturbed xs in energy bin #9 [4.20000e-02, 5.00000e-02] perturbed xs in energy bin #10 [5.00000e-02, 5.80000e-02] perturbed xs in energy bin #11 [5.80000e-02, 6.70000e-02] perturbed xs in energy bin #12 [6.70000e-02, 8.00000e-02] perturbed xs in energy bin #13 [8.00000e-02, 1.00000e-01] perturbed xs in energy bin #14 [1.00000e-01, 1.40000e-01] perturbed xs in energy bin #15 [1.40000e-01, 1.80000e-01] perturbed xs in energy bin #16 [1.80000e-01, 2.20000e-01] perturbed xs in energy bin #17 [2.20000e-01, 2.50000e-01] perturbed xs in energy bin #18 [2.50000e-01, 2.80000e-01] perturbed xs in energy bin #19 [2.80000e-01, 3.00000e-01] perturbed xs in energy bin #20 [3.00000e-01, 3.20000e-01] perturbed xs in energy bin #21 [3.20000e-01, 3.50000e-01] perturbed xs in energy bin #22 [3.50000e-01, 4.00000e-01] perturbed xs in energy bin #23 [4.00000e-01, 5.00000e-01] perturbed xs in energy bin #24 [5.00000e-01, 6.25000e-01] perturbed xs in energy bin #25 [6.25000e-01, 7.80000e-01] perturbed xs in energy bin #26 [7.80000e-01, 8.50000e-01] perturbed xs in energy bin #27 [8.50000e-01, 9.10000e-01] perturbed xs in energy bin #28 [9.10000e-01, 9.50000e-01] perturbed xs in energy bin #29 [9.50000e-01, 9.72000e-01] perturbed xs in energy bin #30 [9.72000e-01, 9.96000e-01] perturbed xs in energy bin #31 [9.96000e-01, 1.02000e+00] perturbed xs in energy bin #32 [1.02000e+00, 1.04500e+00] perturbed xs in energy bin #33 [1.04500e+00, 1.07100e+00] perturbed xs in energy bin #34 [1.07100e+00, 1.09700e+00] perturbed xs in energy bin #35 [1.09700e+00, 1.12300e+00] perturbed xs in energy bin #36 [1.12300e+00, 1.15000e+00] perturbed xs in energy bin #37 [1.15000e+00, 1.30000e+00] perturbed xs in energy bin #38 [1.30000e+00, 1.50000e+00] perturbed xs in energy bin #39 [1.50000e+00, 2.10000e+00] perturbed xs in energy bin #40 [2.10000e+00, 2.60000e+00] perturbed xs in energy bin #41 [2.60000e+00, 3.30000e+00] perturbed xs in energy bin #42 [3.30000e+00, 4.00000e+00] perturbed xs in energy bin #43 [4.00000e+00, 9.87700e+00] perturbed xs in energy bin #44 [9.87700e+00, 1.59680e+01] perturbed xs in energy bin #45 [1.59680e+01, 2.77000e+01] perturbed xs in energy bin #46 [2.77000e+01, 4.80520e+01] perturbed xs in energy bin #47 [4.80520e+01, 7.55010e+01] perturbed xs in energy bin #48 [7.55010e+01, 1.48730e+02] perturbed xs in energy bin #49 [1.48730e+02, 3.67260e+02] perturbed xs in energy bin #50 [3.67260e+02, 9.06900e+02] perturbed xs in energy bin #51 [9.06900e+02, 1.42510e+03] perturbed xs in energy bin #52 [1.42510e+03, 2.23950e+03] perturbed xs in energy bin #53 [2.23950e+03, 3.51910e+03] perturbed xs in energy bin #54 [3.51910e+03, 5.53000e+03] perturbed xs in energy bin #55 [5.53000e+03, 9.11800e+03] perturbed xs in energy bin #56 [9.11800e+03, 1.50300e+04] perturbed xs in energy bin #57 [1.50300e+04, 2.47800e+04] perturbed xs in energy bin #58 [2.47800e+04, 4.08500e+04] perturbed xs in energy bin #59 [4.08500e+04, 6.73400e+04] perturbed xs in energy bin #60 [6.73400e+04, 1.11000e+05] perturbed xs in energy bin #61 [1.11000e+05, 1.83000e+05] perturbed xs in energy bin #62 [1.83000e+05, 3.02500e+05] perturbed xs in energy bin #63 [3.02500e+05, 5.00000e+05] perturbed xs in energy bin #64 [5.00000e+05, 8.21000e+05] perturbed xs in energy bin #65 [8.21000e+05, 1.35300e+06] perturbed xs in energy bin #66 [1.35300e+06, 2.23100e+06] perturbed xs in energy bin #67 [2.23100e+06, 3.67900e+06] perturbed xs in energy bin #68 [3.67900e+06, 6.06550e+06] perturbed xs in energy bin #69 [6.06550e+06, 1.00000e+07]
Entire energy range
fig, ax = plt.subplots(figsize=(7, 3.5), dpi=100)
ax = xs.data[(mat,mt)].plot(logx=True, logy=True, color="dodgerblue", linewidth=2, ax=ax)
for xspert in perturbed_xs:
ax = xspert.data[(mat,mt)].plot(logx=True, logy=True, alpha=0.8, linewidth=0.9, ax=ax)
ax.set_ylabel("cross section / b")
ax.set_xlabel("energy / eV")
ax.set_xlim([1e-5, 2e7])
ax.set_ylim([1e-1, 1e4])
fig.tight_layout();
Thermal part
fig, ax = plt.subplots(figsize=(7, 3.5), dpi=100)
ax = xs.data[(mat,mt)].plot(logx=True, logy=True, color="dodgerblue", linewidth=2, ax=ax)
for xspert in perturbed_xs:
ax = xspert.data[(mat,mt)].plot(logx=True, logy=True, alpha=0.8, linewidth=0.9, ax=ax)
ax.set_ylabel("cross section / b")
ax.set_xlabel("energy / eV")
ax.set_xlim([1e-4, 1])
ax.set_ylim([1e1, 1e4])
fig.tight_layout();
Fast part
fig, ax = plt.subplots(figsize=(7, 3.5), dpi=100)
ax = xs.data[(mat,mt)].plot(logx=True, logy=False, color="dodgerblue", linewidth=2, ax=ax)
for xspert in perturbed_xs:
ax = xspert.data[(mat,mt)].plot(logx=True, logy=False, alpha=0.8, linewidth=0.9, ax=ax)
ax.set_ylabel("cross section / b")
ax.set_xlabel("energy / eV")
ax.set_xlim([1e3, 1e6])
ax.set_ylim([1e0, 1e1])
fig.tight_layout();
Finally, we can write each perturbed cross section into a copy of the original PENDF file for further use.
ipert = 0
perturbed_xs[ipert].to_endf6(tape).to_file(f"copy_pert{ipert}.pendf")