Merge pull request #5 from fusion-energy/adding_two_contacting_cubes

Adding two contacting cubes

examples/compare_csg_cad_two_touching_cuboids.py

@@ -0,0 +1,77 @@
+from model_benchmark_zoo import TwoTouchingCuboids
+import openmc
+
+# materials used in both simulations
+mat1 = openmc.Material(name='1')
+mat1.add_nuclide('Fe56', 1)
+mat1.set_density('g/cm3', 1)
+mat2 = openmc.Material(name='2')
+mat2.add_nuclide('Fe56', 1)
+mat2.set_density('g/cm3', 1)
+my_materials = openmc.Materials([mat1, mat2])
+
+# geometry used in both simulations
+common_geometry_object = TwoTouchingCuboids(
+    materials=my_materials, width1=10, width2=4)
+# just writing a CAD step file for visulisation
+common_geometry_object.export_stp_file("TwoTouchingCuboids.stp")
+
+mat1_filter = openmc.MaterialFilter(mat1)
+tally1 = openmc.Tally(name='mat1_flux_tally')
+tally1.filters = [mat1_filter]
+tally1.scores = ['flux']
+
+mat2_filter = openmc.MaterialFilter(mat2)
+tally2 = openmc.Tally(name='mat2_flux_tally')
+tally2.filters = [mat2_filter]
+tally2.scores = ['flux']
+my_tallies = openmc.Tallies([tally1, tally2])
+
+my_settings = openmc.Settings()
+my_settings.batches = 10
+my_settings.inactive = 0
+my_settings.particles = 500
+my_settings.run_mode = 'fixed source'
+
+# Create a DT point source
+my_source = openmc.Source()
+my_source.space = openmc.stats.Point((0, 0, 0))
+my_source.angle = openmc.stats.Isotropic()
+my_source.energy = openmc.stats.Discrete([14e6], [1])
+my_settings.source = my_source
+
+# making openmc.Model with CSG geometry
+csg_model = common_geometry_object.csg_model()
+csg_model.materials = my_materials
+csg_model.tallies = my_tallies
+csg_model.settings = my_settings
+
+output_file_from_csg = csg_model.run()
+
+# extracting the tally result from the CSG simulation
+with openmc.StatePoint(output_file_from_csg) as sp_from_csg:
+    csg_result1 = sp_from_csg.get_tally(name="mat1_flux_tally")
+    csg_result2 = sp_from_csg.get_tally(name="mat2_flux_tally")
+csg_result_1 = f'CSG tally mat 1 mean {csg_result1.mean} std dev {csg_result1.std_dev}'
+csg_result_2 = f'CSG tally mat 2 mean {csg_result2.mean} std dev {csg_result2.std_dev}'
+
+# making openmc.Model with DAGMC geometry and specifying mesh sizes to get a good representation of a TwoTouchingCuboids
+dag_model = common_geometry_object.dagmc_model(min_mesh_size=0.01, max_mesh_size=0.5)
+dag_model.materials = my_materials
+dag_model.tallies = my_tallies
+dag_model.settings = my_settings
+
+output_file_from_cad = dag_model.run()
+
+# extracting the tally result from the DAGMC simulation
+with openmc.StatePoint(output_file_from_cad) as sp_from_cad:
+    cad_result1 = sp_from_cad.get_tally(name="mat1_flux_tally")
+    cad_result2 = sp_from_cad.get_tally(name="mat2_flux_tally")
+cad_result_1 = f'CAD tally mat 1 mean {cad_result1.mean} std dev {cad_result1.std_dev}'
+cad_result_2 = f'CAD tally mat 2 mean {cad_result2.mean} std dev {cad_result2.std_dev}'
+
+# printing both tally results
+print(csg_result_1)
+print(cad_result_1)
+print(csg_result_2)
+print(cad_result_2)

pyproject.toml

