OpenMC Guide
Plotting Results
Why Plot Your Results
Visualization serves two critical purposes in OpenMC: verifying your geometry is correct before running expensive simulations, and understanding your results after the calculation completes. Both save time and prevent errors.
OpenMC provides built-in plotting for geometry verification and integrates well with Python visualization libraries for results analysis.
Geometry Verification
Always plot your 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()This creates a .png file showing your geometry. Check for overlaps, missing materials, and verify your model looks as expected.
Multiple Views
For complex 3D geometries, create multiple slice views to understand the full structure:
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()Plotting Simulation Results
After your simulation completes, use matplotlib and other Python tools to visualize the 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()Energy Spectrum Plots
For energy-dependent tallies, create spectrum plots to understand 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
energy_bins = energy_filter.bins
# Calculate bin centers for plotting
energy_centers = []
for i in range(len(energy_bins) - 1):
center = np.sqrt(energy_bins[i] * energy_bins[i+1]) # Geometric mean
energy_centers.append(center)
# Plot spectrum
plt.figure(figsize=(10, 6))
plt.loglog(energy_centers, df['mean'], '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
plt.figure(figsize=(10, 6))
plt.errorbar(energy_centers, df['mean'], yerr=df['std. dev.'],
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()Spatial Distribution Plots
For spatial tallies, create heat maps to visualize distributions across your geometry:
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_bins = len(mesh.x_grid) - 1
y_bins = len(mesh.y_grid) - 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=[mesh.x_grid[0], mesh.x_grid[-1],
mesh.y_grid[0], mesh.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 = [(mesh.x_grid[i] + mesh.x_grid[i+1])/2 for i in range(x_bins)]
y_centers = [(mesh.y_grid[i] + mesh.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()Quick Plotting Tips
Here are practical tips for effective plotting:
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.keff_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')Remember: if your plots don't make physical sense, check your model setup before running expensive calculations. Visualization is your best debugging tool.