Hello everyone,
I’m currently trying to run an OpenMC program using MGXS to generate multi-group efficient sections. It consists firstly of calculating multi-group cross-sections using a mgxs library, and then using these calculated cross-sections. These are determined for each zone of my system. This latter consists of a ring of material (tungsten) surrounded by a vacuum (divided into two zones (the inner part of the ring and the outer)). I should point out that the outer walls (those of the red part) are of the ‘vacuum’ type.
My goal is to tally the “scatter” reaction rate depending on the scattering angle.
However, strangely, the program stops working (raise RuntimeError(error_msg)) when the number of simulated particles (nb_particles) becomes too large (for example nb_particles=1000000 and nb_batches=10 leads to the bug)… This problem only appears on the multigroup part (the continuous part works fine).
I ran the program in debug mode (as indicated here How to turn on debug). What I obtained is shown in the photo.
I thought the problem must be with the way I was allocating angles in multigroup. So I tested “histogram” and “Legendre” (at different orders) for mgxs_lib.scatter_format, but to no avail. This might also be a problem with the source as I have got a warning "/usr/local/lib/python3.10/dist-packages/openmc/source.py:388: FutureWarning: This class is deprecated in favor of ‘IndependentSource’ "
Do you have any idea why this might be?
Thanks for your help
Best
PS:the program is just after the photos
#les import importants
import openmc
import sys
import os
import matplotlib.pyplot as plt
import numpy as np
import math
def run_multigroupe(energy_groups, liste_SS316, rayon, bord, ordre_remplissage, nb_batches, nb_particles, nb_cutoff):
#beginning of the continuous where cross sections will be computed
if (len(rayon)!=len(ordre_remplissage)):
print("Please make sure rayon and ordre_remplissage are of the same size")
sys.exit()
nom_des_fichiers=["geometry.xml", "settings.xml", "summary.h5", "tallies.out", "tallies.xml"]
for loop in range (len(nom_des_fichiers)):
if os.path.exists(nom_des_fichiers[loop]):
os.remove(nom_des_fichiers[loop])
sphere=[openmc.Sphere(r=rayon[j]) for j in range (len(rayon))]
plan_z_up=openmc.ZPlane(z0=bord, boundary_type = 'vacuum')
plan_z_down=openmc.ZPlane(z0=-bord, boundary_type = 'vacuum')
plan_x_up=openmc.XPlane(x0=bord, boundary_type = 'vacuum')
plan_x_down=openmc.XPlane(x0=-bord, boundary_type = 'vacuum')
plan_y_up=openmc.YPlane(y0=bord, boundary_type = 'vacuum')
plan_y_down=openmc.YPlane(y0=-bord, boundary_type = 'vacuum')
cells=[]
for j in range (len(rayon)):
if (ordre_remplissage[j]!="void"):
cell=openmc.Cell(fill=ordre_remplissage[j])
else:
cell=openmc.Cell()
if (j==0):
region=-sphere[j]
else :
region=-sphere[j] & +sphere[j-1]
cell.region=region
cells.append(cell)
cell_fin=openmc.Cell()
cell_fin.region=+sphere[len(sphere)-1] & -plan_z_up & -plan_x_up & -plan_y_up & +plan_x_down & +plan_y_down & +plan_z_down
cells.append(cell_fin)
universe=openmc.Universe(cells=cells)
plot=universe.plot (width=(900, 900), color_by="cell", colors={cells[0]:"green", cells[1]:"yellow", cells[2]:"red"})
geom=openmc.Geometry(universe)
geom.export_to_xml()
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)
groups = openmc.mgxs.EnergyGroups(energy_groups)
mgxs_lib = openmc.mgxs.Library(geom)
mgxs_lib.energy_groups = groups
mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission',
'nu-scatter matrix', 'multiplicity matrix', 'chi']
mgxs_lib.domain_type = "cell"
mgxs_lib.domains = geom.get_all_cells().values()
mgxs_lib.by_nuclide = False
mgxs_lib.scatter_format ="legendre"
#mgxs_lib.histogram_bins=nombre_angle
mgxs_lib.legendre_order = 1
mgxs_lib.check_library_for_openmc_mgxs()
mgxs_lib.build_library()
tallies=openmc.Tallies()
mgxs_lib.add_to_tallies_file(tallies, merge=True)
tallies.export_to_xml()
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.batches = nb_batches # Number of batches
settings.particles = nb_particles # Number of particles
settings.cutoff = {'energy_neutron': nb_cutoff} # Kill neutrons with energy < 0.001 eV
settings.source = source
settings.export_to_xml()
openmc.run(output="False")
#end of the continuous part#
ce_spfile='calcul_sections_efficaces.h5'
ce_sumfile='summary_sections_efficaces.h5'
os.rename('statepoint.{}.h5'.format(nb_batches), ce_spfile)
os.rename('summary.h5'.format(nb_batches), ce_sumfile)
sp = openmc.StatePoint(ce_spfile, autolink=False)
su = openmc.Summary(ce_sumfile)
sp.link_with_summary(su)
nom_des_fichiers=["geometry.xml", "settings.xml", "summary.h5", "tallies.out", "tallies.xml"]
for loop in range (len(nom_des_fichiers)):
if os.path.exists(nom_des_fichiers[loop]):
os.remove(nom_des_fichiers[loop])
mgxs_lib.load_from_statepoint(sp)
mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=["step{}".format(j) for j in range (len(cells))])
mgxs_file.export_to_hdf5()
materiaux_liste=[]
for loop in range (len(cells)):
materiel=openmc.Material(name=("materiau{}".format(loop)), material_id=10+loop)
materiel.add_macroscopic("step{}".format(loop))
materiaux_liste.append(materiel)
materiaux=openmc.Materials((materiaux_liste[0], materiaux_liste[1], materiaux_liste[2]))
materiaux.cross_sections="mgxs.h5"
materiaux.export_to_xml()
sphere=[openmc.Sphere(r=rayon[j]) for j in range (len(rayon))]
plan_z_up=openmc.ZPlane(z0=bord, boundary_type = 'vacuum')
plan_z_down=openmc.ZPlane(z0=-bord, boundary_type = 'vacuum')
plan_x_up=openmc.XPlane(x0=bord, boundary_type = 'vacuum')
plan_x_down=openmc.XPlane(x0=-bord, boundary_type = 'vacuum')
plan_y_up=openmc.YPlane(y0=bord, boundary_type = 'vacuum')
plan_y_down=openmc.YPlane(y0=-bord, boundary_type = 'vacuum')
cells=[]
for j in range (len(rayon)):
cell=openmc.Cell(fill=materiaux_liste[j])
if (j==0):
region=-sphere[j]
else :
region=-sphere[j] & +sphere[j-1]
cell.region=region
cells.append(cell)
cell_fin=openmc.Cell()
cell_fin.region=+sphere[len(sphere)-1] & -plan_z_up & -plan_x_up & -plan_y_up & +plan_x_down & +plan_y_down & +plan_z_down
cells.append(cell_fin)
universe=openmc.Universe(cells=cells)
#plot=universe.plot (width=(900, 900), color_by="material", colors={FLIBE:"green", SS316:"yellow", tungsten:"red", water:'blue'})
geom=openmc.Geometry(universe)
geom.export_to_xml()
angular_filter=openmc.MuFilter(nombre_angle)
tallies=openmc.Tallies()
tally=openmc.Tally()
tally.name=str(loop)
tally.filters=[angular_filter]
tally.scores=["scatter"]
tallies.append (tally)
tallies.export_to_xml()
source=openmc.Source()
source.space=openmc.stats.Point((0., 0., 0.))
source.angle=openmc.stats.Isotropic()
source.energy = openmc.stats.Uniform(a=energy_groups[len(energy_groups)-2], b=energy_groups[len(energy_groups)-1])
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.energy_mode = 'multi-group'
settings.batches = nb_batches # Number of batches
settings.particles = nb_particles # Number of particles
settings.cutoff = {'energy_neutron': nb_cutoff} # Kill neutrons with energy < 0.001 eV
settings.source = source
settings.export_to_xml()
openmc.run(output="False")
os.rename('statepoint.{}.h5'.format(nb_batches), 'reference.h5')
liste_SS316=[0.07, 2.0, 1.0, 0.05, 0.02, 18.5, 2.5, 13.0, 0] # les fractions en % de la composition du SS316#
#mise en forme des données
liste_SS316=[liste_SS316[j]/100 for j in range (len(liste_SS316))]
liste_SS316.append(1-sum(liste_SS316))
#définition de tous les matériaux#
FLIBE=openmc.Material(name="FLIBE")
FLIBE.add_element("F", 4.)
FLIBE.add_element ("Li", 2.)
FLIBE.add_element ("Be", 1.)
FLIBE.set_density("g/cm3", 1.94)
tungsten=openmc.Material(name="tungsten")
tungsten.add_element ("W", 1.0)
tungsten.set_density("g/cm3", 19.25)
SS316=openmc.Material(name="SS316")
SS316.add_element("C", liste_SS316[0])
SS316.add_element("Mn", liste_SS316[1])
SS316.add_element("Si", liste_SS316[2])
SS316.add_element("P", liste_SS316[3])
SS316.add_element("S", liste_SS316[4])
SS316.add_element("Cr", liste_SS316[5])
SS316.add_element("Mo", liste_SS316[6])
SS316.add_element("Ni", liste_SS316[7])
SS316.add_element("Ti", liste_SS316[8])
SS316.add_element("Fe", liste_SS316[9])
SS316.set_density("g/cm3", 7.99)
water=openmc.Material(name="water")
water.add_element("H", 2.)
water.add_element("O", 1.)
water.set_density("g/cm3", 1.)
water.add_s_alpha_beta('c_H_in_H2O')
materiaux=openmc.Materials([SS316, FLIBE, tungsten, water])
materiaux.export_to_xml()
rayon=[115, 220]
bord=400
ordre_remplissage=["void", tungsten]
nb_batches=10
nb_particles=1000000
nb_cutoff=0.01
nombre_angle=100
energy_groups=[0., 10e6, 14.0e6, 14.07e6]
run_multigroupe(energy_groups, liste_SS316, rayon, bord, ordre_remplissage, nb_batches, nb_particles, nb_cutoff)