Added basic axial reflector definition and fixed driver script

1 files changed, 64 insertions, 52 deletions

diff --git a/hpmr.py b/hpmr.py
index 8a1cd42..c29bd08 100755
--- a/hpmr.py
+++ b/hpmr.py
@@ -9,9 +9,6 @@ import time
 # HP-MR Model in OpenMC
 # C.Moore, <connor.moore@psi.ch>
 
-# For LCLRS
-#openmc.config["cross_sections"] = "/libs/endfb81_official/cross_sections.xml"
-
 ############ COMMAND LINE ARGUMENTS ############ 
 cli_parser = argparse.ArgumentParser(
         prog="HP-RM TRISO Homogenization Study",
@@ -28,9 +25,7 @@ if cli_args.lclrs:
     openmc.config["cross_sections"] = "/libs/endfb81_official/cross_sections.xml"
 
 ############ INTPUT-OUTPUT MAPPING ############ 
-path = sys.argv[0]
-
-with open(path, "rb") as file:
+with open(sys.argv[0], "rb") as file:
     contents = file.read()
 
 filehash = hashlib.md5(contents).hexdigest()
@@ -430,13 +425,13 @@ mat_core_b4c_central.set_density("g/cc",1.25)
 
 ############ GEOMETRY ############ 
 # All measurements are in cm unless specified
-geo_constant_pin_pitch = 2.3
-geo_constant_lattice_pitch = 26.752/3**0.5
-geo_constant_triso_radius = 1.0
-geo_constant_triso_packing = 0.40
-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_const_pin_pitch = 2.3
+geo_const_lattice_pitch = 26.752/3**0.5
+geo_const_triso_radius = 1.0
+geo_const_triso_packing = 0.40
+geo_const_height = 1.0
+geo_const_z_max = openmc.ZPlane(z0=+geo_const_height,boundary_type="reflective")
+geo_const_z_min = openmc.ZPlane(z0=-geo_const_height,boundary_type="reflective")
 
 geo_none_universe = openmc.Universe(name="Empty Universe",cells=[openmc.Cell(name="Empty Cell",fill=None)])
 
@@ -452,14 +447,14 @@ geo_triso_cells = [openmc.Cell(name="TRISO Kernel UCO",fill=mat_triso_uco,region
 geo_triso_universe = openmc.Universe(name="TRISO Kernel Universe",cells=geo_triso_cells)
 
 ## TRISO compact ##
-geo_triso_compact_region = -openmc.ZCylinder(r=geo_constant_triso_radius,x0=0.0,y0=0.0)&+geo_constant_z_min&-geo_constant_z_max
+geo_triso_compact_region = -openmc.ZCylinder(r=geo_const_triso_radius,x0=0.0,y0=0.0)&+geo_const_z_min&-geo_const_z_max
 geo_triso_compact_domain = (10,10,10)
 
 # Generate packing
 geo_triso_compact_packing = openmc.model.pack_spheres(
         radius=geo_triso_spheres[4].r,
         region=geo_triso_compact_region,
-        pf=geo_constant_triso_packing
+        pf=geo_const_triso_packing
         )
 
 # Create the TRISO particles
@@ -468,8 +463,8 @@ geo_triso_compact_trisos = [openmc.model.TRISO(geo_triso_spheres[4].r, geo_triso
 # Support with a lattice for domain decomposition
 geo_triso_compact_lattice = openmc.model.create_triso_lattice(
         trisos=geo_triso_compact_trisos,
-        lower_left=(-geo_constant_triso_radius,-geo_constant_triso_radius,-geo_constant_height),
-        pitch=(geo_constant_triso_radius*2,geo_constant_triso_radius*2,geo_constant_height*2),
+        lower_left=(-geo_const_triso_radius,-geo_const_triso_radius,-geo_const_height),
+        pitch=(geo_const_triso_radius*2,geo_const_triso_radius*2,geo_const_height*2),
         shape=geo_triso_compact_domain,
         background=mat_graphite_matrix
         )
@@ -483,19 +478,19 @@ geo_triso_compact_universe = openmc.Universe(name="TRISO Compact Universe",cells
 ## 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_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_constant_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_constant_height))**0.5:.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_constant_triso_radius**2*2*geo_constant_height
+geo_triso_compact_vwh_v_total = math.pi*geo_const_triso_radius**2*2*geo_const_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_total = math.pi*cli_args.radius**2*2*geo_const_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 = []
@@ -558,9 +553,9 @@ for mat, tsl_list in mat_tsl_saved.items():
     mat._sab = tsl_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_r_out = math.sqrt(4*geo_n_triso*geo_triso_spheres[4].r**3/(6*geo_const_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
+                & +geo_const_z_min & -geo_const_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)
@@ -609,7 +604,7 @@ geo_mod_universe = openmc.Universe(name="Moderator Pin Universe",cells=geo_mod_c
 
 geo_assembly_lattice = openmc.HexLattice()
 geo_assembly_lattice.center=(0,0)
-geo_assembly_lattice.pitch=(geo_constant_pin_pitch,)
+geo_assembly_lattice.pitch=(geo_const_pin_pitch,)
 geo_assembly_lattice.outer=openmc.Universe(cells=[openmc.Cell(fill=mat_graphite_matrix)])
 
 # Create the lattice structure
@@ -629,17 +624,23 @@ geo_assembly_lattice.universes=geo_assembly_lattice_map
 
 # Bound it to a unit assembly
 geo_assembly_region = openmc.model.HexagonalPrism(
-        edge_length=geo_constant_lattice_pitch,
+        edge_length=geo_const_lattice_pitch,
         orientation="y",
         origin=(0.0,0.0),
         boundary_type="transmission" if cli_args.core_level else "reflective"
         )
 
 # 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_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])
 
 ## Core-level Geometry ##
+geo_const_core_z_max = openmc.ZPlane(z0=80)
+geo_const_core_z_min = openmc.ZPlane(z0=-80)
+geo_const_refl_z_max = openmc.ZPlane(z0=100, boundary_type="vacuum")
+geo_const_refl_z_min = openmc.ZPlane(z0=-100, boundary_type="vacuum")
+geo_const_core_cyl = openmc.ZCylinder(r=115, x0=0.0, y0=0.0, boundary_type="vacuum")
+
 # 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])
@@ -652,41 +653,54 @@ geo_cd_universe = openmc.Universe(name="Core Control Drum Universe",cells=[geo_c
 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.pitch=(geo_const_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_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([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
 
+# Cell and universe
+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])
 
-# Reflector lattice map
+
+# Axial Reflector lattice
 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.orientation="x"
+geo_core_rad_refl_lattice.pitch=(geo_const_lattice_pitch*3**0.5,)
 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])
+# 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([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.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,
+                                     region=-geo_const_core_cyl & +geo_const_core_z_max & -geo_const_refl_z_max)
+geo_core_refl_bot_cell = openmc.Cell(name="Bottom Axial Reflector Cell",fill=geo_core_rad_refl_lattice,
+                                     region=-geo_const_core_cyl & -geo_const_core_z_min & +geo_const_refl_z_min)
 
 # 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.lower_left = (-sys.float_info.max/2, -sys.float_info.max/2, -int(160//geo_const_height)/2*geo_const_height)
+geo_core_axial_lattice.pitch=(sys.float_info.max, sys.float_info.max, geo_const_height)
+geo_core_axial_lattice.universes = [[[geo_core_rad_universe]] for _ in range(int(160//geo_const_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])
+geo_core_cell = openmc.Cell(name="Core 3D Cell",fill=geo_core_axial_lattice, region=-geo_const_core_cyl & +geo_const_core_z_min & -geo_const_core_z_max)
+geo_core_universe = openmc.Universe(cells=[geo_core_cell, geo_core_refl_top_cell, geo_core_refl_bot_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)
@@ -696,17 +710,15 @@ settings = openmc.Settings()
 settings.particles = 100000
 settings.batches = 200
 settings.inactive = 50
+#settings.verbosity = 10
 
-# Source sampling
-#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_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_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)
+            lower_left = (-set_dyn_radial, -set_dyn_radial, -set_dyn_axial),
+            upper_right = (set_dyn_radial, set_dyn_radial, set_dyn_axial)
             ))
         
 
@@ -714,8 +726,8 @@ settings.source = set_source_pts
 
 # Mesh for Shannon Entropy
 set_entropy_mesh = openmc.RegularMesh()
-set_entropy_mesh.lower_left=(-geo_constant_lattice_pitch,-geo_constant_lattice_pitch,-geo_constant_height)
-set_entropy_mesh.upper_right=(+geo_constant_lattice_pitch,+geo_constant_lattice_pitch,+geo_constant_height)
+set_entropy_mesh.lower_left=(-set_dyn_radial,-set_dyn_radial,-set_dyn_axial)
+set_entropy_mesh.upper_right=(+set_dyn_radial,+set_dyn_radial,+set_dyn_axial)
 set_entropy_mesh.dimension=(30,30,30)
 
 settings.entropy_mesh=set_entropy_mesh