Maybe I am completely missing something in the way my geometry and materials are set up, but any time I try to simulate photons, there is lots of absorption and scattering through what should be empty space. This changes when I switch to neutrons. Below are two slices of my model, which at the time only has a 2x2x2 mm Argon box. The first is the absorption with photons, and the second is the absorption with neutrons. Can anyone help in explaining what is causing the difference between the two?
That looks interesting, are you able to let us know what density of Argon you have in your material card or perhaps better still are you able to share the model.
By the way I hear there is a package for plotting regular mesh tallies that includes geometry outlines if that helps 
Here is my entire code:
import cadquery as cq
from cadquery import Assembly
from cad_to_dagmc import CadToDagmc
%matplotlib inline
import os
import openmc
import matplotlib.pyplot as plt
import openmc.lib
import numpy as np
from mpl_toolkits.mplot3d import axes3d
from stl_to_h5m import stl_to_h5m
assert(openmc.lib._dagmc_enabled())
Aluminum = openmc.Material(name='aluminum')
Aluminum.add_element('Al', 1, percent_type='ao')
#Aluminum.set_density('g/cm3', 2.7)
FR4 = openmc.Material(name='FR4')
FR4.add_element('Si', .46, percent_type='ao')
FR4.add_element('O', .29, percent_type='ao')
FR4.add_element('C', .13, percent_type='ao')
FR4.add_element('H', .03, percent_type='ao')
FR4.add_element('Br',.09, percent_type='ao')
Silicon = openmc.Material(name='Silicon')
Silicon.add_element('Si',1,percent_type = 'ao')
Derlin = openmc.Material(name='Derlin')
Derlin.add_element('C', .69, percent_type='ao')
Derlin.add_element('H', .12, percent_type='ao')
Derlin.add_element('O', .19, percent_type='ao')
Tungsten = openmc.Material(name = 'Tungsten')
Tungsten.add_element('W', 1, percent_type='ao')
Argon = openmc.Material(name = 'Argon')
Argon.add_element('Ar',1)
mats = openmc.Materials()
mats.append(Aluminum)
mats.append(FR4)
mats.append(Derlin)
mats.append(Tungsten)
mats.append(Argon)
mats.append(Silicon)
mats.export_to_xml()
dag_univ = openmc.DAGMCUniverse(filename = 'ArgonBox.h5m')
# creates an edge of universe boundary at a large radius
vac_surf = openmc.Sphere(r=300, y0 = 100,z0=0, surface_id=9999, boundary_type="vacuum")
# specifies the region as below the universe boundary
region = -vac_surf
# creates a cell from the region and fills the cell with the dagmc geometry
containing_cell = openmc.Cell(cell_id=9999, region=region, fill=dag_univ)
model = openmc.Model()
model.geometry = openmc.Geometry(root=[containing_cell])
model.geometry.export_to_xml()
settings = openmc.Settings()
settings.batches = 50
settings.inactive = 0
settings.particles = 50000
settings.electron_treatment = 'ttb'
settings.run_mode = 'fixed source'
def sample_spherical(npoints, ndim=3):
vec = np.random.randn(ndim, npoints)
vec /= np.linalg.norm(vec, axis=0)
return vec
phi = np.linspace(0, np.pi, 20)
theta = np.linspace(0, 2*np.pi, 40)
radius = 275
x = np.outer(np.sin(theta), np.cos(phi))*radius
y = np.outer(np.sin(theta), np.sin(phi))*radius
z = np.outer(np.cos(theta), np.ones_like(phi))*radius
xi, yi, zi = sample_spherical(10000)*radius
src_locations = []
for i in range(len(xi)):
src_locations.append([xi[i], yi[i]+100,zi[i]])
fig, ax = plt.subplots(1, 1, subplot_kw={'projection':'3d', 'aspect':'auto'})
ax.plot_wireframe(x, y, z, color='k', rstride=1, cstride=1)
ax.scatter(xi, yi, zi, s=100, c='r', zorder=10)
src_e = openmc.stats.Uniform(0,3e5)
angle1 = openmc.stats.Isotropic()
# create source for each location
sources = []
for loc in src_locations:
src_pnt = openmc.stats.Point(xyz=loc)
src = openmc.Source(space=src_pnt,angle = angle1, energy=src_e, particle = 'photon')
sources.append(src)
settings.rel_max_lost_particles = .9999
settings.source = sources
settings.export_to_xml()
mesh = openmc.RegularMesh().from_domain(
model.geometry, # the corners of the mesh are being set automatically to surround the geometry
dimension=[1,200, 200] # (change to cube for vtk file or (100,1,100) for plot)
)
sensor = openmc.MaterialFilter(Silicon)
mat_filter = openmc.MaterialFilter(Aluminum)
tallies = openmc.Tallies()
# Create mesh filter for tally
mesh_filter = openmc.MeshFilter(mesh)
particle_filter = openmc.ParticleFilter('photon')
# Create mesh tally to score tritium production
mesh_tally_1 = openmc.Tally(name='absorption')
mesh_tally_1.filters = [mesh_filter,particle_filter]
mesh_tally_1.scores = ['absorption']
tallies.append(mesh_tally_1)
mesh_tally_2 = openmc.Tally(name='scattering')
mesh_tally_2.filters = [mesh_filter]
mesh_tally_2.scores = ['scatter']
tallies.append(mesh_tally_2)
mesh_tally_3 = openmc.Tally(name='photoelectric')
mesh_tally_3.filters = [mesh_filter,mat_filter]
mesh_tally_3.scores = ['photoelectric']
tallies.append(mesh_tally_3)
energy_bins = np.linspace(0,3e5,10000)
energyfilter = openmc.EnergyFilter(energy_bins)
count_tally = openmc.Tally(name='counts')
count_tally.filters = [energyfilter,mat_filter]
count_tally.scores = ['photoelectric']
tallies.append(count_tally)
tallies.export_to_xml()
model = openmc.model.Model(model.geometry, mats, settings, tallies)
# deletes old files
!rm summary.h5
!rm statepoint.*.h5
# runs the simulation
sp_filename = model.run()
# open the results file
results = openmc.StatePoint(sp_filename)
my_tally = results.get_tally(scores=['photoelectric'])
my_slice = my_tally.get_slice(scores=['photoelectric'])
my_slice.mean.shape = (mesh.dimension[1], mesh.dimension[2]) # setting the resolution to the mesh dimensions
fig = plt.subplot()
plt.show(fig.imshow(my_slice.mean, extent=[-250,250,-250,250]))
sp = openmc.StatePoint('statepoint.50.h5')
tally = sp.get_tally(name='counts')
pulse_height_values = tally.get_values(scores=['photoelectric']).flatten()
energy_bin_centers = energy_bins[1:] + 0.5 * (energy_bins[1] - energy_bins[0])
plt.figure()
plt.loglog(energy_bin_centers/1000, pulse_height_values)
plt.xlabel('Energy [keV]')
plt.ylabel('Counts')
plt.title('Pulse Height Values')
plt.grid(True)
plt.tight_layout()
The Argon Box file is literally just a 2x2x2 cube dagmc geometry made in cubit, with one material group made for Argon.
Thanks for the nice script. I see this is a DAGMC geometry. Also tha ks for let us know this was made with cubit.
It looks like Argon does not have a density set. I guess openmc defaults to sum which might not be what you want. Neutrons might travel further if you set a low density
I can see the photon absorption tally in the input script but i can’t see the neutron absorption tally.
Is it worth plotting the geometry. Sometimes we see units in cad not being converted to cm needed by openmc. This would just check the small squere is inside the sphere geometry. Or if not plotting then using nbconvert to make some stl geometry to view in paraview
I’m wondering if it is worth using the bounded_universe method to make the bounding sphere. This would take care of the potential sizing issue.
I should have specified, but in order to get the differences in the outcomes, I simply change the source and particle filter to neutron, and proceed with the simulation in the exact same way. I have tried to use the bounded_universe method, but cannot seem to run the code with the external source without getting the error that too many particles are lost. It would make sense that the source is beyond the boundary, so all particles are killed by it, but I do not know how to change the size of the boundary without specifying it in the method I have now.