openmc-dev / plotter

Native plotting GUI for model design and verification
MIT License
48 stars 18 forks source link

bug report: plot mesh with EnergyFunctionFilter shows zeros everywhere #106

Closed shimwell closed 1 year ago

shimwell commented 1 year ago

Really liking the plotting capabilities of this package

I think I stubbled across a small bug

When plotting regular mesh tallies that include a EnergyFunctionFilter the displayed result is all zeros.

This would be a useful feature as plotting dose maps could then be performed

shimwell commented 1 year ago

Screenshot from 2023-01-12 16-03-16

shimwell commented 1 year ago

minimal code to reproduce

import openmc

mat_1 = openmc.Material()
mat_1.add_nuclide('Li6', 1)
mat_1.set_density('g/cm3', 1)
my_materials = openmc.Materials([mat_1])

surface_1 = openmc.Sphere(r=100, boundary_type='vacuum')
region_1 = -surface_1
cell_1 = openmc.Cell(region=region_1)
cell_1.fill = mat_1

my_geometry = openmc.Geometry([cell_1])

reg_mesh = openmc.RegularMesh().from_domain(
    my_geometry,
    dimension=[5, 5, 5]
)
reg_mesh_filter = openmc.MeshFilter(reg_mesh)

energy_bins_n, dose_coeffs_n = openmc.data.dose_coefficients(
    particle="neutron", geometry="AP"
)
energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n)
energy_function_filter_n.interpolation == "cubic"

neutron_particle_filter = openmc.ParticleFilter("neutron")

n_reg_mesh_tally = openmc.Tally(name="neutron_dose_on_cell")
# note that the EnergyFunctionFilter is included as a filter
n_reg_mesh_tally.filters = [
    neutron_particle_filter,
    energy_function_filter_n,
    reg_mesh_filter
]
n_reg_mesh_tally.id = 1
n_reg_mesh_tally.scores = ["flux"]
my_tallies = openmc.Tallies([n_reg_mesh_tally])

my_settings = openmc.Settings()

my_settings.batches = 10
my_settings.particles = 100000
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.0, 0.0))

my_settings.source = source

my_model = openmc.Model(settings=my_settings, geometry=my_geometry, tallies=my_tallies, materials=my_materials)

sp_filename = my_model.run()

sp = openmc.StatePoint(sp_filename)

# shows that the tally is not zero
print(sp.tallies[1].mean)
shimwell commented 1 year ago

changing this line

https://github.com/openmc-dev/plotter/blob/4f353ccee8645caf04fb6c5d71cd6136183d50d5/openmc_plotter/plotmodel.py#L695

to

 data = data[np.array(selected_bins)].sum(axis=0) [0]

and this line

https://github.com/openmc-dev/plotter/blob/4f353ccee8645caf04fb6c5d71cd6136183d50d5/openmc_plotter/plotmodel.py#L689

to

if type(tally_filter) in [openmc.MeshFilter, openmc.EnergyFunctionFilter]:

is a temporary fix

letibdesneiges commented 1 year ago

Hi,

I have encountered the same problem when trying to plot doserate data from a flux mesh tally with an EnergyFunctionFilter converting flux to doserate. All values are displayed as zero. The above solution did not help unfortunately.

Attached the statePoint file. statepoint.200.h5.zip

openmc-plotter v. 0.3.1 installed via conda on Debian 11.

letibdesneiges commented 1 year ago

Implementing @shimwell changes to plotmodel.py and changing:

https://github.com/openmc-dev/plotter/blob/eeec020b1663b020efbfd5654e4554ded2d2a18a/openmc_plotter/plotmodel.py#L726-L727

to

# set image data, reverse y-axis
data=data.squeeze()
image_data = data[::-1, ...]

solved it for my specific case.