ci working locally for cad-to-dagmc

.github/workflows/simple_tokamak_cad_to_dagmc.yml

@@ -1,4 +1,4 @@
-name: simple tokamak - cad to openmc
+name: simple tokamak - cad to dagmc
 
 on:
   pull_request:

.github/workflows/simple_tokamak_cad_to_openmc.yml

@@ -1,4 +1,4 @@
-name: simple tokamak - cad to dagmc
+name: simple tokamak - cad to openmc
 
 on:
   pull_request:

README.md

@@ -12,6 +12,7 @@ models (CSG) for comparing neutronics simulations with both geometry types.
 |         <p align="center"><img src="https://raw.githubusercontent.com/fusion-energy/model_benchmark_zoo/main/examples/nestedcylinder.png" width="100"></p>  |     [![nested_cylinder](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/nested_cylinder_cad_to_dagmc.yml/badge.svg)](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/nested_cylinder_cad_to_dagmc.yml) <br> [![nested cylinder - cad to openmc](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/nested_cylinder_cad_to_openmc.yml/badge.svg)](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/nested_cylinder_cad_to_openmc.yml) |
 |         <p align="center"><img src="https://raw.githubusercontent.com/fusion-energy/model_benchmark_zoo/main/examples/circulartorus.png" width="100"></p>  |     [![circular_torus](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/circular_torus_cad_to_dagmc.yml/badge.svg)](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/circular_torus_cad_to_dagmc.yml) <br> [![circular torus - cad to openmc](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/circular_torus_cad_to_openmc.yml/badge.svg)](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/circular_torus_cad_to_openmc.yml) |
 |         <p align="center"><img src="https://raw.githubusercontent.com/fusion-energy/model_benchmark_zoo/main/examples/ellipticaltorus.png" width="100"></p>  |     [![elliptical_torus](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/elliptical_torus_cad_to_dagmc.yml/badge.svg)](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/elliptical_torus_cad_to_dagmc.yml) <br> [![elliptical torus - cad to openmc](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/elliptical_torus_cad_to_openmc.yml/badge.svg)](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/elliptical_torus_cad_to_openmc.yml) |
