OpenMC Guide
Plotting Results
Why Plot Your Results
Visualization serves two purposes in OpenMC: verifying geometry before running simulations, and understanding results after the calculation completes. OpenMC provides built-in plotting for geometry verification and integrates with Python visualization libraries for results analysis.
Geometry Verification
Plot geometry before running simulations. A 2D slice can reveal errors that would cause runtime failures or incorrect results.
python
import openmc
import matplotlib.pyplot as plt
# Create a basic geometry plot
plot = openmc.Plot()
plot.basis = 'xy' # View from above (z-direction)
plot.origin = (0, 0, 0) # Center point
plot.width = (5.0, 5.0) # 5cm × 5cm view area
plot.pixels = (400, 400) # Resolution
plot.color_by = 'material' # Color by material type
# Generate the plot file
plots = openmc.Plots([plot])
plots.export_to_xml()
# Run plotting (creates .png files)
openmc.plot_geometry()Tutorial snippet — no separate file in examples repo
This creates a .png file showing the geometry. Check for overlaps, missing materials, and correct region assignments.
Multiple Views
For 3D geometries, create multiple orthogonal slice views:
python
# Create three orthogonal views
xy_plot = openmc.Plot(name='xy_view')
xy_plot.basis = 'xy'
xy_plot.origin = (0, 0, 0)
xy_plot.width = (10, 10)
xy_plot.pixels = (500, 500)
xy_plot.color_by = 'material'
xz_plot = openmc.Plot(name='xz_view')
xz_plot.basis = 'xz'
xz_plot.origin = (0, 0, 0)
xz_plot.width = (10, 100) # Different scale for height
xz_plot.pixels = (500, 500)
xz_plot.color_by = 'material'
yz_plot = openmc.Plot(name='yz_view')
yz_plot.basis = 'yz'
yz_plot.origin = (0, 0, 0)
yz_plot.width = (10, 100)
yz_plot.pixels = (500, 500)
yz_plot.color_by = 'material'
plots = openmc.Plots([xy_plot, xz_plot, yz_plot])
plots.export_to_xml()
openmc.plot_geometry()Tutorial snippet — no separate file in examples repo
Plotting Simulation Results
Use matplotlib to visualize simulation results:
python
import numpy as np
import matplotlib.pyplot as plt
# Load simulation results
sp = openmc.StatePoint('statepoint.100.h5')
# Get flux tally
flux_tally = sp.get_tally(name='flux')
flux_mean = flux_tally.mean.flatten()
flux_std = flux_tally.std_dev.flatten()
# Basic flux plot
plt.figure(figsize=(10, 6))
plt.errorbar(range(len(flux_mean)), flux_mean, yerr=flux_std,
fmt='o-', capsize=3, markersize=4)
plt.xlabel('Tally Bin')
plt.ylabel('Flux (n/cm²/src)')
plt.title('Neutron Flux Distribution')
plt.grid(True, alpha=0.3)
plt.show()
# Log scale for wide-ranging data
plt.figure(figsize=(10, 6))
plt.semilogy(flux_mean, 'o-')
plt.xlabel('Position')
plt.ylabel('Flux (n/cm²/src)')
plt.title('Flux Distribution (Log Scale)')
plt.grid(True, alpha=0.3)
plt.show()Tutorial snippet — no separate file in examples repo
Energy Spectrum Plots
For energy-dependent tallies, spectrum plots show the neutron energy distribution:
python
# Load energy-dependent flux tally
energy_tally = sp.get_tally(name='energy_flux')
df = energy_tally.get_pandas_dataframe()
# Extract energy filter
energy_filter = energy_tally.filters[0] # Assuming first filter is energy
bin_edges = np.append(energy_filter.bins[:, 0], energy_filter.bins[-1, 1])
# Calculate bin centers and lethargy widths
energy_centers = np.sqrt(bin_edges[:-1] * bin_edges[1:]) # Geometric mean
lethargy_widths = np.log(bin_edges[1:] / bin_edges[:-1])
# Compute flux per unit lethargy
flux_vals = df['mean'].values
flux_per_lethargy = flux_vals / lethargy_widths
# Plot spectrum
plt.figure(figsize=(10, 6))
plt.loglog(energy_centers, flux_per_lethargy, 'o-', linewidth=2)
plt.xlabel('Energy (eV)')
plt.ylabel('Flux per unit lethargy')
plt.title('Neutron Energy Spectrum')
plt.grid(True, alpha=0.3)
plt.show()
# Add error bars if desired
std_per_lethargy = df['std. dev.'].values / lethargy_widths
plt.figure(figsize=(10, 6))
plt.errorbar(energy_centers, flux_per_lethargy, yerr=std_per_lethargy,
fmt='o-', capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Energy (eV)')
plt.ylabel('Flux per unit lethargy')
plt.title('Neutron Energy Spectrum with Uncertainties')
plt.grid(True, alpha=0.3)
plt.show()Tutorial snippet — no separate file in examples repo
Spatial Distribution Plots
Spatial tallies can be visualized as heat maps:
python
# Load spatial mesh tally
mesh_tally = sp.get_tally(name='mesh_flux')
df = mesh_tally.get_pandas_dataframe()
# Reshape data to 2D array (assuming rectangular mesh)
mesh_filter = mesh_tally.filters[0]
mesh = mesh_filter.mesh
x_grid = np.linspace(mesh.lower_left[0], mesh.upper_right[0], mesh.dimension[0] + 1)
y_grid = np.linspace(mesh.lower_left[1], mesh.upper_right[1], mesh.dimension[1] + 1)
x_bins = mesh.dimension[0]
y_bins = mesh.dimension[1]
flux_2d = df['mean'].values.reshape((y_bins, x_bins))
# Create heat map
plt.figure(figsize=(10, 8))
plt.imshow(flux_2d, origin='lower', aspect='equal',
extent=[x_grid[0], x_grid[-1],
y_grid[0], y_grid[-1]])
plt.colorbar(label='Flux (n/cm²/src)')
plt.xlabel('X Position (cm)')
plt.ylabel('Y Position (cm)')
plt.title('Spatial Flux Distribution')
plt.show()
# Contour plot alternative
plt.figure(figsize=(10, 8))
x_centers = [(x_grid[i] + x_grid[i+1])/2 for i in range(x_bins)]
y_centers = [(y_grid[i] + y_grid[i+1])/2 for i in range(y_bins)]
X, Y = np.meshgrid(x_centers, y_centers)
levels = np.logspace(np.log10(flux_2d.min()), np.log10(flux_2d.max()), 20)
plt.contourf(X, Y, flux_2d, levels=levels, extend='both')
plt.colorbar(label='Flux (n/cm²/src)')
plt.xlabel('X Position (cm)')
plt.ylabel('Y Position (cm)')
plt.title('Flux Contours')
plt.show()Tutorial snippet — no separate file in examples repo
Quick Plotting Tips
Utility functions for common plotting tasks:
python
# Always check your geometry first
def quick_geometry_check(geometry):
plot = openmc.Plot()
plot.basis = 'xy'
plot.origin = (0, 0, 0)
plot.width = (20, 20) # Adjust to your model size
plot.pixels = (400, 400)
plot.color_by = 'material'
plot.filename = 'geometry_check'
plots = openmc.Plots([plot])
plots.export_to_xml()
openmc.plot_geometry()
print("Check geometry_check.png before running simulation!")
# Quick result summary
def plot_summary(statepoint_file):
sp = openmc.StatePoint(statepoint_file)
# K-effective plot if eigenvalue calculation
if hasattr(sp, 'keff'):
keff_gen = sp.k_generation
plt.figure(figsize=(10, 6))
plt.plot(keff_gen, 'b-', alpha=0.7)
plt.axhline(sp.keff.nominal_value, color='r', linestyle='--')
plt.xlabel('Generation')
plt.ylabel('k-effective')
plt.title('k-effective Convergence')
plt.grid(True, alpha=0.3)
plt.show()
print(f"k-effective: {sp.keff.nominal_value:.5f} ± {sp.keff.std_dev:.5f}")
# Usage examples
quick_geometry_check(geometry)
plot_summary('statepoint.100.h5')Tutorial snippet — no separate file in examples repo
If plots don't make physical sense, check the model setup before running expensive calculations.