Added support for pin-level calculations (by default)

2 files changed, 22 insertions, 21 deletions

diff --git a/hpmr.py b/hpmr.py
index c29bd08..f05658b 100755
--- a/hpmr.py
+++ b/hpmr.py
@@ -18,7 +18,8 @@ cli_parser = argparse.ArgumentParser(
 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.8",type=float,help="Radius for RPT or RRPT homogenization [cm], default 0.8 cm. In the case of RRPT, this is the inner radius.")
 cli_parser.add_argument("-l","--lclrs",default=False,action="store_true",help="Flag for specifying LCLRS accessible cross-sections, default false.")
-cli_parser.add_argument("-c","--core-level",default=False,action="store_true",help="Whether or not to run a core-level computation, default false.")
+cli_parser.add_argument("-g","--geometry-level",default="pin",choices=["pin","ass","core"],help="Geometry level to run the computation at [pin/ass/core], default pin.")
+cli_parser.add_argument("-e","--execute",default=False,action="store_true",help="Whether or not to automatically run OpenMC, default false.")
 cli_args = cli_parser.parse_args()
 
 if cli_args.lclrs:
@@ -30,11 +31,11 @@ with open(sys.argv[0], "rb") as file:
 
 filehash = hashlib.md5(contents).hexdigest()
 
-print("="*90)
+print("="*110)
 print(f"Run date CET:     {time.ctime()}")
 print(f"Input file MD5:   {filehash}")
 print(f"Input arguments:  {cli_args}")
-print("="*90)
+print("="*110)
 
 ############ MATERIALS ############ 
 # All materials are specified in the report on HEAT PIPE MICROREACTOR MODELING WITH BLUECRAB
@@ -480,12 +481,12 @@ geo_vwh_fracs = []
 geo_rpt_fracs = []
 geo_triso_rpt_region = -openmc.ZCylinder(r=cli_args.radius,x0=0.0,y0=0.0)&+geo_const_z_min&-geo_const_z_max
 geo_n_triso = len(geo_triso_compact_packing)
-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_const_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_const_height))**0.5:.4f}")
-print("    TRISO Stacking calculation complete!\n")
+print("\nTRISO 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_const_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_const_height))**0.5:.4f}")
+print("TRISO Stacking calculation complete!\n")
 
 # Get volume fractions for mixing
 geo_triso_compact_vwh_v_total = math.pi*geo_const_triso_radius**2*2*geo_const_height
@@ -618,16 +619,16 @@ for i in range(0,7):
         ring[::2] = [geo_dyn_universe]*len(ring[::2])
     geo_assembly_lattice_map.append(ring)
 
-geo_assembly_lattice_map[-1] = [geo_pipe_universe]
+geo_assembly_lattice_map[-1] = [geo_pipe_universe if cli_args.geometry_level!="pin" else geo_dyn_universe]
 
 geo_assembly_lattice.universes=geo_assembly_lattice_map
 
 # Bound it to a unit assembly
 geo_assembly_region = openmc.model.HexagonalPrism(
-        edge_length=geo_const_lattice_pitch,
-        orientation="y",
+        edge_length=geo_const_lattice_pitch if cli_args.geometry_level!="pin" else geo_const_pin_pitch/3**0.5,
+        orientation="y" if cli_args.geometry_level !="pin" else "x",
         origin=(0.0,0.0),
-        boundary_type="transmission" if cli_args.core_level else "reflective"
+        boundary_type="transmission" if cli_args.geometry_level=="core" else "reflective"
         )
 
 # Clip the lattice to the space and make a universe
@@ -662,7 +663,7 @@ geo_core_rad_lattice_map.append([geo_cd_universe if i%4==2 else geo_none_univers
 geo_core_rad_lattice_map.append([geo_cd_universe if i%3==0 else geo_assembly_universe for i in range(18)])
 geo_core_rad_lattice_map.append([geo_assembly_universe]*12)
 geo_core_rad_lattice_map.append([geo_assembly_universe]*6)
-geo_core_rad_lattice_map.append([geo_refl_universe])
+geo_core_rad_lattice_map.append([geo_none_universe])
 geo_core_rad_lattice.universes=geo_core_rad_lattice_map
 
 # Cell and universe
@@ -683,7 +684,7 @@ geo_core_rad_refl_lattice_map.append([geo_cd_universe if i%4==2 else geo_none_un
 geo_core_rad_refl_lattice_map.append([geo_cd_universe if i%3==0 else geo_refl_universe for i in range(18)])
 geo_core_rad_refl_lattice_map.append([geo_refl_universe]*12)
 geo_core_rad_refl_lattice_map.append([geo_refl_universe]*6)
-geo_core_rad_refl_lattice_map.append([geo_refl_universe])
+geo_core_rad_refl_lattice_map.append([geo_none_universe])
 geo_core_rad_refl_lattice.universes = geo_core_rad_refl_lattice_map
 
 geo_core_refl_top_cell = openmc.Cell(name="Top Axial Reflector Cell",fill=geo_core_rad_refl_lattice,
@@ -703,7 +704,7 @@ geo_core_universe = openmc.Universe(cells=[geo_core_cell, geo_core_refl_top_cell
 #geo_core_universe = openmc.Universe(cells=[geo_core_refl_bot_cell])
 
 ## Final Geometry Universe
-geometry = openmc.Geometry(geo_core_universe if cli_args.core_level else geo_assembly_universe)
+geometry = openmc.Geometry(geo_core_universe if cli_args.geometry_level=="core" else geo_assembly_universe)
 
 ############ SETTINGS ############ 
 settings = openmc.Settings()
@@ -712,8 +713,8 @@ settings.batches = 200
 settings.inactive = 50
 #settings.verbosity = 10
 
-set_dyn_radial = geo_const_core_cyl.r if cli_args.core_level else geo_const_lattice_pitch
-set_dyn_axial = 180 if cli_args.core_level else geo_const_height
+set_dyn_radial = geo_const_core_cyl.r if cli_args.geometry_level=="core" else geo_const_pin_pitch if cli_args.geometry_level=="pin" else geo_const_lattice_pitch
+set_dyn_axial = 180 if cli_args.geometry_level=="core" else geo_const_height
 
 set_source_pts = openmc.IndependentSource(
         space=openmc.stats.Box(
@@ -741,4 +742,4 @@ geometry.export_to_xml()
 settings.export_to_xml()
 
 # Run!
-openmc.run()
+openmc.run() if cli_args.execute else None
diff --git a/triso_study_driver.sh b/triso_study_driver.sh
index b7404c8..3c9e58a 100755
--- a/triso_study_driver.sh
+++ b/triso_study_driver.sh
@@ -1,9 +1,9 @@
 #!/bin/bash
 
 if [[ $1 == "lcl" ]]; then
-    prog=(openmc-inst01.py hpmr.py -l)
+    prog=(openmc-inst01.py hpmr.py -l -g pin)
 else
-    prog=./hpmr.py
+    prog=(./hpmr.py -g pin)
 fi
 
 mkdir -p results/ results/rpt results/rrpt