openmc-dev / openmc

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

spherical mesh tally bug #2414

Open shimwell opened 1 year ago

shimwell commented 1 year ago

During PR #2397 we think that we found a bug with the way that spherical mesh tallies are recorded. Further details in the PR

shimwell commented 1 year ago

I've been running this script to test the spherical mesh without normalisation


import openmc
import math
import numpy as np
surf1 = openmc.Sphere(r=10, boundary_type='vacuum')

cell1 = openmc.Cell(region=-surf1)

my_geometry = openmc.Geometry([cell1])

my_settings = openmc.Settings()
batches = 2
my_settings.batches = batches
my_settings.inactive = 0
my_settings.particles = 500000
my_settings.run_mode = 'fixed source'

source = openmc.Source()
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14e6], [1])
source.space = openmc.stats.Point((0, 0, 0.1))

my_settings.source = source

mesh = openmc.SphericalMesh()
mesh.r_grid = (0,1,2,3,4,5,6,7,8,9,10)
mesh.theta_grid = np.linspace(0, math.pi, 20)
mesh.phi_grid = np.linspace(0, 2*math.pi, 20)
mesh_filter = openmc.MeshFilter(mesh)

mesh_tally_1 = openmc.Tally(name='flux_on_mesh')
mesh_tally_1.filters = [mesh_filter]
mesh_tally_1.scores = ['flux']  # where X is a wildcard
my_tallies = openmc.Tallies([mesh_tally_1])

model = openmc.Model(my_geometry, openmc.Materials([]), my_settings, my_tallies)
sp_filename = model.run()

statepoint = openmc.StatePoint(sp_filename)

my_mesh_tally = statepoint.get_tally(name='flux_on_mesh')

mesh.write_data_to_vtk(
    filename="sphere_no_norm_fixPR.vtk",
    datasets={"mean": my_mesh_tally.mean},
    volume_normalization=False
)

The sphere_no_norm_fixPR.vtk shows odd results. It looks like I'm getting a stripe of values across one axis so I thought this might be a systematic bug

Screenshot from 2023-03-01 16-38-08 Screenshot from 2023-03-01 16-37-55

shimwell commented 1 year ago

closing as this is a duplicate of #2416

pshriwise commented 1 year ago

I think this is different than #2416. This effect is still apparent after the fix applied in #2417