+|         <p align="center"><img src="https://raw.githubusercontent.com/fusion-energy/model_benchmark_zoo/main/examples/simple_tokamak.png" width="100"></p>  |     [![simple_tokamak](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/simple_tokamak_cad_to_dagmc.yml/badge.svg)](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/simple_tokamak_cad_to_dagmc.yml) <br> [![elliptical torus - cad to openmc](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/simple_tokamak_cad_to_openmc.yml/badge.svg)](https://github.com/fusion-energy/model_benchmark_zoo/actions/workflows/simple_tokamak_cad_to_openmc.yml) |
 
 # Installation prerequisite
 

examples/compare_csg_cad_simpletokamak.py

@@ -2,7 +2,6 @@
 import openmc
 import math
 
-# single material used in both simulations
 mat1 = openmc.Material(name='1')
 mat1.add_nuclide('Fe56', 1)
 mat1.set_density('g/cm3', 1)
@@ -11,12 +10,15 @@
 mat2.add_nuclide('Be9', 1)
 mat2.set_density('g/cm3', 1)
 
-my_materials = openmc.Materials([mat1, mat2])
-
 # geometry used in both simulations
-common_geometry_object = SimpleTokamak(materials=my_materials, radius=500, blanket_thicknesses=100, center_column_thicknesses=50)
+common_geometry_object = SimpleTokamak(
+    radius=500,
+    blanket_thicknesses=100,
+    center_column_thicknesses=150,
+    center_column_extent_beyond_blanket=20,
+)
 # just writing a CAD step file for visulisation
-common_geometry_object.export_stp_file("nestedsphere.stp")
+common_geometry_object.export_stp_file("simpletokamak.stp")
 
 mat1_filter = openmc.MaterialFilter(mat1)
 tally1 = openmc.Tally(name='mat1_flux_tally')
@@ -37,7 +39,7 @@
 my_settings.run_mode = 'fixed source'
 
 # Create a DT ring source
-my_source = openmc.Source()
+my_source = openmc.IndependentSource()
 source_r = common_geometry_object.center_column_thicknesses + (common_geometry_object.radius-common_geometry_object.center_column_thicknesses) /2
 r = openmc.stats.Discrete([source_r], [1])
 phi = openmc.stats.Uniform(0, 2*math.pi)
@@ -47,7 +49,7 @@
 my_settings.source = my_source
 
 # making openmc.Model with CSG geometry
-csg_model = common_geometry_object.csg_model()
+csg_model = common_geometry_object.csg_model(materials=[mat1, mat2])
 csg_model.tallies = my_tallies
 csg_model.settings = my_settings
 
@@ -57,11 +59,20 @@
 with openmc.StatePoint(output_file_from_csg) as sp_from_csg:
     csg_result_mat_1 = sp_from_csg.get_tally(name="mat1_flux_tally")
     csg_result_mat_2 = sp_from_csg.get_tally(name="mat2_flux_tally")
+
 csg_result_mat_1_str = f'CSG tally mean {csg_result_mat_1.mean} std dev {csg_result_mat_1.std_dev}'
 csg_result_mat_2_str = f'CSG tally mean {csg_result_mat_2.mean} std dev {csg_result_mat_2.std_dev}'
 
+common_geometry_object.export_h5m_file_with_cad_to_dagmc(
+    h5m_filename='simpletokamak.h5m',
+    material_tags=['1', '2'],
+    # the small mesh sizes make a large detailed mesh which is needed to get similar answers
+    min_mesh_size=0.01,
+    max_mesh_size=0.5
+)
 # making openmc.Model with DAGMC geometry and specifying mesh sizes to get a good representation of a sphere
-dag_model = common_geometry_object.dagmc_model(min_mesh_size=1., max_mesh_size=50.)
+dag_model = common_geometry_object.dagmc_model(
+    h5m_filename='simpletokamak.h5m', materials=[mat1, mat2])
 dag_model.tallies = my_tallies
 dag_model.settings = my_settings
 
@@ -71,11 +82,15 @@
 with openmc.StatePoint(output_file_from_cad) as sp_from_cad:
     cad_result_mat_1 = sp_from_cad.get_tally(name="mat1_flux_tally")
     cad_result_mat_2 = sp_from_cad.get_tally(name="mat2_flux_tally")
+
 cad_result_mat_1_str = f'CAD tally mean {cad_result_mat_1.mean} std dev {cad_result_mat_1.std_dev}'
 cad_result_mat_2_str = f'CAD tally mean {cad_result_mat_2.mean} std dev {cad_result_mat_2.std_dev}'
 
 # printing both tally results
 print(csg_result_mat_1_str)
 print(cad_result_mat_1_str)
 print(csg_result_mat_2_str)
-print(cad_result_mat_2_str)
+print(cad_result_mat_2_str)
+
+assert math.isclose(cad_result_mat_1.mean, csg_result_mat_1.mean, rel_tol=0.01)
+assert math.isclose(cad_result_mat_2.mean, csg_result_mat_2.mean, rel_tol=0.01)

src/model_benchmark_zoo/simpletokamak.py

@@ -2,29 +2,43 @@
 from .utils import BaseCommonGeometryObject
 
 class SimpleTokamak(BaseCommonGeometryObject):
-    def __init__(self, radius=500, blanket_thicknesses=100, center_column_thicknesses=50):
+    def __init__(
+            self,
+            radius=500,
+            blanket_thicknesses=100,
+            center_column_thicknesses=150,
+            center_column_extent_beyond_blanket=10,
+        ):
         self.radius = radius
         self.blanket_thicknesses = blanket_thicknesses
         self.center_column_thicknesses = center_column_thicknesses
+        self.center_column_extent_beyond_blanket = center_column_extent_beyond_blanket
 
     def csg_model(self, materials):
         import openmc
 
-        surface1 = openmc.Sphere(r=self.radius)
-        surface2 = openmc.Sphere(r=self.radius + self.blanket_thicknesses, boundary_type="vacuum")
-        surface3 = openmc.ZCylinder(r=self.center_column_thicknesses)
+        center_column_height = (self.radius + self.blanket_thicknesses + self.center_column_extent_beyond_blanket)*2
+
+        surface_inner_wall = openmc.Sphere(r=self.radius)
+        surface_outer_wall = openmc.Sphere(r=self.radius + self.blanket_thicknesses)
+        surface_center_cylinder = openmc.ZCylinder(r=self.center_column_thicknesses)
+        surface_top_cy = openmc.ZPlane(z0=center_column_height/2, boundary_type='vacuum')
+        surface_bot_cy = openmc.ZPlane(z0=-center_column_height/2, boundary_type='vacuum')
+        outer_surface = openmc.ZCylinder(r=self.radius + self.blanket_thicknesses, boundary_type='vacuum')
 
-        region1 = -surface1 & +surface3  # plasma
-        region2 = +surface1 & -surface2  # blanket
-        region3 = -surface2 & -surface3  # center column
+        region1 = -surface_inner_wall & +surface_center_cylinder  # plasma
+        region2 = +surface_inner_wall & -surface_outer_wall & +surface_center_cylinder # blanket
+        region3 = -surface_top_cy & +surface_bot_cy & -surface_center_cylinder  # center column
+        region4 = +surface_outer_wall & -surface_top_cy & +surface_bot_cy & +surface_center_cylinder & -outer_surface  # outer vessel
 
-        cell1 = openmc.Cell(region=region1)
-        cell2 = openmc.Cell(region=region2)
+        cell1 = openmc.Cell(region=region1) # plasma
+        cell2 = openmc.Cell(region=region2) # blanket
         cell2.fill = materials[0]
-        cell3 = openmc.Cell(region=region3)
+        cell3 = openmc.Cell(region=region3) # center column
         cell3.fill = materials[1]
+        cell4 = openmc.Cell(region=region4) # outer vessel
 
-        geometry = openmc.Geometry([cell1, cell2, cell3])
+        geometry = openmc.Geometry([cell1, cell2, cell3, cell4])
         materials = openmc.Materials([materials[0], materials[1]])
         model = openmc.Model(geometry=geometry, materials=materials)
         return model
@@ -34,12 +48,16 @@ def cadquery_assembly(self):
 
         assembly = cq.Assembly(name="simpletokamak")
 
-        sphere_envelope = cq.Workplane().sphere(self.radius + self.blanket_thicknesses)
+        # sphere_envelope = cq.Workplane().sphere(self.radius + self.blanket_thicknesses)
 
-        center_column = cq.Workplane("XY").circle(self.radius).extrude(self.radius+self.blanket_thicknesses).intersect(sphere_envelope)
+        center_column_height = self.radius + self.blanket_thicknesses + self.center_column_extent_beyond_blanket
 
         sphere1 = cq.Workplane().sphere(self.radius)
-        sphere2 = cq.Workplane().sphere(self.radius + self.blanket_thicknesses).cut(sphere1).cut(center_column)
+        sphere2 = cq.Workplane().sphere(self.radius + self.blanket_thicknesses)
+
+        center_column = cq.Workplane("XY").circle(self.center_column_thicknesses).extrude(center_column_height, both=True)
+
+        sphere2 = sphere2.cut(sphere1).cut(center_column)
 
         assembly.add(sphere2)
         assembly.add(center_column)

tests/test_cad_to_dagmc/test_csg_cad_simple_tokamak.py

@@ -62,11 +62,11 @@ def test_compare():
     common_geometry_object.export_h5m_file_with_cad_to_dagmc(
         h5m_filename='simpletokamak.h5m',
         material_tags=['1', '2'],
-        min_mesh_size=0.01,
-        max_mesh_size=0.5
+        max_mesh_size=5,
+        min_mesh_size=0.05
     )
     # making openmc.Model with DAGMC geometry and specifying mesh sizes to get a good representation of a sphere
-    dag_model = common_geometry_object.dagmc_model(materials=[mat1, mat2])
+    dag_model = common_geometry_object.dagmc_model(h5m_filename='simpletokamak.h5m', materials=[mat1, mat2])
     dag_model.tallies = my_tallies
     dag_model.settings = my_settings
 
@@ -77,5 +77,5 @@ def test_compare():
         cad_result_mat_1 = sp_from_cad.get_tally(name="mat1_flux_tally")
         cad_result_mat_2 = sp_from_cad.get_tally(name="mat2_flux_tally")
 
-    assert math.isclose(cad_result_mat_1.mean, csg_result_mat_1.mean)
-    assert math.isclose(cad_result_mat_2.mean, csg_result_mat_2.mean)
+    assert math.isclose(cad_result_mat_1.mean, csg_result_mat_1.mean, rel_tol=0.01)
+    assert math.isclose(cad_result_mat_2.mean, csg_result_mat_2.mean, rel_tol=0.01)

tests/test_cad_to_openmc/test_csg_cad_simple_tokamak.py

@@ -64,7 +64,7 @@ def test_compare():
         material_tags=['1', '2'],
     )
     # making openmc.Model with DAGMC geometry and specifying mesh sizes to get a good representation of a sphere
-    dag_model = common_geometry_object.dagmc_model(materials=[mat1, mat2])
+    dag_model = common_geometry_object.dagmc_model(h5m_filename='simpletokamak.h5m', materials=[mat1, mat2])
     dag_model.tallies = my_tallies
     dag_model.settings = my_settings