Added core-level control drums

1 files changed, 36 insertions, 16 deletions

diff --git a/hpmr.py b/hpmr.py
index f05658b..4fb858b 100755
--- a/hpmr.py
+++ b/hpmr.py
@@ -5,6 +5,18 @@ import argparse
 import hashlib
 import sys
 import time
+import string
+
+
+############ HLPER FUNCTIONS ############ 
+
+def rot(base: openmc.Universe, angle: float = 0.0) -> openmc.Cell:
+    # Take the cells from the existing universe and create a new rotated universe
+    new_cell = openmc.Cell(fill=base)
+    new_cell.rotation=(0.0,0.0,angle)
+    new_cell.name = f"Control Drum Cell, Rotated {angle} deg"
+    return openmc.Universe(name=base.name + f" Rotated {angle} deg", cells=[new_cell])
+
 
 # HP-MR Model in OpenMC
 # C.Moore, <connor.moore@psi.ch>
@@ -18,7 +30,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("-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("-g","--geometry-level",default="core",choices=["pin","ass","core"],help="Geometry level to run the computation at [pin/ass/core], default pin.")
+cli_parser.add_argument("-d","--drums-in",default=False,action="store_true",help="Flag for specifying the control drums are rotated in, default false")
 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()
 
@@ -31,11 +44,11 @@ with open(sys.argv[0], "rb") as file:
 
 filehash = hashlib.md5(contents).hexdigest()
 
-print("="*110)
+print("="*52)
 print(f"Run date CET:     {time.ctime()}")
 print(f"Input file MD5:   {filehash}")
-print(f"Input arguments:  {cli_args}")
-print("="*110)
+print(f"Input arguments:  {str(cli_args)[10:-1].replace(',' , '\n'+17*' ').replace('=',' = ')}")
+print("="*52)
 
 ############ MATERIALS ############ 
 # All materials are specified in the report on HEAT PIPE MICROREACTOR MODELING WITH BLUECRAB
@@ -415,6 +428,7 @@ mat_core_b4c_control_drum.add_nuclide("B11",0.04)
 mat_core_b4c_control_drum.add_nuclide("C12",0.20)
 mat_core_b4c_control_drum.temperature = mat_temp_global
 mat_core_b4c_control_drum.set_density("g/cc",2.51)
+mat_list.append(mat_core_b4c_control_drum)
 
 mat_core_b4c_central = openmc.Material(503,name="Core B4C Central Absorber")
 mat_core_b4c_central.add_nuclide("B10",0.76)
@@ -422,12 +436,13 @@ mat_core_b4c_central.add_nuclide("B11",0.04)
 mat_core_b4c_central.add_nuclide("C12",0.20)
 mat_core_b4c_central.temperature = mat_temp_global
 mat_core_b4c_central.set_density("g/cc",1.25)
+mat_list.append(mat_core_b4c_central)
 
 
 ############ GEOMETRY ############ 
 # All measurements are in cm unless specified
 geo_const_pin_pitch = 2.3
-geo_const_lattice_pitch = 26.752/3**0.5
+geo_const_lattice_pitch = 26.752
 geo_const_triso_radius = 1.0
 geo_const_triso_packing = 0.40
 geo_const_height = 1.0
@@ -625,7 +640,7 @@ 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 if cli_args.geometry_level!="pin" else geo_const_pin_pitch/3**0.5,
+        edge_length=geo_const_lattice_pitch/3**0.5 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.geometry_level=="core" else "reflective"
@@ -633,7 +648,7 @@ geo_assembly_region = openmc.model.HexagonalPrism(
 
 # Clip the lattice to the space and make a universe
 geo_assembly_cell = openmc.Cell(name="Unit Assembly Cell",fill=geo_assembly_lattice,region=-geo_assembly_region&+geo_const_z_min&-geo_const_z_max)
-geo_assembly_universe = openmc.Universe(cells=[geo_assembly_cell])
+geo_assembly_universe = openmc.Universe(name="Unit Lattice Assembly Universe",cells=[geo_assembly_cell])
 
 ## Core-level Geometry ##
 geo_const_core_z_max = openmc.ZPlane(z0=80)
@@ -647,20 +662,25 @@ geo_refl_cell = openmc.Cell(name="Core Reflector Cell",fill=mat_core_be_reflecto
 geo_refl_universe = openmc.Universe(name="Core Reflector Universe",cells=[geo_refl_cell])
 
 # Control Drum cell and universe
-geo_cd_cell = openmc.Cell(name="Core Control Drum Cell",fill=mat_core_be_control_drum)
-geo_cd_universe = openmc.Universe(name="Core Control Drum Universe",cells=[geo_cd_cell])
+geo_cd_cyl = [openmc.ZCylinder(r=rad, x0=0.0, y0=0.0) for rad in [12.250, 13.250]]
+geo_cd_pl = openmc.XPlane(x0=88.92-3*geo_const_lattice_pitch)
+geo_cd_abs_cell = openmc.Cell(name="Control Drum Absorber Cell",fill=mat_core_b4c_control_drum,region=+geo_cd_cyl[0]&-geo_cd_cyl[1]&+geo_cd_pl)
+geo_cd_drum_cell = openmc.Cell(name="Core Control Drum Cell",fill=mat_core_be_control_drum,region=~geo_cd_abs_cell.region&-geo_cd_cyl[1])
+geo_cd_outer_cell = openmc.Cell(name="Core Control Drum Outer Cell",fill=None,region=+geo_cd_cyl[1])
+geo_cd_universe = openmc.Universe(name="Core Control Drum Universe",cells=[geo_cd_abs_cell,geo_cd_drum_cell,geo_cd_outer_cell])
 
 # Core 2D radial lattice
 geo_core_rad_lattice=openmc.HexLattice()
 geo_core_rad_lattice.orientation="x"
 geo_core_rad_lattice.center=(0,0)
-geo_core_rad_lattice.pitch=(geo_const_lattice_pitch*3**0.5,)
+geo_core_rad_lattice.pitch=(geo_const_lattice_pitch,)
 geo_core_rad_lattice.outer=geo_none_universe
 
 # Assemble lattice map
+geo_dyn_cd_rot = 180 if cli_args.drums_in else 0
 geo_core_rad_lattice_map=[]
-geo_core_rad_lattice_map.append([geo_cd_universe if i%4==2 else geo_none_universe if i%4==0 else geo_refl_universe for i in range(24)])
-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([rot(geo_cd_universe,geo_dyn_cd_rot-30-(i-2)*15) if i%4==2 else geo_none_universe if i%4==0 else geo_refl_universe for i in range(24)])
+geo_core_rad_lattice_map.append([rot(geo_cd_universe,geo_dyn_cd_rot-i*20) 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_none_universe])
@@ -675,13 +695,13 @@ geo_core_rad_universe = openmc.Universe(cells=[geo_core_rad_cell])
 geo_core_rad_refl_lattice=openmc.HexLattice()
 geo_core_rad_refl_lattice.center=(0,0)
 geo_core_rad_refl_lattice.orientation="x"
-geo_core_rad_refl_lattice.pitch=(geo_const_lattice_pitch*3**0.5,)
+geo_core_rad_refl_lattice.pitch=(geo_const_lattice_pitch,)
 geo_core_rad_refl_lattice.outer=geo_none_universe
 
 # Axial Reflector lattice map
 geo_core_rad_refl_lattice_map = []
-geo_core_rad_refl_lattice_map.append([geo_cd_universe if i%4==2 else geo_none_universe if i%4==0 else geo_refl_universe for i in range(24)])
-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([rot(geo_cd_universe,geo_dyn_cd_rot-30-(i-2)*15) if i%4==2 else geo_none_universe if i%4==0 else geo_refl_universe for i in range(24)])
+geo_core_rad_refl_lattice_map.append([rot(geo_cd_universe,geo_dyn_cd_rot-i*20) 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_none_universe])
@@ -713,7 +733,7 @@ settings.batches = 200
 settings.inactive = 50
 #settings.verbosity = 10
 
-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_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/3**0.5
 set_dyn_axial = 180 if cli_args.geometry_level=="core" else geo_const_height
 
 set_source_pts = openmc.IndependentSource(