-
Notifications
You must be signed in to change notification settings - Fork 2
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Photon Transport Simulation with Two Layers of Material & Unstructured Mesh #27
base: main
Are you sure you want to change the base?
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Glad to see this is all (mostly) working
Co-authored-by: Paul Wilson <[email protected]>
- Reading total strengths directly from Source_Mesh_Reader instead of writing a new list in this script
- Brought the MeshBase object inside the loop over the energy bounds
Assign material library path to new function and delete unused lines
These two files can be imported into a neutron or photon transport problem
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A number of suggestions here. We should meet to discuss more.
I think all of these files should be in the WC_Layers
directory, and the neutron transport file should be updated to use the same geometry and materials scripts.
|
||
# Create geometry | ||
#Spherical shell: | ||
def make_spherical_shell(inner_radius_W, outer_radius_W, inner_radius_C, M_1, M_2): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
docstrings here would help be clear about what everything is doing
Cells = [Void, W_Shell, C_Shell] | ||
geometry = openmc.Geometry(Cells) | ||
geometry.export_to_xml() | ||
return geometry, Void, W_Shell, C_Shell, Cells |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
seems like you are returning more than necessary since Cells
contains W_Shell
and C_Shell
, and geometry
contains Cells
C_Shell = openmc.Cell(fill=M_2, region=S_C_3) | ||
Cells = [Void, W_Shell, C_Shell] | ||
geometry = openmc.Geometry(Cells) | ||
geometry.export_to_xml() |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Probably don't want to export here - separation of concerns points to putting computing
and input/output
in different methods.
@author: Anupama Rajendra | ||
""" | ||
|
||
import openmc |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not clear that we need to import openmc
here, if this method is only ever used when importing into another script that also imports openmc
S_W_1= openmc.Sphere(r=inner_radius_W) #sphere of radius 1000cm | ||
inside_W_sphere_1 = -S_W_1 | ||
outside_W_sphere_1 = +S_W_1 | ||
S_W_2 = openmc.Sphere(r=outer_radius_W, boundary_type='vacuum') | ||
inside_W_sphere_2 = -S_W_2 | ||
outside_W_sphere_2 = +S_W_2 | ||
S_W_3 = outside_W_sphere_1 & inside_W_sphere_2 #filled with specified material |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
S_W_1
and S_W_2
are surfaces, but S_W_3
is a region - better variable names will make this more readable.
It's also not clear that you need to define the intermediate variables inside_...
and outside_...
. If S_W_1
and S_W_2
(and S_C_1
below, etc) were named better, then the OpenMC notation with -
and +
might just be clear
S_W_1= openmc.Sphere(r=inner_radius_W) #sphere of radius 1000cm | |
inside_W_sphere_1 = -S_W_1 | |
outside_W_sphere_1 = +S_W_1 | |
S_W_2 = openmc.Sphere(r=outer_radius_W, boundary_type='vacuum') | |
inside_W_sphere_2 = -S_W_2 | |
outside_W_sphere_2 = +S_W_2 | |
S_W_3 = outside_W_sphere_1 & inside_W_sphere_2 #filled with specified material | |
inner_W_sphere= openmc.Sphere(r=inner_radius_W) #sphere of radius 1000cm | |
outer_W_sphere = openmc.Sphere(r=outer_radius_W, boundary_type='vacuum') | |
W_region = +inner_W_sphere & -outer_W_sphere |
tstt = file['tstt'] | ||
elements = tstt['elements'] | ||
tet_four = elements['Tet4'] | ||
tags = tet_four['tags'] | ||
sd = tags['source_density'][:] | ||
|
||
# Write source density to a separate file for each mesh | ||
with open(f'source_density_{file_index + 1}.txt', 'w') as source: | ||
for tet_element in sd: | ||
source.write(' '.join(map(str, tet_element)) + '\n') | ||
|
||
# Calculate summed (and individual mesh) strengths for each photon group | ||
summed_strengths = [] | ||
strengths_list = [] | ||
for group in range(photon_groups): | ||
total_strengths = np.sum(sd[:, group]) # Sum over all mesh elements | ||
summed_strengths.append(total_strengths) | ||
strengths = sd[:,group] | ||
strengths_list.append(strengths) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I recommend putting this in a method that does everything you need for a single file. That kind of modularity enhances readability. You might also want to separate the functions of reading from a file, processing the data, and writing to files, so methods like extract_source_data
, arrange_data
, write_source_density
summed_strengths = [] | ||
strengths_list = [] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
make these numpy arrays to start with?
total_strengths = np.sum(sd[:, group]) # Sum over all mesh elements | ||
summed_strengths.append(total_strengths) | ||
strengths = sd[:,group] | ||
strengths_list.append(strengths) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
How is strengths_list
different from sd
? sd
is a 2-D numpy array of size (num_elements
, num_groups
). It looks like strengths_list
is a list of lists with num_group
lists, each of length num_elements
.
Also - you don't appear to use any of the variables made in this loop?
flux_spectrum = sp.get_tally(id=2) | ||
mesh=sp.meshes[1] | ||
# Get the reshaped tally data | ||
tally_data_reshaped = flux_spectrum.get_reshaped_data(value='mean') | ||
|
||
flux_sum_en = tally_data_reshaped.sum(axis=(0,1,3,4,5)) | ||
#Vitamin-J energy filter: | ||
e_filter = flux_spectrum.filters[2] | ||
#Lower bounds of the energy bins | ||
e_filter_lower = e_filter.bins[:, 0] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A lot of magic numbers here including tally ID, mesh number, axes for sum, etc... more readable if these are named as variables, and might make sense if they lived in a function that took some of these as arguments.
# Plot flux spectrum | ||
fix, ax = plt.subplots() | ||
ax.loglog(e_filter_lower, flux_sum_en, drawstyle='steps-post') | ||
ax.set_xlabel('Energy [eV]') | ||
ax.set_ylabel('Flux [photon-cm/source]') | ||
ax.grid(True, which='both') | ||
plt.savefig('Photon_flux_vs_energy.png') | ||
plt.show() | ||
|
||
mesh_data = tally_data_reshaped.sum(axis=(0,2,3,4,5)) | ||
vtk = mesh.write_data_to_vtk(filename="Photon_Flux.vtk", datasets={"mean":mesh_data.flatten()}) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
separation of concerns: make different functions/methods for reading data, processing data, writing data, and plotting data.
OpenMC_W = make_W(element_1, OpenMC_SF) | ||
OpenMC_C = make_C(element_2, OpenMC_SF) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OpenMC_W = make_W(element_1, OpenMC_SF) | |
OpenMC_C = make_C(element_2, OpenMC_SF) | |
materials = [] | |
for material_id, element in enumerate(elements): | |
materials.append(make_element(element, material_id+1, OpenMC_SF)) |
This model currently uses the relative source intensities between different photon groups as the source strengths (and samples equally between individual mesh elements for any given energy group). However, it should be possible to add spatial variation of source intensity onto the MeshSpatial distribution.