diff --git a/core_design/openmc_template_HPMR.py b/core_design/openmc_template_HPMR.py index fa1ea57..4964c96 100644 --- a/core_design/openmc_template_HPMR.py +++ b/core_design/openmc_template_HPMR.py @@ -52,9 +52,11 @@ def create_fuel_pin(params, materials_database): fuel_materials.append(materials_database[params['Moderator']]) # creating the fuel pin universe fuel_cells = create_cells(fuel_pin_regions, fuel_materials) + # The fuel region cell (to be used in distribcell tally) + fuel_cell = fuel_cells['fuel_meat'] fuel_pin_universe = openmc.Universe(cells=fuel_cells.values()) - return fuel_pin_universe, fuel_materials + return fuel_pin_universe, fuel_materials, fuel_cell def create_htpipe_pin(params, materials_database): @@ -456,7 +458,7 @@ def build_openmc_model_HPMR(params): # Create the fuel pin - fuel_pin_universe, fuel_materials = create_fuel_pin(params, materials_database) + fuel_pin_universe, fuel_materials, fuel_cell = create_fuel_pin(params, materials_database) # Create the heat pipe htpipe_universe, htpipe_materials = create_htpipe_pin(params, materials_database) @@ -547,6 +549,13 @@ def build_openmc_model_HPMR(params): mgxs_lib.domains = [core] mgxs_lib.build_library() mgxs_lib.add_to_tallies_file(tallies_file, merge=False) + + # Peaking factor tally (pin power) + pin_filter = openmc.DistribcellFilter(fuel_cell) + pin_power = openmc.Tally(name='pin_power_kappa') + pin_power.scores = ['kappa-fission'] + pin_power.filters = [pin_filter] + tallies_file.append(pin_power) tallies_file.export_to_xml() diff --git a/core_design/peaking_factor.py b/core_design/peaking_factor.py new file mode 100644 index 0000000..8b2c5e8 --- /dev/null +++ b/core_design/peaking_factor.py @@ -0,0 +1,103 @@ +import re +import glob +from pathlib import Path + +import openmc +import pandas as pd + + +def natural_sort_key(s: str): + """Natural sort key: n0, n1, ..., n10 instead of n0, n1, n10, n2...""" + return [int(text) if text.isdigit() else text for text in re.split(r'(\d+)', s)] + + +def compute_pin_peaking_factors(current_dir="."): + """ + Compute pin peaking factors for all OpenMC depletion statepoints + + For each file: + Computes per-pin kappa-fission power + Computes peaking factor PF = P_i / P_mean + Prints table: Rod_ID, Peaking_Factor for that depletion step + + Prints: + Per-step PF tables + Final summary: [Step, Max_PF, Rod_ID_Max] + + Returns: + summary : DataFrame with columns [Step, Max_PF, Rod_ID_Max] + per_step_data : dict[step] -> DataFrame [Rod_ID, Peaking_Factor, Step] + """ + + base = Path(current_dir) + + # Find depletion statepoint files: openmc_simulation_n*.h5 + sp_files = glob.glob(str(base / "openmc_simulation_n*.h5")) + sp_files = sorted(sp_files, key=natural_sort_key) + + if not sp_files: + print("\n[PF] No depletion statepoint files found in:", base) + print("[PF] Expected files like 'openmc_simulation_n0.h5', 'openmc_simulation_n1.h5', ...\n") + return pd.DataFrame(), {} + + tally_name = "pin_power_kappa" + results = [] + per_step_data = {} + + print("\n================ PEAKING FACTOR RESULTS ================\n") + + for sp_file in sp_files: + sp_path = Path(sp_file) + basename = sp_path.name + + # Extract raw numeric index from "openmc_simulation_nX.h5" + m = re.search(r"n(\d+)\.h5", basename) + if m: + step_raw = int(m.group(1)) # 0, 1, 2, ..., 15 + step = step_raw + 1 # 1, 2, 3, ..., 16 (shifted numbering) + else: + # Fallback: if pattern doesn't match, just keep the basename + step = basename + + # Load statepoint and tally + sp = openmc.StatePoint(str(sp_path)) + t = sp.get_tally(name=tally_name) + + # Avoid needing summary.h5 for distribcell paths + df = t.get_pandas_dataframe(paths=False) + + # Pick distribcell index column (or fallback to first) + id_col = "distribcell" if "distribcell" in df.columns else df.columns[0] + + # Per-pin power and PF + per_pin = df.groupby(id_col)["mean"].sum() + pf = per_pin / per_pin.mean() + + # Per-step PF + out = pd.DataFrame({ + "Rod_ID": per_pin.index, + "Peaking_Factor": pf.values, + "Step": step + }) + per_step_data[step] = out + + # Print per-step PF table + print(f"--- Peaking factors for depletion step {step} ---") + print(out[["Rod_ID", "Peaking_Factor"]].to_string(index=False)) + print() + + # Collect for summary + results.append({ + "Step": step, + "Max_PF": float(pf.max()), + "Rod_ID_Max": pf.idxmax() + }) + + # Build and print summary + summary = pd.DataFrame(results).sort_values("Step") + + print("========== Peaking Factor Summary ==========") + print(summary.to_string(index=False)) + print("============================================\n") + + return summary, per_step_data diff --git a/core_design/utils.py b/core_design/utils.py index 6757c7c..89e224c 100644 --- a/core_design/utils.py +++ b/core_design/utils.py @@ -7,6 +7,7 @@ import matplotlib.pyplot as plt import matplotlib.patches as mpatches from core_design.correction_factor import corrected_keff_2d +from core_design.peaking_factor import compute_pin_peaking_factors @@ -183,21 +184,32 @@ def openmc_depletion(params, lattice_geometry, settings): depletion_2d_results_file = openmc.deplete.Results("./depletion_results.h5") # Example file path fuel_lifetime_days = corrected_keff_2d(depletion_2d_results_file, params['Active Height'] + 2 * params['Axial Reflector Thickness']) + + try: + pf_summary, pf_per_step = compute_pin_peaking_factors(".") + 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') - return fuel_lifetime_days, mass_U235, mass_U238 + return fuel_lifetime_days, mass_U235, mass_U238, pf_summary def run_depletion_analysis(params): openmc.run() lattice_geometry = openmc.Geometry.from_xml() settings = openmc.Settings.from_xml() - depletion_results = openmc_depletion(params, lattice_geometry, settings) - params['Fuel Lifetime'] = depletion_results[0] # days - params['Mass U235'] = depletion_results[1] # grams - params['Mass U238'] = depletion_results[2] # grams - params['Uranium Mass'] = (params['Mass U235'] + params['Mass U238']) / 1000 # Kg + fuel_lifetime_days, mass_U235, mass_U238, pf_summary = \ + openmc_depletion(params, lattice_geometry, settings) + + params['Fuel Lifetime'] = fuel_lifetime_days # days + params['Mass U235'] = mass_U235 # grams + params['Mass U238'] = mass_U238 # grams + params['Uranium Mass'] = (mass_U235 + mass_U238) / 1000 # kg + params['PF Summary'] = pf_summary.to_dict(orient="list") def monitor_heat_flux(params):