diff --git a/core_design/utils.py b/core_design/utils.py index 2b96e27..2364172 100644 --- a/core_design/utils.py +++ b/core_design/utils.py @@ -11,6 +11,27 @@ import pandas,copy +def _get_mpi_rank(): + try: + from mpi4py import MPI + return MPI.COMM_WORLD.Get_rank() + except ImportError: + return 0 + +def _mpi_barrier(): + try: + from mpi4py import MPI + MPI.COMM_WORLD.Barrier() + except ImportError: + pass + +def _mpi_bcast(obj, root=0): + try: + from mpi4py import MPI + return MPI.COMM_WORLD.bcast(obj, root=root) + except ImportError: + return obj + def circle_area(r): return (np.pi) * r **2 @@ -206,7 +227,6 @@ def create_universe_plot(materials_database, universe, plot_width, num_pixels, f def openmc_depletion(params, lattice_geometry, settings): openmc.config['cross_sections'] = params['cross_sections_xml_location'] - # depletion operator, performing transport simulations, is created using the geometry and settings xml files operator = openmc.deplete.CoupledOperator(openmc.Model(geometry=lattice_geometry, settings=settings), @@ -230,31 +250,47 @@ def openmc_depletion(params, lattice_geometry, settings): integrator.integrate() print("End Depletion") - depletion_2d_results_file = openmc.deplete.Results("./depletion_results.h5") # Example file path + # Only rank 0 handles post-processing, then broadcasts results + rank = _get_mpi_rank() - fuel_lifetime_days, keff_2d_values, keff_2d_values_corrected = corrected_keff_2d(depletion_2d_results_file, params['Active Height'] + 2 * params['Axial Reflector Thickness']) + if rank == 0: + depletion_2d_results_file = openmc.deplete.Results("./depletion_results.h5") + fuel_lifetime_days, keff_2d_values, keff_2d_values_corrected = corrected_keff_2d( + depletion_2d_results_file, + params['Active Height'] + 2 * params['Axial Reflector Thickness'] + ) - # Compute pin peaking factors - try: - pf_summary, pf_per_step = compute_pin_peaking_factors(".") - # Store peaking factor results in params for Excel output - idx_max = pf_summary['Max_PF'].idxmax() - params['Max Peaking Factor'] = pf_summary.loc[idx_max, 'Max_PF'] - params['Step with Max Peaking Factor'] = pf_summary.loc[idx_max, 'Step'] - params['Rod ID with Max Peaking Factor'] = pf_summary.loc[idx_max, 'Rod_ID_Max'] - params['Max Peaking Factors per Step'] = pf_summary['Max_PF'].tolist() - params['PF Summary'] = pf_summary.to_dict(orient='list') + try: + pf_summary, pf_per_step = compute_pin_peaking_factors(".") + idx_max = pf_summary['Max_PF'].idxmax() + params['Max Peaking Factor'] = pf_summary.loc[idx_max, 'Max_PF'] + params['Step with Max Peaking Factor'] = pf_summary.loc[idx_max, 'Step'] + params['Rod ID with Max Peaking Factor'] = pf_summary.loc[idx_max, 'Rod_ID_Max'] + params['Max Peaking Factors per Step'] = pf_summary['Max_PF'].tolist() + params['PF Summary'] = pf_summary.to_dict(orient='list') + except Exception as e: + print("[PF] WARNING: compute_pin_peaking_factors failed:", e) + pf_summary = None + pf_per_step = None + orig_material = depletion_2d_results_file.export_to_materials(0) + mass_U235 = orig_material[0].get_mass('U235') + mass_U238 = orig_material[0].get_mass('U238') + data_k = pandas.DataFrame() + data_k['keff 2D'] = keff_2d_values + data_k['keff 3D (2D corrected)'] = keff_2d_values_corrected + data_k.to_csv('./objectives_keff.csv', index_label='time') - except Exception as e: - print("[PF] WARNING: compute_pin_peaking_factors failed:", e) - pf_summary = None - pf_per_step = None + results = (fuel_lifetime_days, mass_U235, mass_U238, pf_summary, + keff_2d_values, keff_2d_values_corrected) + else: + results = None - orig_material = depletion_2d_results_file.export_to_materials(0) - mass_U235 = orig_material[0].get_mass('U235') - mass_U238 = orig_material[0].get_mass('U238') + # Broadcast results to all ranks + _mpi_barrier() + results = _mpi_bcast(results, root=0) + fuel_lifetime_days, mass_U235, mass_U238, pf_summary, keff_2d_values, keff_2d_values_corrected = results params['keff 2D'] = keff_2d_values params['keff 3D (2D corrected)'] = keff_2d_values_corrected @@ -263,7 +299,7 @@ def openmc_depletion(params, lattice_geometry, settings): def run_depletion_analysis(params): - openmc.run() + # openmc.run() lattice_geometry = openmc.Geometry.from_xml() settings = openmc.Settings.from_xml() fuel_lifetime_days, mass_U235, mass_U238, pf_summary = \ diff --git a/examples/watts_exec_LTMR_mpi.py b/examples/watts_exec_LTMR_mpi.py new file mode 100644 index 0000000..d342e54 --- /dev/null +++ b/examples/watts_exec_LTMR_mpi.py @@ -0,0 +1,344 @@ +# Copyright 2025, Battelle Energy Alliance, LLC, ALL RIGHTS RESERVED + +""" +This script performs a bottom-up cost estimate for a Liquid Metal Thermal Microreactor (LTMR). +OpenMC is used for core design calculations, and other Balance of Plant components are estimated. +Users can modify parameters in the "params" dictionary below. +""" +from mpi4py import MPI +import numpy as np +import watts # Simulation workflows for one or multiple codes +from core_design.openmc_template_LTMR import * +from core_design.pins_arrangement import LTMR_pins_arrangement +from core_design.utils import * +from core_design.drums import * +from reactor_engineering_evaluation.fuel_calcs import fuel_calculations +from reactor_engineering_evaluation.BOP import * +from reactor_engineering_evaluation.vessels_calcs import * +from reactor_engineering_evaluation.tools import * +from cost.cost_estimation import detailed_bottom_up_cost_estimate + +import warnings +warnings.filterwarnings("ignore") + +import time +time_start = time.time() + +params = watts.Parameters() + +def update_params(updates): + params.update(updates) +def _get_mpi_rank(): + try: + from mpi4py import MPI + return MPI.COMM_WORLD.Get_rank() + except ImportError: + return 0 + +def _mpi_barrier(): + try: + from mpi4py import MPI + MPI.COMM_WORLD.Barrier() + except ImportError: + pass + +rank = _get_mpi_rank() +# ************************************************************************************************************************** +# Sec. 0: Settings +# ************************************************************************************************************************** + +update_params({ + 'plotting': "N", # "Y" or "N": Yes or No + 'cross_sections_xml_location': '/hpc-common/data/openmc/endfb-viii.0-hdf5/cross_sections.xml', # on INL HPC + 'simplified_chain_thermal_xml': '/home/garcsamu/OpenMC/TEMA/data/chain_casl_pwr.xml' # on INL HPC +}) + +# ************************************************************************************************************************** +# Sec. 1: Materials +# ************************************************************************************************************************** + +update_params({ + 'reactor type': "LTMR", # LTMR or GCMR + 'TRISO Fueled': "No", + 'Fuel': 'TRIGA_fuel', + 'Enrichment': 0.1975, # Fraction between 0 and 1 + "H_Zr_ratio": 1.6, # Proportion of hydrogen to zirconium atoms + 'U_met_wo': 0.3, # Weight ratio of Uranium to total fuel weight (less than 1) + 'Coolant': 'NaK', + 'Radial Reflector': 'Graphite', + 'Axial Reflector': 'Graphite', + 'Moderator': 'ZrH', + 'Control Drum Absorber': 'B4C_enriched', + 'Control Drum Reflector': 'Graphite', + 'Common Temperature': 600, # Kelvins + 'HX Material': 'SS316' +}) + +# ************************************************************************************************************************** +# Sec. 2: Geometry: Fuel Pins, Moderator Pins, Coolant, Hexagonal Lattice +# ************************************************************************************************************************** + +update_params({ + 'Fuel Pin Materials': ['Zr', None, params['Fuel'], None, 'SS304'], + 'Fuel Pin Radii': [0.28575, 0.3175, 1.5113, 1.5367, 1.5875], # cm + 'Moderator Pin Materials': ['ZrH', 'SS304'], + 'Moderator Pin Inner Radius': 1.5367, # cm + 'Moderator Pin Radii': [1.5367, 1.5875], # [params['Moderator Pin Inner Radius'], params['Fuel Pin Radii'][-1]] + "Pin Gap Distance": 0.1, # cm + 'Pins Arrangement': LTMR_pins_arrangement, + 'Number of Rings per Assembly': 12, # the number of rings can be 12 or lower as long as the heat flux criteria is not violated + 'Radial Reflector Thickness': 14, # cm +}) + +params['Lattice Radius'] = calculate_lattice_radius(params) +params['Active Height'] = 78.4 # Or it is 2 * params['Lattice Radius'] +params['Axial Reflector Thickness'] = params['Radial Reflector Thickness'] # cm +params['Fuel Pin Count'] = calculate_pins_in_assembly(params, "FUEL") +params['Moderator Pin Count'] = calculate_pins_in_assembly(params, "MODERATOR") +params['Moderator Mass'] = calculate_moderator_mass(params) +params['Core Radius'] = params['Lattice Radius'] + params['Radial Reflector Thickness'] + +# ************************************************************************************************************************** +# Sec. 3: Control Drums +# ************************************************************************************************************************** + +update_params({ + 'Drum Radius': 9.016, # or it is 0.23 * params['Lattice Radius'], # cm + 'Drum Absorber Thickness': 1, # cm + 'Drum Height': params['Active Height'] + 2*params['Axial Reflector Thickness'] +}) + +calculate_drums_volumes_and_masses(params) +calculate_reflector_mass_LTMR(params) + +# ************************************************************************************************************************** +# Sec. 4: Overall System +# ************************************************************************************************************************** + +update_params({ + 'Power MWt': 20, # MWt + 'Thermal Efficiency': 0.31, + 'Heat Flux Criteria': 0.9, # MW/m^2 + 'Burnup Steps': [0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0, 120.0, 140.0] # MWd_per_Kg +}) +params['Power MWe'] = params['Power MWt'] * params['Thermal Efficiency'] +params['Heat Flux'] = calculate_heat_flux(params) + +# ************************************************************************************************************************** +# Sec. 5: Running OpenMC +# ************************************************************************************************************************** + +# --- Shutdown Margin (SDM) --- +# When True, an additional OpenMC simulation is run with all control drums rotated +# to the fully inserted (ARI - All Rods In) position. The SDM is then calculated +# as the difference in reactivity (in pcm) between the ARO and ARI configurations. +# A positive SDM means the reactor can be safely shut down with all drums inserted. +# Recommended: True for final design verification; can be set to False to save +# computation time during early design exploration. +params['SD Margin Calc'] = False # True or False + +# --- Isothermal Temperature Coefficient --- +# When True, two additional OpenMC simulations are run: one at 'Common Temperature' +# and one at 'Common Temperature' + 'Temperature Perturbation'. The temperature +# coefficient is then calculated in units of pcm/K. +# A negative coefficient indicates the reactor is self-stabilizing (desired behavior). +# Recommended: True for safety analysis; can be set to False to save computation time. +params['Isothermal Temperature Coefficients'] = False # True or False + +# --- Temperature Perturbation --- +# The temperature step (in Kelvin) used for the isothermal temperature coefficient calculation. +# Must be large enough to produce a keff difference above OpenMC Monte Carlo statistical +# noise, but small enough to stay in the linear reactivity regime. +# Typical range: 50–300 K. 100 K is chosen here as a balance between accuracy and +# avoiding nonlinear effects. +# Units: Kelvin +# This parameter is REQUIRED only when 'Isothermal Temperature Coefficients' is True. +params['Temperature Perturbation'] = 100 # K +heat_flux_monitor = monitor_heat_flux(params) +run_openmc(build_openmc_model_LTMR, heat_flux_monitor, params) +fuel_calculations(params) # calculate the fuel mass and SWU + +if rank ==0: + # ************************************************************************************************************************** + # Sec. 6: Primary Loop + Balance of Plant + # ************************************************************************************************************************** + + update_params({ + 'Secondary HX Mass': 0, + 'Primary Pump': 'Yes', + 'Secondary Pump': 'No', + 'Pump Isentropic Efficiency': 0.8, + 'Primary Loop Inlet Temperature': 430 + 273.15, # K + 'Primary Loop Outlet Temperature': 520 + 273.15, # K + 'Secondary Loop Inlet Temperature': 395 + 273.15, # K + 'Secondary Loop Outlet Temperature': 495 + 273.15, # K, + }) + + params['Primary HX Mass'] = calculate_heat_exchanger_mass(params) # Kg + # Update BoP Parameters + params.update({ + 'BoP Count': 2, # Number of BoP present in plant + 'BoP per loop load fraction': 0.5, # based on assuming that each BoP Handles the total load evenly (1/2) + }) + params['BoP Power kWe'] = 1000 * params['Power MWe'] * params['BoP per loop load fraction'] + # calculate coolant mass flow rate + mass_flow_rate(params) + calculate_primary_pump_mechanical_power(params) + + # ************************************************************************************************************************** + # Sec. 7: Shielding + # ************************************************************************************************************************** + + update_params({ + 'In Vessel Shield Thickness': 10.16, # cm + 'In Vessel Shield Inner Radius': params['Core Radius'], + 'In Vessel Shield Material': 'B4C_natural', + 'Out Of Vessel Shield Thickness': 39.37, # cm + 'Out Of Vessel Shield Material': 'WEP', + 'Out Of Vessel Shield Effective Density Factor': 0.5 # The out of vessel shield is not fully made of the out of vessel material (e.g. WEP) so we use an effective density factor + }) + + params['In Vessel Shield Outer Radius'] = params['Core Radius'] + params['In Vessel Shield Thickness'] + + # ************************************************************************************************************************** + # Sec. 8: Vessels Calculations + # ************************************************************************************************************************** + + update_params({ + 'Vessel Radius': params['Core Radius'] + params['In Vessel Shield Thickness'], + 'Vessel Thickness': 1, # cm + 'Vessel Lower Plenum Height': 42.848 - 40, # cm, based on Reflecting Barrel~RPV Liner (-Reflector Thickness, which is currently missing in CAD), # cm + 'Vessel Upper Plenum Height': 47.152, # cm + 'Vessel Upper Gas Gap': 0, + 'Vessel Bottom Depth': 32.129, + 'Vessel Material': 'stainless_steel', + 'Gap Between Vessel And Guard Vessel': 2, # cm + 'Guard Vessel Thickness': 0.5, # cm + 'Guard Vessel Material': 'stainless_steel', + 'Gap Between Guard Vessel And Cooling Vessel': 5, # cm + 'Cooling Vessel Thickness': 0.5, # cm + 'Cooling Vessel Material': 'stainless_steel', + 'Gap Between Cooling Vessel And Intake Vessel': 3, # cm + 'Intake Vessel Thickness': 0.5, # cm + 'Intake Vessel Material': 'stainless_steel' + }) + + vessels_specs(params) # calculate the volumes and masses of the vessels + calculate_shielding_masses(params) # calculate the masses of the shieldings + + # ************************************************************************************************************************** + # Sec. 9: Operation + # ************************************************************************************************************************** + + update_params({ + 'Operation Mode': "Autonomous", + 'Number of Operators': 2, + 'Levelization Period': 60, # years + 'Refueling Period': 7, + 'Emergency Shutdowns Per Year': 0.2, + 'Startup Duration after Refueling': 2, + 'Startup Duration after Emergency Shutdown': 14, + 'Reactors Monitored Per Operator': 10, + 'Security Staff Per Shift': 1 + }) + ## Calculated based on 1 tanks + ## Density of NaK=855 kg/m3, Volume=8.2402 m3 (standard tank size) + params['Onsite Coolant Inventory'] = 1 * 855 * 8.2402 # kg + params['Replacement Coolant Inventory'] = 0 # assume that NaK does not need to be replaced. + # params['Annual Coolant Supply Frequency'] # LTMR should not require frequent refilling + + total_refueling_period = params['Fuel Lifetime'] + params['Refueling Period'] + params['Startup Duration after Refueling'] # days + total_refueling_period_yr = total_refueling_period/365 + params['A75: Vessel Replacement Period (cycles)'] = np.floor(10/total_refueling_period_yr) # change each 10 years similar to the ATR + params['A75: Core Barrel Replacement Period (cycles)'] = np.floor(10/total_refueling_period_yr) + params['A75: Reflector Replacement Period (cycles)'] = np.floor(10/total_refueling_period_yr) + params['A75: Drum Replacement Period (cycles)'] = np.floor(10/total_refueling_period_yr) + params['Mainenance to Direct Cost Ratio'] = 0.015 + # A78: Annualized Decommisioning Cost + params['A78: CAPEX to Decommissioning Cost Ratio'] = 0.15 + + # ************************************************************************************************************************** + # Sec. 10: Buildings & Economic Parameters + # ************************************************************************************************************************** + + update_params({ + 'Land Area': 18, # acres + 'Escalation Year': 2024, + + 'Excavation Volume': 412.605, # m^3 + 'Reactor Building Slab Roof Volume': (9750*6502.4*1500)/1e9, # m^3 + 'Reactor Building Basement Volume': (9750*6502.4*1500)/1e9, # m^3 + 'Reactor Building Exterior Walls Volume': ((2*9750*3500*1500)+(3502.4*3500*(1500+750)))/1e9, # m^3 + 'Reactor Building Superstructure Area': ((2*3500*3500)+(2*7500*3500))/1e6, # m^2 + + # Connected to the Reactor Building (contains steel liner) + 'Integrated Heat Exchanger Building Slab Roof Volume': 0, # m^3 + 'Integrated Heat Exchanger Building Basement Volume': 0, # m^3 + 'Integrated Heat Exchanger Building Exterior Walls Volume': 0, # m^3 + 'Integrated Heat Exchanger Building Superstructure Area': 0, # m^2 + + # Assumed to be High 40' CONEX Container with 20 cm wall thickness (including conex wall) + 'Turbine Building Slab Roof Volume': (12192*2438*200)/1e9, # m^3 + 'Turbine Building Basement Volume': (12192*2438*200)/1e9, # m^3 + 'Turbine Building Exterior Walls Volume': ((12192*2496*200)+(2038*2496*200))*2/1e9, # m^3 + + # Assumed to be High 40' CONEX Container with 20 cm wall thickness (including conex wall) + 'Control Building Slab Roof Volume': (12192*2438*200)/1e9, # m^3 + 'Control Building Basement Volume': (12192*2438*200)/1e9, # m^3 + 'Control Building Exterior Walls Volume': ((12192*2496*200)+(2038*2496*200))*2/1e9, # m^3 + + # Manipulator Building + 'Manipulator Building Slab Roof Volume': (4876.8*2438.4*400)/1e9, # m^3 + 'Manipulator Building Basement Volume': (4876.8*2438.4*1500)/1e9, # m^3 + 'Manipulator Building Exterior Walls Volume': ((4876.8*4445*400)+(2038.4*4445*400*2))/1e9, # m^3 + + 'Refueling Building Slab Roof Volume': 0, # m^3 + 'Refueling Building Basement Volume': 0, # m^3 + 'Refueling Building Exterior Walls Volume': 0, # m^3 + + 'Spent Fuel Building Slab Roof Volume': 0, # m^3 + 'Spent Fuel Building Basement Volume': 0, # m^3 + 'Spent Fuel Building Exterior Walls Volume': 0, # m^3 + + 'Emergency Building Slab Roof Volume': 0, # m^3 + 'Emergency Building Basement Volume': 0, # m^3 + 'Emergency Building Exterior Walls Volume': 0, # m^3 + + # Building to host operational spares (CO2, He, filters, etc.) + 'Storage Building Slab Roof Volume': (8400*3500*400)/1e9, # m^3 + 'Storage Building Basement Volume': (8400*3500*400)/1e9, # m^3 + 'Storage Building Exterior Walls Volume': ((8400*2700*400)+(3100*2700*400*2))/1e9, # m^3 + + 'Radwaste Building Slab Roof Volume': 0, # m^3 + 'Radwaste Building Basement Volume': 0, # m^3 + 'Radwaste Building Exterior Walls Volume': 0, # m^3, + + 'Interest Rate': 0.07, + 'Construction Duration': 12, # months + 'Debt To Equity Ratio': 1, + 'Annual Return': 0.0475, # Annual return on decommissioning costs + 'NOAK Unit Number': 100 + }) + + # --- ITC (Investment Tax Credit) --- + # The ITC is a one-time credit applied to the Overnight Capital Cost (OCC) of the plant. + # Under the IRA (Inflation Reduction Act), advanced nuclear facilities placed in service + # after Dec 31, 2024 may qualify for the Clean Electricity ITC (Section 48E). + # The ITC level depends on whether the project meets certain requirements: + # - Base rate (no prevailing wage): 6% of OCC + # - With prevailing wage + apprenticeship requirements: 30% of OCC + # - With prevailing wage + domestic content bonus: 40% of OCC + # - With prevailing wage + domestic content + energy community bonus: 50% of OCC + # Typical values: 0.06, 0.30, 0.40, 0.50 + # Note: ITC and PTC are mutually exclusive — only one can be selected per project. + # To disable ITC, remove or comment out this parameter. + params['ITC credit level'] = 0.30 # fraction — assumes prevailing wage requirements are met + + # ************************************************************************************************************************** + # Sec. 11: Post Processing + # ************************************************************************************************************************** + params['Number of Samples'] = 100 # Accounting for cost uncertainties + # Estimate costs using the cost database file and save the output to an Excel file + estimate = detailed_bottom_up_cost_estimate('/home/garcsamu/OpenMC/MOUSE/cost/Cost_Database.xlsx', params, "/home/garcsamu/OpenMC/MOUSE/examples/output_LTMR_parallel.xlsx") + elapsed_time = (time.time() - time_start) / 60 # Calculate execution time + print('Execution time:', np.round(elapsed_time, 1), 'minutes') \ No newline at end of file diff --git a/myjob_mpi.sh b/myjob_mpi.sh new file mode 100644 index 0000000..00576f0 --- /dev/null +++ b/myjob_mpi.sh @@ -0,0 +1,16 @@ +#!/bin/bash +#SBATCH --job-name=mpi_example +#SBATCH --partition=short +#SBATCH --time=5:58:00 +#SBATCH --nodes=1 +#SBATCH --ntasks=2 +#SBATCH --cpus-per-task=25 +#SBATCH --wckey=edu_res +#SBATCH --error=serial.err.%J +#SBATCH --output=serial.txt.%J + +source ~/.bashrc +conda activate snap + +cd examples/ +mpirun -np 2 python -u watts_exec_LTMR_mpi.py \ No newline at end of file