Photon Transport Simulation with Two Layers of Material & Unstructured Mesh #27
Conversation
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
TwoLayers_Geometry.py
Outdated
|
|
||
| # 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.
docstrings here would help be clear about what everything is doing
TwoLayers_Geometry.py
Outdated
| @author: Anupama Rajendra | ||
| """ | ||
|
|
||
| import openmc |
There was a problem hiding this comment.
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
There was a problem hiding this comment.
I might be missing something here, but OpenMC_PhotonTransport attempts to access openmc inside TwoLayers_Geometry even when only make_spherical_shells is imported
Source_Mesh_Reader.py
Outdated
| 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.
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
Photon_TallytoVtk.py
Outdated
| 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.
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.
Co-authored-by: Paul Wilson <[email protected]>
- Created dictionary with material densities - Created OpenMC Material and Materials objects with separate function
TwoLayers_Geometry.py
Outdated
|
|
||
| # 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.
This makes multiple shells, so maybe change the name.
Could also generalize this to take an inner radius plus layers:
| def make_spherical_shell(inner_radius_W, outer_radius_W, inner_radius_C, M_1, M_2): | |
| def make_spherical_shells(inner_radius_W, layers): | |
| ''' | |
| Builds a set of nested spherical shells, with an inner_radius and | |
| multiple layers each with a material definition and a thickness | |
| inputs | |
| ===== | |
| inner_radius : the radius of the inner surface of the first shell | |
| layers : list of tuples where each tuple is (material, thickness) | |
| ''' |
There was a problem hiding this comment.
I think materials will need to be an input as this function is also assigning materials to the geometry. I'm wondering if this can be done with one set of radii/thickness inputs and the materials from TwoLayers_Materials.
| def make_spherical_shell(inner_radius_W, outer_radius_W, inner_radius_C, M_1, M_2): | |
| def make_spherical_shell(radii, materials): | |
| ''' | |
| Takes a list of materials and radii (outermost to innermost) | |
| and creates an OpenMC Geometry object. | |
| element_list = list of elements (str) in order from radially outermost to innermost | |
| radii = list of the radii (outermost to innermost) that serve as the boundaries of the | |
| spherical shells; must be in the same order as the specified materials | |
| ''' | |
| cells = [] | |
| for index, material in enumerate(materials): | |
| outer_sphere = openmc.Sphere(r = radii[index]) # add boundary_type here | |
| inner_sphere = openmc.Sphere(r = radii[index + 1]) | |
| element_region = +inner_sphere & -outer_sphere | |
| element_shell = openmc.Cell(fill = materials[index], region = element_region) | |
| cells.append(element_shell) | |
| if material == materials[-1]: | |
| void = openmc.Cell(fill = None, region = -inner_sphere) | |
| cells.append(void) | |
| geometry = openmc.Geometry(cells) | |
| return geometry | |
There was a problem hiding this comment.
I agree that the materials are necessary. My version had them bound in the layers. This kind of data structure can lead to clarity because the data gets bound together in instructive ways. An example of calling it might be:
inner_radius = 10.0
# material_a and material_b are OpenMC Materials
layers = [(material_a, 5.0), (material_b, 4.0)]
geometry = make_spherical_shells(inner_radius, layers)
The function may be:
def make_spherical_shells(inner_radius, layers):
inner_sphere = openmc.Sphere(inner_radius)
cells = [openmc.Cell(region=-inner_sphere)]
for (material, thickness) in layers:
outer_radius = inner_radius + thickness
outer_sphere = openmc.Sphere(outer_radius)
cells.append(openmc.Cell(fill = material, region = inner_sphere & -outer_sphere))
inner_radius = outer_radius
inner_sphere = outer_sphere
cells.append(openmc.Cell(region = outer_sphere)
geometry = openmc.Geometry(cells)
return geometry
(There may be minor syntax errors, but this is the basic idea)
Make dictionary of densities with ALARA element library Co-authored-by: Paul Wilson <[email protected]>
Redefine inner and outer W spheres Co-authored-by: Paul Wilson <[email protected]>
- Redefined make_spherical_shells to take layers and inner_radius as input - Deleted reference to TwoLayers_Materials to maintain generality
There was a problem hiding this comment.
Some requests for this file:
- add argparse for filename (with defaults)
- create separate functions to:
- read data from statepoint
- plot data (optoinal via argparse)
- write data to file (filename set by argparse)
- all of these called from a main function
WC_Layers/OpenMC_PhotonTransport.py
Outdated
| parser = argparse.ArgumentParser(description="Specify required inputs: file path to ALARA Element Library, element name, inner radius [cm], outer_radius [cm]") | ||
|
|
||
| #Required user inputs: | ||
| parser.add_argument('--filepath', type=str, required=True) | ||
| parser.add_argument('--element_1', type=str, required=True) | ||
| parser.add_argument('--element_2', type=str, required=True) | ||
| #Shell radii are left as user inputs, but the mesh is specific to W_inner_radius = 1000, W_outer_radius = 1005, C_inner_radius = 995 | ||
| parser.add_argument('--W_inner_radius', type=float, required=True) | ||
| parser.add_argument('--W_outer_radius', type=float, required=True) | ||
| parser.add_argument('--C_inner_radius', type=float, required=True) | ||
|
|
||
| args = parser.parse_args() |
WC_Layers/OpenMC_PhotonTransport.py
Outdated
| from TwoLayers_Materials import * | ||
| from TwoLayers_Geometry import * | ||
|
|
||
| Mesh_File = 'OpenMC_Mesh.h5m' |
There was a problem hiding this comment.
set this with argparse, possibly with default
WC_Layers/OpenMC_PhotonTransport.py
Outdated
| fp = args.filepath | ||
| E_1 = args.element_1 | ||
| E_2 = args.element_2 | ||
| R_W_1 = args.W_inner_radius | ||
| R_W_2 = args.W_outer_radius | ||
| R_C_1 = args.C_inner_radius |
WC_Layers/OpenMC_PhotonTransport.py
Outdated
| return Source_List, Particle_Filter, Total_Mesh, Mesh_Dist | ||
|
|
||
| # Define tallies | ||
| def tallies(W_Shell, C_Shell, Particle_Filter, Total_Mesh): |
There was a problem hiding this comment.
given the flexibility of the method to make cells, this may make more sense if it took a list of cells
|
|
||
| def save_source_density(sd_list, sd_filename): | ||
| ''' | ||
| Saves source density data as a text file. |
There was a problem hiding this comment.
It looks like you open multiple files. Can you add more to the docstring about the difference between files and the structure of the files?
- Modified loop over bounds in make_source (no longer uses enumerate)
Co-authored-by: Paul Wilson <[email protected]>
…sport_Inputs.yaml
| outputs: | ||
| talls: OpenMC Tallies object | ||
| ''' | ||
| particle_filter = openmc.ParticleFilter('photon') |
There was a problem hiding this comment.
You define this particle filter, which makes sense for a photon transport problem, but all the rest of your tally definition seems like it is meant for neutrons, especially the choice of scores, energy filter
There was a problem hiding this comment.
Yes, tally 1 originated from the neutron transport model. Perhaps a heating/absorption photon tally (with the appropriate variable names) would provide more relevant information?
There was a problem hiding this comment.
For R2S problems, it's customary to do an equivalent dose tally.
| import openmc | ||
| import matplotlib.pyplot as plt | ||
|
|
||
| def read_statepoint(sp_filename, flux_spectrum_tally_id, mesh_number, summed_axes_energy_filter, energy_filter_index): |
There was a problem hiding this comment.
With the addiction of a docstring here, I think some better variable names may emerge:
flux_spectrum=>flux_spectrum_tallytally_data_reshaped=>flux_spectrum_meanflux_sum_energy=>total_flux
| mesh=sp.meshes[mesh_number] | ||
| # Return tally data condensed into 1 dimension per filter | ||
| tally_data_reshaped = flux_spectrum.get_reshaped_data(value='mean') | ||
| flux_sum_energy = tally_data_reshaped.sum(axis=summed_axes_energy_filter) |
There was a problem hiding this comment.
It's not clear to me what axes this data is being summed over. The data is originally a flux spectrum of G groups for each of M mesh elements. I assumed that this sum was summing over G, to produce a total flux for each of the M mesh elements. However, when you use this in the plot below, it seems like it may be summing over M to produce a single flux spectrum of G groups. If so, then we need to make sure that the volume normalization is being done correctly. (Full disclosure: I'm not sure when/how volume normalization of tallies happens in OpenMC. If this was MCNP-like, then this operation would be wrong because the results would already be volume normalized per mesh element, and summing would require a volume-weighting.)
There was a problem hiding this comment.
For the purpose of plotting the data, the summation is happening over M so that the flux is divided into G groups.
From digging into r2s.py for a bit, I believe that the outputted photon source mesh info is normalized by volume. That said, I don't think there is any division by volume automatically happening on the OpenMC side, so ALARA/r2s may be outputting the wrong units if it is expecting values from MCNP. However, I do have the volume normalization set to False within MeshSpatial of the PhotonTransport code, so there is no remultiplication by volume of the source strengths. I think this is still worth looking into further.
…Mesh_Reader_Inputs.yaml
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.