Finished RPT support and added RRPT support

1 files changed, 21 insertions, 2 deletions

diff --git a/hpmr.py b/hpmr.py
index a11b066..3864214 100755
--- a/hpmr.py
+++ b/hpmr.py
@@ -207,7 +207,6 @@ geo_triso_compact_vwh_v_matrix = geo_triso_compact_vwh_v_total - geo_n_triso * 4
 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_v_layers = []
 geo_temp_prev_volume = 0.0
 
@@ -217,6 +216,7 @@ for layer in geo_triso_spheres:
     geo_v_layers.append(geo_temp_shell_volume*geo_n_triso)
     geo_temp_prev_volume=geo_temp_cumulative_volume
 
+geo_kernel_fracs = [v / sum(geo_v_layers) for v in geo_v_layers]
 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])]
 
@@ -227,6 +227,14 @@ for mat in mat_list:
         mat_tsl_saved[mat]=list(mat._sab)
         mat._sab=[]
 
+# TRISO only VWH
+mat_kernel_vwh = openmc.Material.mix_materials(
+        materials=[mat_triso_uco, mat_triso_buffer, mat_triso_pyc1, mat_triso_sic, mat_triso_pyc2],
+        fracs=geo_kernel_fracs,
+        percent_type="vo",
+        name="VWH kernel"
+        )
+
 # Simple VWH
 mat_triso_vwh = openmc.Material.mix_materials(
         materials=[mat_triso_uco, mat_triso_buffer, mat_triso_pyc1, mat_triso_sic, mat_triso_pyc2, mat_graphite_matrix],
@@ -248,14 +256,21 @@ for mat, tsl_list in mat_tsl_saved.items():
     mat._sab = tsl_list
 
 # Add TSL for new materials
+mat_kernel_vwh.add_s_alpha_beta("c_Graphite")
 mat_triso_vwh.add_s_alpha_beta("c_Graphite")
 mat_triso_rpt.add_s_alpha_beta("c_Graphite")
 
 # Add new materials to global list
+mat_list.append(mat_kernel_vwh)
 mat_list.append(mat_triso_vwh)
 mat_list.append(mat_triso_rpt)
 materials = openmc.Materials(mat_list)
 
+# RRPT transformation
+geo_rrpt_r_out = math.sqrt(4*geo_n_triso*geo_triso_spheres[4].r**3/(6*geo_constant_height) + cli_args.radius**2)
+geo_rrpt_region = +openmc.ZCylinder(r=cli_args.radius,x0=0.0,y0=0.0) & -openmc.ZCylinder(r=geo_rrpt_r_out,x0=0.0,y0=0.0) \
+                & +geo_constant_z_min & -geo_constant_z_max
+
 # Cast into cells and a universe
 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)
@@ -265,12 +280,16 @@ geo_triso_rpt_fuel_cell = openmc.Cell(name="TRISO RPT Fuel Cell",fill=mat_triso_
 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])
 
+geo_triso_rrpt_fuel_cell = openmc.Cell(name="TRISO RRPT Fuel Cell",fill=mat_kernel_vwh,region=geo_rrpt_region)
+geo_triso_rrpt_surrouding_cell = openmc.Cell(name="TRISO RRPT Surrounding Cell",fill=mat_graphite_matrix,region=~geo_rrpt_region)
+geo_triso_rrpt_universe = openmc.Universe(name="TRISO RRPT Universe",cells=[geo_triso_rrpt_fuel_cell,geo_triso_rrpt_surrouding_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
+        "rrpt":geo_triso_rrpt_universe
         }
 
 geo_dyn_universe = geo_homo_map[cli_args.technique]