search

openmc-dev / openmc

OpenMC Monte Carlo Code
https://docs.openmc.org
Other
720 stars 460 forks source link

Diff_burnable_mats fails on hex lattices #2841

Open gridley opened 5 months ago

gridley commented 5 months ago

Bug Description

Somehow, the differentiated burnable mats are not making it into the materials file when using CoupledOperator. I have been able to see this error ONLY on hex lattices. When I try the examples/lattice/simple problem, there seems to not be a problem.

Steps to Reproduce

Below is a script which can be placed in the lattices/hexagonal directory, based on depleting that example using the pincell_depletion chain, that exhibits the following error. The error specifically is:

 Reading settings XML file...
 Reading cross sections XML file...
 Reading materials XML file...
 Reading geometry XML file...
 ERROR: Could not find material 4 specified on cell 101

Interestingly, telling the model to differentiate the materials beforehand seems to work, whereas the legacy approach does not. Perhaps we should just deprecate the latter, and have the CoupledOperator call the model's volume splitting method. It would be nice to understand what's going wrong here, though. It is 100% for sure a python thing.

# WORKS:
# model.differentiate_depletable_mats("divide equally")
# op = openmc.deplete.CoupledOperator(model, chain_file)# diff_burnable_mats=True)

# BROKEN:
op = openmc.deplete.CoupledOperator(model, chain_file, diff_burnable_mats=True)

import openmc
import openmc.deplete

###############################################################################
#                      Simulation Input File Parameters
###############################################################################

# OpenMC simulation parameters
batches = 20
inactive = 10
particles = 10000

###############################################################################
#                 Exporting to OpenMC materials.xml File
###############################################################################

# Instantiate some Materials and register the appropriate Nuclides
fuel = openmc.Material(material_id=1, name='fuel')
fuel.set_density('g/cc', 4.5)
fuel.add_nuclide('U235', 1.)
fuel.volume = 69.420

moderator = openmc.Material(material_id=2, name='moderator')
moderator.set_density('g/cc', 1.0)
moderator.add_element('H', 2.)
moderator.add_element('O', 1.)
moderator.add_s_alpha_beta('c_H_in_H2O')

iron = openmc.Material(material_id=3, name='iron')
iron.set_density('g/cc', 7.9)
iron.add_element('Fe', 1.)

# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([moderator, fuel, iron])

###############################################################################
#                 Exporting to OpenMC geometry.xml file
###############################################################################

# Instantiate Surfaces
left = openmc.XPlane(surface_id=1, x0=-3, name='left')
right = openmc.XPlane(surface_id=2, x0=3, name='right')
bottom = openmc.YPlane(surface_id=3, y0=-4, name='bottom')
top = openmc.YPlane(surface_id=4, y0=4, name='top')
fuel_surf = openmc.ZCylinder(surface_id=5, x0=0, y0=0, r=0.4)

left.boundary_type = 'vacuum'
right.boundary_type = 'vacuum'
top.boundary_type = 'vacuum'
bottom.boundary_type = 'vacuum'

# Instantiate Cells
cell1 = openmc.Cell(cell_id=1, name='Cell 1')
cell2 = openmc.Cell(cell_id=101, name='cell 2')
cell3 = openmc.Cell(cell_id=102, name='cell 3')
cell4 = openmc.Cell(cell_id=500, name='cell 4')
cell5 = openmc.Cell(cell_id=600, name='cell 5')
cell6 = openmc.Cell(cell_id=601, name='cell 6')

# Use surface half-spaces to define regions
cell1.region = +left & -right & +bottom & -top
cell2.region = -fuel_surf
cell3.region = +fuel_surf
cell5.region = -fuel_surf
cell6.region = +fuel_surf

# Register Materials with Cells
cell2.fill = fuel
cell3.fill = moderator
cell4.fill = moderator
cell5.fill = iron
cell6.fill = moderator

# Instantiate Universe
univ1 = openmc.Universe(universe_id=1)
univ2 = openmc.Universe(universe_id=3)
univ3 = openmc.Universe(universe_id=4)
root = openmc.Universe(universe_id=0, name='root universe')

# Register Cells with Universe
univ1.add_cells([cell2, cell3])
univ2.add_cells([cell4])
univ3.add_cells([cell5, cell6])
root.add_cell(cell1)

