Add Pin Peaking Factor Calculations for HPMR

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()
 
 

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

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):