Hi all,
I’ve been trying to use to produce multigroup cross sections using the MGXS module, and compare them to cross-sections obtained from tallying the multigroup flux and running GROUPR.
As a minimal example, I am considering a sphere of NaClPuCl3 with a radius of r=35cm and a spherical reflector of SiC with a radius of r=200cm. I am trying to compute the multigroup cross-section of the (n,p) reaction of Cl35 using the ECCO-1968, ECCO-33 and 1G energy meshes. I am using JEFF-3.3 cross-sections as provided on the OpenMC website.
Here is the complete script for running the case:
import openmc
import numpy as np
import pandas as pd
import argparse
openmc.config["cross_sections"] = "/data/nuclear_data/endfb8/HDF5/cross_sections.xml"
ecco33 = np.array(
[
1.000010e-05,
1.000000e-01,
5.400000e-01,
4.000000e00,
8.315287e00,
1.370959e01,
2.260329e01,
4.016900e01,
6.790405e01,
9.166088e01,
1.486254e02,
3.043248e02,
4.539993e02,
7.485183e02,
1.234098e03,
2.034684e03,
3.354626e03,
5.530844e03,
9.118820e03,
1.503439e04,
2.478752e04,
4.086771e04,
6.737947e04,
1.110900e05,
1.831564e05,
3.019738e05,
4.978707e05,
8.208500e05,
1.353353e06,
2.231302e06,
3.678794e06,
6.065307e06,
1.000000e07,
1.964033e07,
]
)
GROUP_STRUCTURES = { "ECCO-1968": openmc.mgxs.GROUP_STRUCTURES["ECCO-1968"],
"ECCO-33": ecco33,
"1G": np.array([1e-5, 2e7])
}
def build_model():
pu = openmc.Material(name="Pu")
pu.add_components(
{"Pu238": 1, "Pu239": 50, "Pu240": 35, "Pu241": 7, "Pu242": 7},
percent_type="wo",
)
cl = openmc.Material(name="Cl")
cl.add_element("Cl", percent=1)
na = openmc.Material(name="Na")
na.add_element("Na", percent=1)
fuel = openmc.Material.mix_materials(
[pu, cl, na], [1 / 6, 4 / 6, 1/6], percent_type="ao")
fuel.set_density("g/cm3", 3)
sic = openmc.Material()
sic.add_elements_from_formula("SiC")
boron = openmc.Material()
boron.add_element("B", percent=1)
ppm = 1e-6
refl = openmc.Material.mix_materials([sic, boron], [1 - 1500*ppm, 1500*ppm])
refl.set_density("g/cm3", 3)
water = openmc.Material()
water.add_elements_from_formula("H2O")
water.set_density("g/cm3", 1)
materials = openmc.Materials([fuel, refl])
sphere = openmc.Sphere(r=35, boundary_type="transmission")
reflector = openmc.Sphere(r=200, boundary_type="vacuum")
cells = [
openmc.Cell(name="Fuel", region=-sphere, fill=fuel),
openmc.Cell(name="Reflector", region=-reflector & +sphere, fill=refl),
]
root = openmc.Universe(cells=cells)
geometry = openmc.Geometry(root=root)
source = openmc.IndependentSource(space=openmc.stats.Point((25, 0, 0)))
settings = openmc.Settings()
settings.source = source
settings.batches = 100
settings.inactive = 10
settings.particles = 1_000_000
settings.ptables = False
settings.output = {'tallies': False}
settings.temperature = {"method": "interpolation"}
settings.photon_transport = True
plot = openmc.Plot()
plot.basis = "yz"
plot.pixels = (1000, 1000)
plot.origin = (0., 0., 0.)
plot.width = (100, 100)
plot.color_by = "material"
plots = openmc.Plots(plots=[plot])
tallies = openmc.Tallies()
for gname, g in GROUP_STRUCTURES.items():
tally = openmc.Tally(name=f"Spectrum {gname}")
tally.scores = ["flux"]
tally.filters = [
openmc.MaterialFilter([fuel]),
openmc.EnergyFilter(GROUP_STRUCTURES[gname]),
]
tallies.append(tally)
xs = {
"ECCO-1968": {},
"ECCO-33": {},
"1G": {}
}
for gname, g in GROUP_STRUCTURES.items():
groups = openmc.mgxs.EnergyGroups(group_edges=g)
for r in ["(n,gamma)", "(n,p)"]:
xs[gname][r] = openmc.mgxs.ArbitraryXS(
r,
energy_groups=groups,
domain=fuel,
by_nuclide=True,
)
xs[gname]["fission"] = openmc.mgxs.FissionXS(
energy_groups=groups,
domain=fuel,
by_nuclide=True,
)
for g in xs:
for r in xs[g]:
tallies+= xs[g][r].tallies.values()
model = openmc.Model(
materials=materials, geometry=geometry, settings=settings, tallies=tallies, plots=plots
)
return model, xs
parser = argparse.ArgumentParser(prog='SphericalModel')
parser.add_argument("resdir")
parser.add_argument("--photons", "-p", action="store_true")
if __name__ == "__main__":
args = parser.parse_args()
model, xs = build_model()
model.plot_geometry(cwd=args.resdir)
model.settings.photon_transport = args.photons
model.run(cwd=args.resdir,
mpi_args=["mpirun", "-n", "4"],
export_model_xml=False)
batches = model.settings.batches
with openmc.StatePoint(f"{args.resdir}/statepoint.{batches}.h5") as sp:
for gname, energies in GROUP_STRUCTURES.items():
tally = sp.get_tally(name=f"Spectrum {gname}")
avg = tally.mean.ravel()
unc = tally.std_dev.ravel()
spectrum_data = {
"Emin": np.array([]),
"Emax": np.array([]),
"Flux": np.array([]),
"Flux std. dev.": np.array([]),
"Flux/du": np.array([]),
"Flux/du std. dev.": np.array([]),
}
lethargy = -np.log(20e6 / energies)
spectrum_data["Emin"] = energies[:-1]
spectrum_data["Emax"] = energies[1:]
spectrum_data["Flux"] = avg
spectrum_data["Flux std. dev."] = unc
spectrum_data["Flux/du"] = avg / np.diff(lethargy)
spectrum_data["Flux/du std. dev."] = unc / np.diff(lethargy)
df = pd.DataFrame(data=spectrum_data)
df.to_csv(f"{args.resdir}/spectrum-{gname}.csv", index=False)
for groupstruct in GROUP_STRUCTURES:
for r in ["(n,gamma)", "(n,p)", "fission"]:
xs[groupstruct][r].load_from_statepoint(sp)
xs[groupstruct][r].build_hdf5_store(directory=args.resdir,
filename=f"{groupstruct}.h5",
xs_type="micro")
Here is the script for running NJOY and processing the results in a pdf:
import pandas as pd
import sandy
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
import h5py
import openmc
import openmc.mgxs
from pathlib import Path
name_to_energy = {
"ECCO-33": sandy.energy_grids.ECCO33,
"ECCO-1968": openmc.mgxs.GROUP_STRUCTURES["ECCO-1968"],
"1G": np.array([1e-5, 2e7]),
}
ignnames = {
2: "CSEWG-239",
19: "ECCO-33",
20: "ECCO-1968",
}
MTnames = {
102: "(n,gamma)",
103: "(n,p)",
"(n,gamma)": 102,
"(n,p)": 103,
}
def gendf(results="results"):
df = pd.read_csv(f"{results}/spectrum-ECCO-1968.csv")
df["DE"] = df["Emax"] - df["Emin"]
tapes = {}
tape = sandy.get_endf6_file("jeff_33", "xs", 170350)
spectrum = np.array([df["Emin"], df["Flux"]/df["DE"]]).T.flatten()
tapes["1G"] = tape.get_gendf(temperature=294, groupr_kws={"spectrum": spectrum, "ek": [1e-5, 2e7]})
Path(f"{results}/gendf/1G/").mkdir(parents=True, exist_ok=True)
tapes["1G"].to_file(f"{results}/gendf/1G/Cl35.gendf")
for ign in [2, 19, 20]:
tapes[ign] = tape.get_gendf(temperature=294, groupr_kws={"spectrum": spectrum, "ign": ign})
Path(f"{results}/gendf/{ignnames[ign]}/").mkdir(parents=True, exist_ok=True)
tapes[ign].to_file(f"{results}/gendf/{ignnames[ign]}/Cl35.gendf")
def fullplot(gname, mt, xlim=(8e-6, 1e8), results="results"):
mtname = MTnames[mt]
fig, axes = plt.subplots(2, 1, figsize=(11.69,8.27), height_ratios=[3, 1], sharex=True)
ax = axes[0]
with h5py.File("/data/nuclear_data/jeff33/HDF5/Cl35.h5") as f:
energy = f["Cl35/energy/294K"][...]
xs = f[f"Cl35/reactions/reaction_{mt}/294K/xs"][...]
ax.plot(energy, xs, label="Continuous")
with h5py.File(f"{results}/{gname}.h5") as f:
xs = f[f"material/4/{mtname}/Cl35/average"][...][::-1]
std = f[f"material/4/{mtname}/Cl35/std. dev."][...][::-1]
energy = name_to_energy[gname]
if len(xs) == 1:
mgxs_toplot = [xs[0], xs[0]]
ax.plot(energy, mgxs_toplot, drawstyle="steps-post", label=f"{gname} MGXS")
else:
mgxs_toplot = xs
ax.plot(energy[:-1], xs, drawstyle="steps-post", label=f"{gname} MGXS")
ax.fill_between(energy[:-1], xs - std, xs + std, step="post", color="tab:orange", alpha=0.3)
gendf = sandy.gendf.Gendf.from_file(f"{results}/gendf/{gname}/Cl35.gendf")
xs = gendf.get_xs(mt=[mt]).data
Eleft = np.array([i.left for i in xs.index])
Eright = np.array([i.right for i in xs.index])
njoy_toplot = xs[1725, mt]
if len(njoy_toplot) == 1:
njoy_toplot = [njoy_toplot.values[0], njoy_toplot.values[0]]
ax.plot([Eleft, Eright], njoy_toplot, drawstyle="steps-post", label=f"{gname} NJOY")
else:
ax.plot(Eleft, njoy_toplot, drawstyle="steps-post", label=f"{gname} NJOY")
njoy_toplot = np.array(njoy_toplot)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(*xlim)
# ax.set_ylim(1e-5, 1)
ax.set_ylabel("Cross-section [barn]")
ax.legend(loc="lower left")
ax1 = ax.twinx()
df = pd.read_csv(f"{results}/spectrum-ECCO-1968.csv")
df["DE"] = df["Emax"] - df["Emin"]
flux = np.array(df["Flux"]/df["DE"])
ax1.plot(df["Emin"], flux, ds="steps-post", color="tab:red", label="Spectrum")
top = (df["Flux"] + 3*df["Flux std. dev."])/df["DE"]
bot = (df["Flux"] - 3*df["Flux std. dev."])/df["DE"]
ax1.fill_between(df["Emin"], top, bot, step="post", color="tab:red", alpha=0.3)
ax1.legend()
ax1.set_ylim(0.9*np.min(flux), 1.1*np.max(flux))
ax1.set_yscale("log")
ax1.set_ylabel(r"Flux [$\mathrm{neutron}\cdot \mathrm{cm} \cdot \mathrm{cm}^{-3} \cdot \mathrm{s}^{-1} \cdot \mathrm{sp}^{-1}$]")
ax.grid(which="both")
ax.set_title(f"{mtname} cross-section of Cl35")
ax = axes[1]
ax.set_xlabel("Energy [eV]")
ax.set_ylabel("Relative difference [%]")
xs = np.array(xs[1725, 103])
val = abs(njoy_toplot-mgxs_toplot)/mgxs_toplot * 100
if len(Eleft) == 1:
ax.plot([Eleft, Eright], val, color="k")
else:
ax.plot(Eleft, val, ds="steps-post", color="k")
ax.set_xscale("log")
ax.grid()
fig.tight_layout()
return fig, axes
if __name__=="__main__":
for results in ["photons", "nophotons"]:
gendf(results)
with PdfPages(f"{results}.pdf") as pdf:
fig, ax = fullplot("ECCO-1968", mt=103, results=results)
pdf.savefig(fig)
fig, ax = fullplot("ECCO-1968", mt=103, xlim=(1e2, 1e6), results=results)
pdf.savefig(fig)
fig, ax = fullplot("ECCO-1968", mt=103, xlim=(1.5e4, 2e4), results=results)
pdf.savefig(fig)
fig, ax = fullplot("ECCO-33", mt=103, results=results)
pdf.savefig(fig)
fig, ax = fullplot("1G", mt=103, results=results)
pdf.savefig(fig)
Without photon transport
When photon transport is disabled, I get the following result:
Beside low energies when the tally uncertainty is too large, and some resonnances which are not correctly meshed by the energy discretization, the cross-sections produced by MGXS and by NJOY are in good agreement. The one group cross-section (in which I am most interested for back-of-the-napkin calculations) show less than a percent of discrepancy:
I find this quite satisfactory!
With photon transport
When enabling coupled photon-neutron transport however I get the following:
For energies > 1e4 eV, the discrepancies between MGXS and NJOY are much larger, and we can also see additionnal “positive” resonnances in the flux. The 1 group calculation results in a relative difference of 70% percent between the two methods making the values somewhat dubious.
I was under the impression that photon transport in OpenMC was mostly aimed at gamma source calculations and accurate heat deposition, and was surprised it had such a large effect on the neutron flux and the MGXS results. Does anyone have an explanation for this?
Best regards,
Nicolas
