Question for mesh tallys and cells, very different values

Hello,

Using OpenMC I’ve got some strange values on some cases.

I’ve found that OpenMC is able to give me Mesh tally values non consistant on the same space for spherical or regular Mesh

===================> TALLY 8: DOSE G AVANT 1M <====================

Mesh Index (1, 1, 1)
Energy Function f([1.0e+04, …, 1.0e+10]) = [2.9e-02, …, 2.8e+02]
Particle: photon
Total Material
Flux 13151.5 +/- 61.4526
==================> TALLY 9: DOSE G AVANT 1M R <===================

Mesh Index (1, 1, 1)
Energy Function f([1.0e+04, …, 1.0e+10]) = [2.9e-02, …, 2.8e+02]
Particle: photon
Total Material
Flux 7.81904e+06 +/- 1591.38

Both tally have the same volume (around 64k)
statepoint.tally.h5 (78.3 KB)
tally_exemple.ipynb (1.3 MB)

But when I’ve modelised the sphere in geometry, values were consitent for the 3 tallys…

===================> TALLY 11: DOSE G AVANT 1M <===================

Mesh Index (1, 1, 1)
Energy Function f([1.0e+04, …, 1.0e+10]) = [2.9e-02, …, 2.8e+02]
Particle: photon
Total Material
Flux 8.00043e+06 +/- 16269.6
==================> TALLY 12: DOSE G AVANT 1M R <==================

Mesh Index (1, 1, 1)
Energy Function f([1.0e+04, …, 1.0e+10]) = [2.9e-02, …, 2.8e+02]
Particle: photon
Total Material
Flux 7.82251e+06 +/- 15634.9
================> TALLY 13: DOSE G AVANT 1M CELL <=================

Cell 6
Energy Function f([1.0e+04, …, 1.0e+10]) = [2.9e-02, …, 2.8e+02]
Particle: photon
Total Material
Flux 8.00513e+06 +/- 16268.7

statepoint.tally_and_vol.h5 (84.5 KB)
tally_and_vol.ipynb (144.9 KB)

I don’t see were is my mistake or if there is some thing to avoid using tallys ?

I’ve done an orher test, with an other Spherical tally just 0.01 cm away from the cell.

The one on the cell is still adequate with the cell, the one not coincident is giving back the odd results.

Is this an identificated bug ?

Dear Thomas,

I have tried to understand your problem and test case. Please find here my thoughts.

In your fixed source simulation, you’ve created a spherical mesh and an isotropic point source located at the origin of the reference frame. The source emits particles in air, in a medium that extends almost to infinity. The spherical mesh is very simple and represents a sphere of radius of 25cm as a single shell centered at a distance of 1m from the source.

Within the spherical mesh, I’ve set two different flux tallies: the first one uses a collision estimator, while the second one uses a tracklength estimator.

The estimate of the flux obtained with the collision estimator is fine, while the one obtained with the tracklength estimator seems unphysical (orders of magnitude lower).

If the source is placed inside the sphere, both estimates are good.

If I use a Cartesian mesh instead of a spherical one, thus using a cube as a detector zone instead of a sphere, both estimates still match.

If I use a spherical volume with the same dimensions as the spherical mesh, but created with a CellFilter, as the detector zone, both estimates still match.
An interesting fact is that, like you said, if the spherical volume created with a CellFilter matches with the spherical mesh, the tracklength and collision estimates match in the tallies set within the spherical mesh too.

I will paste my code in the comment below, since as a new user I can’t upload the file.

We are afraid that there could be a problem with the calculation of the chord length in the case of the spherical mesh. We would like to have the developers’ opinion on this matter.

Thank you very much in advance for your kind attention and support.

import openmc
import numpy as np
import pandas as pd

import os
os.environ[‘OPENMC_CROSS_SECTIONS’] = ‘’
openmc.config[‘cross_sections’] = ‘’

Materials

air=openmc.Material(name=‘air’)
air.add_element(‘N’,80)
air.add_element(‘O’,20)
air.set_density(‘g/cm3’,1.3e-3)

materials = openmc.Materials([air])
materials.export_to_xml()

Geometry

Limiting Planes

Xm = openmc.XPlane( x0=-300, name=‘X-’, boundary_type=‘vacuum’)
Xp = openmc.XPlane( x0=+300, name=‘X+’, boundary_type=‘vacuum’)
Ym = openmc.YPlane( y0=-300, name=‘Y-’, boundary_type=‘vacuum’)
Yp = openmc.YPlane( y0=+300, name=‘Y+’, boundary_type=‘vacuum’)
Bottom= openmc.ZPlane(z0= 0, boundary_type=‘vacuum’)
Top= openmc.ZPlane(z0= 300, boundary_type=‘vacuum’)

Create root Universe

root_universe = openmc.Universe(universe_id=0, name=‘root universe’)

Create a cell to be used as a spherical detector

detector_sphere = openmc.Sphere(r = 25, x0 = 100, y0 = 0, z0 = 150) #Uncomment to simulate the case in which the source is outside the sphere
#detector_sphere = openmc.Sphere(r = 25, x0 = 0, y0 = 0, z0 = 150) #Uncomment to simulate the case in which the source is inside the sphere
#detector_sphere = openmc.Sphere(r = 25, x0 = 0, y0 = 100, z0 = 150) #Uncomment to simulate the effect of superimposing MeshFilter and CellFilter
detector_cell = openmc.Cell(fill = air, region = -detector_sphere)

