Hi everyone,
I am posting on this group because I am currently working on OpenMC Multigroup. I am studying multigroup and continuous energy discrepancies in a tungsten wall. More specifically, I am focusing on the values of the reaction rates (absorption and scatter). I have tried all the different scattering formats available (histogram and Legendre) and found a certain behavior with histogram.
Indeed, with this particular model, the MG scatter reaction rate seems to be weighted by a factor (see figure), i.e. R_(scatter,MG, g)/R_(scatter,CE, g)=constant. This constant is approximately the same for each energy group. I should add that this is not true for Legendre (see figure) nor for absorption reaction rates (regardless of scatter format). I also tried replacing the material with carbon and iron and found the same behavior.
I tried changing the number of bins in histograms and found that the median and mean of R_(scatter,MG, g)/R_(scatter,CE, g) over all energy groups increase (almost linearly) with the number of bins (see figure). Is this normal? Does this excessive MG scatter reaction rate also interfere with the overall MG results ? I am attaching the code I used to plot the first two figures, maybe I am making a mistake in it…
One last question: Are the MG reaction rates also tallied using the Track-length Estimator technique (as the continuous energy results) ?
I appreciate your help
Best
import numpy
import openmc
import os
from matplotlib import pyplot as plt
#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)
cell_filter=openmc.CellFilter(tungsten_cell)
tallies=openmc.Tallies()
tally=openmc.Tally()
tally.name="1"
tally.filters=[energy_filter, cell_filter]
tally.scores=["scatter", "absorption"]
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
#definition of the scatter format (chose either histogram or legendre)
mgxs_lib.scatter_format ="histogram"
mgxs_lib.histogram_bins=1000
#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)
cell_filter=openmc.CellFilter(tungsten_cell)
tallies=openmc.Tallies()
tally=openmc.Tally()
tally.name="1"
tally.filters=[energy_filter, cell_filter]
tally.scores=["scatter", "absorption"]
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")
res={}
dres={}
tally=sp.get_tally(name="1")
mean=tally.get_reshaped_data("mean")
error=tally.get_reshaped_data("std_dev")
reaction={0:"scatter", 1:"absorption"}
dictionary_continuous={}
for i in range(2):
ans=[]
for j in range (len(energy_groups)-1):
ans.append(mean[j][0][0][i])
dictionary_continuous[reaction[i]]=ans
#reading of the multigroup run values
sp = openmc.StatePoint("multigroup.h5")
res={}
dres={}
tally=sp.get_tally(name="1")
mean=tally.get_reshaped_data("mean")
error=tally.get_reshaped_data("std_dev")
reaction={0:"scatter", 1:"absorption"}
dictionary_multigroup={}
ddictionary_multigroup={}
for i in range(2):
ans=[]
dans=[]
for j in range (len(energy_groups)-1):
ans.append(mean[j][0][0][i])
dans.append(error[j][0][0][i])
dictionary_multigroup[reaction[i]]=ans
ddictionary_multigroup[reaction[i]]=dans
#plot part
fig, ax=plt.subplots(1, 2, figsize=(10, 6))
fig.suptitle("Comparison of the scatter and absorption reaction rates with a continuous and multigroup run \n for Legendre 10")
ax[0].step (energy_groups, [0]+dictionary_continuous["scatter"], label="Continuous")
ax[0].step (energy_groups, [0]+dictionary_multigroup["scatter"], label="Multigroup")
ax[0].set_xscale ("log")
ax[0].set_yscale ("log")
ax[0].set_xlim(energy_groups[0], energy_groups[-1])
ax[0].set_title("Scatter")
ax[0].legend()
ax[1].step (energy_groups, [0]+dictionary_continuous["absorption"], label="Continuous")
ax[1].step (energy_groups, [0]+dictionary_multigroup["absorption"], label="Multigroup")
ax[1].set_xscale ("log")
ax[1].set_yscale ("log")
ax[1].set_xlim(energy_groups[0], energy_groups[-1])
ax[1].set_title("Absorption")
ax[1].legend()
fig.tight_layout()