Fixed geometry issue by moving entire reactor down so that boundary is not on z=0

1 files changed, 68 insertions, 61 deletions

diff --git a/krusty.py b/krusty.py
index f010689..6321512 100755
--- a/krusty.py
+++ b/krusty.py
@@ -147,77 +147,69 @@ materials.export_to_xml()
 
 ############ GEOMETRY ############ 
 geo_hp_pitch = 5.2
-geo_z_abs_bot = openmc.ZPlane(z0=0.0,boundary_type="vacuum")
-geo_z_core_bot = openmc.ZPlane(z0=10.16)
-geo_z_core_top = openmc.ZPlane(z0=35.16)
-geo_z_abs_top = openmc.ZPlane(z0=45.32,boundary_type="vacuum")
+geo_z_abs_bot = openmc.ZPlane(z0=-10.16,boundary_type="vacuum")
+geo_z_core_bot = openmc.ZPlane(z0=0.00)
+geo_z_core_top = openmc.ZPlane(z0=25)
+geo_z_abs_top = openmc.ZPlane(z0=35.16,boundary_type="vacuum")
 geo_z_region = +geo_z_abs_bot & -geo_z_abs_top
 
-## Central Air Cylinder ##
-geo_air_cyl = openmc.ZCylinder(r=1.9939, x0=0.0, y0=0.0)
-geo_air_region = -geo_air_cyl & +geo_z_core_bot
-geo_air_cell = openmc.Cell(name="Inner Core Air", region=geo_air_region, fill=mat_air)
-
-## Fuel ##
-geo_fuel_cyl = openmc.ZCylinder(r=5.4991, x0=0.0, y0=0.0)
-geo_fuel_region = -geo_fuel_cyl & +geo_air_cyl & -geo_z_core_top & +geo_z_core_bot
-geo_fuel_cell = openmc.Cell(name="Fuel Block", region=geo_fuel_region, fill=mat_fuel)
-
-## BeO Inner Reflectors ##
-
-# Bottom Reflector
-geo_beo_lower_region = -geo_fuel_cyl & -geo_z_core_bot
-geo_beo_lower_cell = openmc.Cell(name="Lower BeO Reflector", region=geo_beo_lower_region, fill=mat_beo)
-
-# Top reflector
-geo_beo_upper_region = +geo_air_cyl & -geo_fuel_cyl & +geo_z_core_top
-geo_beo_upper_cell = openmc.Cell(name="Upper BeO Reflector", region=geo_beo_upper_region, fill=mat_beo)
-
-## Radial Layers ##
-# Air, SS316, Air, BeO
-geo_rad_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [12.7, 13.3, 14.11, 38.40]]
-geo_rad_cyl[3].boundary_type="vacuum"
-geo_rad_cells = [openmc.Cell(name="Intermediate Core Air",region=+geo_fuel_cyl&-geo_rad_cyl[0],fill=mat_air),
-                 openmc.Cell(name="Core SS316 Cell",region=+geo_rad_cyl[0]&-geo_rad_cyl[1],fill=mat_ss316),
-                 openmc.Cell(name="Outer Core Air",region=+geo_rad_cyl[1]&-geo_rad_cyl[2],fill=mat_air),
-                 openmc.Cell(name="Outer BeO Reflector",region=+geo_rad_cyl[2],fill=mat_beo)]
-
-# Export simplified universe without heat pipes
-geo_rad_universe = openmc.Universe(name="Simplified Radial Universe",cells=[geo_air_cell,geo_fuel_cell,geo_beo_lower_cell,geo_beo_upper_cell]+geo_rad_cells)
-
 ## Heat Pipes ##
 geo_hp_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [1.092, 1.270]]
-geo_hp_unit_cell = [openmc.Cell(name="Heat Pipe Coolant Cell",region=-geo_hp_cyl[0],fill=mat_hp_na),
-                openmc.Cell(name="Heat Pipe Wall Cell",region=+geo_hp_cyl[0],fill=mat_haynes230)]
-                #openmc.Cell(name="Heat Pipe Surrounding Cell",region=+geo_hp_cyl[1],fill=mat_air)]
-
-geo_hp_unit_universe = openmc.Universe(name="Heat Pipe Universe",cells=geo_hp_unit_cell)
 
 # Pattern around the core