Create Main Region

main_cell = openmc.Cell(name=‘source et fill’)
main_cell.fill = air
main_cell.region = -Top & +Bottom & -Yp & +Xm & +Ym & -Xp & +detector_sphere

root_universe.add_cells([main_cell,detector_cell])

Create Geometry and set root Universe

geometry = openmc.Geometry(root_universe)

Export to “geometry.xml”

geometry.export_to_xml()

cubic_mesh = openmc.RegularMesh()
cubic_mesh.dimension = [1, 1, 1]
cubic_mesh.lower_left = [-20, 80, 130] #Uncomment to simulate the case in which the source is outside the sphere
cubic_mesh.upper_right = [+20, 120, 170] #Uncomment to simulate the case in which the source is outside the sphere
#cubic_mesh.lower_left = [-20, -20, 130] #Uncomment to simulate the case in which the source is inside the sphere
#cubic_mesh.upper_right = [+20, 20, 170] #Uncomment to simulate the case in which the source is inside the sphere#
cubic_mesh_filter = openmc.MeshFilter(cubic_mesh)

spherical_mesh = openmc.SphericalMesh()
spherical_mesh.origin = [0,100,150] #Uncomment to simulate the case in which the source is outside the sphere
#spherical_mesh.origin = [0,0,150] #Uncomment to simulate the case in which the source is inside the sphere
spherical_mesh.r_grid = [0,25.]
spherical_mesh_filter = openmc.MeshFilter(spherical_mesh)

cellfilter = openmc.CellFilter(detector_cell)

Instantiate an empty Tallies object

tallies = openmc.Tallies()

Flux tally - Cube Tracklength

cube_track = openmc.Tally(name=‘Cube Tracklength’)
cube_track.filters = [cubic_mesh_filter]
cube_track.scores = [‘flux’]
cube_track.estimator = ‘tracklength’
tallies.append(cube_track)

Flux tally - Sphere Tracklength

sphere_track = openmc.Tally(name=‘Sphere Tracklength’)
sphere_track.filters = [spherical_mesh_filter]
sphere_track.scores = [‘flux’]
sphere_track.estimator = ‘tracklength’
tallies.append(sphere_track)

Flux tally - Cube Collision

cube_coll = openmc.Tally(name=‘Cube Collision’)
cube_coll.filters = [cubic_mesh_filter]
cube_coll.scores = [‘flux’]
cube_coll.estimator = ‘collision’
tallies.append(cube_coll)

Flux tally - Sphere Collision

sphere_coll = openmc.Tally(name=‘Sphere Collision’)
sphere_coll.filters = [spherical_mesh_filter]
sphere_coll.scores = [“flux”]
sphere_coll.estimator = ‘collision’
tallies.append(sphere_coll)
tallies.export_to_xml()

Flux tally - Spherical Cell Tracklength

cell_track = openmc.Tally(name=‘Cell Tracklength’)
cell_track.filters = [cellfilter]
cell_track.scores = [‘flux’]
cell_track.estimator = ‘tracklength’
tallies.append(cell_track)

Flux tally - Spherical Cell Collision

cell_coll = openmc.Tally(name=‘Cell Collision’)
cell_coll.filters = [cellfilter]
cell_coll.scores = [“flux”]
cell_coll.estimator = ‘collision’
tallies.append(cell_coll)
tallies.export_to_xml()

Source

src = openmc.Source()
src.space=openmc.stats.Point((0.0000, 0.0000, 150.))
src.energy = openmc.stats.Discrete([1.3e6],1)
src.strength = 1e10
src.particle =‘photon’
src.angle = openmc.stats.Isotropic()

Settings

settings = openmc.Settings()
settings.source = src
settings.seed = 2
settings.run_mode=‘fixed source’
settings.batches = 1000
settings.inactive = 20
settings.particles = 4000
settings.export_to_xml()

Run the model

openmc.run()

post_processing

with openmc.StatePoint(“statepoint.1000.h5”) as sp:

cube_track = sp.get_tally(name='Cube Tracklength')
sphere_track = sp.get_tally(name='Sphere Tracklength')
cube_coll = sp.get_tally(name='Cube Collision')
sphere_coll = sp.get_tally(name='Sphere Collision')
cell_track = sp.get_tally(name='Cell Tracklength')
cell_coll = sp.get_tally(name='Cell Collision')

sp.close()

print(“Tracklength Estimator - Sphere: " + str(sphere_track.mean[0,0,0]))
print(“Collision Estimator - Sphere: " + str(sphere_coll.mean[0,0,0]))
print(”\n\n”)
print(“Tracklength Estimator - Cube: " + str(cube_track.mean[0,0,0]))
print(“Collision Estimator - Cube: " + str(cube_coll.mean[0,0,0]))
print(”\n\n”)
print("Tracklength Estimator - Spherical Cell: " + str(cell_track.mean[0,0,0]))
print("Collision Estimator - Spherical Cell: " + str(cell_coll.mean[0,0,0]))

Is there any news on this topic ?

Seems to me to be a bug on this feature…