From 34893ab9a0bfe3866a3936ccee8fa6cc9b107e20 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 08:12:08 -0500 Subject: [PATCH 01/27] Add files via upload --- WC_Layers/OpenMC_R2S.py | 284 ++++++++++++++++++++++++++++++++++++++++ WC_Layers/R2S.yaml | 92 +++++++++++++ 2 files changed, 376 insertions(+) create mode 100644 WC_Layers/OpenMC_R2S.py create mode 100644 WC_Layers/R2S.yaml diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py new file mode 100644 index 0000000..a673edb --- /dev/null +++ b/WC_Layers/OpenMC_R2S.py @@ -0,0 +1,284 @@ +import openmc +import openmc.deplete +from pathlib import Path +import argparse +import yaml +import numpy as np +import h5py + +def alara_element_densities(alara_fp): + ''' + Creates a dictionary where keys = element names (str) and values = element density (float) + inputs: + alara_filepath : path to ALARA element library + ''' + with open(alara_fp) 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 + +def make_materials(elements, density_dict): + ''' + Creates an OpenMC Materials object using user-specified elements + inputs: + elements: iterable of element names (str) + density_dict: dictionary with keys = element names (str) and values = element density (float) + ''' + 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 + +def make_spherical_shells(inner_radius, layers, outer_boundary_type): + ''' + Creates a set of concentric spherical shells, each with its own material & inner/outer radius. + inputs: + inner_radius: the radius of the innermost spherical shell + layers: iterable of tuples of OpenMC Material object and its respective thickness (float) + ''' + 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)) + inner_radius = outer_radius + inner_sphere = outer_sphere + material.volume = 4.0/3.0 * np.pi * ((inner_radius + thickness)**3 - (inner_radius)**3) + inner_radius = inner_radius + thickness + outer_sphere.boundary_type = outer_boundary_type + cells.append(openmc.Cell(fill = None, region = +outer_sphere)) + geometry = openmc.Geometry(cells) + return geometry + +def make_neutron_source(energy): + point_source = openmc.stats.Point(xyz=(0.0, 0.0, 0.0)) + energy_dist = openmc.stats.Discrete(energy, 1.0) + neutron_source = openmc.Source(space = point_source, energy = energy_dist, strength = 1.0) + return neutron_source + +def make_settings(neutron_source, total_batches, inactive_batches, num_particles, run_mode): + sets = openmc.Settings() + sets.batches = total_batches + sets.inactive = inactive_batches + sets.particles = num_particles + sets.source = neutron_source + sets.run_mode = run_mode + return sets + +#Only executed if external geometry and materials are imported +def make_depletion_volumes(neutron_model, mesh_file): + materials = neutron_model.materials + mesh_file = Path(mesh_file).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=25000000) + #returns dict-like object that maps mat ids to array of volumes equal to # of mesh elements + + total_volumes = {} + for mat_id, volumes in mat_vols.items(): + total_volumes[mat_id] = np.sum(volumes) + for material in materials: + material.volume = total_volumes[material.id] + return neutron_model + +def deplete_wc(neutron_model, mesh_file, chain_file, timesteps, source_rates, norm_mode, timestep_units): + materials = neutron_model.materials + mesh_file = Path(mesh_file).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + model_nuclide_names = [] + for material in materials: + material.depletable = True + material_nuclides = material.nuclides + for material_nuclide in material_nuclides: + model_nuclide_names.append(material_nuclide.name) + + activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=7000000) + activation_mats_object = openmc.Materials(activation_mats) + activation_mats_object.export_to_xml("Activation_Materials.xml") + + depletion_chain = openmc.deplete.Chain.from_xml(Path(chain_file).resolve()) + + nuclide_names = [] + for nuclide in depletion_chain.nuclides: + if nuclide.name in model_nuclide_names: + nuclide_names.append(nuclide.name) + + fluxes, micros = openmc.deplete.get_microxs_and_flux(neutron_model, unstructured_mesh, nuclide_names, + depletion_chain.reactions, + openmc.mgxs.GROUP_STRUCTURES["VITAMIN-J-175"]) + operator = openmc.deplete.IndependentOperator(openmc.Materials(activation_mats), fluxes, micros, chain_file=chain_file, normalization_mode = norm_mode) + copy_timesteps = list(timesteps) + copy_source_rates = list(source_rates) + integrator_list = [] + for timestep in timesteps: + while copy_source_rates[-1] == 0: + integrator = openmc.deplete.PredictorIntegrator(operator, copy_timesteps, source_rates = copy_source_rates, timestep_units = timestep_units) + integrator_list.append(integrator) + copy_timesteps.pop() + copy_source_rates.pop() + + #Letting integrator_list go from the fewest number of decay intervals to the most + integrator_list.reverse() + for int_index, integrator in enumerate(integrator_list): + integrator.integrate(path=f"depletion_results_decay_set_{int_index}.h5") + return activation_mats, unstructured_mesh, neutron_model + +def make_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs): + sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] + for int_index in range(num_cooling_steps): + results = openmc.deplete.Results(f"depletion_results_decay_set_{int_index}.h5") + photon_sources = np.empty(sd_data.shape[0], dtype=object) + activated_mats_list = set(results[-1].index_mat.keys()) + + for mat_index, mat in enumerate(activation_mats) : + if str(mat.id) in activated_mats_list : + mat = results[-1].get_material(str(mat.id)) + energy = mat.get_decay_photon_energy() + if energy == None: + photon_source = openmc.IndependentSource( + energy = energy, + particle = 'photon', + strength = 0.0) + else: + photon_source = openmc.IndependentSource( + energy = energy, + particle = 'photon', + strength = energy.integral()) + + else: + photon_source = openmc.IndependentSource( + energy = openmc.stats.Discrete(0, 1.0), + particle = 'photon', + strength = 0.0) + + photon_sources[mat_index] = photon_source + photon_model = neutron_model + photon_model.settings.source = openmc.MeshSource(unstructured_mesh, photon_sources) + photon_model.settings.export_to_xml(f'settings_{int_index}.xml') + + return photon_model + +def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): + dose_energy, dose = openmc.data.dose_coefficients('photon', geometry=coeff_geom) + dose_filter = openmc.EnergyFunctionFilter(dose_energy, dose) + spherical_mesh = openmc.SphericalMesh(np.arange(0, 3800, 10), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") + spherical_mesh_filter = openmc.MeshFilter(spherical_mesh) + #Could also use an unstructured mesh filter here - spherical mesh filter for visualization purposes + particle_filter = openmc.ParticleFilter('photon') + energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-42") + + flux_tally = openmc.Tally(tally_id=1, name="photon flux tally") + flux_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter] + flux_tally.scores = ['flux'] + + dose_tally = openmc.Tally(tally_id=2, name="dose tally") + dose_tally.filters = [spherical_mesh_filter, dose_filter, energy_filter_flux, particle_filter] + dose_tally.scores = ['flux'] + + photon_model.tallies = [flux_tally, dose_tally] + #Change boundary condition: + list(photon_model.geometry.get_all_surfaces().keys())[-1].boundary_type="white" + + for int_index in range(num_cooling_steps): + #Reassign settings (which contains source) for each decay step + photon_model.settings = openmc.Settings.from_xml(f'settings_{int_index}.xml') + photon_model.export_to_model_xml(path=f'photon_model_{int_index}.xml') + sp_path = photon_model.run(f'photon_model_{int_index}.xml') + sp_path.rename(f'statepoint_openmc_{int_index}.h5') + +#---------------------------------------------------------------------------------- +#Define variables and execute all functions: + +def parse_args(): + parser = argparse.ArgumentParser() + parser.add_argument('--OpenMC_WC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs for WC_Neutron_Transport") + parser.add_argument('--ext_mat_geom', default = False, help="Specify whether materials and geometry come from external model") + parser.add_argument('--photon_transport', default = False, help="Run only photon transport (no depletion)") + args = parser.parse_args() + return args + +def read_yaml(args): + with open(args.OpenMC_WC_YAML, 'r') as transport_file: + inputs = yaml.safe_load(transport_file) + return inputs + +def create_materials_obj(inputs): + densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) + materials = make_materials(inputs['mat_info']['element_list'], + densities) + return materials + +def create_geometry_obj(materials, inputs): + geom_info = inputs['geom_info'] + layers = zip(materials, geom_info['thicknesses']) + geometry = make_spherical_shells(geom_info['inner_radius'], + layers, + geom_info['outer_boundary_type']) + return geometry + +def create_neutron_model(inputs, materials, geometry): + settings_info = inputs['settings_info'] + neutron_source = make_neutron_source(inputs['particle_energy']) + settings = make_settings(neutron_source, + settings_info['total_batches'], + settings_info['inactive_batches'], + settings_info['num_particles'], + settings_info['run_mode']) + neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = settings) + neutron_model.export_to_model_xml("neutron_model.xml") + return neutron_model + +def run_depletion(inputs, neutron_model): + #geom_info = inputs['geom_info'] + dep_params = inputs['dep_params'] + + activation_mats, unstructured_mesh, integrator, neutron_model = deplete_wc(neutron_model, + inputs['filename_dict']['mesh_file'], + dep_params['chain_file'], + dep_params['timesteps'], + dep_params['source_rates'], + dep_params['norm_mode'], + dep_params['timestep_units']) + return activation_mats, unstructured_mesh, integrator, neutron_model + +def main(): + args = parse_args() + inputs = read_yaml(args) + + openmc.config['chain_file'] = inputs['dep_params']['chain_file'] + if args.ext_mat_geom == True : #Import materials and geometry from external model + ext_model = openmc.model.Model.from_model_xml(inputs['filename_dict']['ext_model']) + materials = ext_model.materials + geometry = ext_model.geometry + else: + materials = create_materials_obj(inputs) + geometry = create_geometry_obj(materials, inputs) + + neutron_model = create_neutron_model(inputs, materials, geometry) + + #Set to True to run photon transport only (no depletion): + if args.photon_transport == True: + activation_mats = openmc.Materials.from_xml("Activation_Materials.xml") + mesh_file = Path(inputs['filename_dict']['mesh_file']).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + else: + if args.ext_mat_geom == True : + neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) + activation_mats, unstructured_mesh, neutron_model = run_depletion(inputs, neutron_model) + + num_cooling_steps = (inputs['dep_params']['source_rates']).count(0) + photon_model = make_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) + photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/WC_Layers/R2S.yaml b/WC_Layers/R2S.yaml new file mode 100644 index 0000000..9e57c67 --- /dev/null +++ b/WC_Layers/R2S.yaml @@ -0,0 +1,92 @@ +filename_dict : + elelib_fp : elelib.std + mesh_file : rbt_mesh_w_blocks_materials_no_graveyard.h5 + sp_filename : statepoint.10.h5 + figure_filename : Photon_flux_vs_energy + vtk_filename : Photon_Flux.vtk + photon_tally_figname : photon_tally + ext_model : HCPB_model_photons.xml + +mat_info : + element_list : #add in order of radially innermost to outermost + - W + - C + +geom_info : + inner_radius : 995 + thicknesses : #corresponds to each of the elements in element_list + - 5 + - 5 + outer_boundary_type : vacuum + +particle_energy : 14.0E+06 + +settings_info : + total_batches : 10 + inactive_batches : 1 + num_particles : 10000 + run_mode : fixed source + +source_meshes : + - source_mesh_hcpb_1.h5m + - source_mesh_hcpb_2.h5m +sd_filename : source_density + +source_info : + phtn_e_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 + +file_indices: + source_mesh_index : 1 + flux_spectrum_tally_id : 2 + photon_tally_id : 1 + mesh_number : 1 + energy_filter_index : 2 + +tally_info : + tallied_elements : #change according to desired tally region/material + - W + - C + +axes_without_energy_bins : (0, 1, 3, 4, 5, 6) +axes_without_mesh : (0,2,3,4,5,6) +axes_without_dose_filter_bin : (0,1,2,3,4,5) +coeff_geom : 'AP' + +#Depletion parameters for OpenMC-only R2S +dep_params : + chain_file : chain_endfb71_sfr.xml + timesteps : #sequence of timesteps following the beginning of operation (not cumulative) + - 3.0E+8 + - 86400 + - 2.6E+6 + source_rates : + - 1.0E+18 + - 0 + - 0 + norm_mode : source-rate + timestep_units : s \ No newline at end of file From 58e2ab4d2b4302f34e9380de3c8128c31b22850a Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 08:12:48 -0500 Subject: [PATCH 02/27] Delete WC_Layers/R2S.yaml --- WC_Layers/R2S.yaml | 92 ---------------------------------------------- 1 file changed, 92 deletions(-) delete mode 100644 WC_Layers/R2S.yaml diff --git a/WC_Layers/R2S.yaml b/WC_Layers/R2S.yaml deleted file mode 100644 index 9e57c67..0000000 --- a/WC_Layers/R2S.yaml +++ /dev/null @@ -1,92 +0,0 @@ -filename_dict : - elelib_fp : elelib.std - mesh_file : rbt_mesh_w_blocks_materials_no_graveyard.h5 - sp_filename : statepoint.10.h5 - figure_filename : Photon_flux_vs_energy - vtk_filename : Photon_Flux.vtk - photon_tally_figname : photon_tally - ext_model : HCPB_model_photons.xml - -mat_info : - element_list : #add in order of radially innermost to outermost - - W - - C - -geom_info : - inner_radius : 995 - thicknesses : #corresponds to each of the elements in element_list - - 5 - - 5 - outer_boundary_type : vacuum - -particle_energy : 14.0E+06 - -settings_info : - total_batches : 10 - inactive_batches : 1 - num_particles : 10000 - run_mode : fixed source - -source_meshes : - - source_mesh_hcpb_1.h5m - - source_mesh_hcpb_2.h5m -sd_filename : source_density - -source_info : - phtn_e_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 - -file_indices: - source_mesh_index : 1 - flux_spectrum_tally_id : 2 - photon_tally_id : 1 - mesh_number : 1 - energy_filter_index : 2 - -tally_info : - tallied_elements : #change according to desired tally region/material - - W - - C - -axes_without_energy_bins : (0, 1, 3, 4, 5, 6) -axes_without_mesh : (0,2,3,4,5,6) -axes_without_dose_filter_bin : (0,1,2,3,4,5) -coeff_geom : 'AP' - -#Depletion parameters for OpenMC-only R2S -dep_params : - chain_file : chain_endfb71_sfr.xml - timesteps : #sequence of timesteps following the beginning of operation (not cumulative) - - 3.0E+8 - - 86400 - - 2.6E+6 - source_rates : - - 1.0E+18 - - 0 - - 0 - norm_mode : source-rate - timestep_units : s \ No newline at end of file From f6648479f8031694a5d9dccf7d9631b5153c7f84 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 08:13:01 -0500 Subject: [PATCH 03/27] Delete WC_Layers/OpenMC_R2S.py --- WC_Layers/OpenMC_R2S.py | 284 ---------------------------------------- 1 file changed, 284 deletions(-) delete mode 100644 WC_Layers/OpenMC_R2S.py diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py deleted file mode 100644 index a673edb..0000000 --- a/WC_Layers/OpenMC_R2S.py +++ /dev/null @@ -1,284 +0,0 @@ -import openmc -import openmc.deplete -from pathlib import Path -import argparse -import yaml -import numpy as np -import h5py - -def alara_element_densities(alara_fp): - ''' - Creates a dictionary where keys = element names (str) and values = element density (float) - inputs: - alara_filepath : path to ALARA element library - ''' - with open(alara_fp) 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 - -def make_materials(elements, density_dict): - ''' - Creates an OpenMC Materials object using user-specified elements - inputs: - elements: iterable of element names (str) - density_dict: dictionary with keys = element names (str) and values = element density (float) - ''' - 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 - -def make_spherical_shells(inner_radius, layers, outer_boundary_type): - ''' - Creates a set of concentric spherical shells, each with its own material & inner/outer radius. - inputs: - inner_radius: the radius of the innermost spherical shell - layers: iterable of tuples of OpenMC Material object and its respective thickness (float) - ''' - 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)) - inner_radius = outer_radius - inner_sphere = outer_sphere - material.volume = 4.0/3.0 * np.pi * ((inner_radius + thickness)**3 - (inner_radius)**3) - inner_radius = inner_radius + thickness - outer_sphere.boundary_type = outer_boundary_type - cells.append(openmc.Cell(fill = None, region = +outer_sphere)) - geometry = openmc.Geometry(cells) - return geometry - -def make_neutron_source(energy): - point_source = openmc.stats.Point(xyz=(0.0, 0.0, 0.0)) - energy_dist = openmc.stats.Discrete(energy, 1.0) - neutron_source = openmc.Source(space = point_source, energy = energy_dist, strength = 1.0) - return neutron_source - -def make_settings(neutron_source, total_batches, inactive_batches, num_particles, run_mode): - sets = openmc.Settings() - sets.batches = total_batches - sets.inactive = inactive_batches - sets.particles = num_particles - sets.source = neutron_source - sets.run_mode = run_mode - return sets - -#Only executed if external geometry and materials are imported -def make_depletion_volumes(neutron_model, mesh_file): - materials = neutron_model.materials - mesh_file = Path(mesh_file).resolve() - unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=25000000) - #returns dict-like object that maps mat ids to array of volumes equal to # of mesh elements - - total_volumes = {} - for mat_id, volumes in mat_vols.items(): - total_volumes[mat_id] = np.sum(volumes) - for material in materials: - material.volume = total_volumes[material.id] - return neutron_model - -def deplete_wc(neutron_model, mesh_file, chain_file, timesteps, source_rates, norm_mode, timestep_units): - materials = neutron_model.materials - mesh_file = Path(mesh_file).resolve() - unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - model_nuclide_names = [] - for material in materials: - material.depletable = True - material_nuclides = material.nuclides - for material_nuclide in material_nuclides: - model_nuclide_names.append(material_nuclide.name) - - activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=7000000) - activation_mats_object = openmc.Materials(activation_mats) - activation_mats_object.export_to_xml("Activation_Materials.xml") - - depletion_chain = openmc.deplete.Chain.from_xml(Path(chain_file).resolve()) - - nuclide_names = [] - for nuclide in depletion_chain.nuclides: - if nuclide.name in model_nuclide_names: - nuclide_names.append(nuclide.name) - - fluxes, micros = openmc.deplete.get_microxs_and_flux(neutron_model, unstructured_mesh, nuclide_names, - depletion_chain.reactions, - openmc.mgxs.GROUP_STRUCTURES["VITAMIN-J-175"]) - operator = openmc.deplete.IndependentOperator(openmc.Materials(activation_mats), fluxes, micros, chain_file=chain_file, normalization_mode = norm_mode) - copy_timesteps = list(timesteps) - copy_source_rates = list(source_rates) - integrator_list = [] - for timestep in timesteps: - while copy_source_rates[-1] == 0: - integrator = openmc.deplete.PredictorIntegrator(operator, copy_timesteps, source_rates = copy_source_rates, timestep_units = timestep_units) - integrator_list.append(integrator) - copy_timesteps.pop() - copy_source_rates.pop() - - #Letting integrator_list go from the fewest number of decay intervals to the most - integrator_list.reverse() - for int_index, integrator in enumerate(integrator_list): - integrator.integrate(path=f"depletion_results_decay_set_{int_index}.h5") - return activation_mats, unstructured_mesh, neutron_model - -def make_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs): - sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] - for int_index in range(num_cooling_steps): - results = openmc.deplete.Results(f"depletion_results_decay_set_{int_index}.h5") - photon_sources = np.empty(sd_data.shape[0], dtype=object) - activated_mats_list = set(results[-1].index_mat.keys()) - - for mat_index, mat in enumerate(activation_mats) : - if str(mat.id) in activated_mats_list : - mat = results[-1].get_material(str(mat.id)) - energy = mat.get_decay_photon_energy() - if energy == None: - photon_source = openmc.IndependentSource( - energy = energy, - particle = 'photon', - strength = 0.0) - else: - photon_source = openmc.IndependentSource( - energy = energy, - particle = 'photon', - strength = energy.integral()) - - else: - photon_source = openmc.IndependentSource( - energy = openmc.stats.Discrete(0, 1.0), - particle = 'photon', - strength = 0.0) - - photon_sources[mat_index] = photon_source - photon_model = neutron_model - photon_model.settings.source = openmc.MeshSource(unstructured_mesh, photon_sources) - photon_model.settings.export_to_xml(f'settings_{int_index}.xml') - - return photon_model - -def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): - dose_energy, dose = openmc.data.dose_coefficients('photon', geometry=coeff_geom) - dose_filter = openmc.EnergyFunctionFilter(dose_energy, dose) - spherical_mesh = openmc.SphericalMesh(np.arange(0, 3800, 10), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") - spherical_mesh_filter = openmc.MeshFilter(spherical_mesh) - #Could also use an unstructured mesh filter here - spherical mesh filter for visualization purposes - particle_filter = openmc.ParticleFilter('photon') - energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-42") - - flux_tally = openmc.Tally(tally_id=1, name="photon flux tally") - flux_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter] - flux_tally.scores = ['flux'] - - dose_tally = openmc.Tally(tally_id=2, name="dose tally") - dose_tally.filters = [spherical_mesh_filter, dose_filter, energy_filter_flux, particle_filter] - dose_tally.scores = ['flux'] - - photon_model.tallies = [flux_tally, dose_tally] - #Change boundary condition: - list(photon_model.geometry.get_all_surfaces().keys())[-1].boundary_type="white" - - for int_index in range(num_cooling_steps): - #Reassign settings (which contains source) for each decay step - photon_model.settings = openmc.Settings.from_xml(f'settings_{int_index}.xml') - photon_model.export_to_model_xml(path=f'photon_model_{int_index}.xml') - sp_path = photon_model.run(f'photon_model_{int_index}.xml') - sp_path.rename(f'statepoint_openmc_{int_index}.h5') - -#---------------------------------------------------------------------------------- -#Define variables and execute all functions: - -def parse_args(): - parser = argparse.ArgumentParser() - parser.add_argument('--OpenMC_WC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs for WC_Neutron_Transport") - parser.add_argument('--ext_mat_geom', default = False, help="Specify whether materials and geometry come from external model") - parser.add_argument('--photon_transport', default = False, help="Run only photon transport (no depletion)") - args = parser.parse_args() - return args - -def read_yaml(args): - with open(args.OpenMC_WC_YAML, 'r') as transport_file: - inputs = yaml.safe_load(transport_file) - return inputs - -def create_materials_obj(inputs): - densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) - materials = make_materials(inputs['mat_info']['element_list'], - densities) - return materials - -def create_geometry_obj(materials, inputs): - geom_info = inputs['geom_info'] - layers = zip(materials, geom_info['thicknesses']) - geometry = make_spherical_shells(geom_info['inner_radius'], - layers, - geom_info['outer_boundary_type']) - return geometry - -def create_neutron_model(inputs, materials, geometry): - settings_info = inputs['settings_info'] - neutron_source = make_neutron_source(inputs['particle_energy']) - settings = make_settings(neutron_source, - settings_info['total_batches'], - settings_info['inactive_batches'], - settings_info['num_particles'], - settings_info['run_mode']) - neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = settings) - neutron_model.export_to_model_xml("neutron_model.xml") - return neutron_model - -def run_depletion(inputs, neutron_model): - #geom_info = inputs['geom_info'] - dep_params = inputs['dep_params'] - - activation_mats, unstructured_mesh, integrator, neutron_model = deplete_wc(neutron_model, - inputs['filename_dict']['mesh_file'], - dep_params['chain_file'], - dep_params['timesteps'], - dep_params['source_rates'], - dep_params['norm_mode'], - dep_params['timestep_units']) - return activation_mats, unstructured_mesh, integrator, neutron_model - -def main(): - args = parse_args() - inputs = read_yaml(args) - - openmc.config['chain_file'] = inputs['dep_params']['chain_file'] - if args.ext_mat_geom == True : #Import materials and geometry from external model - ext_model = openmc.model.Model.from_model_xml(inputs['filename_dict']['ext_model']) - materials = ext_model.materials - geometry = ext_model.geometry - else: - materials = create_materials_obj(inputs) - geometry = create_geometry_obj(materials, inputs) - - neutron_model = create_neutron_model(inputs, materials, geometry) - - #Set to True to run photon transport only (no depletion): - if args.photon_transport == True: - activation_mats = openmc.Materials.from_xml("Activation_Materials.xml") - mesh_file = Path(inputs['filename_dict']['mesh_file']).resolve() - unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - else: - if args.ext_mat_geom == True : - neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) - activation_mats, unstructured_mesh, neutron_model = run_depletion(inputs, neutron_model) - - num_cooling_steps = (inputs['dep_params']['source_rates']).count(0) - photon_model = make_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) - photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) - -if __name__ == "__main__": - main() \ No newline at end of file From 127c9b9b52bf009c219ee625c0faeaacf90be26e Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 08:14:25 -0500 Subject: [PATCH 04/27] Add files via upload --- WC_Layers/OpenMC_R2S.py | 284 ++++++++++++++++++++++++++++++++++++++++ WC_Layers/R2S.yaml | 92 +++++++++++++ 2 files changed, 376 insertions(+) create mode 100644 WC_Layers/OpenMC_R2S.py create mode 100644 WC_Layers/R2S.yaml diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py new file mode 100644 index 0000000..a673edb --- /dev/null +++ b/WC_Layers/OpenMC_R2S.py @@ -0,0 +1,284 @@ +import openmc +import openmc.deplete +from pathlib import Path +import argparse +import yaml +import numpy as np +import h5py + +def alara_element_densities(alara_fp): + ''' + Creates a dictionary where keys = element names (str) and values = element density (float) + inputs: + alara_filepath : path to ALARA element library + ''' + with open(alara_fp) 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 + +def make_materials(elements, density_dict): + ''' + Creates an OpenMC Materials object using user-specified elements + inputs: + elements: iterable of element names (str) + density_dict: dictionary with keys = element names (str) and values = element density (float) + ''' + 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 + +def make_spherical_shells(inner_radius, layers, outer_boundary_type): + ''' + Creates a set of concentric spherical shells, each with its own material & inner/outer radius. + inputs: + inner_radius: the radius of the innermost spherical shell + layers: iterable of tuples of OpenMC Material object and its respective thickness (float) + ''' + 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)) + inner_radius = outer_radius + inner_sphere = outer_sphere + material.volume = 4.0/3.0 * np.pi * ((inner_radius + thickness)**3 - (inner_radius)**3) + inner_radius = inner_radius + thickness + outer_sphere.boundary_type = outer_boundary_type + cells.append(openmc.Cell(fill = None, region = +outer_sphere)) + geometry = openmc.Geometry(cells) + return geometry + +def make_neutron_source(energy): + point_source = openmc.stats.Point(xyz=(0.0, 0.0, 0.0)) + energy_dist = openmc.stats.Discrete(energy, 1.0) + neutron_source = openmc.Source(space = point_source, energy = energy_dist, strength = 1.0) + return neutron_source + +def make_settings(neutron_source, total_batches, inactive_batches, num_particles, run_mode): + sets = openmc.Settings() + sets.batches = total_batches + sets.inactive = inactive_batches + sets.particles = num_particles + sets.source = neutron_source + sets.run_mode = run_mode + return sets + +#Only executed if external geometry and materials are imported +def make_depletion_volumes(neutron_model, mesh_file): + materials = neutron_model.materials + mesh_file = Path(mesh_file).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=25000000) + #returns dict-like object that maps mat ids to array of volumes equal to # of mesh elements + + total_volumes = {} + for mat_id, volumes in mat_vols.items(): + total_volumes[mat_id] = np.sum(volumes) + for material in materials: + material.volume = total_volumes[material.id] + return neutron_model + +def deplete_wc(neutron_model, mesh_file, chain_file, timesteps, source_rates, norm_mode, timestep_units): + materials = neutron_model.materials + mesh_file = Path(mesh_file).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + model_nuclide_names = [] + for material in materials: + material.depletable = True + material_nuclides = material.nuclides + for material_nuclide in material_nuclides: + model_nuclide_names.append(material_nuclide.name) + + activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=7000000) + activation_mats_object = openmc.Materials(activation_mats) + activation_mats_object.export_to_xml("Activation_Materials.xml") + + depletion_chain = openmc.deplete.Chain.from_xml(Path(chain_file).resolve()) + + nuclide_names = [] + for nuclide in depletion_chain.nuclides: + if nuclide.name in model_nuclide_names: + nuclide_names.append(nuclide.name) + + fluxes, micros = openmc.deplete.get_microxs_and_flux(neutron_model, unstructured_mesh, nuclide_names, + depletion_chain.reactions, + openmc.mgxs.GROUP_STRUCTURES["VITAMIN-J-175"]) + operator = openmc.deplete.IndependentOperator(openmc.Materials(activation_mats), fluxes, micros, chain_file=chain_file, normalization_mode = norm_mode) + copy_timesteps = list(timesteps) + copy_source_rates = list(source_rates) + integrator_list = [] + for timestep in timesteps: + while copy_source_rates[-1] == 0: + integrator = openmc.deplete.PredictorIntegrator(operator, copy_timesteps, source_rates = copy_source_rates, timestep_units = timestep_units) + integrator_list.append(integrator) + copy_timesteps.pop() + copy_source_rates.pop() + + #Letting integrator_list go from the fewest number of decay intervals to the most + integrator_list.reverse() + for int_index, integrator in enumerate(integrator_list): + integrator.integrate(path=f"depletion_results_decay_set_{int_index}.h5") + return activation_mats, unstructured_mesh, neutron_model + +def make_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs): + sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] + for int_index in range(num_cooling_steps): + results = openmc.deplete.Results(f"depletion_results_decay_set_{int_index}.h5") + photon_sources = np.empty(sd_data.shape[0], dtype=object) + activated_mats_list = set(results[-1].index_mat.keys()) + + for mat_index, mat in enumerate(activation_mats) : + if str(mat.id) in activated_mats_list : + mat = results[-1].get_material(str(mat.id)) + energy = mat.get_decay_photon_energy() + if energy == None: + photon_source = openmc.IndependentSource( + energy = energy, + particle = 'photon', + strength = 0.0) + else: + photon_source = openmc.IndependentSource( + energy = energy, + particle = 'photon', + strength = energy.integral()) + + else: + photon_source = openmc.IndependentSource( + energy = openmc.stats.Discrete(0, 1.0), + particle = 'photon', + strength = 0.0) + + photon_sources[mat_index] = photon_source + photon_model = neutron_model + photon_model.settings.source = openmc.MeshSource(unstructured_mesh, photon_sources) + photon_model.settings.export_to_xml(f'settings_{int_index}.xml') + + return photon_model + +def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): + dose_energy, dose = openmc.data.dose_coefficients('photon', geometry=coeff_geom) + dose_filter = openmc.EnergyFunctionFilter(dose_energy, dose) + spherical_mesh = openmc.SphericalMesh(np.arange(0, 3800, 10), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") + spherical_mesh_filter = openmc.MeshFilter(spherical_mesh) + #Could also use an unstructured mesh filter here - spherical mesh filter for visualization purposes + particle_filter = openmc.ParticleFilter('photon') + energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-42") + + flux_tally = openmc.Tally(tally_id=1, name="photon flux tally") + flux_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter] + flux_tally.scores = ['flux'] + + dose_tally = openmc.Tally(tally_id=2, name="dose tally") + dose_tally.filters = [spherical_mesh_filter, dose_filter, energy_filter_flux, particle_filter] + dose_tally.scores = ['flux'] + + photon_model.tallies = [flux_tally, dose_tally] + #Change boundary condition: + list(photon_model.geometry.get_all_surfaces().keys())[-1].boundary_type="white" + + for int_index in range(num_cooling_steps): + #Reassign settings (which contains source) for each decay step + photon_model.settings = openmc.Settings.from_xml(f'settings_{int_index}.xml') + photon_model.export_to_model_xml(path=f'photon_model_{int_index}.xml') + sp_path = photon_model.run(f'photon_model_{int_index}.xml') + sp_path.rename(f'statepoint_openmc_{int_index}.h5') + +#---------------------------------------------------------------------------------- +#Define variables and execute all functions: + +def parse_args(): + parser = argparse.ArgumentParser() + parser.add_argument('--OpenMC_WC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs for WC_Neutron_Transport") + parser.add_argument('--ext_mat_geom', default = False, help="Specify whether materials and geometry come from external model") + parser.add_argument('--photon_transport', default = False, help="Run only photon transport (no depletion)") + args = parser.parse_args() + return args + +def read_yaml(args): + with open(args.OpenMC_WC_YAML, 'r') as transport_file: + inputs = yaml.safe_load(transport_file) + return inputs + +def create_materials_obj(inputs): + densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) + materials = make_materials(inputs['mat_info']['element_list'], + densities) + return materials + +def create_geometry_obj(materials, inputs): + geom_info = inputs['geom_info'] + layers = zip(materials, geom_info['thicknesses']) + geometry = make_spherical_shells(geom_info['inner_radius'], + layers, + geom_info['outer_boundary_type']) + return geometry + +def create_neutron_model(inputs, materials, geometry): + settings_info = inputs['settings_info'] + neutron_source = make_neutron_source(inputs['particle_energy']) + settings = make_settings(neutron_source, + settings_info['total_batches'], + settings_info['inactive_batches'], + settings_info['num_particles'], + settings_info['run_mode']) + neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = settings) + neutron_model.export_to_model_xml("neutron_model.xml") + return neutron_model + +def run_depletion(inputs, neutron_model): + #geom_info = inputs['geom_info'] + dep_params = inputs['dep_params'] + + activation_mats, unstructured_mesh, integrator, neutron_model = deplete_wc(neutron_model, + inputs['filename_dict']['mesh_file'], + dep_params['chain_file'], + dep_params['timesteps'], + dep_params['source_rates'], + dep_params['norm_mode'], + dep_params['timestep_units']) + return activation_mats, unstructured_mesh, integrator, neutron_model + +def main(): + args = parse_args() + inputs = read_yaml(args) + + openmc.config['chain_file'] = inputs['dep_params']['chain_file'] + if args.ext_mat_geom == True : #Import materials and geometry from external model + ext_model = openmc.model.Model.from_model_xml(inputs['filename_dict']['ext_model']) + materials = ext_model.materials + geometry = ext_model.geometry + else: + materials = create_materials_obj(inputs) + geometry = create_geometry_obj(materials, inputs) + + neutron_model = create_neutron_model(inputs, materials, geometry) + + #Set to True to run photon transport only (no depletion): + if args.photon_transport == True: + activation_mats = openmc.Materials.from_xml("Activation_Materials.xml") + mesh_file = Path(inputs['filename_dict']['mesh_file']).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + else: + if args.ext_mat_geom == True : + neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) + activation_mats, unstructured_mesh, neutron_model = run_depletion(inputs, neutron_model) + + num_cooling_steps = (inputs['dep_params']['source_rates']).count(0) + photon_model = make_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) + photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/WC_Layers/R2S.yaml b/WC_Layers/R2S.yaml new file mode 100644 index 0000000..9e57c67 --- /dev/null +++ b/WC_Layers/R2S.yaml @@ -0,0 +1,92 @@ +filename_dict : + elelib_fp : elelib.std + mesh_file : rbt_mesh_w_blocks_materials_no_graveyard.h5 + sp_filename : statepoint.10.h5 + figure_filename : Photon_flux_vs_energy + vtk_filename : Photon_Flux.vtk + photon_tally_figname : photon_tally + ext_model : HCPB_model_photons.xml + +mat_info : + element_list : #add in order of radially innermost to outermost + - W + - C + +geom_info : + inner_radius : 995 + thicknesses : #corresponds to each of the elements in element_list + - 5 + - 5 + outer_boundary_type : vacuum + +particle_energy : 14.0E+06 + +settings_info : + total_batches : 10 + inactive_batches : 1 + num_particles : 10000 + run_mode : fixed source + +source_meshes : + - source_mesh_hcpb_1.h5m + - source_mesh_hcpb_2.h5m +sd_filename : source_density + +source_info : + phtn_e_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 + +file_indices: + source_mesh_index : 1 + flux_spectrum_tally_id : 2 + photon_tally_id : 1 + mesh_number : 1 + energy_filter_index : 2 + +tally_info : + tallied_elements : #change according to desired tally region/material + - W + - C + +axes_without_energy_bins : (0, 1, 3, 4, 5, 6) +axes_without_mesh : (0,2,3,4,5,6) +axes_without_dose_filter_bin : (0,1,2,3,4,5) +coeff_geom : 'AP' + +#Depletion parameters for OpenMC-only R2S +dep_params : + chain_file : chain_endfb71_sfr.xml + timesteps : #sequence of timesteps following the beginning of operation (not cumulative) + - 3.0E+8 + - 86400 + - 2.6E+6 + source_rates : + - 1.0E+18 + - 0 + - 0 + norm_mode : source-rate + timestep_units : s \ No newline at end of file From f1fe7e58f2fce4af3c9a3cc304e747762680553f Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 08:30:18 -0500 Subject: [PATCH 05/27] Delete WC_Layers/OpenMC_ALARA_WC.yaml --- WC_Layers/OpenMC_ALARA_WC.yaml | 77 ---------------------------------- 1 file changed, 77 deletions(-) delete mode 100644 WC_Layers/OpenMC_ALARA_WC.yaml diff --git a/WC_Layers/OpenMC_ALARA_WC.yaml b/WC_Layers/OpenMC_ALARA_WC.yaml deleted file mode 100644 index 039c911..0000000 --- a/WC_Layers/OpenMC_ALARA_WC.yaml +++ /dev/null @@ -1,77 +0,0 @@ -filename_dict : - elelib_fp : elelib.std - mesh_file : Mesh.h5 - sp_filename : statepoint.10.h5 - figure_filename : Photon_flux_vs_energy - vtk_filename : Photon_Flux.vtk - photon_tally_figname : photon_tally - -mat_info : - element_list : #add in order of radially innermost to outermost - - W - - C - -geom_info : - inner_radius : 995 - thicknesses : #corresponds to each of the elements in element_list - - 5 - - 5 - outer_boundary_type : vacuum - -particle_energy : 14.0E+06 - -settings_info : - total_batches : 10 - inactive_batches : 1 - num_particles : 10000 - run_mode : fixed source - -source_meshes : - - source_mesh_1.h5m - - source_mesh_2.h5m -sd_filename : source_density - -source_info : - phtn_e_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 - -file_indices: - source_mesh_index : 0 - flux_spectrum_tally_id : 2 - photon_tally_id : 1 - mesh_number : 1 - energy_filter_index : 2 - -tally_info : - tallied_elements : #change according to desired tally region/material - - W - - C - -axes_without_energy_bins : (0, 1, 3, 4, 5, 6) -axes_without_mesh : (0,2,3,4,5,6) -axes_without_dose_filter_bin : (0,1,2,3,4,5) -coeff_geom : 'AP' From 3362dd8b8da2298bf5ca17a13fd3ce73f65bd830 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 21:59:17 -0500 Subject: [PATCH 06/27] Delete WC_Layers/OpenMC_R2S.py Replacing with updated version that can also perform PyNE R2S --- WC_Layers/OpenMC_R2S.py | 284 ---------------------------------------- 1 file changed, 284 deletions(-) delete mode 100644 WC_Layers/OpenMC_R2S.py diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py deleted file mode 100644 index a673edb..0000000 --- a/WC_Layers/OpenMC_R2S.py +++ /dev/null @@ -1,284 +0,0 @@ -import openmc -import openmc.deplete -from pathlib import Path -import argparse -import yaml -import numpy as np -import h5py - -def alara_element_densities(alara_fp): - ''' - Creates a dictionary where keys = element names (str) and values = element density (float) - inputs: - alara_filepath : path to ALARA element library - ''' - with open(alara_fp) 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 - -def make_materials(elements, density_dict): - ''' - Creates an OpenMC Materials object using user-specified elements - inputs: - elements: iterable of element names (str) - density_dict: dictionary with keys = element names (str) and values = element density (float) - ''' - 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 - -def make_spherical_shells(inner_radius, layers, outer_boundary_type): - ''' - Creates a set of concentric spherical shells, each with its own material & inner/outer radius. - inputs: - inner_radius: the radius of the innermost spherical shell - layers: iterable of tuples of OpenMC Material object and its respective thickness (float) - ''' - 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)) - inner_radius = outer_radius - inner_sphere = outer_sphere - material.volume = 4.0/3.0 * np.pi * ((inner_radius + thickness)**3 - (inner_radius)**3) - inner_radius = inner_radius + thickness - outer_sphere.boundary_type = outer_boundary_type - cells.append(openmc.Cell(fill = None, region = +outer_sphere)) - geometry = openmc.Geometry(cells) - return geometry - -def make_neutron_source(energy): - point_source = openmc.stats.Point(xyz=(0.0, 0.0, 0.0)) - energy_dist = openmc.stats.Discrete(energy, 1.0) - neutron_source = openmc.Source(space = point_source, energy = energy_dist, strength = 1.0) - return neutron_source - -def make_settings(neutron_source, total_batches, inactive_batches, num_particles, run_mode): - sets = openmc.Settings() - sets.batches = total_batches - sets.inactive = inactive_batches - sets.particles = num_particles - sets.source = neutron_source - sets.run_mode = run_mode - return sets - -#Only executed if external geometry and materials are imported -def make_depletion_volumes(neutron_model, mesh_file): - materials = neutron_model.materials - mesh_file = Path(mesh_file).resolve() - unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=25000000) - #returns dict-like object that maps mat ids to array of volumes equal to # of mesh elements - - total_volumes = {} - for mat_id, volumes in mat_vols.items(): - total_volumes[mat_id] = np.sum(volumes) - for material in materials: - material.volume = total_volumes[material.id] - return neutron_model - -def deplete_wc(neutron_model, mesh_file, chain_file, timesteps, source_rates, norm_mode, timestep_units): - materials = neutron_model.materials - mesh_file = Path(mesh_file).resolve() - unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - model_nuclide_names = [] - for material in materials: - material.depletable = True - material_nuclides = material.nuclides - for material_nuclide in material_nuclides: - model_nuclide_names.append(material_nuclide.name) - - activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=7000000) - activation_mats_object = openmc.Materials(activation_mats) - activation_mats_object.export_to_xml("Activation_Materials.xml") - - depletion_chain = openmc.deplete.Chain.from_xml(Path(chain_file).resolve()) - - nuclide_names = [] - for nuclide in depletion_chain.nuclides: - if nuclide.name in model_nuclide_names: - nuclide_names.append(nuclide.name) - - fluxes, micros = openmc.deplete.get_microxs_and_flux(neutron_model, unstructured_mesh, nuclide_names, - depletion_chain.reactions, - openmc.mgxs.GROUP_STRUCTURES["VITAMIN-J-175"]) - operator = openmc.deplete.IndependentOperator(openmc.Materials(activation_mats), fluxes, micros, chain_file=chain_file, normalization_mode = norm_mode) - copy_timesteps = list(timesteps) - copy_source_rates = list(source_rates) - integrator_list = [] - for timestep in timesteps: - while copy_source_rates[-1] == 0: - integrator = openmc.deplete.PredictorIntegrator(operator, copy_timesteps, source_rates = copy_source_rates, timestep_units = timestep_units) - integrator_list.append(integrator) - copy_timesteps.pop() - copy_source_rates.pop() - - #Letting integrator_list go from the fewest number of decay intervals to the most - integrator_list.reverse() - for int_index, integrator in enumerate(integrator_list): - integrator.integrate(path=f"depletion_results_decay_set_{int_index}.h5") - return activation_mats, unstructured_mesh, neutron_model - -def make_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs): - sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] - for int_index in range(num_cooling_steps): - results = openmc.deplete.Results(f"depletion_results_decay_set_{int_index}.h5") - photon_sources = np.empty(sd_data.shape[0], dtype=object) - activated_mats_list = set(results[-1].index_mat.keys()) - - for mat_index, mat in enumerate(activation_mats) : - if str(mat.id) in activated_mats_list : - mat = results[-1].get_material(str(mat.id)) - energy = mat.get_decay_photon_energy() - if energy == None: - photon_source = openmc.IndependentSource( - energy = energy, - particle = 'photon', - strength = 0.0) - else: - photon_source = openmc.IndependentSource( - energy = energy, - particle = 'photon', - strength = energy.integral()) - - else: - photon_source = openmc.IndependentSource( - energy = openmc.stats.Discrete(0, 1.0), - particle = 'photon', - strength = 0.0) - - photon_sources[mat_index] = photon_source - photon_model = neutron_model - photon_model.settings.source = openmc.MeshSource(unstructured_mesh, photon_sources) - photon_model.settings.export_to_xml(f'settings_{int_index}.xml') - - return photon_model - -def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): - dose_energy, dose = openmc.data.dose_coefficients('photon', geometry=coeff_geom) - dose_filter = openmc.EnergyFunctionFilter(dose_energy, dose) - spherical_mesh = openmc.SphericalMesh(np.arange(0, 3800, 10), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") - spherical_mesh_filter = openmc.MeshFilter(spherical_mesh) - #Could also use an unstructured mesh filter here - spherical mesh filter for visualization purposes - particle_filter = openmc.ParticleFilter('photon') - energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-42") - - flux_tally = openmc.Tally(tally_id=1, name="photon flux tally") - flux_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter] - flux_tally.scores = ['flux'] - - dose_tally = openmc.Tally(tally_id=2, name="dose tally") - dose_tally.filters = [spherical_mesh_filter, dose_filter, energy_filter_flux, particle_filter] - dose_tally.scores = ['flux'] - - photon_model.tallies = [flux_tally, dose_tally] - #Change boundary condition: - list(photon_model.geometry.get_all_surfaces().keys())[-1].boundary_type="white" - - for int_index in range(num_cooling_steps): - #Reassign settings (which contains source) for each decay step - photon_model.settings = openmc.Settings.from_xml(f'settings_{int_index}.xml') - photon_model.export_to_model_xml(path=f'photon_model_{int_index}.xml') - sp_path = photon_model.run(f'photon_model_{int_index}.xml') - sp_path.rename(f'statepoint_openmc_{int_index}.h5') - -#---------------------------------------------------------------------------------- -#Define variables and execute all functions: - -def parse_args(): - parser = argparse.ArgumentParser() - parser.add_argument('--OpenMC_WC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs for WC_Neutron_Transport") - parser.add_argument('--ext_mat_geom', default = False, help="Specify whether materials and geometry come from external model") - parser.add_argument('--photon_transport', default = False, help="Run only photon transport (no depletion)") - args = parser.parse_args() - return args - -def read_yaml(args): - with open(args.OpenMC_WC_YAML, 'r') as transport_file: - inputs = yaml.safe_load(transport_file) - return inputs - -def create_materials_obj(inputs): - densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) - materials = make_materials(inputs['mat_info']['element_list'], - densities) - return materials - -def create_geometry_obj(materials, inputs): - geom_info = inputs['geom_info'] - layers = zip(materials, geom_info['thicknesses']) - geometry = make_spherical_shells(geom_info['inner_radius'], - layers, - geom_info['outer_boundary_type']) - return geometry - -def create_neutron_model(inputs, materials, geometry): - settings_info = inputs['settings_info'] - neutron_source = make_neutron_source(inputs['particle_energy']) - settings = make_settings(neutron_source, - settings_info['total_batches'], - settings_info['inactive_batches'], - settings_info['num_particles'], - settings_info['run_mode']) - neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = settings) - neutron_model.export_to_model_xml("neutron_model.xml") - return neutron_model - -def run_depletion(inputs, neutron_model): - #geom_info = inputs['geom_info'] - dep_params = inputs['dep_params'] - - activation_mats, unstructured_mesh, integrator, neutron_model = deplete_wc(neutron_model, - inputs['filename_dict']['mesh_file'], - dep_params['chain_file'], - dep_params['timesteps'], - dep_params['source_rates'], - dep_params['norm_mode'], - dep_params['timestep_units']) - return activation_mats, unstructured_mesh, integrator, neutron_model - -def main(): - args = parse_args() - inputs = read_yaml(args) - - openmc.config['chain_file'] = inputs['dep_params']['chain_file'] - if args.ext_mat_geom == True : #Import materials and geometry from external model - ext_model = openmc.model.Model.from_model_xml(inputs['filename_dict']['ext_model']) - materials = ext_model.materials - geometry = ext_model.geometry - else: - materials = create_materials_obj(inputs) - geometry = create_geometry_obj(materials, inputs) - - neutron_model = create_neutron_model(inputs, materials, geometry) - - #Set to True to run photon transport only (no depletion): - if args.photon_transport == True: - activation_mats = openmc.Materials.from_xml("Activation_Materials.xml") - mesh_file = Path(inputs['filename_dict']['mesh_file']).resolve() - unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - else: - if args.ext_mat_geom == True : - neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) - activation_mats, unstructured_mesh, neutron_model = run_depletion(inputs, neutron_model) - - num_cooling_steps = (inputs['dep_params']['source_rates']).count(0) - photon_model = make_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) - photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) - -if __name__ == "__main__": - main() \ No newline at end of file From 7b8ab735a1e351d19d37f5505ec1ee91b04d66a0 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 21:59:57 -0500 Subject: [PATCH 07/27] Add files via upload --- OpenMC_R2S.py | 378 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 378 insertions(+) create mode 100644 OpenMC_R2S.py diff --git a/OpenMC_R2S.py b/OpenMC_R2S.py new file mode 100644 index 0000000..a94ee77 --- /dev/null +++ b/OpenMC_R2S.py @@ -0,0 +1,378 @@ +import openmc +import openmc.deplete +from pathlib import Path +import argparse +import yaml +import numpy as np +import h5py + +def alara_element_densities(alara_fp): + ''' + Creates a dictionary where keys = element names (str) and values = element density (float) + inputs: + alara_filepath : path to ALARA element library + ''' + with open(alara_fp) 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 + +def make_materials(elements, density_dict): + ''' + Creates an OpenMC Materials object using user-specified elements + inputs: + elements: iterable of element names (str) + density_dict: dictionary with keys = element names (str) and values = element density (float) + ''' + 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 + +def make_spherical_shells(inner_radius, layers, outer_boundary_type): + ''' + Creates a set of concentric spherical shells, each with its own material & inner/outer radius. + inputs: + inner_radius: the radius of the innermost spherical shell + layers: iterable of tuples of OpenMC Material object and its respective thickness (float) + ''' + 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)) + inner_radius = outer_radius + inner_sphere = outer_sphere + material.volume = 4.0/3.0 * np.pi * ((inner_radius + thickness)**3 - (inner_radius)**3) + inner_radius = inner_radius + thickness + outer_sphere.boundary_type = outer_boundary_type + cells.append(openmc.Cell(fill = None, region = +outer_sphere)) + geometry = openmc.Geometry(cells) + return geometry + +def make_neutron_source(energy): + point_source = openmc.stats.Point(xyz=(0.0, 0.0, 0.0)) + energy_dist = openmc.stats.Discrete(energy, 1.0) + neutron_source = openmc.Source(space = point_source, energy = energy_dist, strength = 1.0) + return neutron_source + +def make_neutron_tallies(mesh_file): + mesh_filter = openmc.MeshFilter(openmc.UnstructuredMesh(mesh_file, library='moab')) + particle_filter = openmc.ParticleFilter('neutron') + energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-175") + + neutron_flux_tally = openmc.Tally(name="Neutron flux spectrum") + neutron_flux_tally.filters = [mesh_filter, energy_filter_flux, particle_filter] + neutron_flux_tally.scores = ['flux'] + + return openmc.Tallies([neutron_flux_tally]) + +def make_settings(neutron_source, total_batches, inactive_batches, num_particles, run_mode): + sets = openmc.Settings() + sets.batches = total_batches + sets.inactive = inactive_batches + sets.particles = num_particles + sets.source = neutron_source + sets.run_mode = run_mode + return sets + +#Only executed if external geometry and materials are imported +def make_depletion_volumes(neutron_model, mesh_file): + materials = neutron_model.materials + mesh_file = Path(mesh_file).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=25000000) + #returns dict-like object that maps mat ids to array of volumes equal to # of mesh elements + + total_volumes = {} + for mat_id, volumes in mat_vols.items(): + total_volumes[mat_id] = np.sum(volumes) + for material in materials: + material.volume = total_volumes[material.id] + return neutron_model + +def deplete_wc(neutron_model, mesh_file, chain_file, timesteps, source_rates, norm_mode, timestep_units): + materials = neutron_model.materials + mesh_file = Path(mesh_file).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + model_nuclide_names = [] + for material in materials: + material.depletable = True + material_nuclides = material.nuclides + for material_nuclide in material_nuclides: + model_nuclide_names.append(material_nuclide.name) + + activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=7000000) + activation_mats_object = openmc.Materials(activation_mats) + activation_mats_object.export_to_xml("Activation_Materials.xml") + + depletion_chain = openmc.deplete.Chain.from_xml(Path(chain_file).resolve()) + + nuclide_names = [] + for nuclide in depletion_chain.nuclides: + if nuclide.name in model_nuclide_names: + nuclide_names.append(nuclide.name) + + fluxes, micros = openmc.deplete.get_microxs_and_flux(neutron_model, unstructured_mesh, nuclide_names, + depletion_chain.reactions, + openmc.mgxs.GROUP_STRUCTURES["VITAMIN-J-175"]) + operator = openmc.deplete.IndependentOperator(openmc.Materials(activation_mats), fluxes, micros, chain_file=chain_file, normalization_mode = norm_mode) + copy_timesteps = list(timesteps) + copy_source_rates = list(source_rates) + integrator_list = [] + for timestep in timesteps: + while copy_source_rates[-1] == 0: + integrator = openmc.deplete.PredictorIntegrator(operator, copy_timesteps, source_rates = copy_source_rates, timestep_units = timestep_units) + integrator_list.append(integrator) + copy_timesteps.pop() + copy_source_rates.pop() + + #Letting integrator_list go from the fewest number of decay intervals to the most + integrator_list.reverse() + for int_index, integrator in enumerate(integrator_list): + integrator.integrate(path=f"depletion_results_decay_set_{int_index}.h5") + return activation_mats, unstructured_mesh, neutron_model + +def extract_source_data(source_mesh_list, num_elements, num_photon_groups): + ''' + Identifies the location of the source density dataset within each mesh file. + + input: + source_mesh_list: iterable of .h5m filenames (str), whose files contain 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(source_mesh_list), num_elements, num_photon_groups)) + for source_index, source_name in enumerate(source_mesh_list): + file = h5py.File(source_name, 'r') + sd_list[source_index,:] = file['tstt']['elements']['Tet4']['tags']['source_density'][:] + return sd_list + +def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_list, neutron_model): + ''' + Creates a list of OpenMC sources, complete with the relevant space and energy distributions + + inputs: + bounds : iterable of photon energy bounds (float) + cells: list of OpenMC Cell objects + mesh_file: .h5/.h5m mesh onto which photon source will be distributed + source_mesh_indices: iterable of indices specifying the photon source from which data is extracted + + output: + all_sources: dictionary of OpenMC independent sources + unstructured_mesh: OpenMC Unstructured Mesh object + photon_model : OpenMC Model object with photon sources + ''' + all_sources = {} + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + for source_mesh_index in source_mesh_indices: + source_list = [] + for index, (lower_bound, upper_bound) in enumerate(zip(bounds[:-1],bounds[1:])): + mesh_dist = openmc.stats.MeshSpatial(unstructured_mesh, strengths=sd_list[source_mesh_index][:,index], volume_normalized=False) + energy_dist = openmc.stats.Uniform(a=lower_bound, b=upper_bound) + source_list.append(openmc.IndependentSource(space=mesh_dist, energy=energy_dist, strength=np.sum(sd_list[source_mesh_index][:, index]), particle='photon', domains=cells)) + all_sources[source_mesh_index] = source_list + photon_model = neutron_model + photon_model.settings.source = all_sources[source_mesh_index] + photon_model.settings.export_to_xml(f'settings_{source_mesh_index}.xml') + return photon_model + +def make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs): + sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] + for int_index in range(num_cooling_steps): + results = openmc.deplete.Results(f"depletion_results_decay_set_{int_index}.h5") + photon_sources = np.empty(sd_data.shape[0], dtype=object) + activated_mats_list = set(results[-1].index_mat.keys()) + + for mat_index, mat in enumerate(activation_mats) : + if str(mat.id) in activated_mats_list : + mat = results[-1].get_material(str(mat.id)) + energy = mat.get_decay_photon_energy() + if energy == None: + photon_source = openmc.IndependentSource( + energy = energy, + particle = 'photon', + strength = 0.0) + else: + photon_source = openmc.IndependentSource( + energy = energy, + particle = 'photon', + strength = energy.integral()) + + else: + photon_source = openmc.IndependentSource( + energy = openmc.stats.Discrete(0, 1.0), + particle = 'photon', + strength = 0.0) + + photon_sources[mat_index] = photon_source + photon_model = neutron_model + photon_model.settings.source = openmc.MeshSource(unstructured_mesh, photon_sources) + photon_model.settings.export_to_xml(f'settings_{int_index}.xml') + + return photon_model + +def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): + dose_energy, dose = openmc.data.dose_coefficients('photon', geometry=coeff_geom) + dose_filter = openmc.EnergyFunctionFilter(dose_energy, dose) + + #Could also use an unstructured mesh filter here - spherical mesh filter for visualization purposes + spherical_mesh = openmc.SphericalMesh(np.arange(0, 3800, 10), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") + spherical_mesh_filter = openmc.MeshFilter(spherical_mesh) + + particle_filter = openmc.ParticleFilter('photon') + energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-42") + + flux_tally = openmc.Tally(tally_id=1, name="photon flux tally") + flux_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter] + flux_tally.scores = ['flux'] + + dose_tally = openmc.Tally(tally_id=2, name="dose tally") + dose_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter, dose_filter] + dose_tally.scores = ['flux'] + + photon_model.tallies = [flux_tally, dose_tally] + #Change boundary condition: + list(photon_model.geometry.get_all_surfaces().keys())[-1].boundary_type="white" + + for int_index in range(num_cooling_steps): + #Reassign settings (which contains source) for each decay step + photon_model.settings = openmc.Settings.from_xml(f'settings_{int_index}.xml') + photon_model.export_to_model_xml(path=f'photon_model_{int_index}.xml') + sp_path = photon_model.run(f'photon_model_{int_index}.xml') + sp_path.rename(f'statepoint_openmc_{int_index}.h5') + +#---------------------------------------------------------------------------------- +#Define variables and execute all functions: + +def parse_args(): + parser = argparse.ArgumentParser() + parser.add_argument('--OpenMC_WC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs for WC_Neutron_Transport") + parser.add_argument('--ext_mat_geom', default = False, help="Specify whether materials and geometry come from external model") + parser.add_argument("--neutron_transport", default=False, help="Create neutron transport model") + parser.add_argument('--pyne_r2s', default = False, help="Specify whether PyNE R2S or OpenMC R2S steps are executed (OpenMC by default)") + parser.add_argument('--photon_transport', default = False, help="Run only photon transport (no depletion)") + args = parser.parse_args() + return args + +def read_yaml(args): + with open(args.OpenMC_WC_YAML, 'r') as transport_file: + inputs = yaml.safe_load(transport_file) + return inputs + +def create_materials_obj(inputs): + densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) + materials = make_materials(inputs['mat_info']['element_list'], + densities) + return materials + +def create_geometry_obj(materials, inputs): + geom_info = inputs['geom_info'] + layers = zip(materials, geom_info['thicknesses']) + geometry = make_spherical_shells(geom_info['inner_radius'], + layers, + geom_info['outer_boundary_type']) + return geometry + +def create_neutron_model(inputs, materials, geometry): + settings_info = inputs['settings_info'] + neutron_source = make_neutron_source(inputs['particle_energy']) + settings = make_settings(neutron_source, + settings_info['total_batches'], + settings_info['inactive_batches'], + settings_info['num_particles'], + settings_info['run_mode']) + neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = settings) + neutron_model.export_to_model_xml("neutron_model.xml") + return neutron_model + +#Convert the output of R2S Step2 to a format suitable for OpenMC photon transport: +def read_source_mesh(inputs): + #Find the size of the first source density dataset (assumed to be the same for all other datasets): + sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] + sd_list = extract_source_data(inputs['source_meshes'], + sd_data.shape[0], + sd_data.shape[1]) + return sd_list + +def create_alara_photon_model(inputs, neutron_model, sd_list): + cells = list(neutron_model.geometry.get_all_cells().values()) + photon_model = make_alara_photon_sources( + inputs['source_info']['phtn_e_bounds'], + cells, + inputs['filename_dict']['mesh_file'], + inputs['file_indices']['source_mesh_indices'], + sd_list, + neutron_model + ) + return photon_model + +def run_depletion(inputs, neutron_model): + dep_params = inputs['dep_params'] + + activation_mats, unstructured_mesh, integrator, neutron_model = deplete_wc(neutron_model, + inputs['filename_dict']['mesh_file'], + dep_params['chain_file'], + dep_params['timesteps'], + dep_params['source_rates'], + dep_params['norm_mode'], + dep_params['timestep_units']) + return activation_mats, unstructured_mesh, integrator, neutron_model + +def main(): + args = parse_args() + inputs = read_yaml(args) + + openmc.config['chain_file'] = inputs['dep_params']['chain_file'] + if args.ext_mat_geom == True : #Import materials and geometry from external model + ext_model = openmc.model.Model.from_model_xml(inputs['filename_dict']['ext_model']) + materials = ext_model.materials + geometry = ext_model.geometry + else: + materials = create_materials_obj(inputs) + geometry = create_geometry_obj(materials, inputs) + + #Settings also assigned here + neutron_model = create_neutron_model(inputs, materials, geometry) + + #Choose between PyNE R2S steps and OpenMC R2S steps + if args.pyne_r2s == True: + if args.neutron_transport == True: + neutron_model.tallies = make_neutron_tallies(inputs['filename_dict']['mesh_file']) + neutron_model.export_to_model_xml("neutron_model.xml") + else: + sd_list = read_source_mesh(inputs) + photon_model = create_alara_photon_model(inputs, neutron_model, sd_list) + num_cooling_steps = len(inputs['file_indices']['source_mesh_indices']) + + else: + # Run OpenMC-only R2S + + # Set to True to run photon transport only (no depletion): + if args.photon_transport == True: + activation_mats = openmc.Materials.from_xml("Activation_Materials.xml") + mesh_file = Path(inputs['filename_dict']['mesh_file']).resolve() + unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + else: + if args.ext_mat_geom == True : + neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) + activation_mats, unstructured_mesh, neutron_model = run_depletion(inputs, neutron_model) + + num_cooling_steps = (inputs['dep_params']['source_rates']).count(0) + photon_model = make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) + + photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) + +if __name__ == "__main__": + main() \ No newline at end of file From a283d38d127e1872faa38c6bb50ee447c60889c2 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 22:00:55 -0500 Subject: [PATCH 08/27] Rename OpenMC_R2S.py to WC_Layers/OpenMC_R2S.py --- OpenMC_R2S.py => WC_Layers/OpenMC_R2S.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename OpenMC_R2S.py => WC_Layers/OpenMC_R2S.py (99%) diff --git a/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py similarity index 99% rename from OpenMC_R2S.py rename to WC_Layers/OpenMC_R2S.py index a94ee77..d06bf38 100644 --- a/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -375,4 +375,4 @@ def main(): photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) if __name__ == "__main__": - main() \ No newline at end of file + main() From 3f9e2e4949ed7ff2a6df164a112271144681349f Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 29 May 2025 22:02:24 -0500 Subject: [PATCH 09/27] Changing source_mesh_index into a list of indices --- WC_Layers/R2S.yaml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/WC_Layers/R2S.yaml b/WC_Layers/R2S.yaml index 9e57c67..59260b1 100644 --- a/WC_Layers/R2S.yaml +++ b/WC_Layers/R2S.yaml @@ -61,7 +61,9 @@ source_info : - 2.00e+7 file_indices: - source_mesh_index : 1 + source_mesh_indices: + - 0 + - 1 flux_spectrum_tally_id : 2 photon_tally_id : 1 mesh_number : 1 @@ -89,4 +91,4 @@ dep_params : - 0 - 0 norm_mode : source-rate - timestep_units : s \ No newline at end of file + timestep_units : s From edf0b509b275a3895be0dcb91daf5a81bf1e536d Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 26 Jun 2025 14:31:52 -0500 Subject: [PATCH 10/27] Fix geometry in make_spherical_shells() --- WC_Layers/OpenMC_R2S.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index d06bf38..a59e474 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -52,10 +52,9 @@ def make_spherical_shells(inner_radius, layers, outer_boundary_type): outer_radius = inner_radius + thickness outer_sphere = openmc.Sphere(r = outer_radius) cells.append(openmc.Cell(fill = material, region = +inner_sphere & -outer_sphere)) + material.volume = 4.0/3.0 * np.pi * ((outer_radius)**3 - (inner_radius)**3) inner_radius = outer_radius inner_sphere = outer_sphere - material.volume = 4.0/3.0 * np.pi * ((inner_radius + thickness)**3 - (inner_radius)**3) - inner_radius = inner_radius + thickness outer_sphere.boundary_type = outer_boundary_type cells.append(openmc.Cell(fill = None, region = +outer_sphere)) geometry = openmc.Geometry(cells) From 876ff49d99735c35f74fcb548453feb082ddb435 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 10 Jul 2025 10:10:41 -0500 Subject: [PATCH 11/27] Remove manual void assignment to region external to geometry --- WC_Layers/OpenMC_R2S.py | 1 - 1 file changed, 1 deletion(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index a59e474..1afa3fb 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -56,7 +56,6 @@ def make_spherical_shells(inner_radius, layers, outer_boundary_type): inner_radius = outer_radius inner_sphere = outer_sphere outer_sphere.boundary_type = outer_boundary_type - cells.append(openmc.Cell(fill = None, region = +outer_sphere)) geometry = openmc.Geometry(cells) return geometry From 112524370db20373a1ca322217ce60be04caaab0 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 10 Jul 2025 11:06:24 -0500 Subject: [PATCH 12/27] Run OpenMC neutron model --- WC_Layers/OpenMC_R2S.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 1afa3fb..b719a65 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -348,7 +348,9 @@ def main(): if args.pyne_r2s == True: if args.neutron_transport == True: neutron_model.tallies = make_neutron_tallies(inputs['filename_dict']['mesh_file']) - neutron_model.export_to_model_xml("neutron_model.xml") + neutron_model.export_to_model_xml('neutron_model.xml') + neutron_model_sp = neutron_model.run('neutron_model.xml') + neutron_model_sp.rename('neutron_model.statepoint.h5') else: sd_list = read_source_mesh(inputs) photon_model = create_alara_photon_model(inputs, neutron_model, sd_list) From 26fb5f9b2f5cacd03388f76598ab3ca989f2ecbf Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 10 Jul 2025 14:32:08 -0500 Subject: [PATCH 13/27] Delete run_depletion() and rename deplete_wc to deplete_model --- WC_Layers/OpenMC_R2S.py | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index b719a65..0bb5769 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -100,7 +100,7 @@ def make_depletion_volumes(neutron_model, mesh_file): material.volume = total_volumes[material.id] return neutron_model -def deplete_wc(neutron_model, mesh_file, chain_file, timesteps, source_rates, norm_mode, timestep_units): +def deplete_model(neutron_model, mesh_file, chain_file, timesteps, source_rates, norm_mode, timestep_units): materials = neutron_model.materials mesh_file = Path(mesh_file).resolve() unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') @@ -315,22 +315,11 @@ def create_alara_photon_model(inputs, neutron_model, sd_list): neutron_model ) return photon_model - -def run_depletion(inputs, neutron_model): - dep_params = inputs['dep_params'] - - activation_mats, unstructured_mesh, integrator, neutron_model = deplete_wc(neutron_model, - inputs['filename_dict']['mesh_file'], - dep_params['chain_file'], - dep_params['timesteps'], - dep_params['source_rates'], - dep_params['norm_mode'], - dep_params['timestep_units']) - return activation_mats, unstructured_mesh, integrator, neutron_model def main(): args = parse_args() inputs = read_yaml(args) + dep_params = inputs['dep_params'] openmc.config['chain_file'] = inputs['dep_params']['chain_file'] if args.ext_mat_geom == True : #Import materials and geometry from external model @@ -367,9 +356,15 @@ def main(): else: if args.ext_mat_geom == True : neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) - activation_mats, unstructured_mesh, neutron_model = run_depletion(inputs, neutron_model) + activation_mats, unstructured_mesh, neutron_model = deplete_model(neutron_model, + inputs['filename_dict']['mesh_file'], + dep_params['chain_file'], + dep_params['timesteps'], + dep_params['source_rates'], + dep_params['norm_mode'], + dep_params['timestep_units']) - num_cooling_steps = (inputs['dep_params']['source_rates']).count(0) + num_cooling_steps = (dep_params['source_rates']).count(0) photon_model = make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) From 870d22cde232fc0dc2638118eacf94b5dd079290 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Thu, 10 Jul 2025 15:03:21 -0500 Subject: [PATCH 14/27] Add indentation for clarity --- WC_Layers/OpenMC_R2S.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 0bb5769..ee5b2b7 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -199,20 +199,20 @@ def make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_ energy = mat.get_decay_photon_energy() if energy == None: photon_source = openmc.IndependentSource( - energy = energy, - particle = 'photon', - strength = 0.0) + energy = energy, + particle = 'photon', + strength = 0.0) else: photon_source = openmc.IndependentSource( - energy = energy, - particle = 'photon', - strength = energy.integral()) + energy = energy, + particle = 'photon', + strength = energy.integral()) else: photon_source = openmc.IndependentSource( - energy = openmc.stats.Discrete(0, 1.0), - particle = 'photon', - strength = 0.0) + energy = openmc.stats.Discrete(0, 1.0), + particle = 'photon', + strength = 0.0) photon_sources[mat_index] = photon_source photon_model = neutron_model From f7758eff1618ecb41acac9e2486a99f6d3b71165 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Mon, 14 Jul 2025 11:56:09 -0500 Subject: [PATCH 15/27] Update OpenMC_R2S.py - Generalize make_settings() to take source as input instead of neutron_source - Call make_settings() twice to assign neutron settings and photon settings separately - Rearrange some .export_to_model_xml() calls --- WC_Layers/OpenMC_R2S.py | 52 ++++++++++++++++++++++++----------------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index ee5b2b7..a2c11a7 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -76,12 +76,12 @@ def make_neutron_tallies(mesh_file): return openmc.Tallies([neutron_flux_tally]) -def make_settings(neutron_source, total_batches, inactive_batches, num_particles, run_mode): +def make_settings(source, total_batches, inactive_batches, num_particles, run_mode): sets = openmc.Settings() sets.batches = total_batches sets.inactive = inactive_batches sets.particles = num_particles - sets.source = neutron_source + sets.source = source sets.run_mode = run_mode return sets @@ -157,7 +157,7 @@ def extract_source_data(source_mesh_list, num_elements, num_photon_groups): sd_list[source_index,:] = file['tstt']['elements']['Tet4']['tags']['source_density'][:] return sd_list -def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_list, neutron_model): +def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_list): ''' Creates a list of OpenMC sources, complete with the relevant space and energy distributions @@ -169,8 +169,6 @@ def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_ output: all_sources: dictionary of OpenMC independent sources - unstructured_mesh: OpenMC Unstructured Mesh object - photon_model : OpenMC Model object with photon sources ''' all_sources = {} unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') @@ -181,10 +179,7 @@ def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_ energy_dist = openmc.stats.Uniform(a=lower_bound, b=upper_bound) source_list.append(openmc.IndependentSource(space=mesh_dist, energy=energy_dist, strength=np.sum(sd_list[source_mesh_index][:, index]), particle='photon', domains=cells)) all_sources[source_mesh_index] = source_list - photon_model = neutron_model - photon_model.settings.source = all_sources[source_mesh_index] - photon_model.settings.export_to_xml(f'settings_{source_mesh_index}.xml') - return photon_model + return all_sources def make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs): sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] @@ -284,15 +279,17 @@ def create_geometry_obj(materials, inputs): return geometry def create_neutron_model(inputs, materials, geometry): - settings_info = inputs['settings_info'] + neutron_settings_info = inputs['neutron_settings_info'] neutron_source = make_neutron_source(inputs['particle_energy']) - settings = make_settings(neutron_source, - settings_info['total_batches'], - settings_info['inactive_batches'], - settings_info['num_particles'], - settings_info['run_mode']) - neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = settings) - neutron_model.export_to_model_xml("neutron_model.xml") + neutron_settings = make_settings(neutron_source, + neutron_settings_info['total_batches'], + neutron_settings_info['inactive_batches'], + neutron_settings_info['num_particles'], + neutron_settings_info['run_mode']) + neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = neutron_settings) + neutron_model.export_to_model_xml('neutron_model.xml') + neutron_model_sp = neutron_model.run('neutron_model.xml') + neutron_model_sp.rename('neutron_model.statepoint.h5') return neutron_model #Convert the output of R2S Step2 to a format suitable for OpenMC photon transport: @@ -305,17 +302,28 @@ def read_source_mesh(inputs): return sd_list def create_alara_photon_model(inputs, neutron_model, sd_list): - cells = list(neutron_model.geometry.get_all_cells().values()) - photon_model = make_alara_photon_sources( + ''' + photon_model: neutron_model to which this function adds photon-relevant settings + ''' + cells = list(photon_model.geometry.get_all_cells().values()) + all_sources = make_alara_photon_sources( inputs['source_info']['phtn_e_bounds'], cells, inputs['filename_dict']['mesh_file'], inputs['file_indices']['source_mesh_indices'], - sd_list, - neutron_model + sd_list, ) + photon_settings_info = inputs['photon_settings_info'] + for source_mesh_index, photon_source in enumerate(list(all_sources.values())): + photon_settings = make_settings(photon_source, + photon_settings_info['total_batches'], + photon_settings_info['inactive_batches'], + photon_settings_info['num_particles'], + photon_settings_info['run_mode']) + photon_model.settings = photon_settings + photon_model.export_to_model_xml('alara_photon_model_{source_mesh_index}.xml') return photon_model - + def main(): args = parse_args() inputs = read_yaml(args) From f37756b4bbff6e77ccdf8a6c253862cf5418201b Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Mon, 14 Jul 2025 12:11:37 -0500 Subject: [PATCH 16/27] Replace export_to_model_xml with settings.export_to_xml() --- WC_Layers/OpenMC_R2S.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index a2c11a7..71ecd9f 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -321,7 +321,7 @@ def create_alara_photon_model(inputs, neutron_model, sd_list): photon_settings_info['num_particles'], photon_settings_info['run_mode']) photon_model.settings = photon_settings - photon_model.export_to_model_xml('alara_photon_model_{source_mesh_index}.xml') + photon_model.settings.export_to_xml(f'settings_{source_mesh_index}.xml') return photon_model def main(): From 3ba8c4dc22663675b8b739b2588b0edeca88cb2f Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 16 Jul 2025 10:36:00 -0500 Subject: [PATCH 17/27] Remove create_materials_obj() and create_geometry_obj() --- WC_Layers/OpenMC_R2S.py | 26 +++++++++----------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 71ecd9f..65976b1 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -263,21 +263,7 @@ def read_yaml(args): with open(args.OpenMC_WC_YAML, 'r') as transport_file: inputs = yaml.safe_load(transport_file) return inputs - -def create_materials_obj(inputs): - densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) - materials = make_materials(inputs['mat_info']['element_list'], - densities) - return materials - -def create_geometry_obj(materials, inputs): - geom_info = inputs['geom_info'] - layers = zip(materials, geom_info['thicknesses']) - geometry = make_spherical_shells(geom_info['inner_radius'], - layers, - geom_info['outer_boundary_type']) - return geometry - + def create_neutron_model(inputs, materials, geometry): neutron_settings_info = inputs['neutron_settings_info'] neutron_source = make_neutron_source(inputs['particle_energy']) @@ -335,8 +321,14 @@ def main(): materials = ext_model.materials geometry = ext_model.geometry else: - materials = create_materials_obj(inputs) - geometry = create_geometry_obj(materials, inputs) + densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) + materials = make_materials(inputs['mat_info']['element_list'], + densities) + geom_info = inputs['geom_info'] + layers = zip(materials, geom_info['thicknesses']) + geometry = make_spherical_shells(geom_info['inner_radius'], + layers, + geom_info['outer_boundary_type']) #Settings also assigned here neutron_model = create_neutron_model(inputs, materials, geometry) From 28fd523e7f7b56f661f78a1824bda54b99f1be7f Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 16 Jul 2025 10:50:50 -0500 Subject: [PATCH 18/27] Combine read_source_mesh() with extract_source_data() --- WC_Layers/OpenMC_R2S.py | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 65976b1..07c74b3 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -142,20 +142,21 @@ def deplete_model(neutron_model, mesh_file, chain_file, timesteps, source_rates, integrator.integrate(path=f"depletion_results_decay_set_{int_index}.h5") return activation_mats, unstructured_mesh, neutron_model -def extract_source_data(source_mesh_list, num_elements, num_photon_groups): +def extract_source_data(inputs): ''' Identifies the location of the source density dataset within each mesh file. - - input: - source_mesh_list: iterable of .h5m filenames (str), whose files contain 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 - ''' + ''' + source_mesh_list = inputs['source_meshes'] + sd_data = h5py.File(source_mesh_list[0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] + num_elements = sd_data.shape[0] + num_photon_groups = sd_data.shape[1] sd_list = np.ndarray((len(source_mesh_list), num_elements, num_photon_groups)) for source_index, source_name in enumerate(source_mesh_list): file = h5py.File(source_name, 'r') sd_list[source_index,:] = file['tstt']['elements']['Tet4']['tags']['source_density'][:] - return sd_list + return sd_list def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_list): ''' @@ -279,13 +280,6 @@ def create_neutron_model(inputs, materials, geometry): return neutron_model #Convert the output of R2S Step2 to a format suitable for OpenMC photon transport: -def read_source_mesh(inputs): - #Find the size of the first source density dataset (assumed to be the same for all other datasets): - sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] - sd_list = extract_source_data(inputs['source_meshes'], - sd_data.shape[0], - sd_data.shape[1]) - return sd_list def create_alara_photon_model(inputs, neutron_model, sd_list): ''' @@ -341,7 +335,7 @@ def main(): neutron_model_sp = neutron_model.run('neutron_model.xml') neutron_model_sp.rename('neutron_model.statepoint.h5') else: - sd_list = read_source_mesh(inputs) + sd_list = extract_source_data(inputs) photon_model = create_alara_photon_model(inputs, neutron_model, sd_list) num_cooling_steps = len(inputs['file_indices']['source_mesh_indices']) From 58ef07506aca81d6d1ed715b74e4ab46e98fdba4 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 16 Jul 2025 11:34:01 -0500 Subject: [PATCH 19/27] Modify input to read_yaml() --- WC_Layers/OpenMC_R2S.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 07c74b3..f5a8e81 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -252,7 +252,7 @@ def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): def parse_args(): parser = argparse.ArgumentParser() - parser.add_argument('--OpenMC_WC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs for WC_Neutron_Transport") + parser.add_argument('--OpenMC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs") parser.add_argument('--ext_mat_geom', default = False, help="Specify whether materials and geometry come from external model") parser.add_argument("--neutron_transport", default=False, help="Create neutron transport model") parser.add_argument('--pyne_r2s', default = False, help="Specify whether PyNE R2S or OpenMC R2S steps are executed (OpenMC by default)") @@ -260,8 +260,12 @@ def parse_args(): args = parser.parse_args() return args -def read_yaml(args): - with open(args.OpenMC_WC_YAML, 'r') as transport_file: +def read_yaml(yaml_arg): + ''' + input: + yaml_arg : output of parse_args() corresponding to args.OpenMC_YAML + ''' + with open(yaml_arg, 'r') as transport_file: inputs = yaml.safe_load(transport_file) return inputs @@ -306,7 +310,7 @@ def create_alara_photon_model(inputs, neutron_model, sd_list): def main(): args = parse_args() - inputs = read_yaml(args) + inputs = read_yaml(args.OpenMC_YAML) dep_params = inputs['dep_params'] openmc.config['chain_file'] = inputs['dep_params']['chain_file'] From f1e36e556bb119c83b3dbcf39ce2ba18f200f871 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 16 Jul 2025 14:26:28 -0500 Subject: [PATCH 20/27] Revise execution logic for PyNE and OpenMC workflows - Added new functions import_ext_model(), make_native_model(), run_pyne_r2s(), run_openmc_r2s() to simplify logic executed in main() - Modify argparse options - Remove args.photon_transport condition for OpenMC R2S (always run all R2S steps) --- WC_Layers/OpenMC_R2S.py | 139 +++++++++++++++++++++++----------------- 1 file changed, 81 insertions(+), 58 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index f5a8e81..77e2cda 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -253,10 +253,12 @@ def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): def parse_args(): parser = argparse.ArgumentParser() parser.add_argument('--OpenMC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs") - parser.add_argument('--ext_mat_geom', default = False, help="Specify whether materials and geometry come from external model") - parser.add_argument("--neutron_transport", default=False, help="Create neutron transport model") - parser.add_argument('--pyne_r2s', default = False, help="Specify whether PyNE R2S or OpenMC R2S steps are executed (OpenMC by default)") - parser.add_argument('--photon_transport', default = False, help="Run only photon transport (no depletion)") + parser.add_argument('--ext_mat_geom', default = True, help="Specify whether materials and geometry come from external model") + parser.add_argument('--pyne_r2s', default = False, help="Choose to run pyne r2s steps") + parser.add_argument('--openmc_r2s', default = False, help="Choose to run openmc r2s steps") + if args.pyne_r2s and args.openmc_r2s: + parser.error("Cannot run PyNE and OpenMC workflows at the same time") + parser.add_argument("--pyne_neutron_transport", default=True, help="If True, run only neutron transport on PyNE model. If False, run only photon transport on PyNE model") args = parser.parse_args() return args @@ -278,9 +280,6 @@ def create_neutron_model(inputs, materials, geometry): neutron_settings_info['num_particles'], neutron_settings_info['run_mode']) neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = neutron_settings) - neutron_model.export_to_model_xml('neutron_model.xml') - neutron_model_sp = neutron_model.run('neutron_model.xml') - neutron_model_sp.rename('neutron_model.statepoint.h5') return neutron_model #Convert the output of R2S Step2 to a format suitable for OpenMC photon transport: @@ -308,64 +307,88 @@ def create_alara_photon_model(inputs, neutron_model, sd_list): photon_model.settings.export_to_xml(f'settings_{source_mesh_index}.xml') return photon_model -def main(): - args = parse_args() - inputs = read_yaml(args.OpenMC_YAML) - dep_params = inputs['dep_params'] +def import_ext_model(inputs): + ''' + Import openmc materials and geometry objects from some external openmc model + ''' + ext_model = openmc.model.Model.from_model_xml(inputs['filename_dict']['ext_model']) + materials = ext_model.materials + geometry = ext_model.geometry + + #Settings also assigned here + neutron_model = create_neutron_model(inputs, materials, geometry) - openmc.config['chain_file'] = inputs['dep_params']['chain_file'] - if args.ext_mat_geom == True : #Import materials and geometry from external model - ext_model = openmc.model.Model.from_model_xml(inputs['filename_dict']['ext_model']) - materials = ext_model.materials - geometry = ext_model.geometry - else: - densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) - materials = make_materials(inputs['mat_info']['element_list'], - densities) - geom_info = inputs['geom_info'] - layers = zip(materials, geom_info['thicknesses']) - geometry = make_spherical_shells(geom_info['inner_radius'], - layers, - geom_info['outer_boundary_type']) + return neutron_model +def make_native_model(inputs): + ''' + Run make_materials() and make_spherical_shells() + ''' + densities = alara_element_densities(inputs['filename_dict']['elelib_fp']) + materials = make_materials(inputs['mat_info']['element_list'], + densities) + geom_info = inputs['geom_info'] + layers = zip(materials, geom_info['thicknesses']) + geometry = make_spherical_shells(geom_info['inner_radius'], + layers, + geom_info['outer_boundary_type']) + #Settings also assigned here neutron_model = create_neutron_model(inputs, materials, geometry) - #Choose between PyNE R2S steps and OpenMC R2S steps - if args.pyne_r2s == True: - if args.neutron_transport == True: - neutron_model.tallies = make_neutron_tallies(inputs['filename_dict']['mesh_file']) - neutron_model.export_to_model_xml('neutron_model.xml') - neutron_model_sp = neutron_model.run('neutron_model.xml') - neutron_model_sp.rename('neutron_model.statepoint.h5') - else: - sd_list = extract_source_data(inputs) - photon_model = create_alara_photon_model(inputs, neutron_model, sd_list) - num_cooling_steps = len(inputs['file_indices']['source_mesh_indices']) - - else: - # Run OpenMC-only R2S - - # Set to True to run photon transport only (no depletion): - if args.photon_transport == True: - activation_mats = openmc.Materials.from_xml("Activation_Materials.xml") - mesh_file = Path(inputs['filename_dict']['mesh_file']).resolve() - unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - else: - if args.ext_mat_geom == True : - neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) - activation_mats, unstructured_mesh, neutron_model = deplete_model(neutron_model, - inputs['filename_dict']['mesh_file'], - dep_params['chain_file'], - dep_params['timesteps'], - dep_params['source_rates'], - dep_params['norm_mode'], - dep_params['timestep_units']) - - num_cooling_steps = (dep_params['source_rates']).count(0) - photon_model = make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) + return neutron_model + +def run_pyne_r2s(inputs, args): + if args.ext_geom == True: #Import materials and geometry from external model + materials, geometry = import_ext_model(inputs) + if args.ext_geom == False: #Run make_materials() and make_spherical_shells() + materials, geometry = make_native_model(inputs) + + #Settings also assigned here + neutron_model = create_neutron_model(inputs, materials, geometry) + + if args.pyne_neutron_transport == True: + neutron_model.tallies = make_neutron_tallies(inputs['filename_dict']['mesh_file']) + neutron_model.export_to_model_xml('neutron_model_pyne.xml') + neutron_model_sp = neutron_model.run('neutron_model_pyne.xml') + neutron_model_sp.rename('neutron_model_pyne.statepoint.h5') + + if args.pyne_neutron_transport == False: + sd_list = extract_source_data(inputs) + photon_model = create_alara_photon_model(inputs, neutron_model, sd_list) + num_cooling_steps = len(inputs['file_indices']['source_mesh_indices']) + photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) + +def run_openmc_r2s(inputs, args): + dep_params = inputs['dep_params'] + if args.ext_geom == True: #Import materials and geometry from external model + neutron_model = import_ext_model(inputs) + neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) + if args.ext_geom == False: #Run make_materials() and make_spherical_shells() + neutron_model = make_native_model(inputs) + + activation_mats, unstructured_mesh, neutron_model = deplete_model(neutron_model, + inputs['filename_dict']['mesh_file'], + dep_params['chain_file'], + dep_params['timesteps'], + dep_params['source_rates'], + dep_params['norm_mode'], + dep_params['timestep_units']) + num_cooling_steps = (dep_params['source_rates']).count(0) + photon_model = make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) + +def main(): + args = parse_args() + inputs = read_yaml(args.OpenMC_YAML) + + openmc.config['chain_file'] = inputs['dep_params']['chain_file'] + + if pyne_r2s == True: + run_pyne_r2s(inputs, args) + if openmc_r2s == True: + run_openmc_r2s(inputs, args) if __name__ == "__main__": main() From 076afa3969433b453db56549076dd28b3d68dbea Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 16 Jul 2025 14:56:49 -0500 Subject: [PATCH 21/27] Change discrete energy dist to None type --- WC_Layers/OpenMC_R2S.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 77e2cda..53b5cf8 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -206,7 +206,7 @@ def make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_ else: photon_source = openmc.IndependentSource( - energy = openmc.stats.Discrete(0, 1.0), + energy = None, particle = 'photon', strength = 0.0) From 8dceacfc8051ae1bf83504842d5611960ded4eb5 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 16 Jul 2025 20:46:27 -0500 Subject: [PATCH 22/27] Fix function calls --- WC_Layers/OpenMC_R2S.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 53b5cf8..e377672 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -340,12 +340,9 @@ def make_native_model(inputs): def run_pyne_r2s(inputs, args): if args.ext_geom == True: #Import materials and geometry from external model - materials, geometry = import_ext_model(inputs) + neutron_model = import_ext_model(inputs) if args.ext_geom == False: #Run make_materials() and make_spherical_shells() - materials, geometry = make_native_model(inputs) - - #Settings also assigned here - neutron_model = create_neutron_model(inputs, materials, geometry) + neutron_model = make_native_model(inputs) if args.pyne_neutron_transport == True: neutron_model.tallies = make_neutron_tallies(inputs['filename_dict']['mesh_file']) From 37bdc270ebbd401510f11fbf3f6b83dc78157fb3 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 16 Jul 2025 21:10:06 -0500 Subject: [PATCH 23/27] Fix argparse-related calls --- WC_Layers/OpenMC_R2S.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index e377672..7e94eb4 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -253,13 +253,13 @@ def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): def parse_args(): parser = argparse.ArgumentParser() parser.add_argument('--OpenMC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs") - parser.add_argument('--ext_mat_geom', default = True, help="Specify whether materials and geometry come from external model") + parser.add_argument('--ext_model', default = True, help="Specify whether materials and geometry come from external model") parser.add_argument('--pyne_r2s', default = False, help="Choose to run pyne r2s steps") parser.add_argument('--openmc_r2s', default = False, help="Choose to run openmc r2s steps") - if args.pyne_r2s and args.openmc_r2s: - parser.error("Cannot run PyNE and OpenMC workflows at the same time") parser.add_argument("--pyne_neutron_transport", default=True, help="If True, run only neutron transport on PyNE model. If False, run only photon transport on PyNE model") args = parser.parse_args() + if args.pyne_r2s and args.openmc_r2s: + parser.error("Cannot run PyNE and OpenMC workflows at the same time") return args def read_yaml(yaml_arg): @@ -339,9 +339,9 @@ def make_native_model(inputs): return neutron_model def run_pyne_r2s(inputs, args): - if args.ext_geom == True: #Import materials and geometry from external model + if args.ext_model == True: #Import materials and geometry from external model neutron_model = import_ext_model(inputs) - if args.ext_geom == False: #Run make_materials() and make_spherical_shells() + if args.ext_model == False: #Run make_materials() and make_spherical_shells() neutron_model = make_native_model(inputs) if args.pyne_neutron_transport == True: @@ -358,10 +358,10 @@ def run_pyne_r2s(inputs, args): def run_openmc_r2s(inputs, args): dep_params = inputs['dep_params'] - if args.ext_geom == True: #Import materials and geometry from external model + if args.ext_model == True: #Import materials and geometry from external model neutron_model = import_ext_model(inputs) neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) - if args.ext_geom == False: #Run make_materials() and make_spherical_shells() + if args.ext_model == False: #Run make_materials() and make_spherical_shells() neutron_model = make_native_model(inputs) activation_mats, unstructured_mesh, neutron_model = deplete_model(neutron_model, @@ -382,9 +382,9 @@ def main(): openmc.config['chain_file'] = inputs['dep_params']['chain_file'] - if pyne_r2s == True: + if args.pyne_r2s == True: run_pyne_r2s(inputs, args) - if openmc_r2s == True: + if args.openmc_r2s == True: run_openmc_r2s(inputs, args) if __name__ == "__main__": From cbd504307be2d8c110f27a6db1a58a98a9e3b77a Mon Sep 17 00:00:00 2001 From: anu1217 Date: Wed, 16 Jul 2025 22:21:25 -0500 Subject: [PATCH 24/27] Change how settings are assigned to openmc r2s & modify geometry to allow flux beyond initial geometry --- WC_Layers/OpenMC_R2S.py | 69 ++++++++++++++++++++++++++--------------- 1 file changed, 44 insertions(+), 25 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 7e94eb4..90f95ac 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -56,6 +56,7 @@ def make_spherical_shells(inner_radius, layers, outer_boundary_type): inner_radius = outer_radius inner_sphere = outer_sphere outer_sphere.boundary_type = outer_boundary_type + #cells.append(openmc.Cell(fill = None, region = +outer_sphere)) geometry = openmc.Geometry(cells) return geometry @@ -86,10 +87,12 @@ def make_settings(source, total_batches, inactive_batches, num_particles, run_mo return sets #Only executed if external geometry and materials are imported +#Is this even needed? def make_depletion_volumes(neutron_model, mesh_file): materials = neutron_model.materials mesh_file = Path(mesh_file).resolve() unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') + #magic number mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=25000000) #returns dict-like object that maps mat ids to array of volumes equal to # of mesh elements @@ -110,7 +113,7 @@ def deplete_model(neutron_model, mesh_file, chain_file, timesteps, source_rates, material_nuclides = material.nuclides for material_nuclide in material_nuclides: model_nuclide_names.append(material_nuclide.name) - + #magic number activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=7000000) activation_mats_object = openmc.Materials(activation_mats) activation_mats_object.export_to_xml("Activation_Materials.xml") @@ -156,7 +159,7 @@ def extract_source_data(inputs): for source_index, source_name in enumerate(source_mesh_list): file = h5py.File(source_name, 'r') sd_list[source_index,:] = file['tstt']['elements']['Tet4']['tags']['source_density'][:] - return sd_list + return sd_list def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_list): ''' @@ -184,6 +187,7 @@ def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_ def make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs): sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] + all_photon_sources = np.empty(num_cooling_steps, dtype=object) for int_index in range(num_cooling_steps): results = openmc.deplete.Results(f"depletion_results_decay_set_{int_index}.h5") photon_sources = np.empty(sd_data.shape[0], dtype=object) @@ -211,18 +215,14 @@ def make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_ strength = 0.0) photon_sources[mat_index] = photon_source - photon_model = neutron_model - photon_model.settings.source = openmc.MeshSource(unstructured_mesh, photon_sources) - photon_model.settings.export_to_xml(f'settings_{int_index}.xml') - - return photon_model + all_photon_sources[int_index] = openmc.MeshSource(unstructured_mesh, photon_sources) + return all_photon_sources def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): - dose_energy, dose = openmc.data.dose_coefficients('photon', geometry=coeff_geom) + dose_energy, dose = openmc.data.dose_coefficients('photon', geometry=coeff_geom) #[pSv cm^2] dose_filter = openmc.EnergyFunctionFilter(dose_energy, dose) - #Could also use an unstructured mesh filter here - spherical mesh filter for visualization purposes - spherical_mesh = openmc.SphericalMesh(np.arange(0, 3800, 10), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") + spherical_mesh = openmc.SphericalMesh(np.arange(0, 1505, 5), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") spherical_mesh_filter = openmc.MeshFilter(spherical_mesh) particle_filter = openmc.ParticleFilter('photon') @@ -232,17 +232,25 @@ def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): flux_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter] flux_tally.scores = ['flux'] - dose_tally = openmc.Tally(tally_id=2, name="dose tally") + dose_tally = openmc.Tally(tally_id=2, name="dose tally") #[particle-cm/src] dose_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter, dose_filter] dose_tally.scores = ['flux'] photon_model.tallies = [flux_tally, dose_tally] - #Change boundary condition: - list(photon_model.geometry.get_all_surfaces().keys())[-1].boundary_type="white" + #Change boundary condition to obtain flux/dose beyond original geometry: + list(photon_model.geometry.get_all_surfaces().values())[-1].boundary_type="transmission" + + tally_sphere = openmc.Sphere(r = 3900) + tally_sphere.boundary_type = "reflective" + tally_cell = openmc.Cell(fill = None, region = +list(photon_model.geometry.get_all_surfaces().values())[-1] & -tally_sphere) + photon_model.geometry.root_universe.add_cell(tally_cell) for int_index in range(num_cooling_steps): #Reassign settings (which contains source) for each decay step photon_model.settings = openmc.Settings.from_xml(f'settings_{int_index}.xml') + for source in photon_model.settings.source: + if "domain_ids" in source.constraints: + source.constraints['domain_ids'].append(source.constraints['domain_ids'][-1]+1) photon_model.export_to_model_xml(path=f'photon_model_{int_index}.xml') sp_path = photon_model.run(f'photon_model_{int_index}.xml') sp_path.rename(f'statepoint_openmc_{int_index}.h5') @@ -253,8 +261,8 @@ def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): def parse_args(): parser = argparse.ArgumentParser() parser.add_argument('--OpenMC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs") - parser.add_argument('--ext_model', default = True, help="Specify whether materials and geometry come from external model") - parser.add_argument('--pyne_r2s', default = False, help="Choose to run pyne r2s steps") + parser.add_argument('--ext_model', default = False, help="Specify whether materials and geometry come from external model") + parser.add_argument('--pyne_r2s', default = True, help="Choose to run pyne r2s steps") parser.add_argument('--openmc_r2s', default = False, help="Choose to run openmc r2s steps") parser.add_argument("--pyne_neutron_transport", default=True, help="If True, run only neutron transport on PyNE model. If False, run only photon transport on PyNE model") args = parser.parse_args() @@ -270,7 +278,7 @@ def read_yaml(yaml_arg): with open(yaml_arg, 'r') as transport_file: inputs = yaml.safe_load(transport_file) return inputs - + def create_neutron_model(inputs, materials, geometry): neutron_settings_info = inputs['neutron_settings_info'] neutron_source = make_neutron_source(inputs['particle_energy']) @@ -282,9 +290,7 @@ def create_neutron_model(inputs, materials, geometry): neutron_model = openmc.model.Model(geometry = geometry, materials = materials, settings = neutron_settings) return neutron_model -#Convert the output of R2S Step2 to a format suitable for OpenMC photon transport: - -def create_alara_photon_model(inputs, neutron_model, sd_list): +def create_alara_photon_model(inputs, photon_model, sd_list): ''' photon_model: neutron_model to which this function adds photon-relevant settings ''' @@ -306,7 +312,19 @@ def create_alara_photon_model(inputs, neutron_model, sd_list): photon_model.settings = photon_settings photon_model.settings.export_to_xml(f'settings_{source_mesh_index}.xml') return photon_model - + +def create_openmc_photon_model(inputs, photon_model, all_photon_sources): + photon_settings_info = inputs['photon_settings_info'] + for mesh_source_index, mesh_source in enumerate(all_photon_sources): + photon_settings = make_settings(mesh_source, + photon_settings_info['total_batches'], + photon_settings_info['inactive_batches'], + photon_settings_info['num_particles'], + photon_settings_info['run_mode']) + photon_model.settings = photon_settings + photon_model.settings.export_to_xml(f'settings_{mesh_source_index}.xml') + return photon_model + def import_ext_model(inputs): ''' Import openmc materials and geometry objects from some external openmc model @@ -342,7 +360,7 @@ def run_pyne_r2s(inputs, args): if args.ext_model == True: #Import materials and geometry from external model neutron_model = import_ext_model(inputs) if args.ext_model == False: #Run make_materials() and make_spherical_shells() - neutron_model = make_native_model(inputs) + neutron_model = make_native_model(inputs) if args.pyne_neutron_transport == True: neutron_model.tallies = make_neutron_tallies(inputs['filename_dict']['mesh_file']) @@ -373,9 +391,10 @@ def run_openmc_r2s(inputs, args): dep_params['timestep_units']) num_cooling_steps = (dep_params['source_rates']).count(0) - photon_model = make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) + all_photon_sources = make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) + photon_model = create_openmc_photon_model(inputs, neutron_model, all_photon_sources) photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) - + def main(): args = parse_args() inputs = read_yaml(args.OpenMC_YAML) @@ -386,6 +405,6 @@ def main(): run_pyne_r2s(inputs, args) if args.openmc_r2s == True: run_openmc_r2s(inputs, args) - + if __name__ == "__main__": - main() + main() \ No newline at end of file From 7269960e02dbf5d2ef9580ac2735d4f50feb7c87 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 16 Jul 2025 22:26:53 -0500 Subject: [PATCH 25/27] Remove commented lines --- WC_Layers/OpenMC_R2S.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index 90f95ac..ac4e211 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -87,12 +87,10 @@ def make_settings(source, total_batches, inactive_batches, num_particles, run_mo return sets #Only executed if external geometry and materials are imported -#Is this even needed? def make_depletion_volumes(neutron_model, mesh_file): materials = neutron_model.materials mesh_file = Path(mesh_file).resolve() unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - #magic number mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=25000000) #returns dict-like object that maps mat ids to array of volumes equal to # of mesh elements @@ -113,7 +111,6 @@ def deplete_model(neutron_model, mesh_file, chain_file, timesteps, source_rates, material_nuclides = material.nuclides for material_nuclide in material_nuclides: model_nuclide_names.append(material_nuclide.name) - #magic number activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=7000000) activation_mats_object = openmc.Materials(activation_mats) activation_mats_object.export_to_xml("Activation_Materials.xml") @@ -407,4 +404,4 @@ def main(): run_openmc_r2s(inputs, args) if __name__ == "__main__": - main() \ No newline at end of file + main() From 65863e0a7ee3adeeb3e0f3844b6aed04cdd6d77d Mon Sep 17 00:00:00 2001 From: anu1217 Date: Wed, 27 Aug 2025 17:30:33 -0500 Subject: [PATCH 26/27] Fix magic numbers & make_openmc_photon_sources(), add dummy neutron to query unstructured mesh --- WC_Layers/OpenMC_R2S.py | 135 +++++++++++++++++++++++++++------------- 1 file changed, 91 insertions(+), 44 deletions(-) diff --git a/WC_Layers/OpenMC_R2S.py b/WC_Layers/OpenMC_R2S.py index ac4e211..0081666 100644 --- a/WC_Layers/OpenMC_R2S.py +++ b/WC_Layers/OpenMC_R2S.py @@ -56,7 +56,7 @@ def make_spherical_shells(inner_radius, layers, outer_boundary_type): inner_radius = outer_radius inner_sphere = outer_sphere outer_sphere.boundary_type = outer_boundary_type - #cells.append(openmc.Cell(fill = None, region = +outer_sphere)) + cells.append(openmc.Cell(fill = None, region = +outer_sphere)) geometry = openmc.Geometry(cells) return geometry @@ -77,6 +77,17 @@ def make_neutron_tallies(mesh_file): return openmc.Tallies([neutron_flux_tally]) +def run_dummy_neutron(neutron_model, mesh_file): + ''' + Run a dummy particle to create a statepoint file from which the unstructured mesh can be queried. + ''' + neutron_model.tallies = make_neutron_tallies(mesh_file) + neutron_model.settings.particles = 1 + neutron_model.settings.batches = 1 + neutron_model.export_to_model_xml("dummy_neutron_model.xml") + neutron_model_sp = neutron_model.run("dummy_neutron_model.xml") + neutron_model_sp.rename("dummy_neutron_sp.h5") + def make_settings(source, total_batches, inactive_batches, num_particles, run_mode): sets = openmc.Settings() sets.batches = total_batches @@ -87,11 +98,11 @@ def make_settings(source, total_batches, inactive_batches, num_particles, run_mo return sets #Only executed if external geometry and materials are imported -def make_depletion_volumes(neutron_model, mesh_file): +def make_depletion_volumes(neutron_model, mesh_file, num_samples): materials = neutron_model.materials mesh_file = Path(mesh_file).resolve() unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') - mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=25000000) + mat_vols = unstructured_mesh.material_volumes(neutron_model, n_samples=num_samples) #returns dict-like object that maps mat ids to array of volumes equal to # of mesh elements total_volumes = {} @@ -101,7 +112,7 @@ def make_depletion_volumes(neutron_model, mesh_file): material.volume = total_volumes[material.id] return neutron_model -def deplete_model(neutron_model, mesh_file, chain_file, timesteps, source_rates, norm_mode, timestep_units): +def deplete_model(neutron_model, mesh_file, num_samples, chain_file, timesteps, source_rates, norm_mode, timestep_units): materials = neutron_model.materials mesh_file = Path(mesh_file).resolve() unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') @@ -111,7 +122,7 @@ def deplete_model(neutron_model, mesh_file, chain_file, timesteps, source_rates, material_nuclides = material.nuclides for material_nuclide in material_nuclides: model_nuclide_names.append(material_nuclide.name) - activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=7000000) + activation_mats = unstructured_mesh.get_homogenized_materials(neutron_model, n_samples=num_samples) activation_mats_object = openmc.Materials(activation_mats) activation_mats_object.export_to_xml("Activation_Materials.xml") @@ -174,52 +185,68 @@ def make_alara_photon_sources(bounds, cells, mesh_file, source_mesh_indices, sd_ all_sources = {} unstructured_mesh = openmc.UnstructuredMesh(mesh_file, library='moab') for source_mesh_index in source_mesh_indices: + summed_strengths_over_energy = np.sum(sd_list[source_mesh_index], axis=1) + with openmc.StatePoint("dummy_neutron_sp.h5") as sp: + vtk_mesh = list(sp.meshes.values())[-1] + vtk_mesh.write_data_to_vtk(filename=f'pyne_r2s_photon_sources_{source_mesh_index}.vtk', datasets={"pyne photon source strengths":summed_strengths_over_energy}) source_list = [] for index, (lower_bound, upper_bound) in enumerate(zip(bounds[:-1],bounds[1:])): mesh_dist = openmc.stats.MeshSpatial(unstructured_mesh, strengths=sd_list[source_mesh_index][:,index], volume_normalized=False) energy_dist = openmc.stats.Uniform(a=lower_bound, b=upper_bound) - source_list.append(openmc.IndependentSource(space=mesh_dist, energy=energy_dist, strength=np.sum(sd_list[source_mesh_index][:, index]), particle='photon', domains=cells)) + source_list.append(openmc.IndependentSource(space=mesh_dist, energy=energy_dist, strength=np.sum(sd_list[source_mesh_index][:, index]), particle='photon')) all_sources[source_mesh_index] = source_list return all_sources -def make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs): +def make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, inputs): sd_data = h5py.File(inputs['source_meshes'][0], 'r')['tstt']['elements']['Tet4']['tags']['source_density'][:] all_photon_sources = np.empty(num_cooling_steps, dtype=object) + for int_index in range(num_cooling_steps): results = openmc.deplete.Results(f"depletion_results_decay_set_{int_index}.h5") - photon_sources = np.empty(sd_data.shape[0], dtype=object) + photon_sources = np.empty(sd_data.shape[0], dtype=object) + photon_source_strengths = np.empty(sd_data.shape[0], dtype=object) activated_mats_list = set(results[-1].index_mat.keys()) - for mat_index, mat in enumerate(activation_mats) : - if str(mat.id) in activated_mats_list : - mat = results[-1].get_material(str(mat.id)) - energy = mat.get_decay_photon_energy() - if energy == None: - photon_source = openmc.IndependentSource( - energy = energy, - particle = 'photon', - strength = 0.0) - else: - photon_source = openmc.IndependentSource( - energy = energy, - particle = 'photon', - strength = energy.integral()) - - else: - photon_source = openmc.IndependentSource( - energy = None, - particle = 'photon', - strength = 0.0) - + def get_photon_energies_strengths(mat, last_results): + mat_id = str(mat.id) + if mat_id not in activated_mats_list: + return None, 0.0 + + decay_mat = last_results.get_material(mat_id) + energy = decay_mat.get_decay_photon_energy() + if energy == None: + return None, 0.0 + + return energy, energy.integral() + + for mat_index, mat in enumerate(activation_mats): + energy, strength = get_photon_energies_strengths(mat, results[-1]) + + photon_source = openmc.IndependentSource( + energy=energy, + particle='photon', + strength=strength + ) photon_sources[mat_index] = photon_source + photon_source_strengths[mat_index] = strength + all_photon_sources[int_index] = openmc.MeshSource(unstructured_mesh, photon_sources) + + with openmc.StatePoint("dummy_neutron_sp.h5") as sp: + vtk_mesh = list(sp.meshes.values())[-1] + vtk_mesh.write_data_to_vtk( + filename=f'openmc_r2s_photon_sources_{int_index}.vtk', + datasets={"source strengths": photon_source_strengths} + ) + return all_photon_sources - -def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): + +def make_photon_tallies(photon_model, coeff_geom): dose_energy, dose = openmc.data.dose_coefficients('photon', geometry=coeff_geom) #[pSv cm^2] dose_filter = openmc.EnergyFunctionFilter(dose_energy, dose) - spherical_mesh = openmc.SphericalMesh(np.arange(0, 1505, 5), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") + outermost_radius = list(photon_model.geometry.get_all_surfaces().values())[-1].r * 1.5 + spherical_mesh = openmc.SphericalMesh(np.arange(0, 5*round(outermost_radius / 5), 5), origin = (0.0, 0.0, 0.0), mesh_id=2, name="spherical_mesh") spherical_mesh_filter = openmc.MeshFilter(spherical_mesh) particle_filter = openmc.ParticleFilter('photon') @@ -233,24 +260,36 @@ def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): dose_tally.filters = [spherical_mesh_filter, energy_filter_flux, particle_filter, dose_filter] dose_tally.scores = ['flux'] - photon_model.tallies = [flux_tally, dose_tally] - #Change boundary condition to obtain flux/dose beyond original geometry: + photon_tallies = openmc.Tallies([flux_tally, dose_tally]) + + return photon_tallies + +def expand_model_geometry(photon_model): + ''' + Change boundary conditions and add new cell to obtain flux/dose beyond original geometry + ''' + outer_cell = list(photon_model.geometry.get_all_cells().values())[-1] + photon_model.geometry.root_universe.remove_cell(outer_cell) list(photon_model.geometry.get_all_surfaces().values())[-1].boundary_type="transmission" - tally_sphere = openmc.Sphere(r = 3900) + outermost_radius = list(photon_model.geometry.get_all_surfaces().values())[-1].r * 1.5 + tally_sphere = openmc.Sphere(r = 5*round(outermost_radius / 5) - 5) tally_sphere.boundary_type = "reflective" tally_cell = openmc.Cell(fill = None, region = +list(photon_model.geometry.get_all_surfaces().values())[-1] & -tally_sphere) photon_model.geometry.root_universe.add_cell(tally_cell) + return photon_model + +def run_photon_transport(photon_model, photon_tallies, num_cooling_steps): + photon_model.tallies = photon_tallies for int_index in range(num_cooling_steps): #Reassign settings (which contains source) for each decay step photon_model.settings = openmc.Settings.from_xml(f'settings_{int_index}.xml') for source in photon_model.settings.source: - if "domain_ids" in source.constraints: - source.constraints['domain_ids'].append(source.constraints['domain_ids'][-1]+1) + source.constraints['domains'] = photon_model.geometry.get_all_cells().values() photon_model.export_to_model_xml(path=f'photon_model_{int_index}.xml') sp_path = photon_model.run(f'photon_model_{int_index}.xml') - sp_path.rename(f'statepoint_openmc_{int_index}.h5') + sp_path.rename(f'statepoint_{int_index}.h5') #---------------------------------------------------------------------------------- #Define variables and execute all functions: @@ -258,7 +297,7 @@ def make_photon_tallies(coeff_geom, photon_model, num_cooling_steps): def parse_args(): parser = argparse.ArgumentParser() parser.add_argument('--OpenMC_YAML', default = "R2S.yaml", help="Path (str) to YAML containing inputs") - parser.add_argument('--ext_model', default = False, help="Specify whether materials and geometry come from external model") + parser.add_argument('--ext_model', default = True, help="Specify whether materials and geometry come from external model") parser.add_argument('--pyne_r2s', default = True, help="Choose to run pyne r2s steps") parser.add_argument('--openmc_r2s', default = False, help="Choose to run openmc r2s steps") parser.add_argument("--pyne_neutron_transport", default=True, help="If True, run only neutron transport on PyNE model. If False, run only photon transport on PyNE model") @@ -367,20 +406,25 @@ def run_pyne_r2s(inputs, args): if args.pyne_neutron_transport == False: sd_list = extract_source_data(inputs) + + run_dummy_neutron(neutron_model, inputs['filename_dict']['mesh_file']) photon_model = create_alara_photon_model(inputs, neutron_model, sd_list) num_cooling_steps = len(inputs['file_indices']['source_mesh_indices']) - photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) + photon_tallies = make_photon_tallies(photon_model, inputs['coeff_geom']) + photon_model = expand_model_geometry(photon_model) + run_photon_transport(photon_model, photon_tallies, num_cooling_steps) def run_openmc_r2s(inputs, args): dep_params = inputs['dep_params'] if args.ext_model == True: #Import materials and geometry from external model neutron_model = import_ext_model(inputs) - neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file']) + neutron_model = make_depletion_volumes(neutron_model, inputs['filename_dict']['mesh_file'], inputs['dep_params']['num_samples']) if args.ext_model == False: #Run make_materials() and make_spherical_shells() neutron_model = make_native_model(inputs) activation_mats, unstructured_mesh, neutron_model = deplete_model(neutron_model, inputs['filename_dict']['mesh_file'], + dep_params['num_samples'], dep_params['chain_file'], dep_params['timesteps'], dep_params['source_rates'], @@ -388,9 +432,12 @@ def run_openmc_r2s(inputs, args): dep_params['timestep_units']) num_cooling_steps = (dep_params['source_rates']).count(0) - all_photon_sources = make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, neutron_model, inputs) + run_dummy_neutron(neutron_model, inputs['filename_dict']['mesh_file']) + all_photon_sources = make_openmc_photon_sources(num_cooling_steps, activation_mats, unstructured_mesh, inputs) photon_model = create_openmc_photon_model(inputs, neutron_model, all_photon_sources) - photon_tallies = make_photon_tallies(inputs['coeff_geom'], photon_model, num_cooling_steps) + photon_tallies = make_photon_tallies(photon_model, inputs['coeff_geom']) + photon_model = expand_model_geometry(photon_model) + run_photon_transport(photon_model, photon_tallies, num_cooling_steps) def main(): args = parse_args() @@ -404,4 +451,4 @@ def main(): run_openmc_r2s(inputs, args) if __name__ == "__main__": - main() + main() \ No newline at end of file From eb10ff9aac4e28b5c5e65cd0c34986fd889e74d9 Mon Sep 17 00:00:00 2001 From: Anupama Rajendra <113371601+anu1217@users.noreply.github.com> Date: Wed, 27 Aug 2025 17:36:55 -0500 Subject: [PATCH 27/27] Add ray firing parameter --- WC_Layers/R2S.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/WC_Layers/R2S.yaml b/WC_Layers/R2S.yaml index 59260b1..97d6387 100644 --- a/WC_Layers/R2S.yaml +++ b/WC_Layers/R2S.yaml @@ -81,6 +81,7 @@ coeff_geom : 'AP' #Depletion parameters for OpenMC-only R2S dep_params : + num_samples : 1000000 chain_file : chain_endfb71_sfr.xml timesteps : #sequence of timesteps following the beginning of operation (not cumulative) - 3.0E+8 @@ -92,3 +93,4 @@ dep_params : - 0 norm_mode : source-rate timestep_units : s +