Skip to content

Commit

Permalink
Merge pull request #36 from fusion-energy/adding_toka_2
Browse files Browse the repository at this point in the history
Adding a simple tokamak model
  • Loading branch information
shimwell committed Mar 21, 2024
2 parents aea8690 + b21a7a5 commit 580db0c
Show file tree
Hide file tree
Showing 10 changed files with 383 additions and 3 deletions.
30 changes: 30 additions & 0 deletions .github/workflows/simple_tokamak_cad_to_dagmc.yml
Original file line number Diff line number Diff line change
@@ -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
30 changes: 30 additions & 0 deletions .github/workflows/simple_tokamak_cad_to_openmc.yml
Original file line number Diff line number Diff line change
@@ -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
3 changes: 2 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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.
Expand Down
2 changes: 1 addition & 1 deletion examples/compare_csg_cad_ellipticaltorus.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
96 changes: 96 additions & 0 deletions examples/compare_csg_cad_simpletokamak.py
Original file line number Diff line number Diff line change
@@ -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)
Binary file added examples/simpletokamak.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
3 changes: 2 additions & 1 deletion src/model_benchmark_zoo/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,4 +6,5 @@
from .cylinder import *
from .circulartorus import *
from .nestedcylinder import *
from .ellipticaltorus import *
from .ellipticaltorus import *
from .simpletokamak import *
62 changes: 62 additions & 0 deletions src/model_benchmark_zoo/simpletokamak.py
Original file line number Diff line number Diff line change
@@ -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
81 changes: 81 additions & 0 deletions tests/test_cad_to_dagmc/test_csg_cad_simple_tokamak.py
Original file line number Diff line number Diff line change
@@ -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)
Loading

0 comments on commit 580db0c

Please sign in to comment.