# Instantiate a Lattice
lattice = openmc.HexLattice(lattice_id=5)
lattice.center = [0., 0., 0.]
lattice.pitch = [1., 2.]
lattice.universes = \
    [ [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ],
      [ [univ2] + [univ1]*11, [univ2] + [univ1]*5, [univ1] ],
      [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ] ]
lattice.outer = univ2

# Fill Cell with the Lattice
cell1.fill = lattice

# Instantiate a Geometry, register the root Universe, and export to XML
geometry = openmc.Geometry(root)

###############################################################################
#                   Exporting to OpenMC settings.xml file
###############################################################################

# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles

# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-1, -1, -1, 1, 1, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.IndependentSource(space=uniform_dist)

settings_file.keff_trigger = {'type' : 'std_dev', 'threshold' : 5E-4}
settings_file.trigger_active = True
settings_file.trigger_max_batches = 100

###############################################################################
#                   Exporting to OpenMC plots.xml file
###############################################################################

plot_xy = openmc.Plot(plot_id=1)
plot_xy.filename = 'plot_xy'
plot_xy.origin = [0, 0, 0]
plot_xy.width = [6, 6]
plot_xy.pixels = [400, 400]
plot_xy.color_by = 'material'

plot_yz = openmc.Plot(plot_id=2)
plot_yz.filename = 'plot_yz'
plot_yz.basis = 'yz'
plot_yz.origin = [0, 0, 0]
plot_yz.width = [8, 8]
plot_yz.pixels = [400, 400]
plot_yz.color_by = 'material'

# Instantiate a Plots collection, add plots, and export to XML
plot_file = openmc.Plots((plot_xy, plot_yz))

###############################################################################
#                   Exporting to OpenMC tallies.xml File
###############################################################################

# Instantiate a distribcell Tally
tally = openmc.Tally(tally_id=1)
tally.filters = [openmc.DistribcellFilter(cell2)]
tally.scores = ['total']

# Instantiate a Tallies collection and export to XML
tallies_file = openmc.Tallies([tally])

model = openmc.Model(tallies=tallies_file,
                     geometry=geometry,
                     materials=materials_file,
                     plots=plot_file,
                     settings=settings_file)
chain_file = '../../pincell_depletion/chain_simple.xml'

# WORKS:
# model.differentiate_depletable_mats("divide equally")
# op = openmc.deplete.CoupledOperator(model, chain_file)# diff_burnable_mats=True)

# BROKEN:
op = openmc.deplete.CoupledOperator(model, chain_file, diff_burnable_mats=True)

# Perform simulation using the predictor algorithm
time_steps = [1.0, 1.0, 1.0, 1.0, 1.0]  # days
power = 174  # W/cm, for 2D simulations only (use W for 3D)
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d')
integrator.integrate()

Environment

arch linux, latest commit

shimwell commented 5 months ago

Just to mention there is an open PR for improving the diff burnable mat over here might be worth checking to see if this changes things.

lewisgross1296 commented 5 months ago

Could this be related to #2798?

icmeyer commented 3 months ago

Ran into this issue independently. Here was my input (only have a broken one!)

import openmc
import openmc.deplete

def run_test():
    # Generate materials
    fuel = openmc.Material(1, "uo2")
    fuel.add_element('U', 1.0, enrichment=3.0)
    fuel.add_element('O', 2.0)
    fuel.set_density('g/cm3', 10.37)
    fuel.temperature = 900
    fuel.volume = 1.0

    materials = openmc.Materials([fuel])

    # Geometry
    fuel_cell = openmc.Cell(fill=fuel) # change to mat_dict['uo2']
    fuel_universe = openmc.Universe(cells=[fuel_cell])

    y0 = openmc.YPlane(0, boundary_type='reflective')
    y1 = openmc.YPlane(1)
    y2 = openmc.YPlane(2, boundary_type='reflective')

    x0 = openmc.XPlane(0, boundary_type='reflective')
    x1 = openmc.XPlane(1)
    x2 = openmc.XPlane(2, boundary_type='reflective')

    fuel1 = openmc.Cell(fill=fuel_universe, region= (+x0 & -x1 & +y0 & -y1))
    fuel2 = openmc.Cell(fill=fuel_universe, region= (+x0 & -x1 & +y1 & -y2))
    fuel3 = openmc.Cell(fill=fuel_universe, region= (+x1 & -x2 & +y0 & -y1))
    fuel4 = openmc.Cell(fill=fuel_universe, region= (+x1 & -x2 & +y1 & -y2))

    universe = openmc.Universe(cells=[fuel1, fuel2, fuel3, fuel4])
    geometry = openmc.Geometry(universe)

    settings = openmc.Settings()
    settings.run_mode = 'eigenvalue'
    settings.particles = 1000
    settings.batches = 100
    settings.inactive = 10

    model = openmc.model.Model(geometry=geometry, 
                               materials=materials,
                               settings=settings)

    chain_file = 'chain_endfb71_pwr.xml'
    chain = openmc.deplete.Chain.from_xml(chain_file)
    operator = openmc.deplete.CoupledOperator(model, chain_file,
                                              diff_burnable_mats=True)
    power = 14e6  # W
    time_steps = [1] * 30 # days
    integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power,
                                                    timestep_units='d')
    integrator.integrate()

