Photon Transport Simulation with Two Layers of Material & Unstructured Mesh

WC_Layers/OpenMC_PhotonTransport.py

@@ -0,0 +1,142 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Aug  7 10:07:33 2024
+
+@author: Anupama Rajendra
+"""
+
+import openmc
+import argparse
+import numpy as np
+from Source_Mesh_Reader import extract_source_data, Files
+from TwoLayers_Materials import *
+from TwoLayers_Geometry import *
+
+# These will be moved to a yaml file
+openmc_mesh_file = 'OpenMC_Mesh.h5m'
+mesh_index = 0
+esd = extract_source_data(Files)
+
+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()
+
+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
+
+# fp = '../../ALARA/data/elelib.std'
+# E_1 = 'W'
+# E_2 = 'C'
+# R_W_1 = 1000
+# R_W_2 = 1005
+# R_C_1 = 995
+bounds = np.array([0, 
+                 1.00e+4, 
+                 2.00e+4, 
+                 5.00e+4, 
+                 1.00e+5, 
+                 2.00e+5, 
+                 3.00e+5, 
+                 4.00e+5, 
+                 6.00e+5, 
+                 8.00e+5, 
+                 1.00e+6, 
+                 1.22e+6, 
+                 1.44e+6, 
+                 1.66e+6, 
+                 2.00e+6, 
+                 2.50e+6, 
+                 3.00e+6, 
+                 4.00e+6, 
+                 5.00e+6, 
+                 6.50e+6, 
+                 8.00e+6, 
+                 1.00e+7, 
+                 1.20e+7, 
+                 1.40e+7, 
+                 2.00e+7])
+
+def make_source(cells, mesh_file):
+      '''
+    Creates a list of OpenMC sources, complete with the relevant space and energy distributions
+
+    inputs:
+        cells: list of OpenMC Cell objects
+        mesh_file: .h5/.h5m mesh onto which photon source will be distributed
+
+    output:
+        source_list: list of OpenMC independent sources
+        total_mesh: OpenMC Unstructured Mesh object
+    '''
+    source_list = []
+    total_mesh = openmc.UnstructuredMesh(mesh_file, library='moab')
+    for index in range(len(bounds) - 1):
+        mesh_dist = openmc.stats.MeshSpatial(total_mesh, strengths=esd[mesh_index][:,index], volume_normalized=False)
+        energy_dist = openmc.stats.Uniform(a=bounds[index], b=bounds[index + 1])
+        source_list.append(openmc.IndependentSource(space=mesh_dist, energy=energy_dist, strength=np.sum(esd[mesh_index][:, index]), particle='photon', domains=cells))
+    return source_list, total_mesh
+
+def tallies(total_mesh, tallied_cells):
+      '''
+    Creates tallies and assigns energy, spatial, and particle filters.
+
+    inputs: 
+        total_mesh: OpenMC unstructured mesh object
+        tallied_cells: OpenMC Cell/iterable of OpenMC Cell objects/iterable of Cell ID #
+
+    outputs:
+        talls: OpenMC Tallies object
+    '''
+    particle_filter = openmc.ParticleFilter('photon')
+    total_filter = openmc.MeshFilter(total_mesh)
+
+    neutron_tally = openmc.Tally(tally_id=1, name="Neutron tally")
+    neutron_tally.scores = ['flux', 'elastic', 'absorption']
+
+    # Implementing filter for neutron tally through shells with material
+    cell_filter = openmc.CellFilter(tallied_cells)
+    neutron_tally.filters = [cell_filter, total_filter]
+
+    # Vitamin-J energy filter created to assign to the flux tally.
+    energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-42")
+
+    spectrum_tally = openmc.Tally(tally_id=2, name="Flux spectrum")
+    spectrum_tally.filters = [cell_filter, total_filter, energy_filter_flux, particle_filter]
+    spectrum_tally.scores = ['flux']
+
+    talls = openmc.Tallies([neutron_tally, spectrum_tally])
+    return talls
+
+def settings(source_list):
+      '''
+    Creates an OpenMC Settings object
+
+    inputs:
+        source_list: iterable of OpenMC SourceBase objects
+    outputs:
+        sets: OpenMC Settings object
+    '''
+    sets = openmc.Settings()
+    sets.batches = 10
+    sets.inactive = 1
+    sets.particles = 100000
+    sets.source = source_list
+    sets.run_mode = 'fixed source'
+    return sets
+
+def create_openmc_model(geom_object, mats_object, sets_object, talls_object):
+    model = openmc.model.Model(geometry = geom_object, materials = mats_object, settings = sets_object, tallies = talls_object)
+    return model

WC_Layers/Photon_TallytoVtk.py

@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct  3 23:18:20 2024
+
+@author: Anupama Rajendra
+"""
+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):
+    with openmc.StatePoint(sp_filename) as sp:
+        flux_spectrum = sp.get_tally(id=flux_spectrum_tally_id) 
+        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)   
+        #Vitamin-J energy filter:
+        e_filter = flux_spectrum.filters[energy_filter_index]
+        #Lower bounds of the energy bins
+        e_filter_lower = e_filter.bins[:, 0]
+    return e_filter_lower, flux_sum_energy, tally_data_reshaped, mesh
+
+def plot_flux_data(e_filter_lower, flux_sum_energy, figure_filename):
+    # Plot flux spectrum
+    fix, ax = plt.subplots()
+    ax.loglog(e_filter_lower, flux_sum_energy, drawstyle='steps-post')
+    ax.set_xlabel('Energy [eV]')
+    ax.set_ylabel('Flux [photon-cm/source]')
+    ax.grid(True, which='both')
+    plt.savefig(figure_filename)
+    plt.show()
+
+def write_to_vtk(tally_data_reshaped, summed_axes_mesh_filter, vtk_filename, mesh):
+    mesh_data = tally_data_reshaped.sum(axis=summed_axes_mesh_filter)
+    mesh.write_data_to_vtk(filename=vtk_filename, datasets={"mean":mesh_data.flatten()})

