OpenMC Guide
Data Processing in OpenMC
What you'll learn
After the first pin · 5 / 1115 min read
- Convert per-source-particle tallies into absolute reaction rates at a given power.
- Reshape a mesh tally back into the geometry it was scored on.
- Collapse a spectrum into thermal, epithermal, and fast groups.
- Check uncertainties before quoting any number.
Before you start
Getting Started with Results
Statepoint files contain all tally results from a simulation. OpenMC provides methods to extract this data in formats suitable for analysis.
Basic Result Extraction
python
import openmc
import numpy as np
import matplotlib.pyplot as plt
# Load simulation results
sp = openmc.StatePoint('statepoint.100.h5')
# Quick access to k-effective
keff = sp.keff
print(f"k-effective: {keff.nominal_value:.5f} ± {keff.std_dev:.5f}")
# Get a specific tally by name
flux_tally = sp.get_tally(name='fuel_flux')
# Access the raw mean and standard deviation arrays
mean_values = flux_tally.mean.flatten()
std_dev_values = flux_tally.std_dev.flatten()
print(f"Average flux: {np.mean(mean_values):.3e}")
print(f"Maximum flux: {np.max(mean_values):.3e}")
print(f"Relative uncertainty: {np.mean(std_dev_values/mean_values)*100:.2f}%")Tutorial snippet — no separate file in examples repo
Working with Different Tally Types
Different tally types require different processing approaches.
Cell-Based Tallies
python
# Cell tallies are the simplest - one value per cell per score.
# Substitute your own tally name; 'fuel_cells' here is a stand-in for
# whatever you called the CellFilter tally when you built the model.
sp = openmc.StatePoint('statepoint.100.h5')
cell_tally = sp.get_tally(name='fuel_cells')
# Extract data by score
df = cell_tally.get_pandas_dataframe()
# Filter by specific scores
flux_data = df[df['score'] == 'flux']
fission_data = df[df['score'] == 'fission']
# Calculate reaction rates (tallies are per source particle)
total_fission_rate = fission_data['mean'].sum()
average_flux = flux_data['mean'].mean()
print(f"Total fission rate: {total_fission_rate:.3e} reactions/source")
print(f"Average flux: {average_flux:.3e} n-cm/source")
# nu-fission / absorption in the tallied region. Without an EnergyFilter this
# is integrated over the whole spectrum, so it is NOT the four-factor eta -
# that one is defined on the thermal group alone.
absorption_data = df[df['score'] == 'absorption']
nu_fission_data = df[df['score'] == 'nu-fission']
ratio = nu_fission_data['mean'].sum() / absorption_data['mean'].sum()
print(f"nu-fission / absorption: {ratio:.3f}")Tutorial snippet — no separate file in examples repo
All OpenMC tallies are reported per source particle. To convert to absolute physical units, multiply by the source strength in neutrons per second. For a fission reactor operating at power P, the source strength is S = P·ν / (keff·Q), where ν is the average neutrons per fission and Q ≈ 200 MeV is the recoverable energy per fission.
Mesh Tallies
python
# Mesh tallies give spatial distributions. 'power_mesh' stands in for
# your own MeshFilter tally name.
mesh_tally = sp.get_tally(name='power_mesh')
mesh_filter = mesh_tally.find_filter(openmc.MeshFilter)
mesh = mesh_filter.mesh
# Get 3D power data (kappa-fission = energy deposited per fission)
power_data = mesh_tally.get_values(scores=['kappa-fission'])
shape = mesh.dimension
power_3d = power_data.reshape(shape)
# Find hot spots
max_power = np.max(power_3d)
max_location = np.unravel_index(np.argmax(power_3d), power_3d.shape)
print(f"Peak power: {max_power:.3e} at mesh indices {max_location}")
# Calculate power peaking factors
avg_power = np.mean(power_3d[power_3d > 0]) # Exclude zero regions
peaking_factor = max_power / avg_power
print(f"Power peaking factor: {peaking_factor:.2f}")
# Extract 2D slice for plotting
midplane = power_3d[:, :, power_3d.shape[2]//2]
plt.figure(figsize=(8, 6))
plt.imshow(midplane, origin='lower', cmap='hot')
plt.colorbar(label='Power Density')
plt.title('Power Distribution (Midplane)')
plt.show()Tutorial snippet — no separate file in examples repo
Energy Spectrum Analysis
When plotting neutron spectra, raw energy-binned flux is misleading because wider bins collect more flux regardless of spectral shape. Dividing by the lethargy width of each bin (Δu = ln(Ehigh/Elow)) removes this bin-width artifact and reveals the true spectral shape. On a log-energy axis, lethargy-normalized flux makes the thermal Maxwellian peak and the fast fission spectrum peak clearly visible — without this normalization, both features are obscured by the logarithmic bin spacing.
python
# Energy tallies show neutron spectra. 'spectrum' stands in for your
# own EnergyFilter tally name.
spectrum_tally = sp.get_tally(name='spectrum')
energy_filter = spectrum_tally.find_filter(openmc.EnergyFilter)
bin_edges = np.append(energy_filter.bins[:, 0], energy_filter.bins[-1, 1])
# Get flux spectrum
df = spectrum_tally.get_pandas_dataframe()
flux_spectrum = df[df['score'] == 'flux']['mean'].values
# Calculate energy group fractions
total_flux = np.sum(flux_spectrum)
thermal_idx = np.searchsorted(bin_edges, 0.625) # 0.625 eV cutoff
fast_idx = np.searchsorted(bin_edges, 100e3) # 100 keV cutoff
thermal_fraction = np.sum(flux_spectrum[:thermal_idx]) / total_flux
epithermal_fraction = np.sum(flux_spectrum[thermal_idx:fast_idx]) / total_flux
fast_fraction = np.sum(flux_spectrum[fast_idx:]) / total_flux
print(f"Thermal flux fraction: {thermal_fraction:.3f}")
print(f"Epithermal flux fraction: {epithermal_fraction:.3f}")
print(f"Fast flux fraction: {fast_fraction:.3f}")
# Plot spectrum (flux per unit lethargy)
energy_mid = np.sqrt(bin_edges[:-1] * bin_edges[1:])
lethargy_widths = np.log(bin_edges[1:] / bin_edges[:-1])
flux_per_lethargy = flux_spectrum / lethargy_widths
plt.figure(figsize=(10, 6))
plt.loglog(energy_mid, flux_per_lethargy, 'b-', linewidth=2)
plt.xlabel('Energy (eV)')
plt.ylabel('Flux per unit lethargy')
plt.title('Neutron Energy Spectrum')
plt.grid(True, which='both', alpha=0.3)
plt.show()Tutorial snippet — no separate file in examples repo
Statistical Analysis and Uncertainty
Monte Carlo results carry statistical uncertainty. Always assess convergence and report confidence intervals.
Checking Convergence
python
# Load statepoint to check convergence
sp = openmc.StatePoint('statepoint.100.h5')
# K-effective convergence
keff_generations = sp.k_generation
active_batches = len(keff_generations)
# Check if k-effective has converged
recent_keff = keff_generations[-20:] # Last 20 batches
keff_std = np.std(recent_keff)
keff_mean = np.mean(recent_keff)
relative_std = keff_std / keff_mean
print(f"K-effective: {keff_mean:.5f}")
print(f"Recent standard deviation: {keff_std:.5f}")
print(f"Relative standard deviation: {relative_std*100:.3f}%")
if relative_std < 0.001: # Less than 0.1%
print("K-effective appears converged")
else:
print("Consider running more batches")
# Plot k-effective vs generation
plt.figure(figsize=(10, 6))
plt.plot(range(1, active_batches + 1), keff_generations, 'b-', alpha=0.7)
plt.axhline(keff_mean, color='r', linestyle='--', label=f'Mean: {keff_mean:.5f}')
plt.xlabel('Generation')
plt.ylabel('k-effective')
plt.title('K-effective Convergence')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()Tutorial snippet — no separate file in examples repo
Tally Uncertainty Analysis
python
# Analyze tally uncertainties
tally = sp.get_tally(name='fuel_flux')
df = tally.get_pandas_dataframe()
# Calculate relative uncertainties
flux_data = df[df['score'] == 'flux']
mean_flux = flux_data['mean'].values
std_flux = flux_data['std. dev.'].values
rel_error = std_flux / mean_flux
# Check for acceptable uncertainties
acceptable_threshold = 0.05 # 5%
high_uncertainty = rel_error > acceptable_threshold
print(f"Average relative uncertainty: {np.mean(rel_error)*100:.2f}%")
print(f"Maximum relative uncertainty: {np.max(rel_error)*100:.2f}%")
print(f"Tallies with >5% uncertainty: {np.sum(high_uncertainty)}")
# Calculate 95% confidence intervals
confidence_level = 1.96 # 95% confidence for normal distribution
lower_bound = mean_flux - confidence_level * std_flux
upper_bound = mean_flux + confidence_level * std_flux
print(f"\nSample flux result with 95% confidence:")
idx = 0 # First tally bin
print(f"Flux: {mean_flux[idx]:.3e} [{lower_bound[idx]:.3e}, {upper_bound[idx]:.3e}]")Tutorial snippet — no separate file in examples repo
Automation and Batch Processing
Parameter studies and design optimization require automated processing of multiple simulations.
Simple Batch Analysis
python
import os
import pandas as pd
def analyze_statepoint(filepath, tally_name='fuel_flux'):
"""Extract key results from a single statepoint file."""
try:
sp = openmc.StatePoint(filepath)
# Get k-effective
keff = sp.keff.nominal_value
keff_std = sp.keff.std_dev
# Get flux tally
tally = sp.get_tally(name=tally_name)
df = tally.get_pandas_dataframe()
flux_data = df[df['score'] == 'flux']
# Calculate summary statistics
total_flux = flux_data['mean'].sum()
max_flux = flux_data['mean'].max()
avg_uncertainty = (flux_data['std. dev.'] / flux_data['mean']).mean()
return {
'filepath': filepath,
'keff': keff,
'keff_std': keff_std,
'total_flux': total_flux,
'max_flux': max_flux,
'avg_uncertainty': avg_uncertainty,
'status': 'success'
}
except Exception as e:
return {
'filepath': filepath,
'status': 'failed',
'error': str(e)
}
# Process multiple cases
case_dirs = ['case_1', 'case_2', 'case_3', 'case_4']
results = []
for case_dir in case_dirs:
statepoint_path = f"{case_dir}/statepoint.100.h5"
if os.path.exists(statepoint_path):
result = analyze_statepoint(statepoint_path)
result['case'] = case_dir
results.append(result)
if result['status'] == 'success':
print(f"Processed {case_dir}: k-eff = {result['keff']:.5f} ± {result['keff_std']:.5f}")
else:
print(f"{case_dir} failed: {result['error']}")
else:
print(f"Skipping {case_dir}: file not found")
# Create summary DataFrame
df_results = pd.DataFrame(results)
successful_cases = df_results[df_results['status'] == 'success']
# Save results
successful_cases.to_csv('batch_analysis_results.csv', index=False)
print(f"\nAnalyzed {len(successful_cases)} cases successfully")
print(f"Results saved to batch_analysis_results.csv")Tutorial snippet — no separate file in examples repo
Parameter Study Analysis
python
import os
import numpy as np
# Analyze results as a function of a parameter
# Assumes you have a parameter study with different enrichments
enrichments = [0.03, 0.035, 0.04, 0.045, 0.05]
keff_values = []
keff_errors = []
for enrichment in enrichments:
case_dir = f"enrichment_{enrichment:.3f}"
statepoint_path = f"{case_dir}/statepoint.100.h5"
if os.path.exists(statepoint_path):
sp = openmc.StatePoint(statepoint_path)
keff_values.append(sp.keff.nominal_value)
keff_errors.append(sp.keff.std_dev)
else:
keff_values.append(np.nan)
keff_errors.append(np.nan)
# Convert to numpy arrays for analysis
enrichments = np.array(enrichments)
keff_values = np.array(keff_values)
keff_errors = np.array(keff_errors)
# Remove any failed cases
valid_mask = ~np.isnan(keff_values)
enrichments = enrichments[valid_mask]
keff_values = keff_values[valid_mask]
keff_errors = keff_errors[valid_mask]
# Plot results
plt.figure(figsize=(10, 6))
plt.errorbar(enrichments * 100, keff_values, yerr=keff_errors,
marker='o', linewidth=2, capsize=5)
plt.axhline(1.0, color='r', linestyle='--', alpha=0.7, label='Critical')
plt.xlabel('Enrichment (%)')
plt.ylabel('k-effective')
plt.title('Criticality vs. Fuel Enrichment')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()
# Find critical enrichment
critical_idx = np.argmin(np.abs(keff_values - 1.0))
critical_enrichment = enrichments[critical_idx]
print(f"Closest to critical: {critical_enrichment*100:.1f}% enrichment")
print(f"k-effective: {keff_values[critical_idx]:.5f}")Tutorial snippet — no separate file in examples repo
Practical Tips and Best Practices
Data Management
A statepoint file is the only durable record of a run, and it is easy to lose track of which one came from which set of inputs. Name them for the case rather than the batch count, keep a parameter study in its own directory, and back them up — regenerating one means paying the CPU time again. The analysis script that reads them deserves the same treatment: a notebook that produced a figure six months ago is worthless if nobody can tell which statepoint it opened.
File size is a useful diagnostic in its own right. A statepoint much larger than you expected usually means a tally with more bins than you intended — a mesh crossed with an energy grid crossed with several nuclides multiplies out quickly — and that same combinatorial blowup is what makes the run slow. Noticing it in the file listing is cheaper than noticing it in the wall clock.
Common Pitfalls
python
# Always check for convergence before trusting results
sp = openmc.StatePoint('statepoint.100.h5')
# 1. Check if simulation ran to completion
if sp.n_batches < 100:
print(f"Warning: Only {sp.n_batches} batches completed")
# 2. Verify sufficient inactive batches
if sp.n_inactive < 10:
print("Warning: Few inactive batches may affect convergence")
# 3. Check tally uncertainty (name your own tally here)
tally = sp.get_tally(name='important_tally')
df = tally.get_pandas_dataframe()
rel_errors = df['std. dev.'] / df['mean']
high_error_count = (rel_errors > 0.1).sum()
if high_error_count > 0:
print(f"Warning: {high_error_count} tally bins have >10% uncertainty")
# 4. Validate physics reasonableness
if sp.keff.nominal_value < 0.5 or sp.keff.nominal_value > 2.0:
print("Warning: k-effective seems unrealistic")
# 5. Check for divide-by-zero in calculations
mean_values = df['mean'].values
zero_tallies = (mean_values == 0).sum()
if zero_tallies > 0:
print(f"Warning: {zero_tallies} tally bins have zero scores")Tutorial snippet — no separate file in examples repo
Always report uncertainties, check convergence, and validate physical reasonableness before drawing conclusions.
Check yourself
- Convert per-source-particle tallies into absolute reaction rates at a given power?
- Reshape a mesh tally back into the geometry it was scored on?
- Check uncertainties before quoting any number?