Problems running a depletion calculation with neutron-photon interaction, for the proper energy deposition calculation and the subsequent proper power normalization

User Support
1

Hello !

So I installed the newest version from the dev branch of OpenMC, and I have the ENDFB71 libraries (with photon data) given by the OpenMC team. If I run the depletion module with the ‘energy-deposition’ attribute and, at the same time, I let the code know that I want a photon transport calculation as well (set to ‘True’ under settings), it gives me a segmentation fault when the transport calculation just begins (attached is the slurm exit with the error, and it seems to be related with the photon transport).

On the other hand, if I just run the depletion module without the photon transport, the depletion calculation runs just fine (both with or without the ‘energy-deposition’ attribute). I also attached the input file in case you would like to look at it.

This is strange to me because for example, once it crashes, the python input file for depletion still export the xml files (in here, the settings.xml would have the True attribute for photon transport). If I just execute openmc in standalone mode with this xml files, then the code actually is able to run successfully the transport calculation (so the first transport calculation at burnup 0, runs in standalone). Thus, it seems that the problem is with the depletion module wrapping the transport photon-neutron calculation. Any idea what it might be?

Thank you in advance.

Augusto

slurm-57224.out (35 KB)

pin_dep.py (8.09 KB)

2

Update: I think is related to this same issue. Now while depleting, I am trying to run a particle filter (in this case, only the one corresponding to the ‘neutron’ particle) in order to estimate the energy deposition as a function of burnup in the materials of the domain (also, the ‘energy-deposition’ attribute is taken into account). In the end, I add a ‘heating-local’ tally to the exercise (all these is depicted below):

slurm-62553.out (10.5 KB)

3

Hi Augusto,

Thanks for reporting this. Yes, it appears that depletion doesn’t work with photon transport at present. I’m going to look into getting this fixed and will hopefully have a pull request soon.

Best,
Paul

4

Hi Augusto,

Just wanted to let you know that a fix for this issue has now been merged in the ‘develop’ branch of OpenMC. Let us know if you have any further issues.

Best,
Paul

5

Hi Paul,

Thanks a lot for letting me know. I already tried it. In fact, if I deplete with photon transport, as well as allowing the ‘energy-deposition’ variable, it works fine (I can even sense how slower the transport calculation is compared to the case when you only have neutrons :slight_smile: ).

Just one more thing. It actually seems that the problem related to the ParticleFilter while depleting is still there. I don’t even need to activate photon transport to get a seg fault while trying to use a ParticleFilter (meaning that just by neutron transport and not even activating ‘energy-deposition’, if I use a neutron ParticleFilter the code crashes).

What I want is to tally the different energy deposition per particle while depleting (i.e. neutron, photon, electron and positron) in different cell/materials. On the other hand, if I do not use the ParticleFilter, and then just tally ‘heating’ using a MaterialFilter and allowing ‘energy-deposition’ and photon transport, then it is ok; then I could actually get the total energy deposition at different materials using a coupled neutron-photon calculation. As I mentioned to you, I just want to separate such energy deposition contribution per type of particle. It actually seems that either with or without photon transport, the ParticleFilter is only working in transport standalone and not with the depletion module.

Hope there can be a solution for this.

Best regards,

Augusto

6

Hi again (and sorry for disturbing so much),

I have an inquiry while running the depletion pin case with photon transport. It is running; however, there are some warnings in my slurm output that seem strange to me. As I mentioned in my previous post, there is a problem (a seg fault) if I add the ParticleFilter during depletion. But if I don’t add such type of filter and I just ask to tally the total deposited heat, it appears the following warning:

WARNING: Particle filter is not used with photon transport on and flux score.
WARNING: Particle filter is not used with photon transport on and fission
score.

I wonder if actually the energy deposition (the one obtained with ‘heating’ variable) would be properly computed by means of the coupled neutron-photon even if the ParticleFilter is not present.

Another thing; also, during the transport calculation it appears the following (even more than once):

WARNING: Particle 32843 underwent maximum number of events.

I attached both the slurm and the input file that is currently being executed, in case you like to see them. Now, I found strange the “maximum number of events” warning, because I had already depleted the pin only with neutron transport and such warning never appeared (only appears now with coupled neutron-photon).

Thanks.

Best regards,

Augusto

slurm-65577.out (22.6 KB)

pin_dep.py (9.8 KB)

7

Hi Augusto – thanks for your persistence. This is indeed another bug; I’ve just submitted a fix for it and confirmed that it works for your model (thanks for sharing). Hopefully the fix will hit the develop branch soon.

