Cell::set_temperature causes exception when trying to set a valid temperature in multigroup mode with temperature interpolation policy

lewisgross1296 commented 1 year ago

During a simulation with settings.energy_mode = 'multi-group' and settings.temperature[method]=interpolation, the following exception occurs when trying to set the temperature of a cell to a temperature that is both within the bounds of settings.temperature[range] and within the bounds used to create the multi-group mode cross section library.

openmc.exceptions.OpenMCError: Temperature is below minimum temperature at which data is available.

Below is a MWE python script that reproduces the error

import openmc
import openmc.lib
import openmc.mgxs as mgxs
import numpy as np
import h5py

Tmin = 293
T0 = 303 # temp in range that should be set-able
Tmax = 313

# generate one group cross section data
groups = mgxs.EnergyGroups()
groups.group_edges = np.array([0.0, 20.0e6]) #  groups in eV

# create temperature data for cross sections
temps = [float(T) for T in range(Tmin,Tmax+1)]
xsdata = openmc.XSdata('slab_xs', energy_groups=groups, temperatures=temps, num_delayed_groups=0)
xsdata.order = 0

# create macroscopic object and export material with temperature XS library generated below
slab = openmc.Material(1, "slab")
slab.set_density('macro',1.)
slab.add_macroscopic('slab_xs')

materials = openmc.Materials([slab])
materials.cross_sections = 'one_gxs.h5'
materials.export_to_xml()

s = 0.45
f = 1.5
nu = f/(1-s)
for T in range(Tmin,Tmax+1):
    Sig_t = 1/T
    Sig_s = s*Sig_t
    nu_Sig_f = f*Sig_t
    Sig_f = nu_Sig_f / nu
    Sig_a = Sig_f # no non fission absorptions
    xsdata.set_total(np.array([Sig_t]),temperature=T)
    xsdata.set_scatter_matrix(np.array([[[Sig_s]]]),temperature=T)
    xsdata.set_absorption(np.array([Sig_a]),temperature=T)
    xsdata.set_fission(np.array([Sig_f]),temperature=T)
    xsdata.set_nu_fission(np.array([nu_Sig_f]),temperature=T)

# export xsdata
one_g_XS_file = openmc.MGXSLibrary(groups) # initialize the library
one_g_XS_file.add_xsdata(xsdata) # add benchmark XS data
one_g_XS_file.export_to_hdf5('one_gxs.h5') # export library as hdf5 format

# define geometry
H = 100.0 # cm
xplane0 = openmc.XPlane(x0  = 0.0, boundary_type='vacuum')
xplane1 = openmc.XPlane(x0 = H, boundary_type='vacuum')
yplane0 = openmc.YPlane(y0 = 0.0, boundary_type='vacuum')
yplane1 = openmc.YPlane(y0 = H, boundary_type='vacuum')
zplane0 = openmc.ZPlane(z0 = 0.0, boundary_type='vacuum')
zplane1 = openmc.ZPlane(z0 = H, boundary_type='vacuum')
cube_region = +xplane0 & -xplane1 & +yplane0 & -yplane1 & +zplane0 & -zplane1
root_cell = openmc.Cell(cell_id=0, region=cube_region, fill=slab)
geometry = openmc.Geometry([root_cell])
geometry.export_to_xml()

# define settings
settings = openmc.Settings()
batches = 200
inactive = 50
particles = 5000
settings.energy_mode = 'multi-group'
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.temperature = {'default': T0,
                        'method': 'interpolation',
                        'range': (Tmin, Tmax)}
bounds = [0, 0, 0, H, H, H]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.Source(space=uniform_dist)
settings.export_to_xml()

# initialize an openmc simulation
openmc.lib.init()
# C-API call to set (valid) temperature of cell that results in error
openmc.lib.Cell.set_temperature(self=openmc.lib.cells[0],T=T0,set_contained=True)
# openmc.lib.run()
# openmc.lib.finalize()

lewisgross1296 commented 1 year ago

@aprilnovak @gonuke FYI

gonuke commented 1 year ago

Can you try this again with a multi-group data file that is already in the repo to confirm that this is not related to the definition of the data file? It is probably the same, but this will help remove speculation from the data generation.

lewisgross1296 commented 1 year ago

I'm not sure there are any specifically multi-group libraries provided by https://openmc.org/official-data-libraries/. It seems like the documentation recommends generating MG libraries from the continuous energy-libraries available. Here is the intro to MG mode using OpenMC's library generation capabilities, where they show how to create your own MG library from continuous energy data using any user-defined group structure.

I think following this notebook as an example and generating an MG library from a continuous one, for say one isotope with one group would be sufficient to remove the possibility that my custom defined data library is causing this issue/error.

