New logic for homogenizations and added RPT support

1 files changed, 62 insertions, 42 deletions

diff --git a/hpmr.py b/hpmr.py
index 4d1e2f8..a11b066 100755
--- a/hpmr.py
+++ b/hpmr.py
@@ -1,37 +1,35 @@
 #!/bin/python3
 import openmc
 import math
-import argparse as ap
-import hashlib as hl
+import argparse
+import hashlib
 import sys
 import time
 
-def id_file():
-    path = sys.argv[0]
-
-    with open(path, "rb") as file:
-        contents = file.read()
-
-    filehash = hl.md5(contents).hexdigest()
-
-    print("="*60)
-    print(f"Run date:       {time.ctime()}")
-    print(f"Input file MD5: {filehash}")
-    print("="*60)
-
-
 ############ COMMAND LINE ARGUMENTS ############ 
-id_file()
-
-cli_parser = ap.ArgumentParser(
+cli_parser = argparse.ArgumentParser(
         prog="HP-RM TRISO Homogenization Study",
         description="Program to investigate various homogenization techniques for TRISO compacts in the HP-MR benchmark."
         )
 
 cli_parser.add_argument("-t","--technique",default="none",choices=["none","vwh","rpt","rrpt"],help="What homogenization technique to use (if any) [none/vwh/rpt/rrpt], default none")
-cli_parser.add_argument("-r","--radius",default="0.5",type=float,help="Radius for RPT or RRPT homogenization [cm], default 0.5 cm. In the case of RRPT, this is the inner radius.")
+cli_parser.add_argument("-r","--radius",default="0.8",type=float,help="Radius for RPT or RRPT homogenization [cm], default 0.5 cm. In the case of RRPT, this is the inner radius.")
 cli_args = cli_parser.parse_args()
 
+############ INTPUT-OUTPUT MAPPING ############ 
+path = sys.argv[0]
+
+with open(path, "rb") as file:
+    contents = file.read()
+
+filehash = hashlib.md5(contents).hexdigest()
+
+print("="*60)
+print(f"Run date CET:     {time.ctime()}")
+print(f"Input file MD5:   {filehash}")
+print(f"Input arguments:  {cli_args}")
+print("="*60)
+
 ############ MATERIALS ############ 
 # All materials are specified in the report on HEAT PIPE MICROREACTOR MODELING WITH BLUECRAB
 # <https://www.tandfonline.com/doi/full/10.1080/00295639.2024.2375175>
@@ -197,30 +195,30 @@ geo_triso_compact_surrounding_cell = openmc.Cell(name="TRISO Compact Surrounding
 
 geo_triso_compact_universe = openmc.Universe(name="TRISO Compact Universe",cells=[geo_triso_compact_fuel_cell,geo_triso_compact_surrounding_cell])
 
-## VWH TRISO Compact ##
+## Homogenized TRISO Compact ##
 geo_vwh_fracs = []
+geo_rpt_fracs = []
+geo_triso_rpt_region = -openmc.ZCylinder(r=cli_args.radius,x0=0.0,y0=0.0)&+geo_constant_z_min&-geo_constant_z_max
 geo_n_triso = len(geo_triso_compact_packing)
 
 # Get volume fractions for mixing
-geo_triso_compact_v_tot = math.pi*geo_constant_triso_radius**2*2*geo_constant_height
+geo_triso_compact_vwh_v_total = math.pi*geo_constant_triso_radius**2*2*geo_constant_height
+geo_triso_compact_vwh_v_matrix = geo_triso_compact_vwh_v_total - geo_n_triso * 4/3*math.pi*geo_triso_spheres[4].r**3
+geo_triso_compact_rpt_v_total = math.pi*cli_args.radius**2*2*geo_constant_height
+geo_triso_compact_rpt_v_matrix = geo_triso_compact_rpt_v_total - geo_n_triso * 4/3*math.pi*geo_triso_spheres[4].r**3
 
-geo_triso_compact_v_uco = geo_n_triso*4/3*math.pi*geo_triso_spheres[0].r**3
-geo_vwh_fracs.append(geo_triso_compact_v_uco/geo_triso_compact_v_tot)
 
-geo_triso_compact_v_buffer = geo_n_triso*4/3*math.pi*geo_triso_spheres[1].r**3 - geo_triso_compact_v_uco
-geo_vwh_fracs.append(geo_triso_compact_v_buffer/geo_triso_compact_v_tot)
+geo_v_layers = []
+geo_temp_prev_volume = 0.0
 
-geo_triso_compact_v_pyc1 = geo_n_triso*4/3*math.pi*geo_triso_spheres[2].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer)
-geo_vwh_fracs.append(geo_triso_compact_v_pyc1/geo_triso_compact_v_tot)
+for layer in geo_triso_spheres:
+    geo_temp_cumulative_volume = 4/3*math.pi*layer.r**3
+    geo_temp_shell_volume = geo_temp_cumulative_volume-geo_temp_prev_volume
+    geo_v_layers.append(geo_temp_shell_volume*geo_n_triso)
+    geo_temp_prev_volume=geo_temp_cumulative_volume
 
-geo_triso_compact_v_sic = geo_n_triso*4/3*math.pi*geo_triso_spheres[3].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer + geo_triso_compact_v_pyc1)
-geo_vwh_fracs.append(geo_triso_compact_v_sic/geo_triso_compact_v_tot)
-
-geo_triso_compact_v_pyc2 = geo_n_triso*4/3*math.pi*geo_triso_spheres[4].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer + geo_triso_compact_v_pyc1 + geo_triso_compact_v_sic)
-geo_vwh_fracs.append(geo_triso_compact_v_pyc2/geo_triso_compact_v_tot)
-
-geo_triso_compact_v_graphite = geo_triso_compact_v_tot - geo_n_triso*4/3*math.pi*geo_triso_spheres[4].r**3
-geo_vwh_fracs.append(geo_triso_compact_v_graphite/geo_triso_compact_v_tot)
+geo_vwh_fracs = [v / geo_triso_compact_vwh_v_total for v in (geo_v_layers + [geo_triso_compact_vwh_v_matrix])]
+geo_rpt_fracs = [v / geo_triso_compact_rpt_v_total for v in (geo_v_layers + [geo_triso_compact_rpt_v_matrix])]
 
 # Strip TSL temporarily for mixing
 mat_tsl_saved = {}