@@ -0,0 +1,38 @@
+[build-system]
+requires = ["setuptools >= 65.4.0", "setuptools_scm[toml]>=7.0.5"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "model_benchmark_zoo"
+authors = [
+  { name="Jonathan Shimwell", email="[email protected]" },
+]
+license = {file = "LICENSE.txt"}
+description = "Collection of geometries for testing neutronics simulations in CSG and CAD format"
+readme = "README.md"
+requires-python = ">=3.8"
+keywords = ["dagmc", "geometry", "test", "benchmark"]
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "License :: OSI Approved :: MIT License",
+    "Operating System :: OS Independent",
+]
+dependencies = []
+dynamic = ["version"]
+
+
+[tool.setuptools_scm]
+write_to = "src/_version.py"
+
+
+[project.optional-dependencies]
+tests = [
+    "pytest"
+]
+
+[project.urls]
+"Homepage" = "https://github.com/fusion-energy/model_benchmark_zoo"
+"Bug Tracker" = "https://github.com/fusion-energy/model_benchmark_zoo/issues"
+
+[tool.setuptools]
+package-dir = {"" = "src"}

src/model_benchmark_zoo/__init__.py

@@ -1,4 +1,5 @@
 from .cuboid import *
 from .sphere import *
 from .cylinder import *
+from .two_touching_cuboids import *
 from .sphericalshell import *

src/model_benchmark_zoo/two_touching_cuboids.py

@@ -0,0 +1,67 @@
+import cadquery as cq
+import openmc
+from cad_to_dagmc import CadToDagmc
+
+#     *----------*
+#     |          |
+#     |          |----*
+#     |          |    |
+#     |          |    |
+#     |          |----*
+#     |          |
+#     *----------*
+
+class TwoTouchingCuboids:
+    def __init__(self, materials, width1=10, width2=4):
+        self.width1 = width1
+        self.width2 = width2
+        self.materials = materials
+
+    def csg_model(self):
+        surface1 = openmc.ZPlane(z0=self.width1*0.5, boundary_type="vacuum")
+        surface2 = openmc.ZPlane(z0=self.width1*-0.5, boundary_type="vacuum")
+        surface3 = openmc.XPlane(x0=self.width1*0.5, boundary_type="vacuum")
+        surface4 = openmc.XPlane(x0=self.width1*-0.5, boundary_type="vacuum")
+        surface5 = openmc.YPlane(y0=self.width1*0.5)
+        surface6 = openmc.YPlane(y0=self.width1*-0.5, boundary_type="vacuum")
+        surface7 = openmc.YPlane(y0=self.width1*0.5+self.width2, boundary_type="vacuum")
+
+        region1 = -surface1 & +surface2 & -surface3 & +surface4 & -surface5 & +surface6
+        region2 = -surface1 & +surface2 & -surface3 & +surface4 & -surface7 & +surface5
+
+        cell1 = openmc.Cell(region=region1)
+        cell1.fill = self.materials[0]
+
+        cell2 = openmc.Cell(region=region2)
+        cell2.fill = self.materials[1]
+
+        geometry = openmc.Geometry([cell1, cell2])
+        model = openmc.Model(geometry=geometry)
+        return model
+
+    def cadquery_assembly(self):
+        assembly = cq.Assembly(name="TwoTouchingCuboids")
+        cuboid1 = cq.Workplane().box(self.width1, self.width1, self.width1)
+        assembly.add(cuboid1)
+        cuboid2 = cq.Workplane().moveTo(0.5*self.width1+ 0.5*self.width2).box(self.width2, self.width2, self.width2)
+        assembly.add(cuboid2)
+        return assembly
+
+    def export_stp_file(self, filename="TwoTouchingCuboids.step"):
+        self.cadquery_assembly().save(filename, "STEP")
+
+    def dagmc_model(self, filename="TwoTouchingCuboids.h5m", min_mesh_size=0.1, max_mesh_size=100.0):
+        assembly = self.cadquery_assembly()
+        ctd = CadToDagmc()
+        material_tags = [self.materials[0].name, self.materials[1].name]
+        ctd.add_cadquery_object(assembly, material_tags=material_tags)
+        ctd.export_dagmc_h5m_file(
+            filename=filename,
+            min_mesh_size=min_mesh_size,
+            max_mesh_size=max_mesh_size,
+            msh_filename='TwoTouchingCuboids.msh'
+        )
+        universe = openmc.DAGMCUniverse(filename).bounded_universe()
+        geometry = openmc.Geometry(universe)
+        model = openmc.Model(geometry=geometry)
+        return model