Hello everyone,
I am posting to this group because I am currently working on OpenMC Multigroup. I am studying multigroup and continuous energy discrepancies in a tungsten wall. In particular, I am focusing on the energy distribution of the outgoing particle after a scattering reaction.
So I have computed “scatter” and “nu-scatter” with an EnergyFilter and an EnergyoutFilter for both continuous and multigroup modes. But as the “nu-scatter” and “scatter” results are different when running in continuous energy, they are strictly the same in multigroup (see Figure). Is this normal? I have to note that these multigroup results are different when running without an EnergyoutFilter (which is coherent).
I am expecting a difference, such as the “nu-scatter” reaction rate from group g to g’ is equal to the “scatter” reaction rate from group g to g’ multiplied by the the coefficient of the multiplicity matrix associated to this energy transfer.
I’d appreciate your help
Best regards
PS: I am attaching the code that allows me to plot the figure
import numpy
import openmc
import os
from matplotlib import pyplot as plt
import matplotlib.patches as mpatches
openmc.config["cross_sections"]="/endfb-vii.1-hdf5/cross_sections.xml"
#definition of all the mandatory files
#definition of the materials
tungsten=openmc.Material(name="tungsten", material_id=0)
tungsten.add_element ("W", 1.0)
tungsten.set_density("g/cm3", 19.25)
materiaux=openmc.Materials([tungsten])
materiaux.export_to_xml()
#definition of the energy groups
vitamin=openmc.mgxs.GROUP_STRUCTURES["VITAMIN-J-42"]
energy_groups=[]
for loop in range (len(vitamin)):
if not((vitamin[loop]>14.06e6) or (vitamin[loop]<0.01)):
energy_groups.append(vitamin[loop])
energy_groups.append(14.06e6)
energy_groups.append(14.07e6)
if (energy_groups[0]!=0.01):
energy_groups=[0.01]+energy_groups
energy_groups.sort()
#definition of the cells
border=400
plane_z_up=openmc.ZPlane(z0=border, boundary_type = 'vacuum')
plane_z_down=openmc.ZPlane(z0=-border, boundary_type = 'vacuum')
plane_x_up=openmc.XPlane(x0=border, boundary_type = 'vacuum')
plane_x_down=openmc.XPlane(x0=-border, boundary_type = 'vacuum')
plane_y_up=openmc.YPlane(y0=border, boundary_type = 'vacuum')
plane_y_down=openmc.YPlane(y0=-border, boundary_type = 'vacuum')
inner_wall=openmc.XPlane(115)
outer_wall=openmc.XPlane(117)
tungsten_cell=openmc.Cell()
tungsten_cell.id=0
tungsten_cell.fill=tungsten
tungsten_cell.region=+inner_wall & -outer_wall & -plane_z_up & -plane_y_up & +plane_y_down & +plane_z_down
void_cell_left=openmc.Cell()
void_cell_left.region=-inner_wall & +plane_x_down & -plane_z_up & -plane_y_up & +plane_y_down & +plane_z_down
void_cell_right=openmc.Cell()
void_cell_right.region=+outer_wall & -plane_x_up & -plane_z_up & -plane_y_up & +plane_y_down & +plane_z_down
#definition of the universe
universe=openmc.Universe(cells=[tungsten_cell, void_cell_left, void_cell_right])
root_cell=openmc.Cell(name="root cell")
root_cell.fill=universe
root_cell.region=-plane_z_up & -plane_x_up & -plane_y_up & +plane_x_down & +plane_y_down & +plane_z_down
root_universe=openmc.Universe(name="root universe", universe_id=0)
root_universe.add_cell(root_cell)
root_universe.plot(origin=(0., 0., 0.), width=(2*border, 2*border), pixels=(100, 100), color_by='material')
root_universe.plot(origin=(0., 0., 0.), width=(2*border, 2*border), pixels=(100, 100), color_by='cell')
geom=openmc.Geometry(root_universe)
geom.export_to_xml()
#definition of the source
source=openmc.Source()
source.space=openmc.stats.Point((0., 0., 0.))
source.angle=openmc.stats.Isotropic()
source.energy = openmc.stats.Uniform(a=14.06e6, b=14.07e6)
source.particle="neutron"
############################################################################
#continuous energy run
#settings
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.energy_mode='continuous-energy'
settings.batches = 10 # Number of batches
settings.particles = 1000000 # Number of particles
settings.cutoff = {'energy_neutron': 0.01} # Kill neutrons with energy < 0.001 eV
settings.source = source
settings.export_to_xml()
#tallies
energy_filter=openmc.EnergyFilter(energy_groups)
energy_out_filter=openmc.EnergyoutFilter(energy_groups)
cell_filter=openmc.CellFilter(tungsten_cell)
tallies=openmc.Tallies()
tally=openmc.Tally()
tally.name="1"
tally.filters=[energy_filter, energy_out_filter, cell_filter]
tally.scores=["scatter", "nu-scatter"]
tallies.append(tally)
tallies.export_to_xml()
#openmc.run
if os.path.exists("summary.h5"):
os.remove("summary.h5")
if os.path.exists('statepoint.{}.h5'.format(10)):
os.remove('statepoint.{}.h5'.format(10))
openmc.run()
os.rename('statepoint.{}.h5'.format(10), "continuous.h5")
#########################################################################################
#multigroup energy run
#computation of the multigroup cross sections
groups = openmc.mgxs.EnergyGroups(energy_groups)
mgxs_lib = openmc.mgxs.Library(geom)
mgxs_lib.energy_groups = groups
types=['total', 'absorption', 'nu-fission', 'fission',
'nu-scatter matrix', 'multiplicity matrix', 'chi', 'nu-scatter']
mgxs_lib.mgxs_types = types
mgxs_lib.domain_type = "material"
mgxs_lib.domains = geom.get_all_materials().values()
mgxs_lib.by_nuclide = False
mgxs_lib.scatter_format ="legendre"
mgxs_lib.legendre_order=10
mgxs_lib.check_library_for_openmc_mgxs()
mgxs_lib.build_library()
tallies_file = openmc.Tallies()
mgxs_lib.add_to_tallies_file(tallies_file, merge=True)
tallies_file.export_to_xml()
if os.path.exists("summary.h5"):
os.remove("summary.h5")
if os.path.exists('statepoint.{}.h5'.format(10)):
os.remove('statepoint.{}.h5'.format(10))
openmc.run()
os.rename('statepoint.{}.h5'.format(10), "statepoint_ce.h5")
#creation of the multigroup mgxs library
ce_spfile = './statepoint_ce.h5'
if (os.path.exists('statepoint.{}.h5'.format(10))):
os.rename('statepoint.{}.h5'.format(10), ce_spfile)
ce_sumfile = './summary_ce.h5'
if (os.path.exists('summary.h5')):
os.rename('summary.h5', ce_sumfile)
sp = openmc.StatePoint(ce_spfile, autolink=False)
su = openmc.Summary(ce_sumfile)
sp.link_with_summary(su)
mgxs_lib.load_from_statepoint(sp)
mgxs_file = mgxs_lib.create_mg_library(xs_type='macro')
if (os.path.exists('mgxs_multigroupe.h5')):
os.remove('mgxs_multigroupe.h5')
mgxs_file.export_to_hdf5('mgxs_multigroupe.h5')
#definition of the multigroup run settings
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.energy_mode = 'multi-group'
settings.batches = 10 # Number of batches
settings.particles = 1000000 # Number of particles
settings.cutoff = {'energy_neutron': 0.01} # Kill neutrons with energy < 0.001 eV
settings.source = source
settings.export_to_xml()
#use of multigroup in materials.xml
if (os.path.exists("materials.xml")):
os.remove("materials.xml")
tungsten_mg=openmc.Material(material_id=0)
tungsten_mg.add_macroscopic("set1")
materials_file=openmc.Materials([tungsten_mg])
materials_file.cross_sections = 'mgxs_multigroupe.h5'
materials_file.export_to_xml()
#tallies
energy_filter=openmc.EnergyFilter(energy_groups)
energy_out_filter=openmc.EnergyoutFilter(energy_groups)
cell_filter=openmc.CellFilter(tungsten_cell)
tallies=openmc.Tallies()
tally=openmc.Tally()
tally.name="1"
tally.filters=[energy_filter, energy_out_filter, cell_filter]
tally.scores=["scatter", "nu-scatter"]
tallies.append(tally)
tallies.export_to_xml()
#multigroup run
openmc.run()
os.rename('statepoint.{}.h5'.format(10), "multigroup.h5")
######################################################################
#ploting
#reading the continuous run values
sp = openmc.StatePoint("continuous.h5")
ce={}
tally=sp.get_tally(name="1")
mean=tally.get_reshaped_data("mean")
reaction={0:"scatter", 1:"nu-scatter"}
for i in range (2):
liste=[[0 for m in range (len(energy_groups)-1)] for n in range (len(energy_groups)-1)]
for m in range (len(energy_groups)-1):
for n in range (len(energy_groups)-1):
liste[m][n]=mean[m][n][0][0][i]
ce[reaction[i]]=liste
#reading of the multigroup run values
sp = openmc.StatePoint("multigroup.h5")
mg={}
tally=sp.get_tally(name="1")
mean=tally.get_reshaped_data("mean")
reaction={0:"scatter", 1:"nu-scatter"}
for i in range (2):
liste=[[0 for m in range (len(energy_groups)-1)] for n in range (len(energy_groups)-1)]
for m in range (len(energy_groups)-1):
for n in range (len(energy_groups)-1):
liste[m][n]=mean[m][n][0][0][i]
mg[reaction[i]]=liste
#plot part
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 6))
energy_group_origin=0
ax1.step(energy_groups, [0]+ce["scatter"][-1-energy_group_origin], label="Scatter", color="red")
ax1.step(energy_groups, [0]+ce["nu-scatter"][-1-energy_group_origin], label="Nu-scatter", color="blue")
ax1.set_title("Continuous energy")
ax1.set_ylabel('Value')
ax1.set_xlabel('Neutron energy (in eV)')
ax1.set_yscale("log")
ax1.set_xscale("log")
ax1.set_xlim(energy_groups[0], energy_groups[-1])
ax2.step(energy_groups, [0]+mg["scatter"][-1-energy_group_origin], label="Scatter", color="red")
ax2.step(energy_groups, [0]+mg["nu-scatter"][-1-energy_group_origin], label="Nu-scatter", color="blue")
ax2.set_title("Multigroup energy")
ax2.set_ylabel('Value')
ax2.set_xlabel('Neutron energy (in eV)')
ax2.set_yscale("log")
ax2.set_xscale("log")
ax2.set_xlim(energy_groups[0], energy_groups[-1])
legend_labels = {
'Curve 1': mpatches.Patch(color='red', label="Scatter"),
'Curve 2': mpatches.Patch(color='blue', label="Nu-scatter"),
}
custom_legend =fig.legend(loc='upper left', handles=legend_labels.values(), bbox_to_anchor=(0.40, 0), fontsize='x-large')
plt.suptitle('Differences between nu-scatter and scatter reaction rates \n when using continuous or multigroup. \n The plots show the reaction rates from energy group {} to all the other groups'.format(energy_group_origin), fontsize=16, x=0.6)
plt.tight_layout()
plt.show()