WIP 3D core geometry support

1 files changed, 69 insertions, 11 deletions

diff --git a/hpmr.py b/hpmr.py
index 8551f7d..8a1cd42 100755
--- a/hpmr.py
+++ b/hpmr.py
@@ -21,6 +21,7 @@ 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_args = cli_parser.parse_args()
 
 if cli_args.lclrs:
@@ -34,11 +35,11 @@ with open(path, "rb") as file:
 
 filehash = hashlib.md5(contents).hexdigest()
 
-print("="*70)
+print("="*90)
 print(f"Run date CET:     {time.ctime()}")
 print(f"Input file MD5:   {filehash}")
 print(f"Input arguments:  {cli_args}")
-print("="*70)
+print("="*90)
 
 ############ MATERIALS ############ 
 # All materials are specified in the report on HEAT PIPE MICROREACTOR MODELING WITH BLUECRAB
@@ -437,6 +438,7 @@ geo_constant_height = 1.0
 geo_constant_z_max = openmc.ZPlane(z0=+geo_constant_height,boundary_type="reflective")
 geo_constant_z_min = openmc.ZPlane(z0=-geo_constant_height,boundary_type="reflective")
 
+geo_none_universe = openmc.Universe(name="Empty Universe",cells=[openmc.Cell(name="Empty Cell",fill=None)])
 
 ### Inidividual Pin Universes ###
 
@@ -616,11 +618,10 @@ for i in range(0,7):
     if i%2 == 0:
         ring = [geo_pipe_universe]*(36-6*i)
         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_dyn_universe]*len(ring[::2])
-        geo_assembly_lattice_map.append(ring)
+    geo_assembly_lattice_map.append(ring)
 
 geo_assembly_lattice_map[-1] = [geo_pipe_universe]
 
@@ -631,14 +632,64 @@ geo_assembly_region = openmc.model.HexagonalPrism(
         edge_length=geo_constant_lattice_pitch,
         orientation="y",
         origin=(0.0,0.0),
-        boundary_type="reflective"
+        boundary_type="transmission" if cli_args.core_level else "reflective"
         )
 
-# Clip the lattice to the space and make it root
+# 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_constant_z_min&-geo_constant_z_max)
 geo_assembly_universe = openmc.Universe(cells=[geo_assembly_cell])
 
-geometry = openmc.Geometry(geo_assembly_universe)
+## Core-level Geometry ##
+# Reflector region cell and universe
+geo_refl_cell = openmc.Cell(name="Core Reflector Cell",fill=mat_core_be_reflector)
+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])
+
+# 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_constant_lattice_pitch*3**0.5,)
+geo_core_rad_lattice.outer=geo_none_universe
+
+# Assemble lattice map
+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 enumerate(range(24))])
+#geo_core_rad_lattice_map.append([geo_assembly_universe]*18)
+geo_core_rad_lattice_map.append([geo_cd_universe if i%3==0 else geo_assembly_universe for i,_ in enumerate(range(18))])
+#print(f"{type(geo_core_rad_lattice_map[1]) = }")
+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.universes=geo_core_rad_lattice_map
+
+
+# Reflector lattice map
+geo_core_rad_refl_lattice=openmc.HexLattice()
+geo_core_rad_refl_lattice.center=(0,0)
+geo_core_rad_refl_lattice.pitch=(geo_constant_lattice_pitch,)
+geo_core_rad_refl_lattice.outer=geo_none_universe
+geo_core_rad_refl_lattice.universes=[[geo_refl_universe]*n for n in [24,18,12,6,1]]
+
+geo_core_rad_cell=openmc.Cell(name="Core Radial Slice Cell",fill=geo_core_rad_lattice)
+geo_core_rad_universe = openmc.Universe(cells=[geo_core_rad_cell])
+
+# 3D Stacking Lattice
+geo_core_axial_lattice = openmc.RectLattice()
+geo_core_axial_lattice.lower_left = (-sys.float_info.max/2, -sys.float_info.max/2, -int(160//geo_constant_height)/2*geo_constant_height)
+#geo_core_axial_lattice.center=(0.0,0.0,0.0)
+geo_core_axial_lattice.pitch=(sys.float_info.max, sys.float_info.max, geo_constant_height)
+geo_core_axial_lattice.universes = [[[geo_core_rad_universe]] for _ in range(int(160//geo_constant_height))]
+geo_core_axial_lattice.outer=geo_none_universe
+
+geo_core_cell = openmc.Cell(name="Core 3D Cell",fill=geo_core_axial_lattice)
+geo_core_universe = openmc.Universe(cells=[geo_core_cell])
+
+## Final Geometry Universe
+geometry = openmc.Geometry(geo_core_universe if cli_args.core_level else geo_assembly_universe)
 
 ############ SETTINGS ############ 
 settings = openmc.Settings()
@@ -647,10 +698,17 @@ settings.batches = 200
 settings.inactive = 50
 
 # Source sampling
-set_source_pts = []
-for loc in geo_triso_compact_packing:
-    point = openmc.stats.Point(xyz=loc)
-    set_source_pts.append(openmc.IndependentSource(space=point))
+#set_source_pts = []
+#for loc in geo_triso_compact_packing[::100]:
+#    point = openmc.stats.Point(xyz=loc)
+#    set_source_pts.append(openmc.IndependentSource(space=point))
+
+set_source_pts = openmc.IndependentSource(
+        space=openmc.stats.Box(
+            lower_left = (-geo_constant_lattice_pitch, -geo_constant_lattice_pitch, -geo_constant_height),
+            upper_right = (geo_constant_lattice_pitch, geo_constant_lattice_pitch, geo_constant_height)
+            ))
+        
 
 settings.source = set_source_pts