if __name__ == "__main__":
    run_test()

Output:

                 | The OpenMC Monte Carlo Code
       Copyright | 2011-2023 MIT, UChicago Argonne LLC, and contributors
         License | https://docs.openmc.org/en/latest/license.html
         Version | 0.14.1-dev
        Git SHA1 | 3901709141f6e1e85707fa1ce901199df081cc29
       Date/Time | 2024-03-11 10:41:51
  OpenMP Threads | 20

 Reading settings XML file...
 Reading cross sections XML file...
 Reading materials XML file...
 Reading geometry XML file...
 ERROR: Could not find material 2 specified on cell 1

This crashes with ERROR: Could not find material 2 specified on cell 1. This comes from the model object calling https://github.com/openmc-dev/openmc/blob/7ed12788df9cc9c1e6582ab4f634a91fde88be89/openmc/model/model.py#L1025-L1074 which calls determine_paths() for the model.materials and then the CoupledOperator exports the old materials object: https://github.com/openmc-dev/openmc/blob/7ed12788df9cc9c1e6582ab4f634a91fde88be89/openmc/deplete/coupled_operator.py#L395-L408

I tried to change the above function to export the model.materials and that gets the simulation a little further, but there is another bug then with the burnable mats that have been set by the original materials in the OpenMCOperator bas class of CoupledOperator.

Output below:

Traceback (most recent call last):
  File "/home/icmeyer/ff_adsr/3d/mitr/depletion_test/test.py", line 57, in <module>
    run_test()
  File "/home/icmeyer/ff_adsr/3d/mitr/depletion_test/test.py", line 54, in run_test
    integrator.integrate()
  File "/home/icmeyer/openmc/openmc/deplete/abc.py", line 781, in integrate
    n = self.operator.initial_condition()
  File "/home/icmeyer/openmc/openmc/deplete/coupled_operator.py", line 391, in initial_condition
    materials = [openmc.lib.materials[int(i)] for i in self.burnable_mats]
  File "/home/icmeyer/openmc/openmc/deplete/coupled_operator.py", line 391, in <listcomp>
    materials = [openmc.lib.materials[int(i)] for i in self.burnable_mats]
  File "/home/icmeyer/openmc/openmc/lib/material.py", line 282, in __getitem__
    raise KeyError(str(e))
KeyError: 'No material exists with ID=1.'

Could be that my input is not the canoncial way to start a depletion calculation, but the Example documents are pretty sparse on depletion. But anyways, not super familiar with the depletion model and how we expect to be handling the different materials objects the process, will try out the suggested PR and see if it gets me anywhere.

gridley commented 3 months ago

@icmeyer did you ever try taking diff_burnable_mats out of the depletion call, and instead calling model.differentiate_depletable_mats("divide equally") before depleting?

icmeyer commented 3 months ago

Just tried and it works great! Thanks for the tip. So could this be fixed by removing the internal CoupledOperator logic and just adding something in its __init__ like if diff_burnable_mats: model.differential_deplatable_mats(option)?

gridley commented 3 months ago

I agree, that'd be the way to do it. Really, I think we should deprecate the diff_burnable_mats kwarg and force users to use the model method first, since you explicitly have to specify the division technique.

Your approach would be best for backwards-compatibility. Maybe someone else can weigh in on this.