Hi
I’m fairly new to openmc and I’m trying to create a sphere of material that is surrounded by two neutron sources which irradiates the sphere and outputs a list of nuclides that have been activated and show their new isotopes. I’ve attached my current code and a ss of the output. Any help would be great!
import openmc
import matplotlib.pyplot as plt
# Material definition
soil = openmc.Material(name='Iron Rich Soil')
soil.add_element('Fe', 0.2)
soil.add_element('Si', 0.25)
soil.add_element('Al', 0.15)
soil.add_element('O', 0.4)
# Geometry definition
radius = 0.5
sphere = openmc.ZCylinder(r=radius, boundary_type='vacuum')
sphere_region = -sphere
sphere_cell = openmc.Cell(region=sphere_region, fill=soil)
# Universe with sphere
universe = openmc.Universe(cells=[sphere_cell])
# Geometry
geometry = openmc.Geometry(universe)
geometry.export_to_xml()
# Plot the geometry (optional)
universe.plot(width=(3.0, 3.0), basis='xy')
plt.show()
# Define the neutron source
source = openmc.Source()
source.space = openmc.stats.Point((0, 0, 0)) # Source located at the center of the sphere
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Watt(a=1.7, b=1)
source.strength = 1e10 # Example: Increase source strength to 1e10 neutrons
# Settings
settings = openmc.Settings()
settings.batches = 20
settings.inactive = 2
settings.particles = 10000
settings.run_mode = 'fixed source'
settings.source = source
settings.export_to_xml()
# Define a mesh tally to record activation rates
mesh = openmc.RegularMesh()
mesh.dimension = [10, 10, 10]
mesh.lower_left = [-radius, -radius, -radius]
mesh.upper_right = [radius, radius, radius]
# Define the tally
tally = openmc.Tally(name='activation')
tally.filters = [openmc.MeshFilter(mesh)]
tally.scores = ['(n,gamma)'] # Ensure nuclide scores are correctly specified here
tallies = openmc.Tallies([tally])
tallies.export_to_xml()
# Run simulation
openmc.run()
# Load statepoint file
sp = openmc.StatePoint('statepoint.10.h5')
# Get activation tally
activation_tally = sp.get_tally(name='activation')
# Extract and print the results
df = activation_tally.get_pandas_dataframe()
# Print the entire DataFrame to inspect its structure
print(df)
# Find activated nuclides
activated_nuclides = df[df['mean'] > 0]
# Display activated nuclides
print("Activated nuclides:")
print(activated_nuclides[['nuclide', 'mean']])
# Plot activation results (optional)
plt.figure()
plt.hist(df['mean'], bins=50, alpha=0.75)
plt.xlabel('Activation Rate')
plt.ylabel('Count')
plt.title('Activation Rate Distribution')
plt.show()