lewisgross1296 commented 1 year ago

I was able to reproduce the error with multi-group data generated from a continuous-energy simulation using utility from OpenMC's mgxs.Library, XSdata, and MGXSLibrary classes. Here are the scripts below. It is recommended to have a subdirectory, e.g. generator at the same level as the MWE script, since the continuous energy simulation will need XML input that is separate from the MWE script's XML input to generate the error. An example directory structure would be

test/MWE.py
test/generator/mgxs_generator.py

In order to reproduce the error, you would first run python mgxs_generator.py from /test/generator. Here is the mgxs_generator.py:

import openmc
import openmc.lib
import openmc.mgxs as mgxs
import numpy as np
import h5py

# create OpenMC model
model = openmc.Model()

Tmin = 293
T0 = 303 # temp in range that should be set-able
Tmax = 313

# create macroscopic object and export material with temperature XS library generated below
fuel = openmc.Material(material_id=1,name="fuel")
fuel.add_element('U', 1., enrichment=3.0)
fuel.set_density('macro',18.9)
materials = openmc.Materials([fuel])
model.materials = materials

# define geometry
H = 100.0 # cm
xplane0 = openmc.XPlane(x0  = 0.0, boundary_type='vacuum')
xmidplane = openmc.XPlane(x0  = H/2, boundary_type='transmission')
xplane1 = openmc.XPlane(x0 = H, boundary_type='vacuum')
yplane0 = openmc.YPlane(y0 = 0.0, boundary_type='vacuum')
yplane1 = openmc.YPlane(y0 = H, boundary_type='vacuum')
zplane0 = openmc.ZPlane(z0 = 0.0, boundary_type='vacuum')
zplane1 = openmc.ZPlane(z0 = H, boundary_type='vacuum')
cold_region = +xplane0 & -xmidplane & +yplane0 & -yplane1 & +zplane0 & -zplane1
hot_region = +xmidplane & -xplane1 & +yplane0 & -yplane1 & +zplane0 & -zplane1
cold_cell = openmc.Cell(cell_id=0, region=cold_region, fill=fuel)
cold_cell.temperature = Tmin
hot_cell = openmc.Cell(cell_id=1, region=hot_region, fill=fuel)
hot_cell.temperature = Tmax
geometry = openmc.Geometry([cold_cell,hot_cell])
model.geometry = geometry

# define settings
settings = openmc.Settings()
batches = 20
inactive = 10
particles = 5000
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.temperature = {'default': T0,
                        'method': 'interpolation',
                        'range': (Tmin, Tmax)} # good to load all temperatures you could encounter in multiphysics
bounds = [0, 0, 0, H, H, H]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.Source(space=uniform_dist)
settings.output = {'tallies': True,'summary':True}
model.settings = settings

# generate multigroup objects and tallies
groups = mgxs.EnergyGroups()
groups.group_edges = np.array([0.0, 20.0e6]) #  groups in eV

# create mgxs library object
onegxs_lib = openmc.mgxs.Library(geometry) # initialize the library
onegxs_lib.energy_groups = groups
onegxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission',
                        'nu-scatter matrix','scatter matrix', 'multiplicity matrix', 'chi']
# Specify a "cell" domain type for the cross section tally filters
onegxs_lib.domain_type = "cell"
# Specify the cell domains over which to compute multi-group cross sections
onegxs_lib.domains = geometry.get_all_cells().values()
# Set the Legendre order to 3 for P3 scattering
onegxs_lib.legendre_order = 3
# check library and build
onegxs_lib.check_library_for_openmc_mgxs()
onegxs_lib.build_library()
onegxs_lib.add_to_tallies_file(model.tallies)

# run model
sp = model.run()

# load statepoint
statepoint = openmc.StatePoint(sp, autolink=True)

# merge tallies and load from statepoint
onegxs_lib.load_from_statepoint(statepoint)
onegxs_file = onegxs_lib.create_mg_library(xs_type='macro')

# create xsdata object w gropu structure and two temperature library
xsdata = openmc.XSdata('temp_lib',energy_groups = groups,temperatures = np.array([Tmin,Tmax]) )
xsdata.order = 3

# set cold data
xsdata.set_total_mgxs(onegxs_lib.get_mgxs(cold_cell,'total'),temperature=Tmin,subdomain='all')
xsdata.set_absorption_mgxs(onegxs_lib.get_mgxs(cold_cell,'absorption'),temperature=Tmin,subdomain='all')
xsdata.set_fission_mgxs(onegxs_lib.get_mgxs(cold_cell,'fission'),temperature=Tmin,subdomain='all')
xsdata.set_scatter_matrix_mgxs(onegxs_lib.get_mgxs(cold_cell,'scatter matrix'),temperature=Tmin,subdomain='all')
xsdata.set_multiplicity_matrix_mgxs(onegxs_lib.get_mgxs(cold_cell,'multiplicity matrix'),temperature=Tmin,subdomain='all')
xsdata.set_chi_mgxs(onegxs_lib.get_mgxs(cold_cell,'chi'),temperature=Tmin,subdomain='all')

