Open shimwell opened 2 years ago
Test the same 3d model made with 3 different methods
Perhaps this code can be used
# # move to neutronics_workflow # class TestSimulationResultsVsCsg(unittest.TestCase): # """Makes a geometry in the paramak and in CSG geometry, simulates and # compares the results""" # def simulate_cylinder_cask_csg( # self, material, source, height, outer_radius, thickness, batches, particles # ): # """Makes a CSG cask geometry runs a simulation and returns the result""" # mats = openmc.Materials([material]) # outer_cylinder = openmc.ZCylinder(r=outer_radius) # inner_cylinder = openmc.ZCylinder(r=outer_radius - thickness) # inner_top = openmc.ZPlane(z0=height * 0.5) # inner_bottom = openmc.ZPlane(z0=-height * 0.5) # outer_top = openmc.ZPlane(z0=(height * 0.5) + thickness) # outer_bottom = openmc.ZPlane(z0=(-height * 0.5) - thickness) # sphere_1 = openmc.Sphere(r=100, boundary_type="vacuum") # cylinder_region = -outer_cylinder & +inner_cylinder & -inner_top & +inner_bottom # cylinder_cell = openmc.Cell(region=cylinder_region) # cylinder_cell.fill = material # top_cap_region = -outer_top & +inner_top & -outer_cylinder # top_cap_cell = openmc.Cell(region=top_cap_region) # top_cap_cell.fill = material # bottom_cap_region = +outer_bottom & -inner_bottom & -outer_cylinder # bottom_cap_cell = openmc.Cell(region=bottom_cap_region) # bottom_cap_cell.fill = material # inner_void_region = -inner_cylinder & -inner_top & +inner_bottom # inner_void_cell = openmc.Cell(region=inner_void_region) # # sphere 1 region is below -sphere_1 and not (~) in the other regions # sphere_1_region = -sphere_1 # sphere_1_cell = openmc.Cell( # region=sphere_1_region # & ~bottom_cap_region # & ~top_cap_region # & ~cylinder_region # & ~inner_void_region # ) # universe = openmc.Universe( # cells=[ # inner_void_cell, # cylinder_cell, # top_cap_cell, # bottom_cap_cell, # sphere_1_cell, # ] # ) # geom = openmc.Geometry(universe) # # Instantiate a Settings object # sett = openmc.Settings() # sett.batches = batches # sett.particles = particles # sett.inactive = 0 # sett.run_mode = "fixed source" # sett.photon_transport = True # sett.source = source # cell_filter = openmc.CellFilter([cylinder_cell, top_cap_cell, bottom_cap_cell]) # tally = openmc.Tally(name="csg_heating") # tally.filters = [cell_filter] # tally.scores = ["heating"] # tallies = openmc.Tallies() # tallies.append(tally) # model = openmc.model.Model(geom, mats, sett, tallies) # sp_filename = model.run() # # open the results file # results = openmc.StatePoint(sp_filename) # # access the tally using pandas dataframes # tally = results.get_tally(name="csg_heating") # tally_df = tally.get_pandas_dataframe() # return tally_df["mean"].sum() # def simulate_cylinder_cask_cad( # self, material, source, height, outer_radius, thickness, batches, particles # ): # """Makes a CAD cask geometry runs a simulation and returns the result""" # top_cap_cell = paramak.RotateStraightShape( # stp_filename="top_cap_cell.stp", # material_tag="test_mat", # points=[ # (outer_radius, height * 0.5), # (outer_radius, (height * 0.5) + thickness), # (0, (height * 0.5) + thickness), # (0, height * 0.5), # ], # ) # bottom_cap_cell = paramak.RotateStraightShape( # stp_filename="bottom_cap_cell.stp", # material_tag="test_mat", # points=[ # (outer_radius, -height * 0.5), # (outer_radius, (-height * 0.5) - thickness), # (0, (-height * 0.5) - thickness), # (0, -height * 0.5), # ], # ) # cylinder_cell = paramak.CenterColumnShieldCylinder( # height=height, # inner_radius=outer_radius - thickness, # outer_radius=outer_radius, # material_tag="test_mat", # ) # my_geometry = paramak.Reactor( # shapes_and_components=[cylinder_cell, bottom_cap_cell, top_cap_cell], # method="pymoab", # ) # my_model = openmc_dagmc_wrapper.NeutronicsModel( # h5m_filename=my_geometry.export_h5m(), # source=source, # simulation_batches=batches, # simulation_particles_per_batch=particles, # materials={"test_mat": material}, # cell_tallies=["heating"], # ) # my_model.simulate() # # scaled from MeV to eV # return ( # my_model.results["test_mat_heating"]["MeV per source particle"]["result"] # * 1e6 # ) # def test_cylinder_cask(self): # """Runs the same source and material with CAD and CSG geoemtry""" # height = 100 # outer_radius = 50 # thickness = 10 # batches = 10 # particles = 500 # test_material = openmc.Material(name="test_material") # test_material.set_density("g/cm3", 7.75) # test_material.add_element("Fe", 0.95, percent_type="wo") # test_material.add_element("C", 0.05, percent_type="wo") # source = openmc.Source() # source.space = openmc.stats.Point((0, 0, 0)) # source.angle = openmc.stats.Isotropic() # source.energy = openmc.stats.Discrete([14e6], [1.0]) # csg_result = self.simulate_cylinder_cask_csg( # test_material, source, height, outer_radius, thickness, batches, particles # ) # cad_result = self.simulate_cylinder_cask_cad( # test_material, source, height, outer_radius, thickness, batches, particles # ) # assert pytest.approx(csg_result, rel=0.02) == cad_result
Test the same 3d model made with 3 different methods
Perhaps this code can be used