-geo_hp_steps = [math.pi*theta/180 for theta in range(0,360,45)]
-geo_hp_regions = [openmc.ZCylinder(r=geo_hp_cyl[1].r,
+geo_hp_steps = [math.pi*(theta+90)/180 for theta in range(0,360,45)]
+
+# Define inner and outer regions
+geo_hp_inner_regions = [openmc.ZCylinder(r=geo_hp_cyl[0].r,
                                    x0=geo_hp_pitch*math.cos(theta),
                                    y0=geo_hp_pitch*math.sin(theta)
                                    ) for theta in geo_hp_steps]
 
-geo_hp_cells = [openmc.Cell(name=f"Heat Pipe {num}",region=-space,fill=geo_hp_unit_universe) for num,space in enumerate(geo_hp_regions, start=1)]
+geo_hp_outer_regions = [openmc.ZCylinder(r=geo_hp_cyl[1].r,
+                                   x0=geo_hp_pitch*math.cos(theta),
+                                   y0=geo_hp_pitch*math.sin(theta)
+                                   ) for theta in geo_hp_steps]
+
+geo_hp_outer_union = openmc.Union(-reg for reg in geo_hp_outer_regions)
+
+# Fill cells for each pipe
+geo_hp_inner_cells = [openmc.Cell(name=f"Heat Pipe {num} Inner Cell",region=-reg&geo_z_region,fill=mat_hp_na) for num,reg in enumerate(geo_hp_inner_regions)]
+geo_hp_outer_cells = [openmc.Cell(name=f"Heat Pipe {num} Outer Cell",region=-reg_out&+reg_in&geo_z_region,fill=mat_haynes230)
+                      for num,(reg_out,reg_in) in enumerate(zip(geo_hp_outer_regions,geo_hp_inner_regions))]
 
-for cell, region in zip(geo_hp_cells, geo_hp_regions):
-    cell.translation = (region.x0, region.y0, 0.0)
+## Core Air Cylinder ##
+geo_air_inner_cyl = openmc.ZCylinder(r=1.9939, x0=0.0, y0=0.0)
+geo_air_inner_region = +geo_z_core_bot & -geo_z_abs_top & -geo_air_inner_cyl
+geo_air_inner_cell = openmc.Cell(name="Inner Air Cell",region=geo_air_inner_region,fill=mat_air)
 
-geo_hp_universe = openmc.Universe(name="Heat Pipes Universe",cells=geo_hp_cells)
+## Core Fuel Block ##
+geo_fuel_cyl = openmc.ZCylinder(r=5.4991, x0=0.0, y0=0.0)
+geo_fuel_region = +geo_z_core_bot & -geo_z_core_top & -geo_fuel_cyl & +geo_air_inner_cyl & ~geo_hp_outer_union
+#geo_fuel_region = +geo_z_abs_bot & -geo_z_abs_top & -geo_fuel_cyl
+geo_fuel_cell = openmc.Cell(name="Core Fuel Cell",region=geo_fuel_region,fill=mat_fuel)
+
+## BeO Reflectors ##
+
+# Lower Reflector
+geo_refl_lower_region = +geo_z_abs_bot & -geo_z_core_bot & -geo_fuel_cyl & ~geo_hp_outer_union
+geo_refl_lower_cell = openmc.Cell(name="Lower BeO Reflector Cell",region=geo_refl_lower_region,fill=mat_beo)
 
-# Combine universes
-geo_hp_comb_cell = openmc.Cell(name="Combined Heat Pipes Cell",region=openmc.Union(-cyl for cyl in geo_hp_regions),fill=geo_hp_universe)
-geo_rad_comb_cell = openmc.Cell(name="Combined Radial Layers Cell",region=~openmc.Union(-cyl for cyl in geo_hp_regions),fill=geo_rad_universe)
+# Upper Reflector
+geo_refl_upper_region = +geo_z_core_top & -geo_z_abs_top & +geo_air_inner_cyl & -geo_fuel_cyl & ~geo_hp_outer_union
+geo_refl_upper_cell = openmc.Cell(name="Upper BeO Reflector Cell",region=geo_refl_upper_region,fill=mat_beo)
 
-geo_root_universe = openmc.Universe(name="Root Universe",cells=[geo_rad_comb_cell,geo_hp_comb_cell])
-geo_clip_cell = openmc.Cell(name="Clipped Root Universe Cell",region=-geo_rad_cyl[3],fill=geo_root_universe)
-geo_clip_universe = openmc.Universe(name="Clipped Root Universe",cells=[geo_clip_cell])
+## Radial Core Surroundings ##
+geo_rad_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [12.70, 13.30, 14.11, 38.40]]
+geo_rad_cyl[-1].boundary_type="vacuum"
+geo_rad_cells = [openmc.Cell(name="Intermediate Air Cell",region=+geo_fuel_cyl&-geo_rad_cyl[0]&geo_z_region,fill=mat_air),
+                 openmc.Cell(name="Steel Annulus Cell",region=+geo_rad_cyl[0]&-geo_rad_cyl[1]&geo_z_region,fill=mat_ss316),
+                 openmc.Cell(name="Outer Air Cell",region=+geo_rad_cyl[1]&-geo_rad_cyl[2]&geo_z_region,fill=mat_air),
+                 openmc.Cell(name="Outer BeO Reflector Cell",region=+geo_rad_cyl[2]&-geo_rad_cyl[3]&geo_z_region,fill=mat_beo)]
 
-geo_object = openmc.Geometry(geo_clip_universe)
-geo_object.export_to_xml()
+geo_root_universe = openmc.Universe(name="Root Universe",cells=[geo_fuel_cell,geo_air_inner_cell,geo_refl_upper_cell,geo_refl_lower_cell]+geo_rad_cells+geo_hp_inner_cells+geo_hp_outer_cells)
 
+geo_obj = openmc.Geometry(geo_root_universe)
+geo_obj.export_to_xml()
 
 ## Settings ##
 set_obj = openmc.Settings()
@@ -227,16 +219,31 @@ set_obj.inactive = 100
 set_obj.temperature = {"method": "interpolation"}
 
 set_src = openmc.IndependentSource()
+set_src.space = openmc.stats.Point((0.0,0.0,30.0))
 set_src.space = openmc.stats.CylindricalIndependent(
-        r = openmc.stats.Uniform(geo_air_cyl.r,geo_fuel_cyl.r),
+        r = openmc.stats.Uniform(geo_air_inner_cyl.r,geo_fuel_cyl.r-1.0),
+        #r = openmc.stats.Uniform(0,2),
         z = openmc.stats.Uniform(geo_z_core_bot.z0,geo_z_core_top.z0),
         phi = openmc.stats.Uniform(0.0,2*math.pi),
-        origin = (0.0,0.0,geo_z_core_bot.z0)
+        origin = (0.0,0.0,0.0)
         )
-
-set_src.track = [(1,1,i) for i in range(1,20)]
-set_src.max_events = 2000
-
+#set_src.space = openmc.stats.Box(lower_left=(-2,-2,-2),upper_right=(2,2,6))
+
+## DEBUG SOURCE
+
+#for i in range(10000):
+#    sample_r = set_src.space.r.sample()
+#    sample_z = set_src.space.z.sample()
+#    sample_phi = set_src.space.phi.sample()
+#    #print(f"{sample_r = }, {sample_z = }, {sample_phi = }")
+#    x = sample_r*math.cos(sample_phi[0])
+#    y = sample_r*math.sin(sample_phi[0])
+#
+#    cell = geo_obj.find((x,y,sample_z))[-1]
+#    print(f"Particle {i}: ({x}, {y}, {sample_z}) -> Cell: {cell}")
+#
+
+#set_obj.verbosity=10
 set_obj.export_to_xml()
 
 # Run!