I am using openmc for fixed-source fusion-related simulations. I create the CAD in FreeCAD and convert it to DAGMC via cad_to_dagmc, which works great. I am currently running some simple centre-column tokamak simulations with a tokamak source created using openmc_plasma_source.
The simulation runs and everything works as expected except some of the tally visualisations. It seems like the axes have swapped places somehow. The 2D geometry plots are all representative of the CAD. My xy-plane tally visualisation should just be a cross section of a column, but I can only see that if I plot the YZ plane, which is very strange.
The xy geometry plot corresponds to YZ tally plot (YZ_mean = tally_values.mean[x_index, :, :]). The plot itself looks like it should.
Maybe I’ve done something strange in my code but I can’t figure it out. Any help would be appreciated!
I would be keen to see this flipping of the axis as described. Are you able to share either the model or a minimal viable example so I can reproduce
My first thoughts are that if you are visualizing a mesh tally then the results might need transposing or rotating to align the mesh slices to the correct axis.
Create a cylinder around the z-axis in some sort of CAD environment and convert to h5m with cad_to_dagmc (I haven’t tried it with just CSG, but I assume it would give the same answer?)
Create a plasma source around the cylinder.
Plot the geometry on all 2D slices. For example:
xy:
xz:
yz:
Run a quick model with some tally, say heating or flux, on a regular, 3D mesh.
Plot the tally on each plane, this is how I did it:
So the YZ-plot is what you would expect the XY-plot to look like in my mind, especially if you compare them to the 2D geometry plots.
Maybe I’m somehow plotting it wrong?
A side note on this is that maybe it would be better to use a cylindrical mesh for this, but I struggle to index and plot them properly once I’ve done the simulation.
Amazing, thanks for the extra effort of making the example. This will be useful for anyone stumbling across this issue in the future and helped me understand the problem.
Ah yes this looks like a familiar issue. I have some code in the openmc_plot that you are welcome to recycle. It flips the mesh tally slice to the correct axis.
This section of code flattens the tally results, then transposes them rotates the slice. Each axis requred different transposes and different rotations. The resulting mpl_image_slice is then sent to plt.imshow and it appears to be aligned correctly.
I hope to add combined regular mesh tally and DAGMC plotting with mesh tallies so the DAGMC geometry can be see on the mesh to openmc_plot soon if it is something that is useful to the community?
Thank you for that! I’m glad to hear it’s “just” a transposing issue an nothing deeper. I’ll have a look at the axes and will probably recycle some of your example code to transpose.
For me, being able to plot both the DAGMC and the regular mesh on the same plot would be very helpful for sure.
This is for anyone who might come across this in the future.
I’ve been looking at this recently and solved it slightly differently. I just figured out which axes that had been swapped and then selected the corresponding indexes. It’s very inspired by your solution but a bit more basic, Shimwell. It’s a part of a post-processing class where we have the DataFrame of the tally, the mesh and mesh size etc:
def plot_tallies(self):
"""
Plots tallies by score using the tally DataFrame.
"""
for score in self.tallies['mesh_tally'].scores:
mean = self.df['mean'].to_numpy().reshape(self.mesh.dimension)
axes_to_slice = ['X', 'Y', 'Z']
for axis in axes_to_slice:
self.plot_slice(mean=mean, axis_to_slice=axis, score=score)
def plot_slice(self, mean, axis_to_slice, score):
"""
Method to plot a mesh tally slice.
When we reshape the tally results from a flattened array to a 3D array we
get the wrong order of the axes. This method takes that into consideration
by selecting the right indexes and allows to plot the proper planes (eg XY, XZ or YZ)
Args:
values (np array): 3D np array of values of each cell in the mesh (could be mean or std. dev.)
axis (str): which axis to plot
score (str): which score it is
"""
if axis_to_slice == "X":
bb_index = [1, 2]
slice_index = int(self.mesh_dimension/2)
image = mean[:, :, slice_index]
x_label = "Y [cm]"
y_label = "Z [cm]"
end = '_YZ'
elif axis_to_slice == "Y":
bb_index = [0, 2]
slice_index = int(self.mesh_dimension/2)
image = mean[:, slice_index, :]
x_label = "X [cm]"
y_label = "Z [cm]"
end = '_XZ'
elif axis_to_slice == 'Z':
bb_index = [0, 1]
slice_index = int(self.mesh_dimension/2)
image = mean[slice_index, :, :]
x_label = "X [cm]"
y_label = "Y [cm]"
end = '_XY'
else:
raise ValueError('Axis needs to be X, Y or Z')
# Plot and save names
plot_name = 'mesh_' + score + end
save_name = plot_name + '.png'
# Imshow extent
left = self.mesh.lower_left[bb_index[0]]
right = self.mesh.upper_right[bb_index[0]]
bottom = self.mesh.lower_left[bb_index[1]]
top = self.mesh.upper_right[bb_index[1]]
extent = (left, right, bottom, top)
# Plot
plt.title(plot_name)
plt.imshow(image, interpolation='None', cmap='jet', extent=extent, aspect='auto')
plt.xlabel(x_label)
plt.ylabel(y_label)
cbar = plt.colorbar()
cbar.set_label(score + ' per source particle')
plt.savefig(save_name)
plt.close()
Just to complicate things the heat map plot in matplotlib and plotly differ a bit and expect the slice in a different order to each other. But I don’t think we should worry about that as OpenMC has matplotlib as a dependency we could adopt that slicing style.
Perhaps we should think about adding some helper functions to the regular mesh tally in the openmc source code.
@maarten I’ve made a PR to openmc to see if this sort of plotting functionality would be useful in the source code. We would value your input on the PR if you have time. The PR has a method for getting a slice of data and there is also a plotting method
