From 76d2da8238f062964fa5ae80dd36b6e1f37d65b3 Mon Sep 17 00:00:00 2001 From: GrimFe <69858810+GrimFe@users.noreply.github.com> Date: Thu, 29 Aug 2024 16:27:06 +0200 Subject: [PATCH] Edistr to Endf6 (#338) * fist commit * write_mf5 * Edistr.to_endf6() * bugfix * pfns.py -> edistr.py * dependencies + bugfix * test import fix * testfix write_mf35 --- sandy/__init__.py | 2 +- sandy/{pfns.py => edistr.py} | 82 +++++++++++++++++- sandy/sections/mf5.py | 164 ++++++++++++++++++++++++++++++++--- 3 files changed, 232 insertions(+), 16 deletions(-) rename sandy/{pfns.py => edistr.py} (87%) diff --git a/sandy/__init__.py b/sandy/__init__.py index 3a0f491a..6e395c68 100755 --- a/sandy/__init__.py +++ b/sandy/__init__.py @@ -12,7 +12,7 @@ from .gls import * from .libraries import * from .pert import * -from .pfns import * +from .edistr import * from sandy.zam import * from .njoy import * from .sections import * diff --git a/sandy/pfns.py b/sandy/edistr.py similarity index 87% rename from sandy/pfns.py rename to sandy/edistr.py index 5b85195b..ca630fe9 100644 --- a/sandy/pfns.py +++ b/sandy/edistr.py @@ -9,7 +9,9 @@ import pandas as pd import numpy as np -import sandy +from .pert import Pert +from .core.endf6 import Endf6 +from .sections.mf5 import write_mf5 __author__ = "Luca Fiorito" __all__ = [ @@ -141,7 +143,7 @@ def filter_by(self, key, value): condition = self.data[key] == value out = self.data.copy()[condition].reset_index(drop=True) if out.empty: - raise sandy.Error("applied filter returned empty dataframe") + raise NotImplementedError("applied filter returned empty dataframe") return self.__class__(out) def get_table(self, mat, mt, k): @@ -209,6 +211,8 @@ def add_energy_point(self, mat, mt, k, enew): Examples -------- + + >>> import sandy >>> orig = Edistr(minimal_edistrtest) >>> new = orig.add_energy_point(9437, 18, 0, 1.5) >>> new @@ -282,6 +286,8 @@ def add_energy_points(self, mat, mt, k, elist): Examples -------- + + >>> import sandy >>> orig = Edistr(minimal_edistrtest) >>> new = orig.add_energy_points(9437, 18, 0, [1, 1.5, 1.7]) >>> new @@ -318,6 +324,8 @@ def get_integrals(self): Examples -------- + + >>> import sandy >>> Edistr(minimal_edistrtest).get_integrals() MAT MT K EIN INTEGRAL 0 9437 18 0 1.00000e+00 1.00000e+07 @@ -350,6 +358,8 @@ def normalize(self): Examples -------- + + >>> import sandy >>> new = Edistr(minimal_edistrtest).normalize() >>> new MAT MT K EIN EOUT VALUE @@ -416,6 +426,8 @@ def custom_perturbation(self, pert, mat, mt, k, ein_low, ein_high): Examples -------- + + >>> import sandy >>> orig = Edistr(minimal_edistrtest) >>> pert = sandy.Pert([1.3], index=[1e-3]) >>> orig.custom_perturbation(pert, 9437, 18, 0, 1.5, 2.5) @@ -437,7 +449,7 @@ def custom_perturbation(self, pert, mat, mt, k, ein_low, ein_high): for ein, df in data[condition].groupby("EIN"): series = pert.reshape(df.EOUT).data.loc[df.EOUT] # truncate extremes and replace them with boundaries - px = sandy.Pert(series).truncate() + px = Pert(series).truncate() df.VALUE *= px.data.values dfs.append(df) out = pd.concat(dfs) @@ -525,5 +537,67 @@ def from_endf6(cls, endf6): } data.append(dct) if not data: - raise sandy.Error("no tabulated energy distribution was found") + raise NotImplementedError("no tabulated energy distribution was found") return cls(data) + + def to_endf6(self, endf6): + """ + Update cross sections in :obj:`~sandy.core.endf6.Endf6` instance with those available in a + :obj:`~sandy.edistr.Edistr` instance. + + Parameters + ---------- + `endf6` : :obj:`~sandy.core.endf6.Endf6` + ENDF6 object. + + Returns + ------- + endf6new : :obj:`~sandy.core.endf6.Endf6` + ENDF6 object with updated xs. + + Examples + -------- + + >>> import sandy + >>> import numpy as np + >>> endf6 = sandy.get_endf6_file('jeff_33', 'xs', 922350) + >>> ed = sandy.Edistr.from_endf6(endf6) + + >>> np.testing.assert_equal( + ... sandy.read_mf5(ed.to_endf6(endf6), 9228, 18), + ... sandy.read_mf5(endf6, 9228, 18) + ... ) + + >>> s = ed.data.query("MT==18")['VALUE'].copy() + >>> s.loc[0] = 1234567 + >>> ed.data["VALUE"] = s + >>> newEndf6 = ed.to_endf6(endf6) + + >>> assert newEndf6.data.keys() == endf6.data.keys() + + >>> assert sandy.read_mf5(newEndf6, 9228, 18)["PDISTR"][0]["EIN"][1e-5]["EDISTR"][0] == 1234567 + + >>> np.testing.assert_equal( + ... sandy.read_mf5(newEndf6, 9228, 18)["PDISTR"][0]["EIN"][1e-5]["EDISTR"][1:], + ... sandy.read_mf5(endf6, 9228, 18)["PDISTR"][0]["EIN"][1e-5]["EDISTR"][1:] + ... ) + + >>> np.testing.assert_equal( + ... sandy.read_mf5(newEndf6, 9228, 18)["PDISTR"][0]["EIN"][1e-5]["EOUT"], + ... sandy.read_mf5(endf6, 9228, 18)["PDISTR"][0]["EIN"][1e-5]["EOUT"] + ... ) + """ + data = endf6.data.copy() + mf = 5 + for (mat, mt), df in self.data.groupby(["MAT", "MT"]): + # Must read original section to extract info not given in `Xs` + # instance, e.g. QI, QM + if (mat, mf, mt) not in endf6.data.keys(): + continue + sec = endf6.read_section(mat, mf, mt) + for (k, ein), dd in df.groupby(['K', 'EIN']): + # EOUT rewritten in case points were added to the grid + sec["PDISTR"][k]["EIN"][ein]["EOUT"] = dd["EOUT"].values + sec["PDISTR"][k]["EIN"][ein]["EDISTR"] = dd["VALUE"].values + data[mat, mf, mt] = write_mf5(sec) + return Endf6(data) diff --git a/sandy/sections/mf5.py b/sandy/sections/mf5.py index 6a40ea2d..fb4aa98e 100644 --- a/sandy/sections/mf5.py +++ b/sandy/sections/mf5.py @@ -13,7 +13,7 @@ string. MAT, MF, MT and line numbers are also added (each line ends with a `\n`). """ -import sandy +from ..core.records import read_cont, read_tab1, read_tab2, write_cont, write_eol, write_tab1, write_tab2 __author__ = "Luca Fiorito" __all__ = [ @@ -55,8 +55,8 @@ def read_mf5(tape, mat, mt): "MF": mf, "MT": mt, } - i = 0 - C, i = sandy.read_cont(df, i) + + C, i = read_cont(df, 0) # Number of partial energy distributions. There will be one subsection # for each partial distribution. NK = C.N1 @@ -65,9 +65,10 @@ def read_mf5(tape, mat, mt): "AWR": C.C2, } out.update(add) + pdistr = {} for j in range(NK): - tp, i = sandy.read_tab1(df, i) + tp, i = read_tab1(df, i) # Flag specifying the energy distribution law used for a particular # subsection (partial energy distribution) LF = tp.L2 @@ -81,6 +82,7 @@ def read_mf5(tape, mat, mt): "P": tp.y, "LF": LF, } + # General Evaporation Spectrum (LF=5) if LF == 5: """ @@ -92,16 +94,17 @@ def read_mf5(tape, mat, mt): 92-U-240g.jeff33 """ sub["U"] = tp.C1 - T, i = sandy.read_tab1(df, i) + T, i = read_tab1(df, i) sub["NBT_THETA"] = T.NBT sub["INT_THETA"] = T.INT sub["E_THETA"] = T.x sub["THETA"] = T.y - T, i = sandy.read_tab1(df, i) + T, i = read_tab1(df, i) sub["NBT_G"] = T.NBT sub["INT_G"] = T.INT sub["E_G"] = T.x sub["G"] = T.y + # Simple Maxwellian Fission Spectrum (LF=7) / # Evaporation Spectrum (LF=9) elif LF in (7, 9): @@ -110,24 +113,26 @@ def read_mf5(tape, mat, mt): 27-Co-59g.jeff33 """ sub["U"] = tp.C1 - T, i = sandy.read_tab1(df, i) + T, i = read_tab1(df, i) sub["NBT_THETA"] = T.NBT sub["INT_THETA"] = T.INT sub["E_THETA"] = T.x sub["THETA"] = T.y + # Energy-Dependent Watt Spectrum (LF=11) elif LF == 11: sub["U"] = tp.C1 - T, i = sandy.read_tab1(df, i) + T, i = read_tab1(df, i) sub["NBT_A"] = T.NBT sub["INT_A"] = T.INT sub["E_A"] = T.x sub["A"] = T.y - T, i = sandy.read_tab1(df, i) + T, i = read_tab1(df, i) sub["NBT_B"] = T.NBT sub["INT_B"] = T.INT sub["E_B"] = T.x sub["B"] = T.y + # Energy-Dependent Fission Neutron Spectrum (Madland and Nix) (LF=12) elif LF == 12: TM, i = sandy.read_tab1(df, i) @@ -137,15 +142,16 @@ def read_mf5(tape, mat, mt): sub["INT_TM"] = T.INT sub["E_TM"] = T.x sub["TM"] = T.y + # Arbitrary Tabulated Function (LF=1) elif LF == 1: - T2, i = sandy.read_tab2(df, i) + T2, i = read_tab2(df, i) NZ = T2.NZ # number of incident energies for which distr. is given sub["NBT_EIN"] = T2.NBT sub["INT_EIN"] = T2.INT edistr = {} for k in range(NZ): - T1, i = sandy.read_tab1(df, i) + T1, i = read_tab1(df, i) e_in = T1.C2 edistr[e_in] = { "EOUT": T1.x, @@ -155,6 +161,142 @@ def read_mf5(tape, mat, mt): } sub["EIN"] = edistr pdistr[j] = sub + if pdistr: out["PDISTR"] = pdistr return out + +def write_mf5(sec): + """ + Given the content of a MF% section as nested dictionaries, write it + to string. + + Returns + ------- + `str` + Multiline string reproducing the content of a ENDF-6 section. + + Notes + ----- + .. note:: The end-of-line records MAT, MF, MT and line number are added at + the end of each line. + + .. note:: Implementation of energy distribution laws (LF) needed to + process JEFF-3.3 data + + .. important:: The string does not endf with a newline symbol `\n`. + + + Examples + -------- + + >>> import sandy + >>> import numpy as np + >>> ## LF = 7 + >>> tape = sandy.get_endf6_file('jeff_33', 'xs', 942420) + >>> mat, mf, mt = 9446, 5, 18 + >>> sec = tape.read_section(mat, mf, mt) + >>> np.testing.assert_equal( + ... sandy.read_mf5(sandy.Endf6.from_text(write_mf5(sec)), mat, mt), sec + ... ) + + >>> ## LF = 9 + >>> tape = sandy.get_endf6_file('jeff_33', 'xs', 942420) + >>> mat, mf, mt = 9446, 5, 37 + >>> sec = tape.read_section(mat, mf, mt) + >>> np.testing.assert_equal( + ... sandy.read_mf5(sandy.Endf6.from_text(write_mf5(sec)), mat, mt), sec + ... ) + + >>> ## LF = 5 + >>> tape = sandy.get_endf6_file('jeff_33', 'xs', 942420) + >>> mat, mf, mt = 9446, 5, 18 + >>> sec = tape.read_section(mat, mf, mt) + >>> np.testing.assert_equal( + ... sandy.read_mf5(sandy.Endf6.from_text(write_mf5(sec)), mat, mt), sec + ... ) + + >>> ## LF = 1 + >>> tape = sandy.get_endf6_file('jeff_33', 'xs', 932390) + >>> mat, mf, mt = 9352, 5, 18 + >>> sec = tape.read_section(mat, mf, mt) + >>> np.testing.assert_equal( + ... sandy.read_mf5(sandy.Endf6.from_text(write_mf5(sec)), mat, mt), sec + ... ) + """ + lines = write_cont( + sec["ZA"], sec["AWR"], 0, 0, len(sec["PDISTR"].keys()), 0 + ) + + for data in sec["PDISTR"].values(): + + if data["LF"] in (5, 7, 9): + lines += write_tab1( + data["U"], + 0, + 0, + data["LF"], + data["NBT_P"], + data["INT_P"], + data["E_P"], + data["P"] + ) + lines += write_tab1( + 0, + 0, + 0, + 0, + data["NBT_THETA"], + data["INT_THETA"], + data["E_THETA"], + data["THETA"] + ) + + if data["LF"] == 5: + lines += write_tab1( + 0, + 0, + 0, + 0, + data["NBT_G"], + data["INT_G"], + data["E_G"], + data["G"] + ) + + elif data["LF"] == 1: + lines += write_tab1( + 0, + 0, + 0, + data["LF"], + data["NBT_P"], + data["INT_P"], + data["E_P"], + data["P"] + ) + lines += write_tab2( + 0, + 0, + 0, + 0, + len(data["EIN"].keys()), + data["NBT_EIN"], + data["INT_EIN"] + ) + + for ein, data_ in data["EIN"].items(): + lines += write_tab1( + 0, + ein, + 0, + 0, + data_['NBT'], + data_["INT"], + data_["EOUT"], + data_["EDISTR"] + ) + + else: + raise NotImplementedError(f"Writing LF = {data['LF']} not implemented yet.") + return '\n'.join(write_eol(lines, sec["MAT"], mf, sec["MT"]))