From e4cd463d4274a642e8cfe739b39e6ab5721a9066 Mon Sep 17 00:00:00 2001
From: Connor Moore <connor.moore@psi.ch>
Date: Mon, 29 Jun 2026 14:51:42 +0200
Subject: Added maximum RRPT radius calculation, driving script for parametric
 studies

---
 .gitignore            |  1 +
 hpmr.py               |  6 +++++-
 triso_study_driver.sh | 42 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 48 insertions(+), 1 deletion(-)
 create mode 100755 triso_study_driver.sh

diff --git a/.gitignore b/.gitignore
index 305c3a1..0608f97 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,4 @@
 *.out
 *.png
 *.pltvw
+results/*
diff --git a/hpmr.py b/hpmr.py
index 4075538..20ea672 100755
--- a/hpmr.py
+++ b/hpmr.py
@@ -6,6 +6,9 @@ import hashlib
 import sys
 import time
 
+# For LCLRS
+openmc.config["cross_sections"] = "/libs/endfb81_official/cross_sections.xml"
+
 ############ COMMAND LINE ARGUMENTS ############ 
 cli_parser = argparse.ArgumentParser(
         prog="HP-RM TRISO Homogenization Study",
@@ -477,6 +480,7 @@ print("\n    TRISO Stacking Information from Packing Calculation")
 print(f"    -> Total TRISO pebbles packed:    {geo_n_triso}")
 print(f"    -> Total TRISO volume in compact: {geo_n_triso*4/3*math.pi*geo_triso_spheres[4].r**3:.4f}")
 print(f"    -> Minimum RPT cylinder radius:   {(geo_n_triso*4/3*geo_triso_spheres[4].r**3/(2*geo_constant_height))**0.5:.4f}")
+print(f"    -> Maximum RRPT ring radius:      {(1.15**2 - geo_n_triso*4/3*geo_triso_spheres[4].r**3/(2*geo_constant_height))**0.5:.4f}")
 print("    TRISO Stacking calculation complete!\n")
 
 # Get volume fractions for mixing
@@ -660,4 +664,4 @@ geometry.export_to_xml()
 settings.export_to_xml()
 
 # Run!
-openmc.run()
+#openmc.run()
diff --git a/triso_study_driver.sh b/triso_study_driver.sh
new file mode 100755
index 0000000..2fd17a0
--- /dev/null
+++ b/triso_study_driver.sh
@@ -0,0 +1,42 @@
+#!/bin/bash
+
+mkdir -p results/ results/rpt results/rrpt
+touch results/rpt_table results/rrpt_table
+echo "\#" `date` > results/rpt_table
+echo "\#" `date` > results/rrpt_table
+
+# 1. Basic calculation without homogenization
+echo "Running explicit TRISO calculation..."
+./hpmr.py > results/explicit.out
+keff=$(cat results/explicit.out| grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "+" -f 1)
+pm=$(cat results/rpt/explicit.out | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "-" -f 2)
+echo Found $keff +/- $pm for explicit!
+
+# 2. VWH calculation
+echo "Running VWH homogenization calculation..."
+./hpmr.py > results/vwh.out
+keff=$(cat results/vwh.out| grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "+" -f 1)
+pm=$(cat results/rpt/vwh.out | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "-" -f 2)
+echo Found $keff +/- $pm for VWH!
+
+# 3. RPT
+# 0.025 step size -> 21 iterations
+for rad in $(seq 0.65 0.025 1.15); do
+    echo Running RPT with $rad cm...
+    ./hpmr.py -t rpt - $rad > results/rpt/$rad.out
+    keff=$(cat results/rpt/$rad.out | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "+" -f 1)
+    pm=$(cat results/rpt/$rad.out | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "-" -f 2)
+    echo Found $keff +/- $pm for RPT with r=$rad cm...
+    echo $rad $keff $pm >> results/rpt_table
+done
+
+# 4. RRPT
+# 0.025 step size -> 35 iterations
+for rad in $(seq 0.1 0.025 0.95); do
+    echo Running RRPT with $rad cm...
+    ./hpmr.py -t rrpt -r $rad > results/rrpt/$rad.out
+    keff=$(cat results/rrpt/$rad.out | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "+" -f 1)
+    pm=$(cat results/rrpt/$rad.out | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "-" -f 2)
+    echo Found $keff +/- $pm for RRPT with r=$rad cm...
+    echo $rad $keff $pm >> results/rrpt_table
+done
-- 
cgit v1.2.3