Photon Transport Simulation with Two Layers of Material & Unstructured Mesh

WC_Layers/OpenMC_PhotonTransport.py

@@ -0,0 +1,138 @@
+# -*- 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 summed_strengths, strengths_list
+from TwoLayers_Materials import *
+from TwoLayers_Geometry import *
+
+Mesh_File = 'OpenMC_Mesh.h5m'
+
+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 = [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]
+
+#Define source:
+def make_source(W_Shell, C_Shell, Cells):
+    Source_List = []
+    Particle_Filter = openmc.ParticleFilter('photon')
+    Total_Mesh = openmc.UnstructuredMesh(Mesh_File, library='moab')
+    for index, bound in enumerate(bounds[:-1]):
+        Mesh_Dist = openmc.stats.MeshSpatial(Total_Mesh, strengths=strengths_list[index], volume_normalized=False)
+        Energy_Dist = openmc.stats.Uniform(a=bounds[index], b=bounds[index + 1])
+        #Source strengths given by strengths_list created by Source_Mesh_Reader
+        Source = Source_List.append(openmc.IndependentSource(space=Mesh_Dist, energy=Energy_Dist, strength=summed_strengths[index], particle='photon', domains=Cells))
+    return Source_List, Particle_Filter, Total_Mesh, Mesh_Dist
+
+# Define tallies
+def tallies(W_Shell, C_Shell, Particle_Filter, Total_Mesh):
+    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 W shell
+    cell_filter = openmc.CellFilter([W_Shell, C_Shell])
+    neutron_tally.filters = [cell_filter, Total_Filter]
+
+    # Creating a tally to get the flux energy spectrum.
+    # An energy filter is 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")
+    # Implementing energy and cell filters for flux spectrum tally
+    spectrum_tally.filters = [cell_filter, Total_Filter, energy_filter_flux, Particle_Filter]
+    spectrum_tally.scores = ['flux']
+
+    tall = openmc.Tallies([neutron_tally, spectrum_tally])
+    tall.export_to_xml()
+    return tall, neutron_tally, spectrum_tally, Total_Filter, cell_filter
+
+# Assign simulation settings
+def settings(Source_List):
+    sets = openmc.Settings()
+    sets.batches = 10
+    sets.inactive = 1
+    sets.particles = 100000
+    sets.source = Source_List
+    sets.run_mode = 'fixed source'
+    sets.export_to_xml()
+    return sets
+
+def plot_universe(Void, W_Shell, C_Shell, Cells):
+    universe_ss = openmc.Universe(cells=Cells)
+    Plot = universe_ss.plot(width=(2500.0, 2500.0), basis='xz',
+              colors={Void: 'blue', W_Shell: 'red', C_Shell: 'green'}, legend=True)
+    Plot.figure.savefig('Universe')
+    return universe_ss
+
+# Exporting materials, geometry, and tallies to .xml
+def export_to_xml(filepath, element_1, element_2, inner_radius_W, outer_radius_W, inner_radius_C):
+    OpenMC_SF = mat_lib(filepath)
+    materials = []
+    for material_id, element in enumerate(elements):
+         materials.append(make_element(element, material_id+1, OpenMC_SF))
+    OpenMC_Mat = all_mat(OpenMC_W, OpenMC_C)
+    OpenMC_Geometry = make_spherical_shell(inner_radius_W, outer_radius_W, inner_radius_C, OpenMC_W, OpenMC_C)
+    OpenMC_Source = make_source(OpenMC_Geometry[2], OpenMC_Geometry[3], OpenMC_Geometry[4])
+    OpenMC_Settings = settings(OpenMC_Source[0])
+    OpenMC_Tallies = tallies(OpenMC_Geometry[2], OpenMC_Geometry[3], OpenMC_Source[1], OpenMC_Source[2])
+    OpenMC_Universe = plot_universe(OpenMC_Geometry[1], OpenMC_Geometry[2], OpenMC_Geometry[3], OpenMC_Geometry[4])
+    #return OpenMC_Materials, OpenMC_Geometry, OpenMC_Source, OpenMC_Settings, OpenMC_Tallies, OpenMC_Universe
+    return OpenMC_SF, OpenMC_W, OpenMC_C, *OpenMC_Geometry, *OpenMC_Source, OpenMC_Settings, *OpenMC_Tallies
+
+Lib_Lines, M_1, M_2, geometry, Void, W_Shell, C_Shell, Cells, Source_List, Particle_Filter, Total_Mesh, Mesh_Dist, tall, neutron_tally, spectrum_tally, Total_Filter, cell_filter, sets = export_to_xml(fp, E_1, E_2, R_W_1, R_W_2, R_C_1)

WC_Layers/Photon_TallytoVtk.py

@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct  3 23:18:20 2024
+
+@author: Anupama Rajendra
+"""
+import openmc
+import matplotlib.pyplot as plt
+
+with openmc.StatePoint("statepoint.10.h5") as sp:
+    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]
+
+    # 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()})

WC_Layers/Source_Mesh_Reader.py

@@ -0,0 +1,41 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Sep 30 15:16:44 2024
+
+@author: Anupama Rajendra
+"""
+
+import h5py
+import numpy as np
+
+# 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]
+
+# Number of photon groups
+photon_groups = 24
+
+# Process each file
+for file_index, file in enumerate(Files):
+    # Extract data from the current file
+    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)
+

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
+