Not sure if this is an issue with cylindrical meshes, or simply with write_data_to_vtk() (or with the manner in which I extract the data from the statepoint, or something else that I am doing), but I have noticed some strange behavior on cylindrical meshes, and more while making the example for this post. The example script for making the following results is at the end of the post. The model uses an isotropic point source at (750,750,0). OpenMC version 0.14.0 installed via conda.
First, some strange results on a mesh with a large radius and small phi extent, defined as follows.
mesh = openmc.CylindricalMesh(
r_grid=np.linspace(76327, 77077, 66),
z_grid=np.linspace(-540, 540, 86),
phi_grid=np.linspace(0.771, 0.801, 109),
origin=np.array([-53604, -53604, 0])
)
It appears that no mesh elements that are not on the same radial coordinate are scoring. See screenshots
If I change the mesh to cover roughly the same domain, but instead with the origin at (0,0,0), the results are a bit different
mesh = openmc.CylindricalMesh(
r_grid=np.linspace(490, 1585, 66),
z_grid=np.linspace(-570, 560, 86),
phi_grid=np.linspace(0, np.pi/2, 109)
)
The flux profile on a plane on the z axis and the point source looks promising (log scale)
But looking at another angle things are still weird
Removing all the zero score mesh elements:
This second one made me think perhaps my method of extracting data for write_data_to_vtk() is flawed. I use get_reshaped_data() and turn it into a 1D array of appropriate dimension for the mesh, but perhaps that is not correct?
To test this I used the volumes variable on the mesh, which is a 3 dimensional array. I first tried just adding this array as another dataset (turning off volume normalization), and this worked, but produced unexpected results, with a gradient along the z axis rather than in the radial direction.
I also tried the same thing, but using np.flatten() on the volume array, and got an identical result when writing to vtk.
Interestingly, transposing the volume matrix prior to writing to vtk produces the expected results (maybe this is expected behavior?)
Perhaps I need to get 3D matrix of the mesh tally results and transpose that prior to writing the data to vtk as I did with the volumes? (though I wouldn’t expect the behavior I am seeing in the screenshots above if the matrix just needs to be transposed prior to flattening it).
Sorry this turned into a bit of an essay, any help or insight would be appreciated, thanks! hoping this is just user error on my part 
As promised, the example script.
import openmc
import numpy as np
model = openmc.Model()
"""# this mesh has badly behaved results
mesh = openmc.CylindricalMesh(
r_grid=np.linspace(76327, 77077, 66),
z_grid=np.linspace(-540, 540, 86),
phi_grid=np.linspace(0.771, 0.801, 109),
origin=np.array([-53604, -53604, 0])
)
"""
# this mesh is badly behaved but differently I think
mesh = openmc.CylindricalMesh(
r_grid=np.linspace(490, 1585, 66),
z_grid=np.linspace(-570, 560, 86),
phi_grid=np.linspace(0, np.pi/2, 109)
)
source = openmc.IndependentSource()
source.space = openmc.stats.Point([750,750,0])
source.energy = openmc.stats.Discrete([14.1e6], [1])
source.angle = openmc.stats.Isotropic()
settings = openmc.Settings()
settings.particles = 100000
settings.batches = 15
settings.source = source
settings.run_mode = 'fixed source'
model.settings = settings
vac_surf = openmc.Sphere(r=10000, boundary_type = 'vacuum')
vac_cell = openmc.Cell(region=-vac_surf)
geometry = openmc.Geometry([vac_cell])
model.geometry = geometry
mesh_filter = openmc.MeshFilter(mesh)
flux_tally = openmc.Tally(name='flux tally')
flux_tally.filters = [mesh_filter]
flux_tally.scores = ['flux']
model.tallies = [flux_tally]
statepoint = model.run()
with openmc.StatePoint(statepoint) as sp:
flux_tally = sp.get_tally(name='flux tally')
print(sp.meshes)
mesh = sp.meshes[1]
volumes = mesh.volumes
mean = flux_tally.get_reshaped_data('mean').flatten()
mesh.write_data_to_vtk(
filename = 'flux.vtk',
datasets = {
'flux':mean,
'volumes':volumes.T.flatten()
},
volume_normalization = False
)