Heating and Flux score affecting each other

User Support
1

I am running a very simple fixed source simulation:

  • Geometry: spherical shell (r_min = 800cm, r_max = 1050cm)
  • Material: Li17Pb83
  • Source: isotropic at 14.1MeV
  • settings.photon_transport = True

I am trying to compute the radial distribution of the neutron flux and power deposition. To do so I wrote the following tallies

radial_mesh = openmc.SphericalMesh(name='Radial Mesh', r_grid=r_coords)
mesh_filter = openmc.MeshFilter(radial_mesh)
particle_filter = openmc.ParticleFilter(['neutron', 'photon'])

radial_tally = openmc.Tally(name='Radial')
radial_tally.filters = [mesh_filter, particle_filter]
radial_tally.scores = ['flux', 'heating']
tallies.append(radial_tally)

The odd behaviour that I am noticing is that the resulting flux distribution changes if I don’t score the heating, but only the flux, i.e. radial_tally.scores = ['flux']. The changes are not minimal but clearly visible when plotting the radial distribution.

Interestingly, this discrepancy doesn’t appear if I filter over the entire Cell instead of using the mesh_filter.

Intuition leads me to think that the cleanest and correct result is obtained separing the scores in different tallies. However, only if I include the heating score then the average flux (weighted by the volumes) obtained from the radial distribution matches the one from filtering over the entire Cell.

Screenshot 2025-12-23 164310
Screenshot 2025-12-23 1643101943Ă—912 89.7 KB

Thank you in advance!

1 Like
2

The problem turned out to be caused by an incorrect normalization automatically applied by OpenMC.

The Track-Length estimator normalizes the flux computed dividing it by the volume of the filter region. This works for cell filters, because OpenMC knows exactly the volume of the cell. But for a mesh filter, OpenMC does not know the true volume inside each mesh voxel.

Adding the heating score magically fixes the mesh flux because it makes OpenMC switch to the Collision estimator, which normalizes by the actual material interaction rate instead of the volume.

The correct way to compute the radial flux is to use concentric spherical cells and a cell filter instead of the spherical mesh.

3

This is not correct. OpenMC does not automatically do any volume normalization and by default it does not know the volumes of cells. Tallied flux results are always given in units of particle-cm/source particle. It is true that you will get slightly different results from applying different tally estimators, but they should be within statistical uncertainty of one another. If you are finding that the same tally gives statistically different results using tracklength vs collision estimator, I would consider that to be a bug, but for sure any differences you are seeing should not be due to volume normalization.

4

Thanks Paul for your contribution.
You are right OpenMC doesn’t do any volume normalization.
However the discrepancy between the results I get from applying different estimators is way bigger than statistical uncertainty (around 50 times standard deviation).
The only way I managed to compute correctly the radial flux is to use concentric spherical cells and a cell filter instead of the spherical mesh.

5

Understood, thanks Luca. If you’re able to share a minimal working example of the behavior you’re seeing, we can take a look and see what’s going on.

6

Here is the code that I am running

# MATERIAL
blanket_material = openmc.Material(name='Breeder')
blanket_material.set_density('g/cm3', 9.5)
enrichment_fraction = 0.075
blanket_material.add_nuclide('Li6', 17 * enrichment_fraction, 'ao')
blanket_material.add_nuclide('Li7', 17 * (1.0 - enrichment_fraction), 'ao')
blanket_material.add_element('Pb', 83, 'ao')
materials = openmc.Materials([blanket_material])
materials.export_to_xml('./materials.xml')

# SOURCE
source = openmc.IndependentSource()
source.name = 'DT Fusion Neutron Source'
source.space = openmc.stats.Point((0, 0, 0))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14.1e6], [1.0])
source.energy.units = 'eV'

# GEOMETRY
r_min = 800.0
r_max = 1050.0
blanket_volume = (4/3) * np.pi * (r_max**3 - r_min**3)
r_coords = np.linspace(r_min, r_max, 101)
volumes = (4/3) * np.pi * (r_coords[1:]**3 - r_coords[:-1]**3)
R_inner = openmc.Sphere(r=r_min)
R_outer = openmc.Sphere(r=r_max, boundary_type='vacuum')
center_void = openmc.Cell(name='Center Void', fill=None)
center_void.region = -R_inner
blanket_cell = openmc.Cell(name='Blanket Shell', fill=blanket_material)
blanket_cell.region = +R_inner & -R_outer
universe = openmc.Universe(cells=[center_void, blanket_cell])
geometry = openmc.Geometry(universe)
geometry.export_to_xml('./geometry.xml')

# SETTINGS
settings = openmc.Settings()
settings.batches = 10
settings.inactive = 0
settings.particles = 100000
settings.source = source
settings.run_mode = 'fixed source'
settings.photon_transport = True
settings.export_to_xml('./settings.xml')

tallies = openmc.Tallies()
radial_mesh = openmc.SphericalMesh(name='Radial Mesh', r_grid=r_coords)
mesh_filter = openmc.MeshFilter(radial_mesh)
cell_filter = openmc.CellFilter(blanket_cell)
particle_filter = openmc.ParticleFilter(['neutron', 'photon'])
blanket_tally = openmc.Tally(name='Blanket')
radial_tally = openmc.Tally(name='Radial')

blanket_tally.scores = ['flux', 'heating'] # change to only flux
radial_tally.scores = ['flux', 'heating']  # change to only flux

blanket_tally.filters = [cell_filter, particle_filter]
radial_tally.filters = [mesh_filter, particle_filter]
tallies.append(blanket_tally)
tallies.append(radial_tally)
tallies.export_to_xml('./tallies.xml')

And this is the code to read the results

sp = openmc.StatePoint('statepoint.{}.h5'.format(settings.batches))
blanket_result = sp.get_tally(name='Blanket')
flux = blanket_result.get_values(scores=['flux'], value='mean', filters=[openmc.ParticleFilter], filter_bins=[('neutron',)])[0][0][0]
print(f"Neutron Flux: {flux/blanket_volume:.4e}")

radial_result = sp.get_tally(name='Radial')
flux_values = radial_result.get_values(scores=['flux'], value='mean', filters=[openmc.ParticleFilter], filter_bins=[('neutron',)]).flatten()
flux_std_dev  = radial_result.get_values(scores=['flux'], value='std_dev', filters=[openmc.ParticleFilter], filter_bins=[('neutron',)]).flatten()
flux_values /= volumes
flux_std_dev /= volumes
radial_flux_mean = np.sum(flux_values * volumes) / blanket_volume
radial_flux_std_dev = np.sqrt(np.sum((flux_std_dev * volumes)**2)) / blanket_volume
print(f"Radial Flux: {radial_flux_mean:.4e} +- {radial_flux_std_dev:.4e}")

sp.close()
7

I’ve opened PR #3792 to fix this problem.