I am modeling a deuterium-tritium fuel pellet undergoing thermonuclear fusion, at a temperature of ___ and a density of ___. This produces neutrons with energies of ___ MeV uniformly throughout the fuel pellet. I am attempting to get the neutron energy spectrum emerging from this fuel pellet and then graph it. However, it seems like it does not like my geometry. When I reduce the number of particles, or decrease the density or increase the radius, it is fine, but at the proper conditions, it gives me the following output:
Reading settings XML file...
Reading cross sections XML file...
Reading materials XML file...
Reading geometry XML file...
Reading H2 from cross_sections_data/H2.h5
Reading H2 WMP data from cross_sections_data/wmp/001002.h5
Reading H3 from cross_sections_data/H3.h5
Reading H3 WMP data from cross_sections_data/wmp/001003.h5
Minimum neutron data temperature: 250.0 K
Maximum neutron data temperature: 2500.0 K
Reading tallies XML file...
Preparing distributed cell instances...
Writing summary.h5 file...
Maximum neutron transport energy: 20000000.0 eV for H3
===============> FIXED SOURCE TRANSPORT SIMULATION <===============
Simulating batch 1
Simulating batch 2
Simulating batch 3
Simulating batch 4
Simulating batch 5
Simulating batch 6
Simulating batch 7
Simulating batch 8
Simulating batch 9
Simulating batch 10
Simulating batch 11
Simulating batch 12
Simulating batch 13
WARNING: Particle 1232008 underwent maximum number of events.
Simulating batch 14
WARNING: Particle 692464 underwent maximum number of events.
Simulating batch 15
Simulating batch 16
Simulating batch 17
Simulating batch 18
WARNING: Could not find the cell containing particle 505877
Traceback (most recent call last):
File "/home/number1son100/Desktop/TF_Link_Verification/neutron_simulate.py", line 7, in <module>
neutron_methods.run(radius)
File "/home/number1son100/Desktop/TF_Link_Verification/neutron_methods.py", line 274, in run
openmc.run()
File "/home/number1son100/anaconda3/envs/openmc-env/lib/python3.9/site-packages/openmc/executor.py", line 227, in run
_run(args, output, cwd)
File "/home/number1son100/anaconda3/envs/openmc-env/lib/python3.9/site-packages/openmc/executor.py", line 38, in _run
raise RuntimeError(error_msg)
RuntimeError: OpenMC aborted unexpectedly.
Does anybody know why this might be? I believe my geometry is simple and should have no gaps in it at all. This was done when running the neutron_simulate.py script, included below along with the function file neutron_methods.py.
neutron_methods.py:
import openmc
import math
import numpy as np
import matplotlib.pyplot as plt
def expand_temperatures():
D = openmc.data.IncidentNeutron.from_njoy('D.endf', temperatures=[300., 600., 1000.])
"""
material_setup(materials)
Inputs: materials Empty materials object, to be filled by the function
Outputs: None
Purpose: Defines all necessary materials, so they can be referenced in other functions
Saves materials to an xml file, to be used by openmc when it is run
"""
def material_setup(materials):
# Set cross-sections path
materials.cross_sections = "cross_sections_data/cross_sections.xml"
D_T = openmc.Material()
D_T.name = 'D-T'
D_T.add_nuclide('H2', 0.5, 'ao')
D_T.add_nuclide('H3', 0.5, 'ao')
D_T.set_density('g/cm3', 300)
D_T.temperature = 10*1000/(8.617333262*10**(-5))
materials.append(D_T)
# Export Data to File
materials.export_to_xml()
"""
get_material(materials, name)
Inputs: materials Materials object holding all relevant materials
name String indicating the name of the material desired
Outputs: The material object desired
Purpose: Searches the materials object for the material with the name provided
"""
def get_material(materials, name):
for material in materials:
if material.name == name:
return material
break
"""
geometry_setup(radius, materials)
Inputs: radius The radius of the fuel pellet
materials Materials object holding all relevant materials
Outputs: cells A list of cells to be used in setting cell tallies
Purpose: Defines all necessary geometry and regions of space, and saves this data to an xml file
The two regions are the fuel pellet and a thin layer without it that is empty
"""
def geometry_setup(radius, materials):
# Setup arrays to hold objects of import
surfaces = []
cells = []
radii = []
# Define fusion region inside the pellet
surf = openmc.Sphere()
surf.r = radius*1.0
radii.append(surf.r)
surf.name = "Inner Region"
surfaces.append(surf)
shape = -surf
material = get_material(materials, 'D-T')
cell = openmc.Cell(fill=material, region=shape)
cell.name = surf.name
cells.append(cell)
# Set the thickness of the outer empty shell to be one percent of the pellet radius
thickness = radius*0.01
# Define the outer surface of the model
surf = openmc.Sphere()
surf.r = radii[0] + thickness
radii.append(surf.r)
surf.name = "Outer Surface"
surf.boundary_type = 'vacuum'
surfaces.append(surf)
# Generates the outer cell and adds it to the list
shape = +surfaces[0] & -surf
cell = openmc.Cell(region=shape)
cell.name = "Empty Outer Shell"
cells.append(cell)
# Sets the root universe to include all the cells that were just created
root_uni = openmc.Universe()
root_uni.add_cells(cells)
# Export this root universe as a geometry object to an xml file
geom = openmc.Geometry(root_uni)
geom.export_to_xml()
return cells
"""
settings_setup(radius)
Inputs: radius The radius of the fuel pellet
Outputs:
Purpose: Defines the neutron source in the fusion region to emit 14.1 MeV neutrons
isotropically throughout
"""
def settings_setup(radius):
# Ensure thet the simulation is fixed source
settings = openmc.Settings()
settings.run_mode = 'fixed source'
# Set the number of batches and particles used in the simulation
settings.particles = 1500000
settings.batches = 150
settings.generations_per_batch = 1
settings.inactive = 0
# (Hopefully) set up settings for doppler broadening
settings.resonance_scattering['enable'] = True
settings.resonance_scattering['mode'] = 'rvs'
settings.resonance_scattering['energy_min'] = 0
settings.resonance_scattering['energy_max'] = 15e6
# Define the neutron source region, which will emit 14.1 MeV neutrons
r = openmc.stats.Uniform(0.0, radius)
phi = openmc.stats.Uniform(0, 2*math.pi)
theta = openmc.stats.Uniform(0, math.pi)
source = openmc.Source()
source.space = openmc.stats.SphericalIndependent(r, phi, theta)
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14.1e6], [1.0])
settings.source = source
# Temperature Settings to allow for high temperature simulations with low temperature cross-sections
settings.temperature["method"] = "nearest"
settings.temperature["tolerance"] = 1000000000
settings.temperature["range"] = (0.0, 1000000000)
settings.temperature["multipole"] = True
# Set the output folder for results
settings.output = {
'path': 'results',
'tallies':True,
'summary':True
}
# Export the settings to an xml document
settings.export_to_xml()
"""
create_energy_spectrum(range, scale):
Inputs: range Range over which energy spectrum is generated
scale Scale by which the energies are multiplied
Outputs: energies Array of energies with this range and scale in increments of 0.05
Purpose: Creates an array of energies with this range and scale in increments of 0.05
"""
def create_energy_spectrum(range, scale):
# Create the spectrum of energy ranges that will be differentiated when tallying the flux
step = 0.05
last = 0.0
energies = [last]
while energies[len(energies)-1]<range*scale:
next = last+step
energies.append(next*scale)
last = next
return energies
"""
tallies_setup(cells)
Inputs: cells A list of the different cells used, omitting the inner one, so flux can be calculated accross them
Outputs:
Purpose: Sets a tally so that the flux across each radius of interest, and at different energies, will be recorded by openmc
"""
def tallies_setup(cells):
# Generate an energy spectrum between 0 and 20 MeV
energies = create_energy_spectrum(20,10**6)
# Make tallies to keep track of flux at different radii and at different energies, omitting the empty inner region
cell_filter = openmc.CellFilter(cells[-1])
energy_filter = openmc.EnergyFilter(energies)
tally = openmc.Tally()
tally.filters = [cell_filter, energy_filter]
# Specify that only flux should be counted
tally.scores = ['flux']
# Export this tally information to an xml file
tallies = openmc.Tallies([tally])
tallies.export_to_xml()
return energies
"""
setup(radius)
Inputs: radius The radius of the fuel pellet
Outputs:
Purpose: Creates the setup files needed for openmc to run the simulation, based on the configuration data provided
"""
def setup(radius):
# Initialize Materials File
materials = openmc.Materials()
# Generate necessary xml files
material_setup(materials)
cells = geometry_setup(radius, materials)
settings_setup(radius)
energies = tallies_setup(cells)
return energies
"""
process_results(energies)
Inputs: energies Array of energies to use in plotting
Outputs:
Purpose: Processes the flux data, and then prints and plots as requested
"""
def process_results(energies):
# Grabs the flux data from a file
global flux
flux = []
sp = openmc.StatePoint('results/statepoint.150.h5')
# Processes the data to appear in an array that is offset from the energy spectrum
tally = sp.get_tally(scores=['flux'])
tally.mean.shape = (1,len(energies)-1)
values = np.array(tally.mean)
flux = values[0]
# Offsets the energy spectrum to align with the data, which was gathered from energy ranges, not discrete energies
new_energies = []
for i in range(len(energies)):
if i != 0:
new_energies.append(((energies[i]-energies[i-1])/2.0 + energies[i-1])/10**6)
energies = new_energies
# Plot flux as a function of energy, in arbitrary units
plot = plt.semilogy(new_energies, flux, marker='')
plt.xlabel("Energies (MeV)")
plt.ylabel("Flux (Arbitrary Units)")
plt.title("Flux in a D-T fusion fuel pellet")
plt.show()
"""
run(r_i, length, r_frc, elongation, regions, degrees)
Inputs: radius The radius of the fuel pellet
Outputs:
Purpose: Creates the setup files needed and then runs openmc. Once results are generated, prints the desired results to console and then plots them, if asked to.
"""
def run(radius):
# Sets up the xml files for the simulation
energies = setup(radius)
# Runs the simulation
import os
#os.system('openmc')
openmc.run()
# Find the data and plot it
process_results(energies)
neutron_simulate.py:
import neutron_methods
import math
radius = 333.3*10**-6
# Run the simulation
neutron_methods.run(radius)