Best,
Paul

8

Dear Paul,

Thanks for letting me know and for your effort in improving the code. I took the liberty to follow the “fix for it” link with the pull request, and I realized that it has been already incorporated into the dev branch (I think so :wink: ). So I took the liberty to download this new version and re-compiled it, and to execute once more the pin depletion case.

Now the ParticleFilter problems have been fixed; in fact, now I can actually even tally the energy deposition at different burnup points (see just below for instance, the output of the tally at the last burnup point, which is around 27 MWd/KgU). Nevertheless, I have to say that I did not like at all the outputs of OpenMC while depleting with neutron-photon coupling (after the particle energy deposition tally, I will keep explaining my founds).

kinf_pin_edep_comp.jpg
kinf_pin_edep_comp.jpg1920×879 60.5 KB

U235_pin_edep.jpg
U235_pin_edep.jpg1920×879 70.7 KB

Pu239_pin_edep.jpg
Pu239_pin_edep.jpg1920×879 65.5 KB

9

Hi Augusto,

Thanks for reporting this. I have observed the same issue and have some fixes underway that I will be sending a pull request for soon. I’ll let you know when I have a branch ready to test if you want to try out your problem with a proposed fix.

Best,
Paul

10

Hi! Paul,
Im meeting the problem now, And the version Im using is 0.12.2
if the WARNING still be there, what should I do to get the flux distribution and the power distribution correctly ?
(for my case , the power is known, and the run_mode is default)
Thank you
lizhuang

11

If you want the neutron flux distribution, you should add a particle filter that specifies neutrons:

particle_filter = openmc.ParticleFilter(['neutron'])
tally = openmc.Tally()
tally.filters = [particle_filter, ...]
...

For the power distribution in a coupled n-γ problem, you can use the “heating” score:

power_tally = openmc.Tally()
power_tally.scores = ['heating']
12

Hello dr. Romano,

I’m sorry for resurrecting this topic! I’m facing the same issue (segment fault) with depletion considering energy deposition and photon transport. I’m using version 0.13. If you can also let me know when this new branch is ready to test, I will be glad.

Thanks in advance,

Luiz.

13

@luizaldeia The fix that was referenced above has been merged in the code since 2020, so if you are experiencing issues they are likely unrelated to the above. If you’re able to share a model that produces the segfault, I can try to reproduce it.

14

Hi @paulromano, for sure I can.

I’m studying BISON depletion capability and planning to couple BISON + OpenMC using Cardinal in the future. I’m using the following script to reproduce a rod irradiation process:

# BISON depletion capability comparison
# Rep_Na_3 assessment case, fuel pellet mesh with 55 radial divisions.
# BISON solves the depletion for regions with the same radius but different volumes.
# The number of radial divisions in the mesh controls the number of the depletion regions

import openmc
import openmc.deplete
import numpy as np
import math
from matplotlib import pyplot


######################################################################
#                      Operational Parameters
######################################################################

fuel_temp = 978.31
gap_temp = 663.69
water_den = 0.7051
ppm_boron = 500
water_temp = 585.25

######################################################################
#                      Geometric Parameters
######################################################################


Rfo = 0.40956  #pellet outer radius
Rci = 0.41776  #clad inner radius
Rco = 0.47736  #clad outer radius
hpitch = 0.63    #half-pitch between rods
fheight = 1.3589 #fuel height
nr = 55 #number of radial divisions

######################################################################
#                      Chain Parameters
######################################################################

#chain_path = "/storage/home/lca5209/OpenMC/chains/chain_simple.xml"
#chain_path = "/storage/home/lca5209/OpenMC/chains/chain_casl_pwr.xml"
chain_path = "/storage/home/lca5209/OpenMC/chains/chain_endfb71_pwr.xml"

######################################################################
#                      Settings
######################################################################

batches = 100
inactive = 20
particles = 50000
seed=np.random.randint(2147483647)*2+1

######################################################################
#                      Setting the materials
######################################################################
#fuel = 3.0 wt% enriched uranium dioxide

fuel = openmc.Material(material_id=1,name="UO2")

fuel.add_nuclide("U238",0.84223593,'wo')
fuel.add_nuclide("U235",0.03918626,'wo')
fuel.add_nuclide("O16",0.11857781,'wo')
fuel.set_density("g/cc",10.394075)
fuel.temperature = fuel_temp
fuel.volume = (math.pi*Rfo**2)*fheight

mat_list = []

