I’m having trouble getting a uniform, random distribution of spheres packed into my geometry with openmc.model.pack_spheres().
I use the following snippet of code to generate the positions of the sphere centers:
outer_radius = 422.5*1e-4
compact_rad = openmc.ZCylinder(r=0.635, boundary_type='transmission')
min_z = openmc.ZPlane(z0=0, boundary_type='transmission')
max_z = openmc.ZPlane(z0=42, boundary_type='transmission')
region = -compact_rad & +min_z & -max_z
centers = openmc.model.pack_spheres(radius=outer_radius, region=region, pf=0.42)
This code does generate 70733 spheres and satisfies the packing fraction of 42%, however, most of the spheres are towards the bottom of the geometry, as seen in the following plot:
I quickly binned them into 6-cm regions along the z-axis, and found the following:
00 < z < 06 : 15628
06 < z < 12 : 10015
12 < z < 18 : 9471
18 < z < 24 : 9298
24 < z < 30 : 9017
30 < z < 36 : 8747
36 < z < 42 : 8557
It seems the spheres are nearly twice as dense at the bottom of the channel as at the top. In fact, the packing fraction between 0 cm and 6 cm seems to be ~64.9%, which actually exceeds the theoretical maximum described in the tutorial.
Can someone please explain this behavior and help me make a triso distribution that is more uniform along the z-axis?
Thanks
Lu
(pronouns: they/them)
Hi Lu,
I tried running your code snippet and binning along the z-axis, and I see a much more uniform distribution:
bins = range(0, 48, 6)
plt.hist(centers[:,2], bins=bins, edgecolor='k')
plt.xticks(bins)
Would you mind sharing how you are generating your plot/histogram?
Thanks,
Amanda
Hi Amanda,
Thanks for writing back to me, and sorry for my late response. I have tried a few things to untangle this issue but I’m still quite confused.
- I reran my input that generated these triso centers (the exact same code I ran above) and I did get a better, distribution, but it is still flawed. Attached below is the plot generated by openmc. I would share the .py input to make that plot, but the forum tells me that new users cannot share files. I have tried to select the most important block of code, but please let me know if you need to see anything more.
Processing: plot.py…
centers = np.empty(shape=3)
with open("newest-centers.txt") as file:
for line in file:
coords = np.array(line.split(" "))
coords[-1] = coords[-1][:-2]
coords2 = np.empty(shape=0)
for coord in coords:
if coord != '':
coords2 = np.append(coords2, float(coord))
centers = np.append(centers, coords2)
centers = centers[3:]
centers = np.reshape(centers, (int(len(centers)/3),3))
trisos = [openmc.model.TRISO(outer_radius, triso_univ, c) for c in centers]
fuel.volume = 4/3 * np.pi * (radii[0]/100) ** 3 * len(trisos)
trisvol = 4/3 * np.pi * (radii[-1]/100) ** 3 * len(trisos)
compvol = np.pi * compact_rad.r ** 2 * (max_z.z0 - min_z.z0)
centers = np.vstack([triso.center for triso in trisos])
comp = openmc.Cell(region=comp_region)
lower_left, upper_right = comp.region.bounding_box
shape = (3, 3, 3)
pitch = (upper_right - lower_left)/shape
lattice = openmc.model.create_triso_lattice(
trisos, lower_left, pitch, shape, graphite)
comp.fill = lattice
-
I also made histograms out of this distribution, but I am only allowed to post one image at a time. The histograms describe the same distribution with two peaks near the origin as the plot above.
-
Additionally, I noticed that many of the centers in the old distribution (which I accidentally created with 0<z<42) and the new distribution (with -21<z<21) have the same x/y coordinates, and the z coordinate shifted down by 21 cm. I wrote a python script to go line by line and compare, and here’s what I found:
- ntrisos in the new file is 70733
- ntrisos in the old file is 70733
- num trisos exactly 21 cm displaced is 64078
- num trisos any other distance displaced is 6655
- of the trisos exactly 21 cm displaced, 61805 have the same x,y coords, and 2273 have different xy coords from their predecessors
- of the trisos not exactly 21 cm displaced, 6375 have the same x,y coords, and 280 have different xy coords from their predecessors
That’s all I have for now. Please let me know what I can do to get a flatter z-distribution, and if you’d like me to share anything else.
Thanks again
Lu
pronouns: they / them
Hi Lu,
Thanks for the extra information. As a quick sanity check, could you run the following code snippet:
outer_radius = 422.5*1e-4
compact_rad = openmc.ZCylinder(r=0.635, boundary_type='transmission')
min_z = openmc.ZPlane(z0=-21, boundary_type='transmission')
max_z = openmc.ZPlane(z0=21, boundary_type='transmission')
region = -compact_rad & +min_z & -max_z
centers = openmc.model.pack_spheres(radius=outer_radius, region=region, pf=0.42)
hist, bin_edges = np.histogram(centers[:,2], bins=range(-21, 27, 6))
print('histogram:', hist)
print('bin edges:', bin_edges)
and send me the output (histogram values and bin edges)?
Thanks,
Amanda
I seem to have failed the sanity check, as the output here is:
histogram: [10058 10161 10101 10134 10097 10088 10094]
bin edges: [-21 -15 -9 -3 3 9 15 21]
Does this mean I am making a mistake when I write the triso centers into a text file and then read them back into openmc? I don’t want to run the pack_spheres() routine every time I edit the model in some small way, as it takes quite a while to generate all those spheres.
Please let me know if there is a stable way to do this. So far I have been running this code:
from math import pi
import numpy as np
import openmc
import openmc.model
outer_radius = 422.5*1e-4
compact_rad = openmc.ZCylinder(r=0.635, boundary_type='transmission')
min_z = openmc.ZPlane(z0=-21, boundary_type='transmission')
max_z = openmc.ZPlane(z0=21, boundary_type='transmission')
region = -compact_rad & +min_z & -max_z
centers = openmc.model.pack_spheres(radius=outer_radius, region=region, pf=0.42)
for _ in centers:
print(_)
#######################
trisovol = len(centers)*4/3*pi*outer_radius**3
regiovol = compact_rad.r ** 2 * pi * (max_z.z0 - min_z.z0)
print(f'packing fraction is {trisovol/regiovol}')
centers = np.vstack(centers)
print(centers.min(axis=0))
print(centers.max(axis=0))
as python dispersed_triso.py > triso_centers.txt
and then in my build_xml.py file I have:
# import centers from the inp file this block takes ~ 1 min to run
centers = np.empty(shape=3)
with open("triso_centers.txt") as file:
for line in file:
coords = np.array(line.split(" "))
coords[-1] = coords[-1][:-2]
coords2 = np.empty(shape=0)
for coord in coords:
if coord != '':
coords2 = np.append(coords2, float(coord))
centers = np.append(centers, coords2)
centers = centers[3:]
centers = np.reshape(centers, (int(len(centers)/3),3))
trisos = [openmc.model.TRISO(outer_radius, triso_univ, c) for c in centers]
fuel.volume = 4/3 * np.pi * (radii[0]/100) ** 3 * len(trisos)
trisvol = 4/3 * np.pi * (radii[-1]/100) ** 3 * len(trisos)
compvol = np.pi * compact_rad.r ** 2 * (max_z.z0 - min_z.z0)
centers = np.vstack([triso.center for triso in trisos])
print(centers.min(axis=0))
print(centers.max(axis=0))
print(f'there are {len(trisos)} trisos in each compact')
print(f'UCO vol is {fuel.volume}')
print(f'triso vol is {trisvol}')
print(f'fuel compact vol is {compvol}')
print(f'packing fraction is {trisvol/compvol}')
I have to say I’m not too surprised to learn I’ve been doing this wrong, but I don’t know another method. Do you see an error here? I’d like to avoid the lengthy wait associated with openmc.model.pack_spheres(). Please let me know the most effective way to accomplish this.
Thanks
Lu
Hi Lu,
One simple way to write the centers to a text file and later read them back is to use numpy.savetxt() and numpy.loadtxt(). For example, in dispersed_triso.py you can save the centers like this:
centers = openmc.model.pack_spheres(radius=outer_radius, region=region, pf=0.42)
np.savetxt('triso_centers.txt', centers)
$ python dispersed_triso.py
Then the first block in build_xml.py can be replaced with:
centers = np.loadtxt('triso_centers.txt')
Let me know if your results look any better after trying this.
Thanks,
Amanda
Amanda,
Thank you so much for the help! My friends and I have been looking at this for ages; I can’t believe none of us thought of this.
Here is the plot I get now. I computed a histogram as well and confirmed that these are well-distributed along the z-axis.
Thanks again,
Lu
Awesome! Glad I could help and that it was a simple fix.