WC_Layers/Source_Mesh_Reader.py

@@ -0,0 +1,47 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Sep 30 15:16:44 2024
+
+@author: Anupama Rajendra
+"""
+
+import h5py
+import numpy as np
+
+source_density_name = 'source_density'
+
+# Open the HDF5 files
+File_1 = h5py.File("source_mesh_1.h5m", 'r')
+File_2 = h5py.File("source_mesh_2.h5m", 'r')
+Files = [File_1, File_2]
+
+def extract_source_data(file_list):
+    '''
+    Identifies the location of the source density dataset within each mesh file.
+
+    input: list of .h5m mesh files (opened by h5py) containing photon source information
+    output: numpy array of source density data with rows = # of mesh elements and 
+    columns = # number of photon groups, with one array per source mesh
+    '''
+    sd_list = np.ndarray((len(file_list), num_elements, photon_groups))
+
+    for file_num, file in enumerate(file_list):
+        sd_list[file_num,:] = file['tstt']['elements']['Tet4']['tags']['source_density'][:]
+    return sd_list  
+
+def save_source_density(sd_list, sd_filename):
+    '''
+    Saves source density data as a text file.
+
+    sd_list: list of source density datasets from h5m files
+    sd_filename: user-provided filename that appears as a prefix for the text files
+    '''
+    for sd_index, source_density in enumerate(sd_list):
+        with open(f'{sd_filename}_{sd_index + 1}.txt', 'w') as source:
+            for tet_element in source_density:
+                source.write(' '.join(map(str, tet_element)) + '\n')
+
+def extract_save_sd(file_list, sd_filename):
+    data_extractor = extract_source_data(file_list)
+    save_source_density(data_extractor, sd_filename)
+    return data_extractor

WC_Layers/TwoLayers_Geometry.py

@@ -0,0 +1,30 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Oct 11 15:56:57 2024
+
+@author: Anupama Rajendra
+"""
+
+import openmc
+
+# Create geometry
+#Spherical shell:
+def make_spherical_shells(layers, inner_radius):    
+    '''
+    Creates a set of concentric spherical shells, each with its own material & inner/outer radius.
+
+    layers: list of tuples with OpenMC Material name and thickness: (material, thickness)
+    inner_radius: the radius of the innermost spherical shell
+    '''
+    inner_sphere = openmc.Sphere(r = inner_radius)
+    cells = [openmc.Cell(fill = None, region = -inner_sphere)]
+    for (material, thickness) in layers:
+        outer_radius = inner_radius + thickness
+        outer_sphere = openmc.Sphere(r = outer_radius)
+        cells.append(openmc.Cell(fill = material, region = +inner_sphere & -outer_sphere))
+        outer_radius = inner_radius
+        outer_sphere = inner_sphere
+    outer_sphere.boundary_type = 'vacuum'     
+    cells.append(openmc.Cell(fill = None, region = +outer_sphere)) 
+    geometry = openmc.Geometry(cells)    
+    return geometry

WC_Layers/TwoLayers_Materials.py

@@ -0,0 +1,33 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Oct 11 13:49:23 2024
+
+@author: Anupama Rajendra
+"""
+import openmc
+
+#ALARA Element Library to read density for specified element:
+def alara_element_densities(filepath):  
+    with open(""+filepath+"") as ALARA_Lib:
+        libLines = ALARA_Lib.readlines()
+    num_lines = len(libLines)
+    density_dict = {}
+    line_num = 0
+    while line_num < num_lines:
+        element_data = libLines[line_num].strip().split()
+        element_name = element_data[0].lower()
+        density_dict[element_name] = float(element_data[3])
+        line_num += int(element_data[4]) + 1
+    return density_dict
+
+# Create materials & export to XML:
+#Simulating tungsten shell:
+def make_element(elements, density_dict):
+    mats = openmc.Materials([])
+    for element_id, element in enumerate(elements):
+        mat = openmc.Material(material_id=element_id+1, name=element)
+        mat.add_element(element, 1.00)
+        mat.set_density('g/cm3', density_dict.get(element.lower()))
+        mats.append(mat)
+    return mats
+