for i in range(nr):
	mat = fuel.clone()
	mat.volume = (math.pi*((((i+1)/55)*Rfo)**2-(((i)/55)*Rfo)**2))*fheight
	mat_list.append(mat)
	
#gap = helium

gap = openmc.Material(name='He')
gap.set_density('g/cc',0.001598)
gap.add_element('He', 1.0,'wo')
gap.temperature = gap_temp
mat_list.append(gap)

#clad = zircaloy

clad = openmc.Material(name='Zircaloy_4')
clad.set_density('g/cc', 6.56)
clad.add_element('Sn', 0.014  , 'wo')
clad.add_element('Fe', 0.00165, 'wo')
clad.add_element('Cr', 0.001  , 'wo')
clad.add_element('Zr', 0.98335, 'wo')
mat_list.append(clad)

#coolant = borated water

coolant = openmc.Material(name='Borated_water')
coolant.set_density('g/cc', water_den) 
coolant.add_nuclide('B10',ppm_boron*0.199E-06,'wo') 
coolant.add_nuclide('B11',ppm_boron*0.801E-06,'wo') 
coolant.add_element('H',0.1108612,'wo')
coolant.add_element('O',0.8886388,'wo')
coolant.remove_nuclide('O18')
coolant.temperature= water_temp 
coolant.add_s_alpha_beta('c_H_in_H2O')
mat_list.append(coolant)

print(mat_list)


materials = openmc.Materials(mat_list)
materials.export_to_xml()

######################################################################
#                      Building the geometry
######################################################################

# Fuel rings, gap, and cladding

cylinder = []

counter_1 = 0
for i in range(nr):
	c = openmc.ZCylinder(surface_id=(i+1), x0=0, y0=0, r=((i+1)/55)*Rfo, name='Fuel '+ str(i+1) + ' OR')
	cylinder.append(c)
	counter_1 +=1
cylinder_g = openmc.ZCylinder(surface_id=counter_1+1, x0=0, y0=0, r=Rci, name='Gap OR')
cylinder_c = openmc.ZCylinder(surface_id=counter_1+2, x0=0, y0=0, r=Rco, name='Clad OR')
cylinder.append(cylinder_g)
cylinder.append(cylinder_c)

# Box surrounding fuel rod
left = openmc.XPlane(surface_id=counter_1+3, x0=-hpitch, name='left')
right = openmc.XPlane(surface_id=counter_1+4, x0=hpitch, name='right')
bottom = openmc.YPlane(surface_id=counter_1+5, y0=-hpitch, name='bottom')
top = openmc.YPlane(surface_id=counter_1+6, y0=hpitch, name='top')

# Surfaces above/below fuel
fsouth = openmc.ZPlane(surface_id=counter_1+7, z0=-fheight/2.0, name='fsouth') 
fnorth = openmc.ZPlane(surface_id=counter_1+8, z0=fheight/2.0, name='fnorth')

# Boundary Conditions
left.boundary_type = 'reflective'
right.boundary_type = 'reflective'
top.boundary_type = 'reflective'
bottom.boundary_type = 'reflective'
fsouth.boundary_type = 'reflective'
fnorth.boundary_type = 'reflective'

cell = []

counter_2 = 0
for i in range(nr):
	f = openmc.Cell(cell_id=(i+1), name='fuel '+str(i+1))
	cell.append(f)
	counter_2 += 1
gap_cell = openmc.Cell(cell_id=counter_2+1, name='Gap cell')
clad_cell = openmc.Cell(cell_id=counter_2+2, name='Clad cell')
water_cell = openmc.Cell(cell_id=counter_2+3, name='Water cell')
cell.append(gap_cell)
cell.append(clad_cell)
cell.append(water_cell)

for i in range(len(cell)-1):
	if i==0:
		cell[i].region = -cylinder[i] & -fnorth & +fsouth
	elif i==(len(cell)-2):
		cell[i+1].region = +cylinder[i] & +left & -right & +bottom & -top & -fnorth & +fsouth
		cell[i].region = +cylinder[(i-1)] & -cylinder[i] & -fnorth & +fsouth
	else:
		cell[i].region = +cylinder[(i-1)] & -cylinder[i] & -fnorth & +fsouth
		
count_3 = 0
for i in range(nr):
	cell[i].fill = mat_list[i]
	count_3 +=1
cell[count_3].fill = gap
cell[count_3+1].fill = clad
cell[count_3+2].fill = coolant


# Instantiate Universe
root = openmc.Universe(universe_id=0, name='root universe')

# Register Cells with Universe
root.add_cells(cell)

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

