diff --git a/SS_Activation_OpenMC.py b/SS_Activation_OpenMC.py deleted file mode 100644 index 7e2e185..0000000 --- a/SS_Activation_OpenMC.py +++ /dev/null @@ -1,170 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu Feb 8 08:16:12 2024 - -@author: Anupama Rajendra -""" -import openmc -import openmc.deplete -import numpy as np -import matplotlib.pyplot as plt -import matplotlib.colors as mcolors -import pandas as pd -import random - -# Importing Vitamin-J energy group structure: -# This excel file contains the energy bounds of the Vitamin J structure -# Vit_J = pd.read_excel('VitaminJEnergyGroupStructure.xlsx') -# ebounds = Vit_J.iloc[:, 1] - -openmc.config['chain_file'] = 'chain_endfb71_sfr.xml' - -# Create materials & export to XML: -#Simulating tungsten shell: -W = openmc.Material(material_id=1, name='W_Shell') -W.set_density('g/cm3', 19.28) -W.add_element('W', 1.00) -materials = openmc.Materials([W]) -materials.export_to_xml() - -# Create geometry -#Spherical shell: - -R_1 = 1000 -S_1= openmc.Sphere(r=R_1) #sphere of radius 1000cm -inside_sphere_1 = -S_1 -outside_sphere_1 = +S_1 -R_2 = 1005 -S_2 = openmc.Sphere(r=R_2, boundary_type='vacuum') -inside_sphere_2 = -S_2 -outside_sphere_2 = +S_2 -S_3 = outside_sphere_1 & inside_sphere_2 #filled with tungsten - -# Mapping materials to geometry: -Void = openmc.Cell(fill=None, region = inside_sphere_1) -Shell = openmc.Cell(fill=W, region=S_3) -Cells = [Void, Shell] -geometry = openmc.Geometry(Cells) -geometry.export_to_xml() - -# Source distribution: -PointSource = openmc.stats.Point(xyz=(0.0, 0.0, 0.0)) -Prob = openmc.stats.Discrete(14E+06, 1.0) - -# Assign simulation settings -settings = openmc.Settings() -settings.batches = 10 -settings.inactive = 1 -settings.particles = 10000 -settings.source = openmc.Source(space=PointSource, energy=Prob, strength = 1.0, particle = 'neutron') -settings.run_mode = 'fixed source' -settings.export_to_xml() - -# Define tallies -neutron_tally = openmc.Tally(name="Neutron tally") -neutron_tally.scores = ['flux', 'elastic', 'absorption'] -# Implementing filter for neutron tally through W shell -cell_filter = openmc.CellFilter([Shell]) -neutron_tally.filters = [cell_filter] - -# Creating a tally to get the flux energy spectrum. -# An energy filter is created to assign to the flux tally using the Vitamin J structure. -energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-175") - -spectrum_tally = openmc.Tally(name="Flux spectrum") -# Implementing energy and cell filters for flux spectrum tally -spectrum_tally.filters = [energy_filter_flux, cell_filter] -spectrum_tally.scores = ['flux'] - -# Collecting and exporting tallies to .xml -tallies = openmc.Tallies([neutron_tally, spectrum_tally]) -tallies.export_to_xml() - -model = openmc.model.Model(geometry=geometry,settings=settings) -#Depletion calculation -W.depletable = True -W.volume = 4.0/3.0 * np.pi * (R_2**3 - R_1**3) #volume of W wall material -fluxes, micros = openmc.deplete.get_microxs_and_flux(model, Cells) -operator = openmc.deplete.IndependentOperator(materials, fluxes[0:1], micros[0:1],normalization_mode='source-rate') -# operator = openmc.deplete.CoupledOperator(model, normalization_mode='source-rate') -time_steps = [3e8, 86400, 2.6e6] -source_rates = [1E+18, 0, 0] -integrator = openmc.deplete.PredictorIntegrator(operator=operator, timesteps=time_steps, source_rates=source_rates, timestep_units='s') -integrator.integrate() - -#Opening statepoint file to read tallies: -with openmc.StatePoint('statepoint.10.h5') as sp: - fl = sp.get_tally(name="Flux spectrum") - nt = sp.get_tally(name="Neutron tally") - -# Get the neutron energies from the energy filter -energy_filter_fl = fl.filters[0] -energies_fl = energy_filter_fl.bins[:, 0] - -# Get the neutron flux values -flux = fl.get_values(value='mean').ravel() - -#Neutron flux/elastic/absorption tallies: -tal = nt.get_values(value='mean').ravel() -print(tal) - -Flux_Data = np.c_[energies_fl, flux] -#Creating an excel file that stores flux data for each energy bin (used as input for ALARA) -#Dividing by volume to obtain proper units of flux (#/cm^2-s) -Flux_Data = np.c_[energies_fl, flux/W.volume] -#ALARA flux inputs go from high energy to low energy -Flux_Data_ALARA = Flux_Data[::-1] -FD_CSV = pd.DataFrame(Flux_Data_ALARA, columns=['Energy [eV]', 'Flux [n-cm/sp]']) -FD_CSV.to_csv('Neutron_Flux.csv', index=False) - -Tallies_CSV = pd.DataFrame(tal) -#Creating a csv file that stores total tally value data -Tallies_CSV.to_csv('Tally_Values.csv', index=False) - -# Depletion results file -results = openmc.deplete.Results(filename='depletion_results.h5') - -# Stable W nuclides present at beginning of operation (will not be plotted) -stable_nuc = ['W180', 'W182', 'W183', 'W184', 'W186'] - -#Store list of nuclides from last timestep as a Materials object -materials_last = results.export_to_materials(-1) -# Storing depletion data from 1st material -mat_dep = materials_last[0] -# Obtaining the list of nuclides from the results file -nuc_last = mat_dep.get_nuclides() - -# Removing stable W nuclides from list so that they do not appear in the plot -for j in stable_nuc : - nuc_last.remove(j) -print(nuc_last) - -colors = list(mcolors.CSS4_COLORS.keys()) -num_dens= {} -pair_list = {} - -with open(r'Densities_CSV.csv', 'a') as density_file: - - for nuclide in nuc_last: - plot_color = random.choice(colors) - time, num_dens[nuclide] = results.get_atoms('1', nuclide, nuc_units = 'atom/cm3') - print(time, num_dens[nuclide]) - density_file.write(f'{nuclide},') - density_file.write(','.join(map(str, num_dens[nuclide])) + '\n') - plt.plot(time, num_dens[nuclide], marker='.', linestyle='solid', color=plot_color, label=nuclide) - -# Adding labels and title -plt.xlabel('Time after beginning of operation [s]') -plt.xlim(1, sum(time_steps) -#plt.gca().set_ylim(bottom=0) -plt.ylabel('Nuclide density [atoms/cm^3]') -plt.xscale("log") -plt.yscale("log") -plt.title('Plot of number density vs time') - -# Adding a legend -plt.legend() - -plt.savefig('Nuclide_density_OpenMC') -# Display the plot -plt.show() diff --git a/SphericalShell/ALARA_OpenMC_Comparison.py b/SphericalShell/ALARA_OpenMC_Comparison.py new file mode 100644 index 0000000..5dfdba0 --- /dev/null +++ b/SphericalShell/ALARA_OpenMC_Comparison.py @@ -0,0 +1,164 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Jul 10 12:54:10 2024 + +@author: Anupama Rajendra +""" + +import matplotlib.pyplot as plt +import csv + +#Initializing arrays to store data from OpenMC results .csv file: +OpenMC_Nuclides = [] +OpenMC_Day = [] +OpenMC_Month = [] +OpenMC_Shutdown = [] + +with open('Densities_CSV.csv', 'r') as OpenMC_Data: + OpenMC_csv_reader = csv.reader(OpenMC_Data) + for OpenMC_line in OpenMC_csv_reader: + OpenMC_Nuclides.append(OpenMC_line[0]) + OpenMC_Day.append(OpenMC_line[2]) + OpenMC_Month.append(OpenMC_line[3]) + OpenMC_Shutdown.append(OpenMC_line[4]) + +# Reading content from ALARA Output file: +with open('ALARA_Data.txt', 'r') as Data_File: + ALARA_Data = Data_File.readlines() + +#Identifying the part of the ALARA output file that contains the relevant density data: +ALARA_FileBounds = [] +for Density_Index, Density_Content in enumerate(ALARA_Data): + if '==' in Density_Content: + ALARA_FileBounds.append(Density_Index) + #If '==' is found twice, end the loop + if len(ALARA_FileBounds) == 2: + break + +#Process the data in between the two lines that begin with '==' +ALARA_Nuclide_Data = ALARA_Data[ALARA_FileBounds[0] + 1:ALARA_FileBounds[1]] + +# Stable W nuclides present at beginning of operation: +Stable_Nuc = ['w-180', 'w-182', 'w-183', 'w-184', 'w-186'] + +#Initializing arrays to store data from ALARA and data common to both ALARA & OpenMC +ALARA_Nuc_no_W = [] +ALARA_no_W_Day = [] +ALARA_no_W_Month = [] +ALARA_no_W_Shutdown = [] +ALARA_List = [ALARA_no_W_Day, ALARA_no_W_Month, ALARA_no_W_Shutdown] + +Common_Nuclides = [] +Common_Day_Diff = [] +Common_Month_Diff = [] +Common_Shutdown_Diff = [] + +#Nuclide data found in ALARA only +ALARA_Nuclides_Only = [] +ALARA_no_W_Day_Only = [] +ALARA_no_W_Month_Only = [] +ALARA_no_W_Shutdown_Only = [] + +for ALARA_Filtered_Lines in ALARA_Nuclide_Data: + #For any nuclide that is not one of the stable W nuclides: + if not any(ALARA_Filtered_Lines.startswith(Nuclide) for Nuclide in Stable_Nuc): + #Formatting change to make OpenMC and ALARA nuclide formats match up + ALARA_Filtered_Lines = ALARA_Filtered_Lines.replace('-','').capitalize() + #Storing all density information from ALARA output + + #The first part of each line is the nuclide + Nuc_Each_Line = (ALARA_Filtered_Lines.strip().split()[0]) + ALARA_Nuc_no_W.append(Nuc_Each_Line) + ALARA_no_W_Day.append(ALARA_Filtered_Lines.strip().split()[1]) + ALARA_no_W_Month.append(ALARA_Filtered_Lines.strip().split()[2]) + ALARA_no_W_Shutdown.append(ALARA_Filtered_Lines.strip().split()[3]) + #Identifying nuclides from ALARA that are also found from OpenMC data: + if Nuc_Each_Line in OpenMC_Nuclides: + Common_Nuclides.append(Nuc_Each_Line) + else: + #Add to lists that contain nuclides/densities only found in ALARA output + ALARA_Nuclides_Only.append(Nuc_Each_Line) + ALARA_no_W_Day_Only.append(ALARA_Filtered_Lines.split()[1]) + ALARA_no_W_Month_Only.append(ALARA_Filtered_Lines.split()[2]) + ALARA_no_W_Shutdown_Only.append(ALARA_Filtered_Lines.split()[3]) + +# Times after shutdown [s] +time_steps = [0, 86400, 2.6864e+06] + +# Plot data for each isotope found in ALARA +for i, isotope in enumerate(ALARA_Nuc_no_W): + densities = [ALARA_no_W_Day[i], ALARA_no_W_Month[i], ALARA_no_W_Shutdown[i]] + plt.plot(time_steps, densities, marker='o', label=isotope) + +plt.xlabel('Time (seconds)') +plt.ylabel('Number Density (atoms/cm^3)') +plt.title('Number Density [atoms/cm^3] vs. Time After Shutdown [s]') +plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1)) +plt.subplots_adjust(right=0.7) +plt.xscale('log') +plt.yscale('log') +plt.grid(True) +plt.savefig('Nuclide_density_ALARA') +plt.show() + +#Initializing arrays to store ratios of nuclide densities from both codes +Ratio_Day = [] +Ratio_Month = [] +Ratio_Shutdown = [] + +#print("All nuclides found in OpenMC", OpenMC_Nuclides) +print("Nuclides common to ALARA and OpenMC", Common_Nuclides) +#print("Nuclides found only in ALARA:", ALARA_Nuclides_Only) + +#Iterating over all nuclides common to both sets of results +for Common_Name in Common_Nuclides: + #Identifying the location of each common nuclide in the OpenMC nuclide list + Index_OpenMC = OpenMC_Nuclides.index(Common_Name) + Density_OpenMC_Day = float(OpenMC_Day[Index_OpenMC]) + Density_OpenMC_Month = float(OpenMC_Month[Index_OpenMC]) + Density_OpenMC_Shutdown = float(OpenMC_Shutdown[Index_OpenMC]) + + #Identifying the location of each common nuclide in the ALARA nuclide list + Index_ALARA = ALARA_Nuc_no_W.index(Common_Name) + Density_ALARA_Day = float(ALARA_no_W_Day[Index_ALARA]) + Density_ALARA_Month = float(ALARA_no_W_Month[Index_ALARA]) + Density_ALARA_Shutdown = float(ALARA_no_W_Shutdown[Index_ALARA]) + + #Ratios between the values from ALARA and OpenMC: + Ratio_Day.append(Density_OpenMC_Day / Density_ALARA_Day) + Ratio_Month.append(Density_OpenMC_Month / Density_ALARA_Month) + Ratio_Shutdown.append(Density_OpenMC_Shutdown / Density_ALARA_Shutdown) + + #Absolute values of the differences between ALARA and OpenMC + Common_Day_Diff.append(abs(Density_OpenMC_Day - Density_ALARA_Day)) + Common_Month_Diff.append(abs(Density_OpenMC_Month - Density_ALARA_Month)) + Common_Shutdown_Diff.append(abs(Density_OpenMC_Shutdown - Density_ALARA_Shutdown)) + +for j, com_iso in enumerate(Common_Nuclides): + common_diff = [Common_Day_Diff[j], Common_Month_Diff[j], Common_Shutdown_Diff[j]] + common_ratios = [Ratio_Day[j], Ratio_Month[j], Ratio_Shutdown[j]] + #plt.plot(time_steps, common_diff, marker='o', label=com_iso) + plt.plot(time_steps, common_ratios, marker = 'o', label = com_iso) + +# plt.xlabel('Time [seconds]') +# plt.xlim(1, 1E+7) +# #plt.ylim(1E+12, 1E+15) +# plt.ylabel('Number Density Difference [atoms/cm^3]') +# plt.title('Number Density [atoms/cm^3] vs. Time After Shutdown [s]') +# plt.subplots_adjust(right=0.7) +# plt.xscale('log') +# plt.yscale('log') +# plt.grid(True) +#plt.savefig('Nuclide_density_diff') +# plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1)) +# plt.show() + +plt.xlabel('Time [seconds]') +plt.ylabel('Ratio Between OpenMC and ALARA') +plt.title('Number Density Ratio vs. Time After Shutdown [s]') +plt.subplots_adjust(right=0.7) +plt.legend(loc='upper left', bbox_to_anchor=(1.15, 1.1)) +plt.grid(True) +plt.savefig('Nuclide_density_ratio') +plt.show() +plt.close()