@@ -238,23 +236,45 @@ mat_triso_vwh = openmc.Material.mix_materials(
         )
 
 # RPT Homogenization
+mat_triso_rpt = openmc.Material.mix_materials(
+        materials=[mat_triso_uco, mat_triso_buffer, mat_triso_pyc1, mat_triso_sic, mat_triso_pyc2, mat_graphite_matrix],
+        fracs=geo_rpt_fracs,
+        percent_type="vo",
+        name="RPT Homogenized TRISO Compact"
+        )
 
-# Restore TSL afterwards
+# Restore TSL for other materials
 for mat, tsl_list in mat_tsl_saved.items():
     mat._sab = tsl_list
 
+# Add TSL for new materials
 mat_triso_vwh.add_s_alpha_beta("c_Graphite")
+mat_triso_rpt.add_s_alpha_beta("c_Graphite")
 
-# Add to list
+# Add new materials to global list
 mat_list.append(mat_triso_vwh)
+mat_list.append(mat_triso_rpt)
 materials = openmc.Materials(mat_list)
 
 # Cast into cells and a universe
-geo_triso_vwh_fuel_cell = openmc.Cell(name="TRISO VWH Cell",fill=mat_triso_vwh,region=geo_triso_compact_region)
+geo_triso_vwh_fuel_cell = openmc.Cell(name="TRISO VWH Fuel Cell",fill=mat_triso_vwh,region=geo_triso_compact_region)
 geo_triso_vwh_surrounding_cell = openmc.Cell(name="TRISO VWH Surrounding Cell",fill=mat_graphite_matrix,region=~geo_triso_compact_region)
-
 geo_triso_vwh_universe = openmc.Universe(name="TRISO VWH Universe",cells=[geo_triso_vwh_fuel_cell,geo_triso_vwh_surrounding_cell])
 
+geo_triso_rpt_fuel_cell = openmc.Cell(name="TRISO RPT Fuel Cell",fill=mat_triso_rpt,region=geo_triso_rpt_region)
+geo_triso_rpt_surrounding_cell = openmc.Cell(name="TRISO RPT Surrounding Cell",fill=mat_graphite_matrix,region=~geo_triso_rpt_region)
+geo_triso_rpt_universe = openmc.Universe(name="TRISO RPT Universe",cells=[geo_triso_rpt_fuel_cell,geo_triso_rpt_surrounding_cell])
+
+# Dynamically pick the correct one
+geo_homo_map = {
+        "none":geo_triso_compact_universe,
+        "vwh":geo_triso_vwh_universe,
+        "rpt":geo_triso_rpt_universe,
+        "rrpt":None
+        }
+
+geo_dyn_universe = geo_homo_map[cli_args.technique]
+
 ## Heat Pipe ##
 geo_pipe_cylinders = [openmc.ZCylinder(r=rad, x0=0.0, y0=0.0) for rad in [0.80, 0.90, 0.97, 1.05, 1.07]]
 geo_pipe_cells = [openmc.Cell(name="Heat Pipe K Gas",fill=mat_pipe_k_gas,region=-geo_pipe_cylinders[0]),
@@ -287,11 +307,11 @@ geo_assembly_lattice_map = []
 for i in range(0,7):
     if i%2 == 0:
         ring = [geo_pipe_universe]*(36-6*i)
-        ring[1::2] = [geo_triso_vwh_universe]*len(ring[1::2])
+        ring[1::2] = [geo_dyn_universe]*len(ring[1::2])
         geo_assembly_lattice_map.append(ring)
     else:
         ring = [geo_mod_universe]*(36-6*i)
-        ring[::2] = [geo_triso_vwh_universe]*len(ring[::2])
+        ring[::2] = [geo_dyn_universe]*len(ring[::2])
         geo_assembly_lattice_map.append(ring)
 
 geo_assembly_lattice_map[-1] = [geo_pipe_universe]