# set hot data
xsdata.set_total_mgxs(onegxs_lib.get_mgxs(hot_cell,'total'),temperature=Tmax,subdomain='all')
xsdata.set_absorption_mgxs(onegxs_lib.get_mgxs(hot_cell,'absorption'),temperature=Tmax,subdomain='all')
xsdata.set_fission_mgxs(onegxs_lib.get_mgxs(hot_cell,'fission'),temperature=Tmax,subdomain='all')
xsdata.set_scatter_matrix_mgxs(onegxs_lib.get_mgxs(hot_cell,'scatter matrix'),temperature=Tmax,subdomain='all')
xsdata.set_multiplicity_matrix_mgxs(onegxs_lib.get_mgxs(hot_cell,'multiplicity matrix'),temperature=Tmax,subdomain='all')
xsdata.set_chi_mgxs(onegxs_lib.get_mgxs(hot_cell,'chi'),temperature=Tmax,subdomain='all')

# creat MGXSLibrary and export to hdf5
temp_lib = openmc.MGXSLibrary(groups)
temp_lib.add_xsdata(xsdata)
temp_lib.export_to_hdf5(filename='/path/to/test/onegxs.h5')

Note you need to replace /path/to with the path to the test directory. At this point you will now have XML files in /test/generator as well as the file onegxs.h5 in /test. Next, you would go back to /test and run python MWE.py using the following script as MWE.py:

import openmc
import openmc.lib
import openmc.mgxs as mgxs
import numpy as np
import h5py

Tmin = 293
T0 = 303 # temp in range that should be set-able
Tmax = 313

# create macroscopic object and export material with temperature XS library generated below
cold_fuel = openmc.Material(material_id=1,name="cold_fuel")
hot_fuel = openmc.Material(material_id=2,name="hot_fuel")
# default names for the cross section set used in the CE file to generate the data are set1,set2,...
# there are two cells in the model and we computed XS for each
cold_fuel.add_macroscopic('temp_lib')
hot_fuel.add_macroscopic('temp_lib')
materials = openmc.Materials([cold_fuel,hot_fuel])
materials.cross_sections = 'onegxs.h5'
materials.export_to_xml()

# define geometry
H = 100.0 # cm
xplane0 = openmc.XPlane(x0  = 0.0, boundary_type='vacuum')
xmidplane = openmc.XPlane(x0  = H/2, boundary_type='transmission')
xplane1 = openmc.XPlane(x0 = H, boundary_type='vacuum')
yplane0 = openmc.YPlane(y0 = 0.0, boundary_type='vacuum')
yplane1 = openmc.YPlane(y0 = H, boundary_type='vacuum')
zplane0 = openmc.ZPlane(z0 = 0.0, boundary_type='vacuum')
zplane1 = openmc.ZPlane(z0 = H, boundary_type='vacuum')
cold_region = +xplane0 & -xmidplane & +yplane0 & -yplane1 & +zplane0 & -zplane1
hot_region = +xmidplane & -xplane1 & +yplane0 & -yplane1 & +zplane0 & -zplane1
cold_cell = openmc.Cell(cell_id=0, region=cold_region, fill=cold_fuel)
cold_cell.temperature = Tmin
hot_cell = openmc.Cell(cell_id=1, region=hot_region, fill=hot_fuel)
hot_cell.temperature = Tmax
geometry = openmc.Geometry([cold_cell,hot_cell])
geometry.export_to_xml()

# define settings
settings = openmc.Settings()
batches = 20
inactive = 10
particles = 5000
settings.energy_mode = 'multi-group'
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.temperature = {'default': T0,
                        'method': 'interpolation',
                        'range': (Tmin, Tmax)} # good to load all temperatures you could encounter in multiphysics
bounds = [0, 0, 0, H, H, H]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.Source(space=uniform_dist)
settings.export_to_xml()

# initialize an openmc simulation
openmc.lib.init()
# C-API call to set (valid) temperature of cell that results in error
openmc.lib.Cell.set_temperature(self=openmc.lib.cells[0],T=T0,set_contained=True)

It should produce the error

openmc.exceptions.OpenMCError: Temperature is below minimum temperature at which data is available.

openmc-dev / openmc

Cell::set_temperature causes exception when trying to set a valid temperature in multigroup mode with temperature interpolation policy #2296