plot = openmc.Plot.from_geometry(geometry)
plot.pixels = (250, 250)
plot.filename = 'pinplot_55_regions'
plot.width = (hpitch*2, hpitch*2)
plot.pixels = (500, 500)
plot.color_by = 'cell'
plot.colors = {cell[-3]: 'purple',cell[-2]: 'grey',
               cell[-1]: 'blue'}
plots = openmc.Plots([plot])
plots.export_to_xml()
openmc.plot_geometry()

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

# Instantiate a Settings object, set all runtime parameters, and export to XML
settings = openmc.Settings()
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.seed = seed
settings.verbosity = 10
settings.photon_transport = False

# Cross-section temperature information
settings.temperature = {
    'method': 'interpolation',
    'tolerance': 100 ,
    'multipole': True
}

# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-hpitch, -hpitch, -fheight/2, hpitch, hpitch, fheight/2]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.source.Source(space=uniform_dist)

entropy_mesh = openmc.RegularMesh()
entropy_mesh.lower_left = [-hpitch, -hpitch, -1.e50]
entropy_mesh.upper_right = [hpitch, hpitch, 1.e50]
entropy_mesh.dimension = [10, 10, 1]


settings.entropy_mesh = entropy_mesh
settings.export_to_xml()

###############################################################################
#                   Exporting to OpenMC tallies.xml file
###############################################################################

# Instantiate a tally mesh
mesh = openmc.RegularMesh()
mesh.type = 'regular'
mesh.dimension = [1, 1, 100]
mesh.lower_left = [-hpitch, -hpitch, -fheight/2]
mesh.upper_right = [hpitch, hpitch, fheight/2]

# Instantiate some tally Filters
energy_filter = openmc.EnergyFilter([0., 0.625, 20.e6])
energy_filter2 = openmc.EnergyFilter(np.logspace(-3, 7, num=100))
mesh_filter = openmc.MeshFilter(mesh)

# Instantiate the Tally
tally = openmc.Tally(tally_id=1, name='Spatial Tally 1')
tally.filters = [energy_filter, mesh_filter]
tally.scores = ['flux']#, 'fission', 'nu-fission']

# Instantiate the Tally
tallyspec = openmc.Tally(tally_id=2, name='Energy Spectrum Tally 1')
tallyspec.filters = [energy_filter2]
tallyspec.scores = ['flux']#, 'fission', 'nu-fission']

particle_filter = openmc.ParticleFilter(['neutron', 'photon'])
mat_filter = openmc.MaterialFilter(materials)

tallydep = openmc.Tally(tally_id=3, name='Spatial Tally 2')
tallydep.filters = [mat_filter, particle_filter]
tallydep.scores = ['fission', 'heating', 'heating-local']

# Instantiate a Tallies collection and export to XML
tallies = openmc.Tallies([tally, tallyspec,tallydep])
tallies.export_to_xml()



###############################################################################
#                      Depletion settings
###############################################################################

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

#setting the transport operator
operator = openmc.deplete.Operator(model,chain_path,diff_burnable_mats='False',normalization_mode='fission-q',fission_q={"U235": 202.27e6})

#setting the system linear power [W]
power = [45.50,91.01,136.51,182.02,227.52,273.03,324.45,324.45,273.03,273.03,227.52,182.02,136.51,91.01,45.50]
time_steps = [8640.00,8640.00,8640.00,8640.00,8640.00,8640.00,4268160.00,45792000.00,1728000.00,46915200.00,8640.00,8640.00,8640.00,8640.00,8640.00]
print(time_steps)

#depleting usin a first-order predictor algorithm
integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power,timestep_units =  's')

integrator.integrate()

If I set settings.photon_transport = True the segfaul error happens.

Thanks for your time and attention!

Rep_Na_3_55.py (9.5 KB)

15

@luizaldeia This is actually unrelated to depletion. Right now multipole data unfortunately does not work in conjunction with photon transport, and fixing it is not going to be simple. In the meantime, I’ll submit a fix to the code to at least give a better error message rather than just segfaulting. Thanks for bringing this up and sorry for the trouble!

One other thing I should point out in your script:

This argument needs to be a boolean value (True or False with no quotes around it). In Python, any non-empty string is considered true as far as logic expressions go, so if you pass a string with ‘False’ in it, OpenMC thinks it is True.

1 Like
16

Hello @paulromano,

Thank you again for your time, attention, and kind answer. Also, thanks for the advice about the diff_burnable_mats. It was a complete lack of attention on my part. I will fix this!!!

Have a great night!

1 Like