Adding tokamak model

.github/workflows/simple_tokamak_cad_to_dagmc.yml

@@ -0,0 +1,30 @@
+name: simple tokamak - cad to dagmc
+
+on:
+  pull_request:
+    branches:
+      - develop
+      - main
+  push:
+    branches:
+      - main
+
+jobs:
+  testing:
+    runs-on: ubuntu-latest
+    steps:
+
+      - name: checkout actions
+        uses: actions/checkout@v4
+
+      - name: use upstream test composite action
+        uses: ./.github/actions/dependencies_cad_to_dagmc
+
+      - shell: bash
+        env:
+          OPENMC_CROSS_SECTIONS: /home/runner/work/model_benchmark_zoo/model_benchmark_zoo/cross_sections.xml
+        run: |
+          source "${HOME}/conda/etc/profile.d/conda.sh"
+          source "${HOME}/conda/etc/profile.d/mamba.sh"
+          mamba activate
+          pytest tests/test_cad_to_dagmc/test_csg_cad_simple_tokamak.py

.github/workflows/simple_tokamak_cad_to_openmc.yml

@@ -0,0 +1,30 @@
+name: simple tokamak - cad to openmc
+
+on:
+  pull_request:
+    branches:
+      - develop
+      - main
+  push:
+    branches:
+      - main
+
+jobs:
+  testing:
+    runs-on: ubuntu-latest
+    steps:
+
+      - name: checkout actions
+        uses: actions/checkout@v4
+
+      - name: use upstream test composite action
+        uses: ./.github/actions/dependencies_cad_to_openmc
+
+      - shell: bash
+        env:
+          OPENMC_CROSS_SECTIONS: /home/runner/work/model_benchmark_zoo/model_benchmark_zoo/cross_sections.xml
+        run: |
+          source "${HOME}/conda/etc/profile.d/conda.sh"
+          source "${HOME}/conda/etc/profile.d/mamba.sh"
+          mamba activate
+          pytest tests/test_cad_to_openmc/test_csg_cad_simple_tokamak.py

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
 
@@ -41,7 +42,7 @@ mamba activate new_env
 Install the dependencies, if this fails to solve the environment you could also try [installing OpenMC from source](https://docs.openmc.org/en/stable/quickinstall.html) which might be prefered.
 
 ```bash
-mamba install -y -c conda-forge moab>=5.3.0 gmsh python-gmsh "openmc=0.14.0=dagmc*nompi*"
+mamba install -y -c conda-forge gmsh python-gmsh "openmc=0.14.0=dagmc*nompi*"
 ```
 
 CadQuery should then be installed, here is the mamba command and the pip command.

examples/compare_csg_cad_ellipticaltorus.py

@@ -27,7 +27,7 @@
 my_settings.particles = 500
 my_settings.run_mode = 'fixed source'
 
-# Create a DT point source
+# Create a DT ring source
 my_source = openmc.IndependentSource()
 r = openmc.stats.Discrete([major_radius], [1])
 phi = openmc.stats.Uniform(0, 2*np.pi)

examples/compare_csg_cad_simpletokamak.py

@@ -0,0 +1,96 @@
+from model_benchmark_zoo import SimpleTokamak
+import openmc
+import math
+
+mat1 = openmc.Material(name='1')
+mat1.add_nuclide('Fe56', 1)
+mat1.set_density('g/cm3', 1)
+
+mat2 = openmc.Material(name='2')
+mat2.add_nuclide('Be9', 1)
+mat2.set_density('g/cm3', 1)
+
+# geometry used in both simulations
+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("simpletokamak.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 ring 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)
+z = openmc.stats.Discrete([0], [1])
+my_source.space = openmc.stats.CylindricalIndependent(r, phi, z)
+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(materials=[mat1, mat2])
+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_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(
+    h5m_filename='simpletokamak.h5m', materials=[mat1, mat2])
+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_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)
+
+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/__init__.py

@@ -6,4 +6,5 @@
 from .cylinder import *
 from .circulartorus import *
 from .nestedcylinder import *
-from .ellipticaltorus import *
+from .ellipticaltorus import *
+from .simpletokamak import *

src/model_benchmark_zoo/simpletokamak.py

@@ -0,0 +1,62 @@
+
+from .utils import BaseCommonGeometryObject
+
+class SimpleTokamak(BaseCommonGeometryObject):
+    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
+
+        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 = -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) # plasma
+        cell2 = openmc.Cell(region=region2) # blanket
+        cell2.fill = materials[0]
+        cell3 = openmc.Cell(region=region3) # center column
+        cell3.fill = materials[1]
+        cell4 = openmc.Cell(region=region4) # outer vessel
+
+        geometry = openmc.Geometry([cell1, cell2, cell3, cell4])
+        materials = openmc.Materials([materials[0], materials[1]])
+        model = openmc.Model(geometry=geometry, materials=materials)
+        return model
+
+    def cadquery_assembly(self):
+        import cadquery as cq
+
+        assembly = cq.Assembly(name="simpletokamak")
+
+        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)
+
+        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)
+        return assembly

tests/test_cad_to_dagmc/test_csg_cad_simple_tokamak.py

@@ -0,0 +1,81 @@
+from model_benchmark_zoo import SimpleTokamak
+import openmc
+import math
+
+def test_compare():
+    # single material 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('Be9', 1)
+    mat2.set_density('g/cm3', 1)
+
+    # geometry used in both simulations
+    common_geometry_object = SimpleTokamak(
+        radius=500,
+        blanket_thicknesses=100,
+        center_column_thicknesses=50
+    )
+
+    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 ring 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)
+    z = openmc.stats.Discrete([0], [1])
+    my_source.space = openmc.stats.CylindricalIndependent(r, phi, z)
+    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(materials=[mat1, mat2])
+    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_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")
+
+    common_geometry_object.export_h5m_file_with_cad_to_dagmc(
+        h5m_filename='simpletokamak.h5m',
+        material_tags=['1', '2'],
+        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(h5m_filename='simpletokamak.h5m', materials=[mat1, mat2])
+    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_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, rel_tol=0.01)
+    assert math.isclose(cad_result_mat_2.mean, csg_result_mat_2.mean, rel_tol=0.01)