From 3181d7ece7a38b518c05d5bc28613bfe8e677cec Mon Sep 17 00:00:00 2001 From: Jon Shimwell Date: Thu, 26 Sep 2024 23:29:39 +0100 Subject: [PATCH] complete project rewrite --- .github/workflows/ci.yml | 55 + .github/workflows/documentation_update.yml | 63 + .gitignore | 177 +++ README.md | 16 + docs/Makefile | 20 + docs/README.md | 11 + docs/conf.py | 70 ++ docs/index.rst | 19 + docs/install.rst | 101 ++ docs/make.bat | 35 + docs/python_api.rst | 30 + docs/usage.rst | 693 +++++++++++ docs/version_switcher.json | 17 + environment_dev.yml | 20 + examples/example_util_functions.py | 102 ++ examples/plasma.py | 12 + .../spherical_tokamak_from_plasma_minimal.py | 48 + ...rical_tokamak_from_plasma_with_divertor.py | 49 + ...cal_tokamak_from_plasma_with_pf_magnets.py | 64 + ...om_plasma_with_pf_magnets_and_divertors.py | 69 ++ ...cal_tokamak_from_plasma_with_tf_magnets.py | 57 + examples/spherical_tokamak_minimal.py | 57 + examples/tokamak_from_plasma_minimal.py | 40 + examples/tokamak_from_plasma_with_divertor.py | 43 + examples/tokamak_minimal.py | 54 + examples/tokamak_with_pf_magnets.py | 75 ++ pyproject.toml | 55 + src/paramak/__init__.py | 23 + src/paramak/assemblies/__init__.py | 0 src/paramak/assemblies/spherical_tokamak.py | 261 ++++ src/paramak/assemblies/tokamak.py | 288 +++++ src/paramak/utils.py | 336 ++++++ src/paramak/workplanes/__init__.py | 0 .../blanket_constant_thickness_arc_h.py | 52 + src/paramak/workplanes/blanket_from_plasma.py | 303 +++++ .../center_column_shield_cylinder.py | 60 + .../workplanes/constant_thickness_dome.py | 190 +++ src/paramak/workplanes/cutting_wedge.py | 44 + .../workplanes/dished_vacuum_vessel.py | 112 ++ src/paramak/workplanes/plasma_simplified.py | 64 + src/paramak/workplanes/poloidal_field_coil.py | 43 + .../workplanes/poloidal_field_coil_case.py | 66 ++ .../toroidal_field_coil_rectangle.py | 106 ++ src/paramak/workplanes/u_shaped_dome.py | 105 ++ tests/test_assemblies/__init__.py | 0 .../test_assemblies/test_spherical_tokamak.py | 102 ++ tests/test_assemblies/test_utils.py | 102 ++ tests/test_utils.py | 1047 +++++++++++++++++ tests/test_workplanes/__init__.py | 0 .../test_blanket_constant_thickness_arc_h.py | 37 + .../test_blanket_from_plasma.py | 154 +++ .../test_center_column_shield_cylinder.py | 29 + .../test_constant_thickness_dome.py | 9 + .../test_dished_vacuum_vessel.py | 30 + .../test_workplanes/test_plasma_simplified.py | 45 + .../test_poloidal_field_coil.py | 20 + .../test_poloidal_field_coil_case.py | 24 + .../test_toroidal_field_coil_rectangle.py | 24 + tests/test_workplanes/test_u_shaped_dome.py | 17 + tests/test_workplanes/utils.py | 102 ++ 60 files changed, 5847 insertions(+) create mode 100644 .github/workflows/ci.yml create mode 100644 .github/workflows/documentation_update.yml create mode 100644 .gitignore create mode 100644 README.md create mode 100644 docs/Makefile create mode 100644 docs/README.md create mode 100644 docs/conf.py create mode 100644 docs/index.rst create mode 100644 docs/install.rst create mode 100644 docs/make.bat create mode 100644 docs/python_api.rst create mode 100644 docs/usage.rst create mode 100644 docs/version_switcher.json create mode 100644 environment_dev.yml create mode 100644 examples/example_util_functions.py create mode 100644 examples/plasma.py create mode 100644 examples/spherical_tokamak_from_plasma_minimal.py create mode 100644 examples/spherical_tokamak_from_plasma_with_divertor.py create mode 100644 examples/spherical_tokamak_from_plasma_with_pf_magnets.py create mode 100644 examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py create mode 100644 examples/spherical_tokamak_from_plasma_with_tf_magnets.py create mode 100644 examples/spherical_tokamak_minimal.py create mode 100644 examples/tokamak_from_plasma_minimal.py create mode 100644 examples/tokamak_from_plasma_with_divertor.py create mode 100644 examples/tokamak_minimal.py create mode 100644 examples/tokamak_with_pf_magnets.py create mode 100644 pyproject.toml create mode 100644 src/paramak/__init__.py create mode 100644 src/paramak/assemblies/__init__.py create mode 100644 src/paramak/assemblies/spherical_tokamak.py create mode 100644 src/paramak/assemblies/tokamak.py create mode 100644 src/paramak/utils.py create mode 100644 src/paramak/workplanes/__init__.py create mode 100644 src/paramak/workplanes/blanket_constant_thickness_arc_h.py create mode 100644 src/paramak/workplanes/blanket_from_plasma.py create mode 100644 src/paramak/workplanes/center_column_shield_cylinder.py create mode 100644 src/paramak/workplanes/constant_thickness_dome.py create mode 100644 src/paramak/workplanes/cutting_wedge.py create mode 100644 src/paramak/workplanes/dished_vacuum_vessel.py create mode 100644 src/paramak/workplanes/plasma_simplified.py create mode 100644 src/paramak/workplanes/poloidal_field_coil.py create mode 100644 src/paramak/workplanes/poloidal_field_coil_case.py create mode 100644 src/paramak/workplanes/toroidal_field_coil_rectangle.py create mode 100644 src/paramak/workplanes/u_shaped_dome.py create mode 100644 tests/test_assemblies/__init__.py create mode 100644 tests/test_assemblies/test_spherical_tokamak.py create mode 100644 tests/test_assemblies/test_utils.py create mode 100644 tests/test_utils.py create mode 100644 tests/test_workplanes/__init__.py create mode 100644 tests/test_workplanes/test_blanket_constant_thickness_arc_h.py create mode 100644 tests/test_workplanes/test_blanket_from_plasma.py create mode 100644 tests/test_workplanes/test_center_column_shield_cylinder.py create mode 100644 tests/test_workplanes/test_constant_thickness_dome.py create mode 100644 tests/test_workplanes/test_dished_vacuum_vessel.py create mode 100644 tests/test_workplanes/test_plasma_simplified.py create mode 100644 tests/test_workplanes/test_poloidal_field_coil.py create mode 100644 tests/test_workplanes/test_poloidal_field_coil_case.py create mode 100644 tests/test_workplanes/test_toroidal_field_coil_rectangle.py create mode 100644 tests/test_workplanes/test_u_shaped_dome.py create mode 100644 tests/test_workplanes/utils.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 000000000..8e9a9e969 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,55 @@ +name: CI testing + +on: + pull_request: + branches: + - develop + - main + push: + branches: + - main + +jobs: + test: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.11"] + + steps: + - name: Install system packages + run: | + sudo apt-get update -y + sudo apt-get install -y libgl1-mesa-glx libgl1-mesa-dev libglu1-mesa-dev freeglut3-dev libosmesa6 libosmesa6-dev libgles2-mesa-dev libarchive-dev libpangocairo-1.0-0 + + - name: checkout actions + uses: actions/checkout@v4 + + - uses: mamba-org/setup-micromamba@v1 + with: + micromamba-version: '1.5.8-0' # any version from https://github.com/mamba-org/micromamba-releases + environment-file: environment_dev.yml + init-shell: bash + cache-environment: true + post-cleanup: 'all' + + - name: install dependencies run tests + shell: bash -el {0} + run: | + python -m pip install .[tests] + + - name: Test with pytest + shell: bash -el {0} + run: | + pytest -v + python examples/plasma.py + python examples/spherical_tokamak_from_plasma_minimal.py + python examples/spherical_tokamak_from_plasma_with_divertor.py + python examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py + python examples/spherical_tokamak_from_plasma_with_pf_magnets.py + python examples/spherical_tokamak_from_plasma_with_tf_magnets.py + python examples/spherical_tokamak_minimal.py + python examples/tokamak_from_plasma_minimal.py + python examples/tokamak_minimal.py + python examples/tokamak_from_plasma_with_divertor.py + diff --git a/.github/workflows/documentation_update.yml b/.github/workflows/documentation_update.yml new file mode 100644 index 000000000..0de8a67f3 --- /dev/null +++ b/.github/workflows/documentation_update.yml @@ -0,0 +1,63 @@ +name: documentation release + +on: + pull_request: + branches: + - main + push: + branches: + - main + tags: + - '*' + + +permissions: + contents: write + +jobs: + testing: + name: Documentation + runs-on: ubuntu-latest + steps: + - name: checkout actions + uses: actions/checkout@v4 + + - name: Install system packages + run: | + sudo apt-get update -y + sudo apt-get install -y libgl1-mesa-glx libgl1-mesa-dev libglu1-mesa-dev freeglut3-dev libosmesa6 libosmesa6-dev libgles2-mesa-dev libarchive-dev libpangocairo-1.0-0 + + + - uses: mamba-org/setup-micromamba@v1 + with: + micromamba-version: '1.5.8-0' # any version from https://github.com/mamba-org/micromamba-releases + environment-file: environment_dev.yml + init-shell: bash + cache-environment: true + post-cleanup: 'all' + + - name: install package + run: | + pip install --upgrade pip + pip install .[docs] + - name: Sphinx build tagged version + if: startsWith(github.ref, 'refs/tags/') + run: | + sphinx-build docs _build/${{ github.ref_name }} + sphinx-build docs _build + - name: Sphinx build dev version + if: (github.event_name == 'push' || github.event_name == 'pull_request') && !startsWith(github.ref, 'refs/tags/') + run: | + sphinx-build docs _build/dev + - name: Deploy to GitHub Pages + if: github.event_name == 'push' + uses: peaceiris/actions-gh-pages@v4 + with: + publish_branch: gh-pages + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: _build + # with next rlease of actions-gh-pages + # issue to allow force_orphan will be fixed + # https://github.com/peaceiris/actions-gh-pages/issues/455 + # force_orphan: true + keep_files: true diff --git a/.gitignore b/.gitignore new file mode 100644 index 000000000..02307b6f6 --- /dev/null +++ b/.gitignore @@ -0,0 +1,177 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +#pdm.lock +# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it +# in version control. +# https://pdm.fming.dev/latest/usage/project/#working-with-version-control +.pdm.toml +.pdm-python +.pdm-build/ + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ + +src/_version.py + +# openmc files +*.h5 +*.out +*.xml +*.h5m +*.h5 + +# cad files +*.step +*.stp +*.sql +*.html diff --git a/README.md b/README.md new file mode 100644 index 000000000..41cab5676 --- /dev/null +++ b/README.md @@ -0,0 +1,16 @@ + +[![N|Python](https://www.python.org/static/community_logos/python-powered-w-100x40.png)](https://www.python.org) + +[![CI testing](https://github.com/fusion-energy/paramak2/actions/workflows/ci.yml/badge.svg)](https://github.com/fusion-energy/paramak2/actions/workflows/ci.yml) + +[![documentation release](https://github.com/fusion-energy/paramak2/actions/workflows/documentation_update.yml/badge.svg)](https://github.com/fusion-energy/paramak2/actions/workflows/documentation_update.yml) + + +# Paramak + +Paramak python package allows rapid production of 3D CAD models and neutronics +models of fusion reactors. The purpose of Paramak is to provide geometry for +parametric studies. Paramak can create geometry in standard CAD formats such as +STP, STL, BRep, HTML and DAGMC h5m. + +:point_right: Please see the [Documentation](https://paramak.readthedocs.io) for installation, usage and API documentation. diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 000000000..d4bb2cbb9 --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 000000000..ee5d080aa --- /dev/null +++ b/docs/README.md @@ -0,0 +1,11 @@ +Documentation section + +To build the documentation locally + +``` +sphinx-build docs _build +``` +or +``` +make html +``` \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 000000000..aaf6b81ac --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,70 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information + +import os +import re +import sys + +sys.path.insert(0, os.path.abspath("../src")) + + +project = "Paramak" +copyright = "2024, J. Shimwell" +author = "J. Shimwell" + +import paramak + +version = paramak.__version__ +release = paramak.__version__ + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = [ + "sphinx.ext.autodoc", + "sphinx.ext.autosummary", + "sphinx_autodoc_typehints", + "sphinx.ext.coverage", + "sphinx.ext.napoleon", + "sphinx.ext.doctest", + "sphinx.ext.viewcode", + "sphinxcadquery.sphinxcadquery", +] + +templates_path = ["_templates"] +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] + + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +html_theme = "pydata_sphinx_theme" +html_static_path = ["_static"] + + +# TODO add logo +# html_favicon = "favicon.ico" + +# Version match must match the 'version' key in version_switcher.json +pattern = re.compile(r"^[0-9]+\.[0-9]+") +version_match = pattern.search(version) +if version_match: + version_match = version_match.group() +elif "dev" in version: + version_match = "dev" +else: + version_match = version + +html_theme_options = { + "github_url": "https://github.com/fusion-energy/paramak2", + "switcher": { + "json_url": "https://raw.githubusercontent.com/fusion-energy/paramak2/main/docs/version_switcher.json", + "version_match": version_match, + }, + "navbar_start": ["version-switcher", "navbar-icon-links"], +} diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 000000000..f8fa606e5 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,19 @@ + +Paramak documentation +===================== + +**Version**: |version| + +Parameter driven CAD creation for fusion reactors. + +Paramak provides parameter driven creation of Tokamak and Spherical Tokamak CAD models as well as DAGMC compatible neutronics models. + +The style of reaction, sizes of components, plasma shape and number of radial or vertical layers can be specified, + + +.. toctree:: + :maxdepth: 3 + + install + usage + python_api diff --git a/docs/install.rst b/docs/install.rst new file mode 100644 index 000000000..2614ffe0f --- /dev/null +++ b/docs/install.rst @@ -0,0 +1,101 @@ +Install +======= + +Prerequisites +------------- + +To use of Paramak you will need Python 3 installed using Miniconda or Anaconda, or Miniforge + +* `Miniforge `_ recommended as it includes Mamba +* `Miniconda `_ +* `Anaconda `_ + + + +Once you have a version of Mamba or Conda installed then proceed with the Paramak specific steps. + + +Install (mamba) +--------------- + +This is the recommended method as it installs all the dependencies and Mamba is faster and requires less RAM than the pure Conda method. + +Create a new environment (with your preferred python version). + +.. code-block:: bash + + mamba create --name paramak_env python=3.11 + + +Then activate the new environment. + +.. code-block:: bash + + mamba activate paramak_env + + +Then install the Paramak. + +.. code-block:: bash + + mamba install -c conda-forge paramak + +Now you should be ready to import paramak from your new python environment. + +Install (conda) +--------------- + +Create a new environment (with your preferred python version). + +.. code-block:: bash + + conda create --name paramak_env python=3.11 + + +Then activate the new environment. + +.. code-block:: bash + + conda activate paramak_env + +Then install the Paramak. + +.. code-block:: bash + + mamba install -c conda-forge paramak + +Now you should be ready to import paramak from your new python environment. + + + +Developer Installation +---------------------- + +If you want to contribute to the paramak or then you might want to install the +package in a more dynamic manner. + +Download and install MiniConda, create a new python environment and activate the +environment as covered in the installation procedure above. + +Then install CadQuery with Conda, Mamba or pip. + +.. code-block:: bash + + conda install -c conda-forge cadquery + mamba install -c conda-forge cadquery + pip install cadquery + + +Then clone the repository + +.. code-block:: bash + + git clone https://github.com/fusion-energy/paramak.git + +Navigate to the paramak repository and within the terminal install the paramak +package and the dependencies using pip with e -e (developer option). + +.. code-block:: bash + + cd paramak + pip install -e . diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 000000000..32bb24529 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/python_api.rst b/docs/python_api.rst new file mode 100644 index 000000000..d4504f06d --- /dev/null +++ b/docs/python_api.rst @@ -0,0 +1,30 @@ +Python API reference +==================== + + +.. currentmodule:: paramak + +Assemblies +---------- + +.. autofunction:: tokamak + +.. autofunction:: tokamak_from_plasma + +.. autofunction:: spherical_tokamak + +.. autofunction:: spherical_tokamak_from_plasma + +Workplanes +---------- + +.. autofunction:: blanket_from_plasma +.. autofunction:: center_column_shield_cylinder +.. autofunction:: constant_thickness_dome +.. autofunction:: cutting_wedge +.. autofunction:: dished_vacuum_vessel +.. autofunction:: plasma_simplified +.. autofunction:: poloidal_field_coil_case +.. autofunction:: poloidal_field_coil +.. autofunction:: toroidal_field_coil_rectangle +.. autofunction:: u_shaped_dome \ No newline at end of file diff --git a/docs/usage.rst b/docs/usage.rst new file mode 100644 index 000000000..91fd7dfc0 --- /dev/null +++ b/docs/usage.rst @@ -0,0 +1,693 @@ +Usage +===== + +There are two main reactors to choose from, Tokamak and Spherical Tokamak. +These can be built with: + - A radial and vertical build + - A radial build and plasma elongation. + +The former gives the user more control of the size of components allowing reactor blankets to vary both radially and vertically. +The later allows reactors to be built with a minimal number of parameters. +In all cases it is possible to add additional components such as divertors, poloidal and toroidal magnets and any self made geometry as a CadQuery Workplane. +The reactors can be varied in terms of their radial build, vertical build, elongation and triangularity which gives a lot of variability. +These examples show how to make various reactors with and without different components, each example is minimal and aims to show a single feature, you will have to combine examples to make a complete model. + +Tokamak from plasma +------------------- + +- The tokamak_from_plasma function provides a parametric tokamak shaped reactor. +- This is characterized by a continuous blanket that goes around the inboard and outboard sides of the plasma. +- This reactor requires few arguments to create as it keeps the vertical build of the blanket layers the same thickness as the radial build. + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + import paramak + result = paramak.tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=90, + ).toCompound() + + +.. code-block:: python + + import paramak + result = paramak.tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=90, + ) + + +Spherical tokamak from plasma +----------------------------- + +- The spherical_tokamak_from_plasma function provides a parametric tokamak shaped reactor. +- This is characterized by a blanket that only goes around the outboard sides of the plasma. +- This reactor requires few arguments to create as it keeps the vertical build of the blanket layers the same thickness as the radial build. + + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + import paramak + result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 60), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=90, + ).toCompound() + + +.. code-block:: python + + import paramak + result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 60), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=90, + ) + result.save('reactor.step') + + + +Tokamak +------- + +- The spherical_tokamak function provides a parametric tokamak shaped reactor. +- This is characterized by a blanket that only goes around the outboard sides of the plasma. +- This reactor allows for a separate vertical and radial build which allows different thickness layers in the blanket. + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + import paramak + + result = paramak.tokamak( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + triangularity=0.55, + rotation_angle=180, + ).toCompound() + +.. code-block:: python + + import paramak + + result = paramak.tokamak( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + triangularity=0.55, + rotation_angle=180, + ) + + result.save(f"tokamak_minimal.step") + + + +Spherical tokamak +----------------- + +- The spherical_tokamak function provides a parametric tokamak shaped reactor. +- This is characterized by a blanket that only goes around the outboard sides of the plasma. +- This reactor allows for a separate vertical and radial build which allows different thickness layers in the blanket. + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + import paramak + + result = paramak.spherical_tokamak( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 10), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 10), + ("layer", 60), + ("layer", 10), + ] + ], + vertical_build=[ + ("layer", 15), + ("layer", 120), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 30), + ("layer", 15), + ], + rotation_angle=180, + triangularity=0.55, + ).toCompound() + +.. code-block:: python + + import paramak + + result = paramak.spherical_tokamak( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 10), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 10), + ("layer", 60), + ("layer", 10), + ] + ], + vertical_build=[ + ("layer", 15), + ("layer", 120), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 30), + ("layer", 15), + ], + rotation_angle=180, + triangularity=0.55, + ) + + result.save(f"spherical_tokamak_minimal.step") + +Reactor with divertor(s) +------------------------ + +- ll reactors support adding additional radial builds for the lower_divertor and or the upper_divertor. +- This example adds two divertors to a spherical_tokamak_from_plasma reactor but and other reactor would also work. + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ], + [("gap", 75), ("lower_divertor", 100)], # this divertor connects to the center column + [("gap", 120), ("upper_divertor", 100)], # this divertor has some blanket between the center colum and itself + ], + elongation=2, + triangularity=0.55, + rotation_angle=180, + ).toCompound() + + +.. code-block:: python + + result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ], + [("gap", 75), ("lower_divertor", 100)], # this divertor connects to the center column + [("gap", 120), ("upper_divertor", 140)], # this divertor has some blanket between the center colum and itself + ], + elongation=2, + triangularity=0.55, + rotation_angle=180, + ) + result.save('reactor.step') + +Reactor with poloidal field coils +--------------------------------- + +- All reactors support adding a sequence of CadQuery shapes (e.g. workplanes) to the reactor using the add_extra_cut_shapes argument +- This example adds PF coils to a spherical_tokamak_from_plasma reactor but and other reactor would also work. + + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + import paramak + + add_extra_cut_shapes = [] + for case_thickness, height, width, center_point in zip( + [10, 15, 15, 10], [20, 50, 50, 20], [20, 50, 50, 20], + [(500, 300), (560, 100), (560, -100), (500, -300)] + ): + add_extra_cut_shapes.append( + paramak.poloidal_field_coil( + height=height, width=width, center_point=center_point, rotation_angle=270 + ) + ) + add_extra_cut_shapes.append( + paramak.poloidal_field_coil_case( + coil_height=height, + coil_width=width, + casing_thickness=case_thickness, + rotation_angle=270, + center_point=center_point, + ) + ) + + result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=270, + add_extra_cut_shapes=add_extra_cut_shapes, + ).toCompound() + + +.. code-block:: python + + import paramak + + add_extra_cut_shapes = [] + for case_thickness, height, width, center_point in zip( + [10, 15, 15, 10], [20, 50, 50, 20], [20, 50, 50, 20], + [(500, 300), (560, 100), (560, -100), (500, -300)] + ): + add_extra_cut_shapes.append( + paramak.poloidal_field_coil( + height=height, width=width, center_point=center_point, rotation_angle=270 + ) + ) + add_extra_cut_shapes.append( + paramak.poloidal_field_coil_case( + coil_height=height, + coil_width=width, + casing_thickness=case_thickness, + rotation_angle=270, + center_point=center_point, + ) + ) + + result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=270, + add_extra_cut_shapes=add_extra_cut_shapes, + ) + + result.save(f"spherical_tokamak_from_plasma_with_pf_magnets.step") + + +Reactor with toroidal field coils +--------------------------------- + +- In a similar way to adding poloidal field coils one can also add toroidal field coils by making use of the add_extra_cut_shapes argument. +- All reactors support adding a sequence of CadQuery shapes (e.g. workplanes) to the reactor using the add_extra_cut_shapes argument +- This example adds TF coils to a spherical_tokamak_from_plasma reactor but and other reactor would also work. +- Also these are rectangle shaped TF coils but other shapes are also available. + + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + import paramak + + tf = paramak.toroidal_field_coil_rectangle( + horizontal_start_point = (10, 520), + vertical_mid_point = (600, 0), + thickness = 50, + distance = 40, + with_inner_leg = True, + azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180], + ) + + result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 70), + ("layer", 10), + ("layer", 10), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 10), + ("layer", 60), + ("layer", 10), + ], + elongation=2.5, + rotation_angle=180, + triangularity=0.55, + add_extra_cut_shapes=[tf] + ).toCompound() + +.. code-block:: python + + import paramak + + tf = paramak.toroidal_field_coil_rectangle( + horizontal_start_point = (10, 520), + vertical_mid_point = (600, 0), + thickness = 50, + distance = 40, + with_inner_leg = True, + azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180], + ) + + result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 70), + ("layer", 10), + ("layer", 10), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 10), + ("layer", 60), + ("layer", 10), + ], + elongation=2.5, + rotation_angle=180, + triangularity=0.55, + add_extra_cut_shapes=[tf] + ) + + result.save(f"spherical_tokamak_minimal.step") + + +Tokamak with negative triangularity +----------------------------------- + +- The triangularity argument can be set to a negative value to make a plasma with a negative triangularity. +- This example makes a tokamak with a negative but this would work on any reactor. + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + import paramak + + result = paramak.tokamak( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + triangularity=-0.55, + rotation_angle=180, + ).toCompound() + +.. code-block:: python + + import paramak + + result = paramak.tokamak( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + triangularity=-0.55, + rotation_angle=180, + ) + + result.save(f"tokamak_minimal.step") + + + +Spherical tokamak with negative triangularity +--------------------------------------------- + +- The triangularity argument can be set to a negative value to make a plasma with a negative triangularity. +- This example makes a spherical tokamak with a negative but this would work on any reactor. + +.. cadquery:: + :gridsize: 0 + :select: result + :color: #00cd00 + :width: 100% + :height: 600px + + import paramak + + result = paramak.spherical_tokamak( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 40), + ("layer", 60), + ("layer", 10), + ] + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + rotation_angle=180, + triangularity=-0.55, + ).toCompound() + +.. code-block:: python + + import paramak + + result = paramak.spherical_tokamak( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 10), + ("layer", 60), + ("layer", 10), + ] + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + rotation_angle=180, + triangularity=-0.55, + ) + result.save(f"spherical_tokamak_minimal.step") diff --git a/docs/version_switcher.json b/docs/version_switcher.json new file mode 100644 index 000000000..f4cb2ff5a --- /dev/null +++ b/docs/version_switcher.json @@ -0,0 +1,17 @@ +[ + { + "name": "dev", + "version": "dev", + "url": "https://fusion-energy.github.io/paramak2/dev" + }, + { + "name": "0.9.1", + "version": "0.9.1", + "url": "https://fusion-energy.github.io/paramak2/0.9.1" + }, + { + "name": "0.9.0", + "version": "0.9.0", + "url": "https://fusion-energy.github.io/paramak2/0.9.0" + } +] diff --git a/environment_dev.yml b/environment_dev.yml new file mode 100644 index 000000000..65c597573 --- /dev/null +++ b/environment_dev.yml @@ -0,0 +1,20 @@ +name: paramak-dev-environment + +dependencies: + - openmc=0.15.0=dagmc* + - ocp=7.7.2.1 + - cadquery=2.4.0 + - gmsh + - python-gmsh + - numpy<=1.26.4 + - pytest + - pip + - mpmath + - sympy + - cad_to_dagmc + - sphinx + - pydata-sphinx-theme + - pip: + - sphinx_autodoc_typehints + - openmc_data_downloader + - sphinxcadquery diff --git a/examples/example_util_functions.py b/examples/example_util_functions.py new file mode 100644 index 000000000..d9f691c0e --- /dev/null +++ b/examples/example_util_functions.py @@ -0,0 +1,102 @@ +def transport_particles_on_h5m_geometry( + h5m_filename: str, + material_tags: list, + nuclides: list = None, + cross_sections_xml: str = None, +): + """A function for testing the geometry file with particle transport in + DAGMC OpenMC. Requires openmc and either the cross_sections_xml to be + specified. Returns the flux on each volume + + Arg: + h5m_filename: The name of the DAGMC h5m file to test + material_tags: The + nuclides: + cross_sections_xml: + + """ + + import openmc + from openmc.data import NATURAL_ABUNDANCE + + if nuclides is None: + nuclides = list(NATURAL_ABUNDANCE.keys()) + + materials = openmc.Materials() + for i, material_tag in enumerate(set(material_tags)): + # simplified material definitions have been used to keen this example minimal + mat_dag_material_tag = openmc.Material(name=material_tag) + mat_dag_material_tag.add_nuclide(nuclides[i], 1, "ao") + mat_dag_material_tag.set_density("g/cm3", 0.01) + + materials.append(mat_dag_material_tag) + + if cross_sections_xml: + openmc.config["cross_sections"] = cross_sections_xml + + else: + with open("cross_sections.xml", "w") as file: + file.write( + """ + + + + + + """ + ) + + openmc.config["cross_sections"] = "cross_sections.xml" + + dag_univ = openmc.DAGMCUniverse(filename=h5m_filename) + bound_dag_univ = dag_univ.bounded_universe() + geometry = openmc.Geometry(root=bound_dag_univ) + + # initializes a new source object + my_source = openmc.IndependentSource() + + center_of_geometry = ( + (dag_univ.bounding_box[0][0] + dag_univ.bounding_box[1][0]) / 2, + (dag_univ.bounding_box[0][1] + dag_univ.bounding_box[1][1]) / 2, + (dag_univ.bounding_box[0][2] + dag_univ.bounding_box[1][2]) / 2, + ) + # sets the location of the source which is not on a vertex + center_of_geometry_nudged = ( + center_of_geometry[0] + 0.1, + center_of_geometry[1] + 0.1, + center_of_geometry[2] + 0.1, + ) + + my_source.space = openmc.stats.Point(center_of_geometry_nudged) + # sets the direction to isotropic + my_source.angle = openmc.stats.Isotropic() + # sets the energy distribution to 100% 14MeV neutrons + my_source.energy = openmc.stats.Discrete([14e6], [1]) + + # specifies the simulation computational intensity + settings = openmc.Settings() + settings.batches = 10 + settings.particles = 10000 + settings.inactive = 0 + settings.run_mode = "fixed source" + settings.source = my_source + + # adds a tally to record the heat deposited in entire geometry + cell_tally = openmc.Tally(name="flux") + cell_tally.scores = ["flux"] + + # groups the two tallies + tallies = openmc.Tallies([cell_tally]) + + # builds the openmc model + my_model = openmc.Model(materials=materials, geometry=geometry, settings=settings, tallies=tallies) + + # starts the simulation + output_file = my_model.run() + + # loads up the output file from the simulation + statepoint = openmc.StatePoint(output_file) + + my_flux_cell_tally = statepoint.get_tally(name="flux") + + return my_flux_cell_tally.mean.flatten()[0] diff --git a/examples/plasma.py b/examples/plasma.py new file mode 100644 index 000000000..ff3f82f73 --- /dev/null +++ b/examples/plasma.py @@ -0,0 +1,12 @@ +from cadquery import exporters + +import paramak + +solid = paramak.plasma_simplified(rotation_angle=180, triangularity=-0.55) +exporters.export(solid, "plasma_simplified.step") +# solid.exportStep("plasma_simplified.step") + +# solid = paramak.plasma_simplified(rotation_angle=360) + +# exporters.export(solid, "plasma_simplified.step") +# solid.exportStep("plasma_simplified.step") diff --git a/examples/spherical_tokamak_from_plasma_minimal.py b/examples/spherical_tokamak_from_plasma_minimal.py new file mode 100644 index 000000000..a8d12eb3c --- /dev/null +++ b/examples/spherical_tokamak_from_plasma_minimal.py @@ -0,0 +1,48 @@ +from pathlib import Path + +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +my_reactor = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ] + ], + elongation=2, + triangularity=0.55, + rotation_angle=180, +) +my_reactor.save(f"spherical_tokamak_from_plasma_minimal.step") + +# needed to downgrade pip with ... python -m pip install pip==24.0 +# from jupyter_cadquery import show +# view = show(my_reactor) +# view.export_html('3d.html') + + +# my_model = CadToDagmc() +# material_tags = [ +# "mat1" +# ] * 6 +# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +# h5m_filename = "dagmc.h5m" +# flux = transport_particles_on_h5m_geometry( +# h5m_filename=h5m_filename, +# material_tags=material_tags, +# nuclides=["H1"] * len(material_tags), +# cross_sections_xml="tests/cross_sections.xml", +# ) +# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_divertor.py b/examples/spherical_tokamak_from_plasma_with_divertor.py new file mode 100644 index 000000000..d34215a33 --- /dev/null +++ b/examples/spherical_tokamak_from_plasma_with_divertor.py @@ -0,0 +1,49 @@ +from pathlib import Path + +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +my_reactor = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ], + [("gap", 75), ("lower_divertor", 100)], # this divertor connects to the center column + [("gap", 120), ("upper_divertor", 140)], # this divertor has some blanket between the center colum and itself + ], + elongation=2, + triangularity=0.55, + rotation_angle=180, +) +my_reactor.save("spherical_tokamak_from_plasma_with_divertor.step") +print('written spherical_tokamak_from_plasma_with_divertor.step') + +# needed to downgrade pip with ... python -m pip install pip==24.0 +# from jupyter_cadquery import show +# view = show(my_reactor) +# view.export_html("3d.html") + + +# my_model = CadToDagmc() +# material_tags = ["mat1"] * 21 # the two divertors split the 3 blanket layers into 9 and the magnets also splt the blanket. +# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +# h5m_filename = "dagmc.h5m" +# flux = transport_particles_on_h5m_geometry( +# h5m_filename=h5m_filename, +# material_tags=material_tags, +# nuclides=["H1"] * len(material_tags), +# cross_sections_xml="tests/cross_sections.xml", +# ) +# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_pf_magnets.py b/examples/spherical_tokamak_from_plasma_with_pf_magnets.py new file mode 100644 index 000000000..7ccb7f2e8 --- /dev/null +++ b/examples/spherical_tokamak_from_plasma_with_pf_magnets.py @@ -0,0 +1,64 @@ +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +add_extra_cut_shapes = [] +for case_thickness, height, width, center_point in zip( + [10, 15, 15, 10], [20, 50, 50, 20], [20, 50, 50, 20], + [(500, 300), (560, 100), (560, -100), (500, -300)] +): + add_extra_cut_shapes.append( + paramak.poloidal_field_coil( + height=height, width=width, center_point=center_point, rotation_angle=270 + ) + ) + add_extra_cut_shapes.append( + paramak.poloidal_field_coil_case( + coil_height=height, + coil_width=width, + casing_thickness=case_thickness, + rotation_angle=270, + center_point=center_point, + ) + ) + + +my_reactor = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=270, + add_extra_cut_shapes=add_extra_cut_shapes, +) +my_reactor.save(f"spherical_tokamak_from_plasma_with_pf_magnets.step") + +# needed to downgrade pip with ... python -m pip install pip==24.0 +# from jupyter_cadquery import show +# view = show(my_reactor) +# view.export_html('spherical_tokamak_from_plasma.html') + + +# my_model = CadToDagmc() +# material_tags = ["mat1"] * 5 +# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +# h5m_filename = "dagmc.h5m" +# flux = transport_particles_on_h5m_geometry( +# h5m_filename=h5m_filename, +# material_tags=material_tags, +# nuclides=["H1"] * len(material_tags), +# cross_sections_xml="tests/cross_sections.xml", +# ) +# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py b/examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py new file mode 100644 index 000000000..6edd914b4 --- /dev/null +++ b/examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py @@ -0,0 +1,69 @@ +from pathlib import Path + +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +rotation_angle = 180 +poloidal_field_coils = [] +for case_thickness, height, width, center_point in zip( + [10, 15, 15, 10], [20, 50, 50, 20], [20, 50, 50, 20], [(500, 300), (560, 100), (560, -100), (500, -300)] +): + poloidal_field_coils.append( + paramak.poloidal_field_coil( + height=height, width=width, center_point=center_point, rotation_angle=rotation_angle + ) + ) + poloidal_field_coils.append( + paramak.poloidal_field_coil_case( + coil_height=height, + coil_width=width, + casing_thickness=case_thickness, + rotation_angle=rotation_angle, + center_point=center_point, + ) + ) + + +my_reactor = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ], + [("gap", 75), ("lower_divertor", 100)], + ], + elongation=2, + triangularity=0.55, + rotation_angle=rotation_angle, + add_extra_cut_shapes=poloidal_field_coils, +) +my_reactor.save(f"spherical_tokamak_from_plasma_with_pf_magnets_and_divertor.step") + +# needed to downgrade pip with ... python -m pip install pip==24.0 +# from jupyter_cadquery import show +# view = show(my_reactor) +# view.export_html('3d.html') + + +# my_model = CadToDagmc() +# material_tags = ["mat1"] * 5 +# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +# h5m_filename = "dagmc.h5m" +# flux = transport_particles_on_h5m_geometry( +# h5m_filename=h5m_filename, +# material_tags=material_tags, +# nuclides=["H1"] * len(material_tags), +# cross_sections_xml="tests/cross_sections.xml", +# ) +# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_tf_magnets.py b/examples/spherical_tokamak_from_plasma_with_tf_magnets.py new file mode 100644 index 000000000..edfa10ee0 --- /dev/null +++ b/examples/spherical_tokamak_from_plasma_with_tf_magnets.py @@ -0,0 +1,57 @@ +from pathlib import Path + +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +tf = paramak.toroidal_field_coil_rectangle( + horizontal_start_point = (10, 520), + vertical_mid_point = (600, 0), + thickness = 50, + distance = 40, + with_inner_leg = True, + azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180], +) + +result = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 70), + ("layer", 10), + ("layer", 10), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 10), + ("layer", 60), + ("layer", 10), + ], + elongation=2.5, + rotation_angle=180, + triangularity=0.55, + add_extra_cut_shapes=[tf] +) + +result.save(f"spherical_tokamak_minimal.step") + + +# needed to downgrade pip with ... python -m pip install pip==24.0 +from jupyter_cadquery import show + +view = show(result) +view.export_html("3d.html") + + +my_model = CadToDagmc() +material_tags = ["mat1"] * 7 +my_model.add_cadquery_object(cadquery_object=result, material_tags=material_tags) +my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +h5m_filename = "dagmc.h5m" +flux = transport_particles_on_h5m_geometry( + h5m_filename=h5m_filename, + material_tags=material_tags, + nuclides=["H1"] * len(material_tags), + cross_sections_xml="tests/cross_sections.xml", +) +assert flux > 0.0 diff --git a/examples/spherical_tokamak_minimal.py b/examples/spherical_tokamak_minimal.py new file mode 100644 index 000000000..67490b2e3 --- /dev/null +++ b/examples/spherical_tokamak_minimal.py @@ -0,0 +1,57 @@ +from pathlib import Path + +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +my_reactor = paramak.spherical_tokamak( + radial_builds=[ + [ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 40), + ("layer", 60), + ("layer", 10), + ] + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + rotation_angle=180, + triangularity=-0.55, +) +my_reactor.save(f"spherical_tokamak_minimal.step") + +# needed to downgrade pip with ... python -m pip install pip==24.0 +from jupyter_cadquery import show + +view = show(my_reactor) +view.export_html("3d.html") + + +my_model = CadToDagmc() +material_tags = ["mat1"] * 6 +my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +h5m_filename = "dagmc.h5m" +flux = transport_particles_on_h5m_geometry( + h5m_filename=h5m_filename, + material_tags=material_tags, + nuclides=["H1"] * len(material_tags), + cross_sections_xml="tests/cross_sections.xml", +) +assert flux > 0.0 diff --git a/examples/tokamak_from_plasma_minimal.py b/examples/tokamak_from_plasma_minimal.py new file mode 100644 index 000000000..1ecb19e49 --- /dev/null +++ b/examples/tokamak_from_plasma_minimal.py @@ -0,0 +1,40 @@ +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +my_reactor = paramak.tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=180, +) +my_reactor.save(f"tokamak_minimal.step") +print(f"Saved as tokamak_minimal.step") + +# my_model = CadToDagmc() +# material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model +# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +# h5m_filename = "dagmc.h5m" +# flux = transport_particles_on_h5m_geometry( +# h5m_filename=h5m_filename, +# material_tags=material_tags, +# nuclides=["H1"] * len(material_tags), +# cross_sections_xml="tests/cross_sections.xml", +# ) +# assert flux > 0.0 diff --git a/examples/tokamak_from_plasma_with_divertor.py b/examples/tokamak_from_plasma_with_divertor.py new file mode 100644 index 000000000..6ec86e8c3 --- /dev/null +++ b/examples/tokamak_from_plasma_with_divertor.py @@ -0,0 +1,43 @@ +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +my_reactor = paramak.tokamak( + radial_builds=[ + [ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + [("gap", 300), ("lower_divertor", 150)], + ], + elongation=2, + triangularity=0.55, + rotation_angle=180, +) +my_reactor.save(f"tokamak_with_divertor.step") +print(f"Saved as tokamak_with_divertor.step") + +# my_model = CadToDagmc() +# material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model +# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +# h5m_filename = "dagmc.h5m" +# flux = transport_particles_on_h5m_geometry( +# h5m_filename=h5m_filename, +# material_tags=material_tags, +# nuclides=["H1"] * len(material_tags), +# cross_sections_xml="tests/cross_sections.xml", +# ) +# assert flux > 0.0 diff --git a/examples/tokamak_minimal.py b/examples/tokamak_minimal.py new file mode 100644 index 000000000..93758e615 --- /dev/null +++ b/examples/tokamak_minimal.py @@ -0,0 +1,54 @@ +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +my_reactor = paramak.tokamak( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + triangularity=0.55, + rotation_angle=180, +) + +from cadquery.vis import show +show(my_reactor) + +my_reactor.save(f"tokamak_minimal.step") +print(f"Saved as tokamak_minimal.step") + +my_model = CadToDagmc() +material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model +my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +h5m_filename = "dagmc.h5m" +flux = transport_particles_on_h5m_geometry( + h5m_filename=h5m_filename, + material_tags=material_tags, + nuclides=["H1"] * len(material_tags), + cross_sections_xml="tests/cross_sections.xml", +) +assert flux > 0.0 diff --git a/examples/tokamak_with_pf_magnets.py b/examples/tokamak_with_pf_magnets.py new file mode 100644 index 000000000..e954d995a --- /dev/null +++ b/examples/tokamak_with_pf_magnets.py @@ -0,0 +1,75 @@ +from cad_to_dagmc import CadToDagmc +from example_util_functions import transport_particles_on_h5m_geometry + +import paramak + +add_extra_cut_shapes = [] +for case_thickness, height, width, center_point in zip( + [10, 15, 15, 10], [20, 50, 50, 20], [20, 50, 50, 20], + [(700, 300), (800, 100), (800, -100), (700, -300)] +): + add_extra_cut_shapes.append( + paramak.poloidal_field_coil( + height=height, width=width, center_point=center_point, rotation_angle=180 + ) + ) + add_extra_cut_shapes.append( + paramak.poloidal_field_coil_case( + coil_height=height, + coil_width=width, + casing_thickness=case_thickness, + rotation_angle=180, + center_point=center_point, + ) + ) + +my_reactor = paramak.tokamak( + radial_builds=[ + ("gap", 10), + ("layer", 30), + ("layer", 50), + ("layer", 10), + ("layer", 120), + ("layer", 20), + ("gap", 60), + ("plasma", 300), + ("gap", 60), + ("layer", 20), + ("layer", 120), + ("layer", 10), + ], + vertical_build=[ + ("layer", 15), + ("layer", 80), + ("layer", 10), + ("gap", 50), + ("plasma", 700), + ("gap", 60), + ("layer", 10), + ("layer", 40), + ("layer", 15), + ], + triangularity=0.55, + rotation_angle=180, + add_extra_cut_shapes=add_extra_cut_shapes, +) + +from cadquery.vis import show +show(my_reactor) + +my_reactor.save(f"tokamak_minimal.step") +print(f"Saved as tokamak_minimal.step") + +my_model = CadToDagmc() +material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model +my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) + +h5m_filename = "dagmc.h5m" +flux = transport_particles_on_h5m_geometry( + h5m_filename=h5m_filename, + material_tags=material_tags, + nuclides=["H1"] * len(material_tags), + cross_sections_xml="tests/cross_sections.xml", +) +assert flux > 0.0 diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 000000000..e859fd881 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,55 @@ +[build-system] +requires = ["setuptools >= 69.5.1", "setuptools_scm[toml]>=8.1.0"] +build-backend = "setuptools.build_meta" + +[project] +name = "paramak" +dynamic = ["version"] +description = "Create 3D fusion reactor CAD models based on input parameters" +readme = "README.md" +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", +] +authors = [ + { name="The Paramak Development Team" }, +] +license = {file = "LICENSE.txt"} +requires-python = ">=3.8" +keywords = ["python", "geometry", "reactor", "model", "cad", "fusion", "parametric", "dagmc", "openmc"] + +[project.urls] +"Homepage" = "https://github.com/fusion-energy/paramak" +"Bug Tracker" = "https://github.com/fusion-energy/paramak/issues" +"Documentation" = "https://paramak.readthedocs.io" + + +[project.optional-dependencies] +tests = [ + "pytest>=5.4.3", + "pytest-cov>=2.12.1", + "dagmc_h5m_file_inspector>=0.5.0", + "openmc_data_downloader" +] +docs = [ + "sphinx_autodoc_typehints", + "openmc_data_downloader", + "sphinxcadquery", + "sphinx", + "pydata-sphinx-theme", + "cadquery", + "numpy<=1.26.4", + "mpmath", + "sympy", + "scipy" +] + +[tool.black] +line-length = 120 + +[tool.setuptools_scm] +write_to = "src/_version.py" + +[tool.setuptools] +package-dir = {"" = "src"} diff --git a/src/paramak/__init__.py b/src/paramak/__init__.py new file mode 100644 index 000000000..59a2cf4eb --- /dev/null +++ b/src/paramak/__init__.py @@ -0,0 +1,23 @@ +from importlib.metadata import version + +from .assemblies.spherical_tokamak import spherical_tokamak, spherical_tokamak_from_plasma +from .assemblies.tokamak import tokamak, tokamak_from_plasma + +from .workplanes.blanket_constant_thickness_arc_h import blanket_constant_thickness_arc_h +from .workplanes.blanket_from_plasma import blanket_from_plasma +from .workplanes.center_column_shield_cylinder import center_column_shield_cylinder +from .workplanes.constant_thickness_dome import constant_thickness_dome +from .workplanes.cutting_wedge import cutting_wedge +from .workplanes.dished_vacuum_vessel import dished_vacuum_vessel +from .workplanes.plasma_simplified import plasma_simplified +from .workplanes.poloidal_field_coil import poloidal_field_coil +from .workplanes.poloidal_field_coil_case import poloidal_field_coil_case +from .workplanes.toroidal_field_coil_rectangle import toroidal_field_coil_rectangle +from .workplanes.u_shaped_dome import u_shaped_dome + +# from .utils import * + +__version__ = version("paramak") + + +__all__ = ["__version__"] diff --git a/src/paramak/assemblies/__init__.py b/src/paramak/assemblies/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/paramak/assemblies/spherical_tokamak.py b/src/paramak/assemblies/spherical_tokamak.py new file mode 100644 index 000000000..d30f871a8 --- /dev/null +++ b/src/paramak/assemblies/spherical_tokamak.py @@ -0,0 +1,261 @@ +from typing import Optional, Sequence, Tuple, Union + +import cadquery as cq +import numpy as np + +from ..utils import ( + build_divertor_modify_blanket, + extract_radial_builds, + get_plasma_index, + get_plasma_value, + sum_up_to_gap_before_plasma, + sum_up_to_plasma, + sum_before_after_plasma, +) +from ..workplanes.blanket_from_plasma import blanket_from_plasma +from ..workplanes.center_column_shield_cylinder import center_column_shield_cylinder +from ..workplanes.plasma_simplified import plasma_simplified + + +def create_blanket_layers_after_plasma( + radial_build, vertical_build, minor_radius, major_radius, triangularity, elongation, rotation_angle, center_column +): + layers = [] + cumulative_thickness_rb = 0 + cumulative_thickness_uvb = 0 + cumulative_thickness_lvb = 0 + + plasma_index_radial = get_plasma_index(radial_build) + plasma_index_vertical = get_plasma_index(vertical_build) + + for i, item in enumerate(radial_build[plasma_index_radial + 1 :]): + upper_thicknees = vertical_build[plasma_index_vertical + 1 + i][1] + lower_thicknees = vertical_build[plasma_index_vertical - 1 - i][1] + radial_thickness = item[1] + + if item[0] == "gap": + cumulative_thickness_rb += radial_thickness + cumulative_thickness_uvb += upper_thicknees + cumulative_thickness_lvb += lower_thicknees + continue + + layer = blanket_from_plasma( + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + thickness=[ + lower_thicknees, + radial_thickness, + upper_thicknees, + ], + offset_from_plasma=[ + cumulative_thickness_lvb, + cumulative_thickness_rb, + cumulative_thickness_uvb, + ], + start_angle=-90, + stop_angle=90, + rotation_angle=rotation_angle, + color=(0.5, 0.5, 0.5), + name=f"layer_{plasma_index_radial+i+1}", + allow_overlapping_shape=True, + connect_to_center=True, + ) + layer = layer.cut(center_column) + cumulative_thickness_rb += radial_thickness + cumulative_thickness_uvb += upper_thicknees + cumulative_thickness_lvb += lower_thicknees + layers.append(layer) + + return layers + + +def create_center_column_shield_cylinders(radial_build, vertical_build, rotation_angle): + cylinders = [] + total_sum = 0 + layer_count = 0 + + before, _ = sum_before_after_plasma(vertical_build) + center_column_shield_height = sum([item[1] for item in vertical_build]) + + for index, item in enumerate(radial_build): + if item[0] == "plasma": + break + if item[0] == "gap" and radial_build[index + 1][0] == "plasma": + break + if item[0] == "gap": + total_sum += item[1] + continue + + layer_count += 1 + # print('inner_radius', total_sum, 'item thickness', item[1], 'layer_count', layer_count) + cylinder = center_column_shield_cylinder( + inner_radius=total_sum, + thickness=item[1], + name=f"layer_{layer_count}", + rotation_angle=rotation_angle, + height=center_column_shield_height, + reference_point=("lower", -before), + ) + cylinders.append(cylinder) + total_sum += item[1] + + return cylinders + + +def spherical_tokamak_from_plasma( + radial_builds: Union[Sequence[Sequence[Tuple[str, float]]], Sequence[Tuple[str, float]]], + elongation: float = 2.0, + triangularity: float = 0.55, + rotation_angle: float = 180.0, + add_extra_cut_shapes: Sequence[cq.Workplane] = [], +): + """_summary_ + + Args: + + elongation (float, optional): _description_. Defaults to 2.0. + triangularity (float, optional): _description_. Defaults to 0.55. + rotation_angle (Optional[str], optional): _description_. Defaults to 180.0. + add_extra_cut_shapes (Sequence, optional): _description_. Defaults to []. + + Returns: + _type_: _description_ + """ + + plasma_radial_build, _ = extract_radial_builds(radial_builds) + + inner_equatorial_point = sum_up_to_plasma(plasma_radial_build) + plasma_radial_thickness = get_plasma_value(plasma_radial_build) + outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness + + # sets major radius and minor radius from equatorial_points to allow a + # radial build. This helps avoid the plasma overlapping the center + # column and other components + major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 + minor_radius = major_radius - inner_equatorial_point + + # make vertical build from outer radial build + pi = get_plasma_index(plasma_radial_build) + upper_vertical_build = plasma_radial_build[pi:] + + plasma_height = 2 * minor_radius * elongation + # slice opperation reverses the list and removes the last value to avoid two plasmas + vertical_build = upper_vertical_build[::-1][:-1] + [("plasma", plasma_height)] + upper_vertical_build[1:] + + return spherical_tokamak( + radial_builds=radial_builds, + vertical_build=vertical_build, + triangularity=triangularity, + rotation_angle=rotation_angle, + add_extra_cut_shapes=add_extra_cut_shapes, + ) + + +def spherical_tokamak( + radial_builds: Union[Sequence[Sequence[Tuple[str, float]]], Sequence[Tuple[str, float]]], + vertical_build: Sequence[Tuple[str, float]], + triangularity: float = 0.55, + rotation_angle: Optional[str] = 180.0, + add_extra_cut_shapes: Sequence[cq.Workplane] = [], +): + """_summary_ + + Args: + + radial_build + elongation (float, optional): _description_. Defaults to 2.0. + triangularity (float, optional): _description_. Defaults to 0.55. + rotation_angle (Optional[str], optional): _description_. Defaults to 180.0. + add_extra_cut_shapes (Sequence, optional): _description_. Defaults to []. + + Returns: + _type_: _description_ + """ + + plasma_radial_build, divertor_radial_builds = extract_radial_builds(radial_builds) + + inner_equatorial_point = sum_up_to_plasma(plasma_radial_build) + plasma_radial_thickness = get_plasma_value(plasma_radial_build) + plasma_vertical_thickness = get_plasma_value(vertical_build) + outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness + + # sets major radius and minor radius from equatorial_points to allow a + # radial build. This helps avoid the plasma overlapping the center + # column and other components + major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 + minor_radius = major_radius - inner_equatorial_point + + # vertical build + elongation = (plasma_vertical_thickness / 2) / minor_radius + blanket_rear_wall_end_height = sum([item[1] for item in vertical_build]) + + plasma = plasma_simplified( + major_radius=major_radius, + minor_radius=minor_radius, + elongation=elongation, + triangularity=triangularity, + rotation_angle=rotation_angle, + ) + + inner_radial_build = create_center_column_shield_cylinders( + radial_build=plasma_radial_build, + vertical_build=vertical_build, + rotation_angle=rotation_angle, + ) + + blanket_cutting_cylinder = center_column_shield_cylinder( + inner_radius=0, + thickness=sum_up_to_gap_before_plasma(plasma_radial_build), + rotation_angle=360, + height=2 * blanket_rear_wall_end_height, + ) + + blanket_layers = create_blanket_layers_after_plasma( + radial_build=plasma_radial_build, + vertical_build=vertical_build, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + rotation_angle=rotation_angle, + center_column=blanket_cutting_cylinder, + ) + + divertor_layers, blanket_layers = build_divertor_modify_blanket( + blanket_layers, divertor_radial_builds, blanket_rear_wall_end_height, rotation_angle + ) + + my_assembly = cq.Assembly() + + for i, entry in enumerate(add_extra_cut_shapes): + + if isinstance(entry, cq.Workplane): + my_assembly.add(entry, name=f"add_extra_cut_shape_{i+1}") + else: + raise ValueError(f"add_extra_cut_shapes should only contain cadquery Workplanes, not {type(entry)}") + + if len(add_extra_cut_shapes) == 0: + for i, entry in enumerate(inner_radial_build): + my_assembly.add(entry, name=f"inboard_layer_{i+1})") + for i, entry in enumerate(blanket_layers): + my_assembly.add(entry, name=f"outboard_layer_{i+1})") + for i, entry in enumerate(divertor_layers): + my_assembly.add(entry, name=f"{entry.name})") # gets upper or lower name + else: + shapes_and_components = [] + for i, entry in enumerate(inner_radial_build + blanket_layers + divertor_layers): + for cutter in add_extra_cut_shapes: + entry = entry.cut(cutter) + # TODO use something like this to return a list of material tags for the solids in order, as some solids get split into multiple + # for subentry in entry.objects: + # print(i, subentry) + shapes_and_components.append(entry) + + for i, entry in enumerate(shapes_and_components): + my_assembly.add(entry, name=f"layer_{i+1})") #TODO track the names of shapes, even when extra shapes are made due to splitting + + my_assembly.add(plasma, name="plasma") + + return my_assembly diff --git a/src/paramak/assemblies/tokamak.py b/src/paramak/assemblies/tokamak.py new file mode 100644 index 000000000..37115eebb --- /dev/null +++ b/src/paramak/assemblies/tokamak.py @@ -0,0 +1,288 @@ +from typing import Optional, Sequence, Tuple, Union + +import cadquery as cq + +from ..utils import build_divertor_modify_blanket, extract_radial_builds, get_plasma_index, sum_after_plasma +from ..workplanes.blanket_from_plasma import blanket_from_plasma +from ..workplanes.center_column_shield_cylinder import center_column_shield_cylinder +from ..workplanes.plasma_simplified import plasma_simplified +from .spherical_tokamak import get_plasma_value, sum_up_to_plasma + + +def count_cylinder_layers(radial_build): + before_plasma = 0 + after_plasma = 0 + found_plasma = False + + for item in radial_build: + if item[0] == "plasma": + found_plasma = True + elif item[0] == "layer": + if not found_plasma: + before_plasma += 1 + else: + after_plasma += 1 + + return before_plasma - after_plasma + + +def create_center_column_shield_cylinders(radial_build, rotation_angle, center_column_shield_height): + cylinders = [] + total_sum = 0 + layer_count = 0 + + number_of_cylinder_layers = count_cylinder_layers(radial_build) + + for index, item in enumerate(radial_build): + if item[0] == "plasma": + break + + if item[0] == "gap": + total_sum += item[1] + continue + + thickness = item[1] + layer_count += 1 + + if layer_count > number_of_cylinder_layers: + break + + cylinder = center_column_shield_cylinder( + inner_radius=total_sum, + thickness=item[1], + name=f"layer_{layer_count}", + rotation_angle=rotation_angle, + height=center_column_shield_height, + ) + total_sum += thickness + cylinders.append(cylinder) + return cylinders + + +def distance_to_plasma(radial_build, index): + distance = 0 + for item in radial_build[index + 1 :]: + if item[0] == "plasma": + break + distance += item[1] + return distance + + +def create_layers_from_plasma( + radial_build, vertical_build, minor_radius, major_radius, triangularity, elongation, rotation_angle, center_column +): + + plasma_index_rb = get_plasma_index(radial_build) + plasma_index_vb = get_plasma_index(vertical_build) + indexes_from_plamsa_to_end = len(radial_build) - plasma_index_rb + layers = [] + + cumulative_thickness_orb = 0 + cumulative_thickness_irb = 0 + cumulative_thickness_uvb = 0 + cumulative_thickness_lvb = 0 + for index_delta in range(indexes_from_plamsa_to_end): + + if radial_build[plasma_index_rb + index_delta][0] == "plasma": + continue + outer_layer_thickness = radial_build[plasma_index_rb + index_delta][1] + inner_layer_thickness = radial_build[plasma_index_rb - index_delta][1] + upper_layer_thickness = vertical_build[plasma_index_vb - index_delta][1] + lower_layer_thickness = vertical_build[plasma_index_vb + index_delta][1] + + if radial_build[plasma_index_rb + index_delta][0] == "gap": + cumulative_thickness_orb += outer_layer_thickness + cumulative_thickness_irb += inner_layer_thickness + cumulative_thickness_uvb += upper_layer_thickness + cumulative_thickness_lvb += lower_layer_thickness + continue + + # build outer layer + if radial_build[plasma_index_rb + index_delta][0] == "layer": + outer_layer = blanket_from_plasma( + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + thickness=[ + upper_layer_thickness, + outer_layer_thickness, + lower_layer_thickness + ], + offset_from_plasma=[ + cumulative_thickness_uvb, + cumulative_thickness_orb, + cumulative_thickness_lvb + ], + start_angle=90, + stop_angle=-90, + rotation_angle=rotation_angle, + color=(0.5, 0.5, 0.5), + name=f"layer_{index_delta}", + allow_overlapping_shape=True, + ) + inner_layer = blanket_from_plasma( + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + thickness=[ + lower_layer_thickness, + inner_layer_thickness, + upper_layer_thickness, + ], + offset_from_plasma=[ + cumulative_thickness_lvb, + cumulative_thickness_irb, + cumulative_thickness_uvb, + ], + start_angle=-90, + stop_angle=-270, + rotation_angle=rotation_angle, + color=(0.5, 0.5, 0.5), + name=f"layer_{index_delta}", + allow_overlapping_shape=True, + ) + layer = outer_layer.union(inner_layer) + layers.append(layer) + # layers.append(inner_layer) + cumulative_thickness_orb += outer_layer_thickness + cumulative_thickness_irb += inner_layer_thickness + cumulative_thickness_uvb += upper_layer_thickness + cumulative_thickness_lvb += lower_layer_thickness + # build inner layer + + # union layers + + return layers + +def tokamak_from_plasma( + radial_builds: Union[Sequence[Sequence[Tuple[str, float]]], Sequence[Tuple[str, float]]], + elongation: float = 2.0, + triangularity: float = 0.55, + rotation_angle: float = 180.0, + add_extra_cut_shapes: Sequence[cq.Workplane] = [], +): + + plasma_radial_build, _ = extract_radial_builds(radial_builds) + + inner_equatorial_point = sum_up_to_plasma(plasma_radial_build) + plasma_radial_thickness = get_plasma_value(plasma_radial_build) + outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness + + # sets major radius and minor radius from equatorial_points to allow a + # radial build. This helps avoid the plasma overlapping the center + # column and other components + major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 + minor_radius = major_radius - inner_equatorial_point + + # make vertical build from outer radial build + pi = get_plasma_index(plasma_radial_build) + upper_vertical_build = plasma_radial_build[pi:] + + plasma_height = 2 * minor_radius * elongation + # slice opperation reverses the list and removes the last value to avoid two plasmas + vertical_build = upper_vertical_build[::-1][:-1] + [("plasma", plasma_height)] + upper_vertical_build[1:] + + return tokamak( + radial_builds=radial_builds, + vertical_build=vertical_build, + triangularity=triangularity, + rotation_angle=rotation_angle, + add_extra_cut_shapes=add_extra_cut_shapes, + ) + +def tokamak( + radial_builds: Union[Sequence[Sequence[Tuple[str, float]]], Sequence[Tuple[str, float]]], + vertical_build: Sequence[Tuple[str, float]], + triangularity: float = 0.55, + rotation_angle: float = 180.0, + add_extra_cut_shapes: Sequence[cq.Workplane] = [], +): + """ + Creates a tokamak fusion reactor from a radial build and plasma parameters. + + Args: + radial_builds (Sequence[tuple[str, float]]): A list of tuples containing the radial build of the reactor. + elongation (float, optional): The elongation of the plasma. Defaults to 2.0. + triangularity (float, optional): The triangularity of the plasma. Defaults to 0.55. + rotation_angle (float, optional): The rotation angle of the plasma. Defaults to 180.0. + add_extra_cut_shapes (Sequence, optional): A list of extra shapes to cut the reactor with. Defaults to []. + + Returns: + CadQuery.Assembly: A CadQuery Assembly object representing the tokamak fusion reactor. + """ + + plasma_radial_build, divertor_radial_builds = extract_radial_builds(radial_builds) + + inner_equatorial_point = sum_up_to_plasma(plasma_radial_build) + plasma_radial_thickness = get_plasma_value(plasma_radial_build) + plasma_vertical_thickness = get_plasma_value(vertical_build) + outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness + + major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 + minor_radius = major_radius - inner_equatorial_point + + elongation = (plasma_vertical_thickness / 2) / minor_radius + blanket_rear_wall_end_height = sum([item[1] for item in vertical_build]) + + plasma = plasma_simplified( + major_radius=major_radius, + minor_radius=minor_radius, + elongation=elongation, + triangularity=triangularity, + rotation_angle=rotation_angle, + ) + + inner_radial_build = create_center_column_shield_cylinders( + plasma_radial_build, rotation_angle, blanket_rear_wall_end_height + ) + + blanket_layers = create_layers_from_plasma( + radial_build=plasma_radial_build, + vertical_build=vertical_build, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + rotation_angle=rotation_angle, + center_column=inner_radial_build[0], # blanket_cutting_cylinder, + ) + + divertor_layers, blanket_layers = build_divertor_modify_blanket( + blanket_layers, divertor_radial_builds, blanket_rear_wall_end_height, rotation_angle + ) + + my_assembly = cq.Assembly() + + for i, entry in enumerate(add_extra_cut_shapes): + + if isinstance(entry, cq.Workplane): + my_assembly.add(entry, name=f"add_extra_cut_shape_{i+1}") + else: + raise ValueError(f"add_extra_cut_shapes should only contain cadquery Workplanes, not {type(entry)}") + + if len(add_extra_cut_shapes) == 0: + for i, entry in enumerate(inner_radial_build): + my_assembly.add(entry, name=f"inboard_layer_{i+1})") + for i, entry in enumerate(blanket_layers): + my_assembly.add(entry, name=f"outboard_layer_{i+1})") + for i, entry in enumerate(divertor_layers): + my_assembly.add(entry, name=f"{entry.name})") # gets upper or lower name + else: + shapes_and_components = [] + for i, entry in enumerate(inner_radial_build + blanket_layers + divertor_layers): + for cutter in add_extra_cut_shapes: + entry = entry.cut(cutter) + # TODO use something like this to return a list of material tags for the solids in order, as some solids get split into multiple + # for subentry in entry.objects: + # print(i, subentry) + shapes_and_components.append(entry) + + for i, entry in enumerate(shapes_and_components): + my_assembly.add(entry, name=f"layer_{i+1})") #TODO track the names of shapes, even when extra shapes are made due to splitting + + + my_assembly.add(plasma, name="plasma") + + return my_assembly diff --git a/src/paramak/utils.py b/src/paramak/utils.py new file mode 100644 index 000000000..7704e9fc0 --- /dev/null +++ b/src/paramak/utils.py @@ -0,0 +1,336 @@ +import typing + +from cadquery import Workplane + + +def instructions_from_points(points): + # obtains the first two values of the points list + XZ_points = [(p[0], p[1]) for p in points] + + # obtains the last values of the points list + connections = [p[2] for p in points[:-1]] + + current_linetype = connections[0] + current_points_list = [] + instructions = [] + # groups together common connection types + for i, connection in enumerate(connections): + if connection == current_linetype: + current_points_list.append(XZ_points[i]) + else: + current_points_list.append(XZ_points[i]) + instructions.append({current_linetype: current_points_list}) + current_linetype = connection + current_points_list = [XZ_points[i]] + instructions.append({current_linetype: current_points_list}) + + if list(instructions[-1].values())[0][-1] != XZ_points[0]: + keyname = list(instructions[-1].keys())[0] + instructions[-1][keyname].append(XZ_points[0]) + return instructions + + +def create_wire_workplane_from_instructions( + instructions, + plane="XY", + origin=(0, 0, 0), + obj=None, +): + solid = Workplane(plane, origin=origin, obj=obj) # offset=extrusion_offset + + all_spline = all(list(entry.keys())[0] == "spline" for entry in instructions) + if all_spline: + entry_values = [(list(entry.values())[0]) for entry in instructions][0][:-1] + res = solid.spline( + entry_values, makeWire=True, tol=1e-1, periodic=True + ) # period smooths out the connecting joint + return res + + for entry in instructions: + if list(entry.keys())[0] == "spline": + solid = solid.spline(listOfXYTuple=list(entry.values())[0]) + if list(entry.keys())[0] == "straight": + solid = solid.polyline(list(entry.values())[0]) + if list(entry.keys())[0] == "circle": + p0 = list(entry.values())[0][0] + p1 = list(entry.values())[0][1] + p2 = list(entry.values())[0][2] + solid = solid.moveTo(p0[0], p0[1]).threePointArc(p1, p2) + + return solid.close() + + +def create_wire_workplane_from_points(points, plane, origin=(0, 0, 0), obj=None): + instructions = instructions_from_points(points) + + return create_wire_workplane_from_instructions( + instructions, + plane=plane, + origin=origin, + obj=obj, + ) + + +def rotate_solid(angles: typing.Sequence[float], solid: Workplane) -> Workplane: + rotation_axis = { + "X": [(-1, 0, 0), (1, 0, 0)], + "-X": [(1, 0, 0), (-1, 0, 0)], + "Y": [(0, -1, 0), (0, 1, 0)], + "-Y": [(0, 1, 0), (0, -1, 0)], + "Z": [(0, 0, -1), (0, 0, 1)], + "-Z": [(0, 0, 1), (0, 0, -1)], + } + + rotated_solids = [] + # Perform separate rotations for each angle + for angle in angles: + rotated_solids.append(solid.rotate(*rotation_axis["Z"], angle)) + solid = Workplane(solid.plane) + + # Joins the solids together + for i in rotated_solids: + solid = solid.union(i) + return solid + + +def sum_up_to_gap_before_plasma(radial_build): + total_sum = 0 + for i, item in enumerate(radial_build): + if item[0] == "plasma": + return total_sum + if item[0] == "gap" and i + 1 < len(radial_build) and radial_build[i + 1][0] == "plasma": + return total_sum + total_sum += item[1] + return total_sum + + +def sum_up_to_plasma(radial_build): + total_sum = 0 + for item in radial_build: + if item[0] == "plasma": + break + total_sum += item[1] + return total_sum + + +def sum_after_plasma(radial_build): + plasma_found = False + total_sum = 0 + for item in radial_build: + if plasma_found: + total_sum += item[1] + if item[0] == "plasma": + plasma_found = True + return total_sum + + +class ValidationError(Exception): + pass + + +def sum_before_after_plasma(vertical_build): + before_plasma = 0 + after_plasma = 0 + plasma_value = 0 + plasma_found = False + + for item in vertical_build: + if item[0] == "plasma": + plasma_value = item[1] / 2 + plasma_found = True + continue + if not plasma_found: + before_plasma += item[1] + else: + after_plasma += item[1] + + before_plasma += plasma_value + after_plasma += plasma_value + + return before_plasma, after_plasma + + +def create_divertor_envelope(divertor_radial_build, blanket_height, rotation_angle): + divertor_name = is_lower_or_upper_divertor(divertor_radial_build) + if divertor_name == "lower_divertor": + z_value_sigh = -1 + elif divertor_name == "upper_divertor": + z_value_sigh = 1 + + points = [ + (divertor_radial_build[0][1], z_value_sigh * blanket_height, "straight"), + (divertor_radial_build[0][1], 0, "straight"), + (divertor_radial_build[0][1] + divertor_radial_build[1][1], 0, "straight"), + (divertor_radial_build[0][1] + divertor_radial_build[1][1], z_value_sigh * blanket_height, "straight"), + ] + points.append(points[0]) + + wire = create_wire_workplane_from_points(points=points, plane="XZ", origin=(0, 0, 0), obj=None) + + divertor_solid = wire.revolve(rotation_angle) + divertor_solid.name = divertor_name + return divertor_solid + + +def build_divertor_modify_blanket(outer_layers, divertor_radial_builds, blanket_rear_wall_end_height, rotation_angle): + + divertor_layers = [] + for divertor_radial_build in divertor_radial_builds: + + divertor_solid = create_divertor_envelope(divertor_radial_build, blanket_rear_wall_end_height, rotation_angle) + + # finds the intersection of the blanket and divertor rectangle envelope. + outer_layer_envelope = outer_layers[0] + for outer_layer in outer_layers[1:]: + outer_layer_envelope = outer_layer_envelope.union(outer_layer) + + divertor_solid = divertor_solid.intersect(outer_layer_envelope) + # we reapply this name as it appears to get lost and we might need it when adding to assembly + divertor_solid.name = is_lower_or_upper_divertor(divertor_radial_build) + + # cuts the diverter out of the outer layers of the blanket + for i, layer in enumerate(outer_layers): + layer = layer.cut(divertor_solid) + outer_layers[i] = layer + divertor_layers.append(divertor_solid) + return divertor_layers, outer_layers + + +def is_plasma_radial_build(radial_build): + for entry in radial_build: + if entry[0] == "plasma": + return True + return False + + +def extract_radial_builds(radial_build): + # TODO more rubust method of finding if it is a single list of tupes or multiple lists + # only one radial build so it should be a plasma based radial build + divertor_radial_builds = [] + if isinstance(radial_build[0][0], str) and ( + isinstance(radial_build[0][1], float) or isinstance(radial_build[0][1], int) + ): + plasma_radial_build = radial_build + else: + for entry in radial_build: + if is_plasma_radial_build(entry): + # TODO this assumes htere is only one radial build, which needs to e checked + plasma_radial_build = entry + else: + divertor_radial_builds.append(entry) + + validate_plasma_radial_build(plasma_radial_build) + for divertor_radial_build in divertor_radial_builds: + validate_divertor_radial_build(divertor_radial_build) + return plasma_radial_build, divertor_radial_builds + + +def validate_divertor_radial_build(radial_build): + if len(radial_build) != 2: + raise ValidationError( + f'The radial build for the divertor should only contain two entries, for example (("gap",10), ("lower_divertor", 10)) not {radial_build}' + ) + + if len(radial_build[0]) != 2 or len(radial_build[1]) != 2: + raise ValidationError( + 'The radial build for the divertor should only contain tuples of length 2,, for example ("gap",10)' + ) + + if radial_build[1][0] not in {"lower_divertor", "upper_divertor"}: + raise ValidationError( + f'The second entry in the radial build for the divertor should be either "lower_divertor" or "upper_divertor" not {radial_build[1][0]}' + ) + + if radial_build[0][0] != "gap": + raise ValidationError( + f'The first entry in the radial build for the divertor should be a "gap" not {radial_build[0][0]}' + ) + + if not isinstance(radial_build[0][1], (int, float)) or not isinstance(radial_build[1][1], (int, float)): + raise ValidationError( + f"The thickness of the gap and the divertor should both be integers or floats, not {type(radial_build[0][1])} and {type(radial_build[1][1])}" + ) + + if radial_build[0][1] <= 0 or radial_build[1][1] <= 0: + raise ValidationError( + f"The thickness of the gap and the divertor should both be positive values, not {radial_build[0][1]} and {radial_build[1][1]}" + ) + + +def validate_plasma_radial_build(radial_build): + # TODO should end with layer, not gap + valid_strings = {"gap", "layer", "plasma"} + plasma_count = 0 + plasma_index = -1 + for index, item in enumerate(radial_build): + if not (isinstance(item[0], str) and isinstance(item[1], (int, float))): + raise ValidationError(f"Invalid tuple structure at index {index}: {item}") + if item[0] not in valid_strings: + raise ValidationError(f"Invalid string '{item[0]}' at index {index}") + if item[1] <= 0: + raise ValidationError(f"Non-positive value '{item[1]}' at index {index}") + if item[0] == "plasma": + plasma_count += 1 + plasma_index = index + if plasma_count > 1: + raise ValidationError("Multiple 'plasma' entries found") + if plasma_count != 1: + raise ValidationError("'plasma' entry not found or found multiple times") + if plasma_index == 0 or plasma_index == len(radial_build) - 1: + raise ValidationError("'plasma' entry must have at least one entry before and after it") + if radial_build[plasma_index - 1][0] != "gap" or radial_build[plasma_index + 1][0] != "gap": + raise ValidationError("'plasma' entry must be preceded and followed by a 'gap'") + + +def is_lower_or_upper_divertor(radial_build): + for item in radial_build: + if item[0] == "lower_divertor": + return "lower_divertor" + if item[0] == "upper_divertor": + return "upper_divertor" + raise ValidationError("neither upper_divertor or lower_divertor found") + + +def get_plasma_value(radial_build): + for item in radial_build: + if item[0] == "plasma": + return item[1] + raise ValueError("'plasma' entry not found") + + +def get_plasma_index(radial_build): + for i, item in enumerate(radial_build): + if item[0] == "plasma": + return i + raise ValueError("'plasma' entry not found") + + +def get_gap_after_plasma(radial_build): + for index, item in enumerate(radial_build): + if item[0] == "plasma": + if index + 1 < len(radial_build) and radial_build[index + 1][0] == "gap": + return radial_build[index + 1][1] + else: + raise ValueError("'plasma' entry is not followed by a 'gap'") + raise ValueError("'plasma' entry not found") + + +def sum_after_gap_following_plasma(radial_build): + found_plasma = False + found_gap_after_plasma = False + total_sum = 0 + + for item in radial_build: + if found_gap_after_plasma: + total_sum += item[1] + elif found_plasma and item[0] == "gap": + found_gap_after_plasma = True + elif item[0] == "plasma": + found_plasma = True + + if not found_plasma: + raise ValueError("'plasma' entry not found") + if not found_gap_after_plasma: + raise ValueError("'plasma' entry is not followed by a 'gap'") + + return total_sum diff --git a/src/paramak/workplanes/__init__.py b/src/paramak/workplanes/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/paramak/workplanes/blanket_constant_thickness_arc_h.py b/src/paramak/workplanes/blanket_constant_thickness_arc_h.py new file mode 100644 index 000000000..fe8147a70 --- /dev/null +++ b/src/paramak/workplanes/blanket_constant_thickness_arc_h.py @@ -0,0 +1,52 @@ +import typing + +import cadquery as cq + +from ..utils import create_wire_workplane_from_points + + +def blanket_constant_thickness_arc_h( + inner_mid_point: typing.Tuple[float, float], + inner_upper_point: typing.Tuple[float, float], + inner_lower_point: typing.Tuple[float, float], + thickness: float, + rotation_angle=90, + plane="XZ", + origin=(0, 0, 0), + obj=None, + color: typing.Tuple[float, float, float, typing.Optional[float]] = ( + 0.0, + 0.333, + 0.0, + ), + name="blanket_constant_thickness_arc_h", +): + points = [ + (inner_upper_point[0], inner_upper_point[1], "circle"), + (inner_mid_point[0], inner_mid_point[1], "circle"), + (inner_lower_point[0], inner_lower_point[1], "straight"), + ( + inner_lower_point[0] + abs(thickness), + inner_lower_point[1], + "circle", + ), + ( + inner_mid_point[0] + abs(thickness), + inner_mid_point[1], + "circle", + ), + ( + inner_upper_point[0] + abs(thickness), + inner_upper_point[1], + "straight", + ), + ] + + points.append(points[0]) + + wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj) + + solid = wire.revolve(rotation_angle) + solid.name = name + solid.color = cq.Color(*color) + return solid diff --git a/src/paramak/workplanes/blanket_from_plasma.py b/src/paramak/workplanes/blanket_from_plasma.py new file mode 100644 index 000000000..71049e341 --- /dev/null +++ b/src/paramak/workplanes/blanket_from_plasma.py @@ -0,0 +1,303 @@ +import warnings +import typing + +from ..utils import create_wire_workplane_from_points +import mpmath +import numpy as np +import sympy as sp +from scipy.interpolate import interp1d + + +def make_callable(attribute, start_angle, stop_angle): + """This function transforms an attribute (thickness or offset) into a + callable function of theta + """ + # if the attribute is a list, create a interpolated object of the + # values + if isinstance(attribute, (tuple, list)): + if isinstance(attribute[0], (tuple, list)) and isinstance(attribute[1], (tuple, list)) and len(attribute) == 2: + # attribute is a list of 2 lists + if len(attribute[0]) != len(attribute[1]): + raise ValueError( + "The length of angles list must equal \ + the length of values list" + ) + list_of_angles = np.array(attribute[0]) + offset_values = attribute[1] + else: + # no list of angles is given + offset_values = attribute + list_of_angles = np.linspace(start_angle, stop_angle, len(offset_values), endpoint=True) + interpolated_values = interp1d(list_of_angles, offset_values) + + def fun(theta): + if callable(attribute): + return attribute(theta) + elif isinstance(attribute, (tuple, list)): + return interpolated_values(theta) + else: + return attribute + + return fun + + +def find_points( + start_angle, + stop_angle, + offset_from_plasma, + major_radius, + minor_radius, + triangularity, + elongation, + vertical_displacement, + thickness, + connect_to_center, + num_points, + allow_overlapping_shape, + angles=None, +): + # create array of angles theta + if angles is None: + thetas = np.linspace( + start_angle, + stop_angle, + num=num_points, + endpoint=True, + ) + else: + thetas = angles + + # create inner points + inner_offset = make_callable(offset_from_plasma, start_angle, stop_angle) + inner_points, overlapping_shape = create_offset_points( + major_radius=major_radius, + minor_radius=minor_radius, + triangularity=triangularity, + elongation=elongation, + vertical_displacement=vertical_displacement, + thetas=thetas, + offset=inner_offset, + ) + inner_points[-1][2] = "straight" + + # create outer points + thickness = make_callable(thickness, start_angle, stop_angle) + + def outer_offset(theta): + return inner_offset(theta) + thickness(theta) + + outer_points, overlapping_shape = create_offset_points( + major_radius=major_radius, + minor_radius=minor_radius, + triangularity=triangularity, + elongation=elongation, + vertical_displacement=vertical_displacement, + thetas=np.flip(thetas), + offset=outer_offset, + ) + + outer_points[-1][2] = "straight" + + if connect_to_center: + inner_points.append([0, inner_points[-1][1], "straight"]) + inner_points = [[0, inner_points[0][1], "straight"]] + inner_points + outer_points.append([0, outer_points[-1][1], "straight"]) + outer_points = [[0, outer_points[0][1], "straight"]] + outer_points + + # assemble + points = inner_points + outer_points + print(overlapping_shape, allow_overlapping_shape) + if overlapping_shape and allow_overlapping_shape is False: + msg = "blanket_from_plasma: Some points with negative R coordinate have " "been ignored." + warnings.warn(msg, category=UserWarning) + + # input() + return points + + +def create_offset_points( + major_radius: float, + minor_radius: float, + triangularity: float, + elongation: float, + vertical_displacement, + thetas, + offset, +): + """generates a list of points following parametric equations with an + offset + + Args: + thetas (np.array): the angles in degrees. + offset (callable): offset value (cm). offset=0 will follow the + parametric equations. + + Returns: + list: list of points [[R1, Z1, connection1], [R2, Z2, connection2], + ...] + """ + # create sympy objects and derivatives + + overlapping_shape = False + theta_sp = sp.Symbol("theta") + + R_sp, Z_sp = distribution( + major_radius, + minor_radius, + triangularity, + elongation, + vertical_displacement, + theta_sp, + pkg=sp, + ) + R_derivative = sp.diff(R_sp, theta_sp) + Z_derivative = sp.diff(Z_sp, theta_sp) + points = [] + + for theta in thetas: + # get local value of derivatives + val_R_derivative = float(R_derivative.subs("theta", theta)) + val_Z_derivative = float(Z_derivative.subs("theta", theta)) + + # get normal vector components + nx = val_Z_derivative + ny = -val_R_derivative + + # normalise normal vector + normal_vector_norm = (nx**2 + ny**2) ** 0.5 + nx /= normal_vector_norm + ny /= normal_vector_norm + + # calculate outer points + val_R_outer = ( + distribution( + major_radius, + minor_radius, + triangularity, + elongation, + vertical_displacement, + theta, + )[0] + + offset(theta) * nx + ) + val_Z_outer = ( + distribution( + major_radius, + minor_radius, + triangularity, + elongation, + vertical_displacement, + theta, + )[1] + + offset(theta) * ny + ) + if float(val_R_outer) > 0: + points.append([float(val_R_outer), float(val_Z_outer), "spline"]) + else: + overlapping_shape = True + return points, overlapping_shape + + +def distribution(major_radius, minor_radius, triangularity, elongation, vertical_displacement, theta, pkg=np): + """Plasma distribution theta in degrees + + Args: + theta (float or np.array or sp.Symbol): the angle(s) in degrees. + pkg (module, optional): Module to use in the funciton. If sp, as + sympy object will be returned. If np, a np.array or a float + will be returned. Defaults to np. + + Returns: + (float, float) or (sympy.Add, sympy.Mul) or + (numpy.array, numpy.array): The R and Z coordinates of the + point with angle theta + """ + if pkg == np: + theta = np.radians(theta) + else: + theta = mpmath.radians(theta) + R = major_radius + minor_radius * pkg.cos(theta + triangularity * pkg.sin(theta)) + Z = elongation * minor_radius * pkg.sin(theta) + vertical_displacement + return R, Z + + +def blanket_from_plasma( + thickness, + start_angle: float, + stop_angle: float, + minor_radius: float = 150.0, + major_radius: float = 450.0, + triangularity: float = 0.55, + elongation: float = 2.0, + vertical_displacement: float = 0.0, + offset_from_plasma: typing.Union[float, typing.Iterable[float]] = 0.0, + num_points: int = 50, + name: str = "blanket_from_plasma", + color: typing.Tuple[float, float, float, typing.Optional[float]] = ( + 0.333, + 0.0, + 0.0, + ), + rotation_angle: float = 90.0, + plane="XZ", + origin=(0, 0, 0), + obj=None, + allow_overlapping_shape=False, + connect_to_center=False, +): + """A blanket volume created from plasma parameters. In might be nessecary + to increase the num_points when making long but thin geometry with this + component. + + Args: + thickness (float or [float] or callable or [(float), (float)]): + the thickness of the blanket (cm). If the thickness is a float then + this produces a blanket of constant thickness. If the thickness is + a tuple of floats, blanket thickness will vary linearly between the + two values. If thickness is callable, then the blanket thickness + will be a function of poloidal angle (in degrees). If thickness is + a list of two lists (thicknesses and angles) then these will be + used together with linear interpolation. + start_angle: the angle in degrees to start the blanket, measured anti + clockwise from 3 o'clock. + stop_angle: the angle in degrees to stop the blanket, measured anti + clockwise from 3 o'clock. + plasma: If not None, the parameters of the plasma Object will be used. + minor_radius: the minor radius of the plasma (cm). + major_radius: the major radius of the plasma (cm). + triangularity: the triangularity of the plasma. + elongation: the elongation of the plasma. + vertical_displacement: the vertical_displacement of the plasma (cm). + offset_from_plasma: the distance between the plasma and the blanket + (cm). If float, constant offset. If list of floats, offset will + vary linearly between the values. If callable, offset will be a + function of poloidal angle (in degrees). If a list of two lists + (angles and offsets) then these will be used together with linear + interpolation. + num_points: number of points that will describe the shape. + allow_overlapping_shape: allows parameters to create a shape that + overlaps itself. + """ + + points = find_points( + thickness=thickness, + start_angle=start_angle, + stop_angle=stop_angle, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + vertical_displacement=vertical_displacement, + offset_from_plasma=offset_from_plasma, + num_points=num_points, + allow_overlapping_shape=allow_overlapping_shape, + connect_to_center=connect_to_center, + ) + points.append(points[0]) + + wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj) + + solid = wire.revolve(rotation_angle) + solid.name = name + solid.color = color + return solid diff --git a/src/paramak/workplanes/center_column_shield_cylinder.py b/src/paramak/workplanes/center_column_shield_cylinder.py new file mode 100644 index 000000000..ce742fea5 --- /dev/null +++ b/src/paramak/workplanes/center_column_shield_cylinder.py @@ -0,0 +1,60 @@ +import typing + +from ..utils import create_wire_workplane_from_points + + +def center_column_shield_cylinder( + height: float, + inner_radius: float, + thickness: float, + reference_point: tuple = ("center", 0), + name: str = "center_column_shield_cylinder", + color: typing.Tuple[float, float, float, typing.Optional[float]] = ( + 0.0, + 0.333, + 0.0, + ), + rotation_angle=90, + plane="XZ", + origin=(0, 0, 0), + obj=None, +): + """A cylindrical center column shield volume with constant thickness. + + Args: + height: height of the center column shield. + inner_radius: the inner radius of the center column shield. + thickness: the outer radius of the center column shield. + reference_point: the vertical coordinates to build te vessel from and + description of the reference point. Can be either the 'center' + with a numerical value or 'lower' with a numerical value. + """ + + outer_radius = inner_radius + thickness + + if reference_point[0] == "center": + center_height = reference_point[1] + elif reference_point[0] == "lower": + center_height = reference_point[1] + 0.5 * height + else: + raise ValueError('reference_point should be a tuple where the first value is either "center" or "lower"') + + if not isinstance(center_height, (int, float)): + msg = f"center_height should be a float or int. Not a {type(center_height)}" + raise TypeError(msg) + + points = [ + (inner_radius, center_height + height / 2, "straight"), + (outer_radius, center_height + height / 2, "straight"), + (outer_radius, center_height + (-height / 2), "straight"), + (inner_radius, center_height + (-height / 2), "straight"), + ] + + points.append(points[0]) + + wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj) + + solid = wire.revolve(rotation_angle) + solid.name = name + solid.color = color + return solid diff --git a/src/paramak/workplanes/constant_thickness_dome.py b/src/paramak/workplanes/constant_thickness_dome.py new file mode 100644 index 000000000..7e7737ed4 --- /dev/null +++ b/src/paramak/workplanes/constant_thickness_dome.py @@ -0,0 +1,190 @@ +import math +import numbers +import typing + +import cadquery as cq + +from ..utils import create_wire_workplane_from_points +from ..workplanes.cutting_wedge import cutting_wedge + + +def constant_thickness_dome( + thickness: float = 10, + chord_center_height: float = 0, + chord_width: float = 100, + chord_height: float = 20, + upper_or_lower: str = "upper", + name: str = "constant_thickness_dome", + plane="XZ", + color: typing.Tuple[float, float, float, typing.Optional[float]] = ( + 0.0, + 0.333, + 0.0, + ), + origin=(0, 0, 0), + rotation_angle=90, + obj=None, + **kwargs, +): + """A cylindrical vessel volume with constant thickness with a simple dished + head. This style of tank head has no knuckle radius or straight flange. The + dished shape is made from a chord of a circle. + + Arguments: + thickness: the radial thickness of the dome. + chord_center_height: the vertical position of the chord center + chord_width: the width of the chord base + chord_height: the height of the chord which is also distance between + the chord_center_height and the inner surface of the dome + upper_or_lower: Curves the dish with a positive or negative direction + to allow the upper section or lower section of vacuum vessel + domes to be made. + name: the name of the shape, used in the graph legend and as a + filename prefix when exporting. + """ + + if not isinstance(chord_width, numbers.Number): + raise ValueError("ConstantThicknessDome.chord_width must be a float. Not", chord_width) + if chord_width <= 0: + msg = f"ConstantThicknessDome.chord_width must be a positive number above 0. Not {chord_width}" + raise ValueError(msg) + + if not isinstance(chord_height, numbers.Number): + raise ValueError("ConstantThicknessDome.chord_height must be a float. Not", chord_height) + if chord_height <= 0: + msg = f"ConstantThicknessDome.chord_height must be a positive number above 0. Not {chord_height}" + raise ValueError(msg) + + if not isinstance(thickness, numbers.Number): + msg = f"VacuumVessel.thickness must be a float. Not {thickness}" + raise ValueError(msg) + if thickness <= 0: + msg = f"VacuumVessel.thickness must be a positive number above 0. Not {thickness}" + raise ValueError(msg) + + # Note these points are not used in the normal way when constructing + # the solid + # + # 6 - + # | - + # 7 - - + # - - + # - 3 + # - | + # cc 1 -- 2 + # chord center + # + # + # cp + # center point + # + # + # + # + # cc 1 -- 2 + # - | + # - 3 + # - - + # 7 - - + # | - + # 6 - + # far side + + if chord_height * 2 >= chord_width: + msg = "ConstantThicknessDome requires that the chord_width " "is at least 2 times as large as the chord height" + raise ValueError(msg) + + radius_of_sphere = ((math.pow(chord_width, 2)) + (4.0 * math.pow(chord_height, 2))) / (8 * chord_height) + + # TODO set to 0 for now, add ability to shift the center of the chord left and right + chord_center = (0, chord_center_height) + + point_1 = (chord_center[0] + (chord_width / 2), chord_center[1], "straight") + + if upper_or_lower == "upper": + center_point = (chord_center[0], chord_center[1] + chord_height - radius_of_sphere) + inner_tri_angle = math.atan((center_point[1] - chord_center[1]) / (chord_width / 2)) + outer_tri_adj = math.cos(inner_tri_angle) * thickness + point_2 = (point_1[0] + thickness, point_1[1], "straight") + outer_tri_opp = math.sqrt(math.pow(thickness, 2) - math.pow(outer_tri_adj, 2)) + point_7 = (chord_center[0], chord_center[1] + radius_of_sphere, "straight") + point_6 = (chord_center[0], chord_center[1] + radius_of_sphere + thickness, "straight") + far_side = (center_point[0], center_point[1] - (radius_of_sphere + thickness)) + point_3 = (point_2[0], point_2[1] + outer_tri_opp, "straight") + elif upper_or_lower == "lower": + center_point = (chord_center[0], chord_center[1] - chord_height + radius_of_sphere) + inner_tri_angle = math.atan((center_point[1] - chord_center[1]) / (chord_width / 2)) + outer_tri_adj = math.cos(inner_tri_angle) * thickness + point_2 = (point_1[0] + thickness, point_1[1], "straight") + outer_tri_opp = math.sqrt(math.pow(thickness, 2) - math.pow(outer_tri_adj, 2)) + point_7 = (chord_center[0], chord_center[1] - radius_of_sphere, "straight") + point_6 = (chord_center[0], chord_center[1] - (radius_of_sphere + thickness), "straight") + far_side = (center_point[0], center_point[1] + radius_of_sphere + thickness) + point_3 = (point_2[0], point_2[1] - outer_tri_opp, "straight") + else: + msg = f'upper_or_lower should be either "upper" or "lower". Not {upper_or_lower}' + raise ValueError(msg) + + points = [point_1, point_2, point_3, point_6, point_7] + + radius_of_sphere = ((math.pow(chord_width, 2)) + (4.0 * math.pow(chord_height, 2))) / (8 * chord_height) + + # TODO set to 0 for now, add ability to shift the center of the chord left and right + chord_center = (0, chord_center_height) + + # far_side is center on x and + if upper_or_lower == "upper": + center_point = (chord_center[0], chord_center[1] + chord_height - radius_of_sphere) + far_side = (center_point[0], center_point[1] - (radius_of_sphere + thickness)) + elif upper_or_lower == "lower": + center_point = (chord_center[0], chord_center[1] - chord_height + radius_of_sphere) + far_side = (center_point[0], center_point[1] + radius_of_sphere + thickness) + else: + raise ValueError("upper_or_lower argument must be set to either 'upper' or 'lower'") + + big_sphere = cq.Workplane(plane).moveTo(center_point[0], center_point[1]).sphere(radius_of_sphere + thickness) + + small_sphere = cq.Workplane(plane).moveTo(center_point[0], center_point[1]).sphere(radius_of_sphere) + + wire = create_wire_workplane_from_points( + points=( + (chord_center[0], chord_center[1], "straight"), + (chord_center[0] + radius_of_sphere + thickness, chord_center[1], "straight"), + (chord_center[0] + radius_of_sphere + thickness, far_side[1], "straight"), + (chord_center[0], far_side[1], "straight"), + (chord_center[0], chord_center[1], "straight"), + ), + plane=plane, + origin=origin, + obj=obj, + ) + inner_cylinder_cutter = wire.revolve(360) + + wire = create_wire_workplane_from_points( + points=( + (chord_center[0], chord_center[1], "straight"), # cc + (points[1][0], points[1][1], "straight"), # point 2 + (points[2][0], points[2][1], "straight"), # point 3 + (points[2][0] + radius_of_sphere, points[2][1], "straight"), # point 3 wider + (points[2][0] + radius_of_sphere, far_side[1], "straight"), + (far_side[0], far_side[1], "straight"), + (chord_center[0], far_side[1], "straight"), + ), + plane=plane, + origin=origin, + obj=obj, + ) + outer_cylinder_cutter = wire.revolve(360) + + cap = big_sphere.cut(small_sphere) + + height = 2 * (radius_of_sphere + abs(center_point[1]) + thickness) + radius = 2 * (radius_of_sphere + abs(center_point[0]) + thickness) + cutter = cutting_wedge(height=height, radius=radius, rotation_angle=rotation_angle, plane=plane) + + cap = cap.cut(outer_cylinder_cutter).cut(inner_cylinder_cutter) + cap = cap.intersect(cutter) + + cap.name = name + cap.color = cq.Color(*color) + return cap diff --git a/src/paramak/workplanes/cutting_wedge.py b/src/paramak/workplanes/cutting_wedge.py new file mode 100644 index 000000000..8a96c5706 --- /dev/null +++ b/src/paramak/workplanes/cutting_wedge.py @@ -0,0 +1,44 @@ +import typing + +import cadquery as cq + +from ..utils import create_wire_workplane_from_points + + +def cutting_wedge( + height: float, + radius: float, + rotation_angle: float = 180.0, + plane="XZ", + origin=(0, 0, 0), + obj=None, + color: typing.Tuple[float, float, float, typing.Optional[float]] = ( + 0.0, + 0.333, + 0.0, + ), + name="cutting_wedge", +): + """Creates a wedge from height, radius and rotation angle arguments than + can be useful for cutting sector models. + + Args: + height: the vertical (z axis) height of the coil (cm). + radius: the horizontal (x axis) width of the coil (cm). + rotation_angle: Defaults to 180.0. + """ + + points = [ + (0, height / 2, "straight"), + (radius, height / 2, "straight"), + (radius, -height / 2, "straight"), + (0, -height / 2, "straight"), + ] + points.append(points[0]) + + wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj) + + solid = wire.revolve(rotation_angle) + solid.name = name + solid.color = cq.Color(*color) + return solid diff --git a/src/paramak/workplanes/dished_vacuum_vessel.py b/src/paramak/workplanes/dished_vacuum_vessel.py new file mode 100644 index 000000000..5056f4fc5 --- /dev/null +++ b/src/paramak/workplanes/dished_vacuum_vessel.py @@ -0,0 +1,112 @@ +import typing + +from paramak import center_column_shield_cylinder, constant_thickness_dome + + +def dished_vacuum_vessel( + radius: float = 300, + reference_point: tuple = ("center", 0), + dish_height: typing.Tuple[float, float] = (20, 50), + cylinder_height: float = 400, + thickness: float = 15, + rotation_angle: float = 90, + name: str = "dished_vessel", + plane="XZ", +): + """A cylindrical vessel volume with constant thickness with a simple dished + head. This style of tank head has no knuckle radius or straight flange. + + Arguments: + radius: the radius from which the centres of the vessel meets the outer + circumference. + reference_point: the x,z coordinates to build te vessel from. Can be + either the 'center' with a value or 'lower' with a + value. For example + dish_height: the height of the lower and upper dish sections. + cylinder_height: the height of the cylindrical section of the vacuum + vessel. + thickness: the radial thickness of the vessel in cm. + """ + + if not isinstance(radius, (float, int)): + raise ValueError(f"radius must be a number. Not {type(radius)}") + if radius <= 0: + msg = "radius must be a positive number above 0. " f"Not {radius}" + raise ValueError(msg) + + if not isinstance(thickness, (float, int)): + msg = f"VacuumVessel.thickness must be a number. Not {type(thickness)}" + raise ValueError(msg) + if thickness <= 0: + msg = f"VacuumVessel.thickness must be a positive number above 0. Not {value}" + raise ValueError(msg) + + # + # - - + # - + # - - - + # - - + # - - + # - | + # | | + # | | + # | | + # c,p | | + # | | + # | | + # | | + # - | + # - - + # - - + # - - - + # - + # - - + # + + if reference_point[0] == "center": + center_height = reference_point[1] + lower_chord_center_height = reference_point[1] - 0.5 * cylinder_height + upper_chord_center_height = reference_point[1] + 0.5 * cylinder_height + elif reference_point[0] == "lower": + center_height = reference_point[1] + thickness + dish_height[0] + 0.5 * cylinder_height + lower_chord_center_height = reference_point[1] + thickness + dish_height[0] + upper_chord_center_height = reference_point[1] + thickness + dish_height[0] + cylinder_height + else: + raise ValueError('reference_point should be a tuple where the first value is either "center" or "lower"') + + cylinder_section = center_column_shield_cylinder( + height=cylinder_height, + inner_radius=radius - thickness, + thickness=thickness, + reference_point=("center", center_height), + rotation_angle=rotation_angle, + plane=plane, + ) + + upper_dome_section = constant_thickness_dome( + thickness=thickness, + chord_center_height=upper_chord_center_height, + chord_width=(radius - thickness) * 2, + chord_height=dish_height[1], + upper_or_lower="upper", + rotation_angle=rotation_angle, + plane=plane, + ) + + lower_dome_section = constant_thickness_dome( + thickness=thickness, + chord_center_height=lower_chord_center_height, + chord_width=(radius - thickness) * 2, + chord_height=dish_height[0], + upper_or_lower="lower", + rotation_angle=rotation_angle, + plane=plane, + ) + + # union can fail, safer to return 3 workplanes + # solid = cylinder_section.union(upper_dome_section).union(lower_dome_section) + + cylinder_section.name = name + upper_dome_section.name = name + lower_dome_section.name = name + return lower_dome_section, cylinder_section, upper_dome_section diff --git a/src/paramak/workplanes/plasma_simplified.py b/src/paramak/workplanes/plasma_simplified.py new file mode 100644 index 000000000..1ebcf62b4 --- /dev/null +++ b/src/paramak/workplanes/plasma_simplified.py @@ -0,0 +1,64 @@ +import typing + +import numpy as np + +from ..utils import create_wire_workplane_from_points + + +def plasma_simplified( + elongation: float = 2.0, + major_radius: float = 450.0, + minor_radius: float = 150.0, + triangularity: float = 0.55, + vertical_displacement: float = 0.0, + num_points: float = 50, + name: str = "tokamak_plasma", + color: typing.Tuple[float, float, float, typing.Optional[float]] = ( + 0.333, + 0.0, + 0.0, + ), + rotation_angle=90, + plane="XZ", + origin=(0, 0, 0), + obj=None, +): + """Creates a double null tokamak plasma shape that is controlled by 4 + shaping parameters. + + Args: + elongation: the elongation of the plasma. + major_radius: the major radius of the plasma (cm). + minor_radius: the minor radius of the plasma (cm). + triangularity: the triangularity of the plasma. + vertical_displacement: the vertical_displacement of the plasma (cm).. + num_points: number of points to describe the shape. + """ + + # create array of angles theta + theta = np.linspace(0, 2 * np.pi, num=num_points, endpoint=False) + + # parametric equations for plasma + def R(theta): + return major_radius + minor_radius * np.cos(theta + triangularity * np.sin(theta)) + + def Z(theta): + return elongation * minor_radius * np.sin(theta) + vertical_displacement + + points = np.stack((R(theta), Z(theta)), axis=1).tolist() + points.append(points[0]) + for point in points: + point.append("spline") + + wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj) + + # avoids shape with surface on join that can't be meshed for 360 degree plasmas + if rotation_angle >= 360: + solid1 = wire.revolve(180, (1, 0, 0), (1, 1, 0)) + solid2 = solid1.mirror(solid1.faces(">X"), union=True) + solid = solid2.union(solid1) # todo try fuzzy bool tol=0.01 + else: + solid = wire.revolve(rotation_angle) + solid.name = name + solid.color = color + return solid diff --git a/src/paramak/workplanes/poloidal_field_coil.py b/src/paramak/workplanes/poloidal_field_coil.py new file mode 100644 index 000000000..211439e1a --- /dev/null +++ b/src/paramak/workplanes/poloidal_field_coil.py @@ -0,0 +1,43 @@ +import typing + +from ..utils import create_wire_workplane_from_points + + +def poloidal_field_coil( + height: float, + width: float, + center_point: float, + name: str = "poloidal_field_coil", + color: typing.Tuple[float, float, float, typing.Optional[float]] = ( + 0.0, + 0.333, + 0.0, + ), + rotation_angle=90, + plane="XZ", + origin=(0, 0, 0), + obj=None, +): + """A rectangular poloidal field coil. + + Args: + height: the vertical (z axis) height of the coil. + width: the horizontal (x axis) width of the coil. + center_point: the center of the coil (x,z) values. + """ + + points = [ + (center_point[0] + width / 2.0, center_point[1] + height / 2.0, "straight"), # upper right + (center_point[0] + width / 2.0, center_point[1] - height / 2.0, "straight"), # lower right + (center_point[0] - width / 2.0, center_point[1] - height / 2.0, "straight"), # lower left + (center_point[0] - width / 2.0, center_point[1] + height / 2.0, "straight"), + ] + + points.append(points[0]) + + wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj) + + solid = wire.revolve(rotation_angle) + solid.name = name + solid.color = color + return solid diff --git a/src/paramak/workplanes/poloidal_field_coil_case.py b/src/paramak/workplanes/poloidal_field_coil_case.py new file mode 100644 index 000000000..626c66187 --- /dev/null +++ b/src/paramak/workplanes/poloidal_field_coil_case.py @@ -0,0 +1,66 @@ +import typing + +from ..utils import create_wire_workplane_from_points + + +def poloidal_field_coil_case( + coil_height: float, + coil_width: float, + casing_thickness: typing.Tuple[float, float], + center_point: typing.Tuple[float, float], + name: str = "poloidal_field_coil_case", + color: typing.Tuple[float, float, float, typing.Optional[float]] = (1.0, 1.0, 0.498), + rotation_angle=90, + plane="XZ", + origin=(0, 0, 0), + obj=None, +): + """Constant thickness casing for a rectangular poloidal coil. + + Args: + coil_height: the vertical (z axis) height of the coil (cm). + coil_width: the horizontal (x axis) width of the coil (cm). + center_point: the center of the coil (x,z) values (cm). + casing_thickness: the thickness of the coil casing (cm). + """ + + inner_points = [ + (center_point[0] + coil_width / 2.0, center_point[1] + coil_height / 2.0, "straight"), # upper right + (center_point[0] + coil_width / 2.0, center_point[1] - coil_height / 2.0, "straight"), # lower right + (center_point[0] - coil_width / 2.0, center_point[1] - coil_height / 2.0, "straight"), # lower left + (center_point[0] - coil_width / 2.0, center_point[1] + coil_height / 2.0, "straight"), # upper left + ] + inner_points.append(inner_points[0]) + + outer_points = [ + ( + center_point[0] + (casing_thickness + coil_width / 2.0), + center_point[1] + (casing_thickness + coil_height / 2.0), + "straight", + ), + ( + center_point[0] + (casing_thickness + coil_width / 2.0), + center_point[1] - (casing_thickness + coil_height / 2.0), + "straight", + ), + ( + center_point[0] - (casing_thickness + coil_width / 2.0), + center_point[1] - (casing_thickness + coil_height / 2.0), + "straight", + ), + ( + center_point[0] - (casing_thickness + coil_width / 2.0), + center_point[1] + (casing_thickness + coil_height / 2.0), + "straight", + ), + ] + outer_points.append(outer_points[0]) + + inner_wire = create_wire_workplane_from_points(points=inner_points, plane=plane, origin=origin, obj=obj) + outer_wire = create_wire_workplane_from_points(points=outer_points, plane=plane, origin=origin, obj=obj) + + inner_solid = inner_wire.revolve(rotation_angle) + solid = outer_wire.revolve(rotation_angle).cut(inner_solid) + solid.name = name + solid.color = color + return solid diff --git a/src/paramak/workplanes/toroidal_field_coil_rectangle.py b/src/paramak/workplanes/toroidal_field_coil_rectangle.py new file mode 100644 index 000000000..61c46977c --- /dev/null +++ b/src/paramak/workplanes/toroidal_field_coil_rectangle.py @@ -0,0 +1,106 @@ +import typing + +from ..utils import create_wire_workplane_from_points, rotate_solid + + +def toroidal_field_coil_rectangle( + horizontal_start_point: typing.Tuple[float, float] = (20, 200), + vertical_mid_point: typing.Tuple[float, float] = (350, 0), + thickness: float = 30, + distance: float = 20, + name: str = "toroidal_field_coil", + with_inner_leg: bool = True, + azimuthal_placement_angles: typing.Sequence[float] = [0, 90, 180], + vertical_displacement: float = 0.0, + color: typing.Tuple[float, float, float, typing.Optional[float]] = (0.0, 0.0, 1.0), + plane: str = "XZ", + origin: typing.Tuple[float, float, float] = (0.0, 0.0, 0.0), + obj=None, +): + """Creates a rectangular shaped toroidal field coil. + + Args: + horizontal_start_point: the (x,z) coordinates of the inner upper + point (cm). + vertical_mid_point: the (x,z) coordinates of the mid point of the + vertical section (cm). + thickness: the thickness of the toroidal field coil. + distance: the extrusion distance. + number_of_coils: the number of tf coils. This changes by the + azimuth_placement_angle dividing up 360 degrees by the number of + coils. + with_inner_leg: include the inner tf leg. Defaults to True. + azimuth_start_angle: The azimuth angle to for the first TF coil which + offsets the placement of coils around the azimuthal angle + """ + + if horizontal_start_point[0] >= vertical_mid_point[0]: + raise ValueError( + "horizontal_start_point x should be smaller than the \ + vertical_mid_point x value" + ) + if vertical_mid_point[1] >= horizontal_start_point[1]: + raise ValueError( + "vertical_mid_point y value should be smaller than the \ + horizontal_start_point y value" + ) + + points = [ + horizontal_start_point, # connection point + ( + horizontal_start_point[0] + thickness, + horizontal_start_point[1], + ), + (vertical_mid_point[0], horizontal_start_point[1]), + (vertical_mid_point[0], -horizontal_start_point[1]), + # connection point + ( + horizontal_start_point[0] + thickness, + -horizontal_start_point[1], + ), + # connection point + (horizontal_start_point[0], -horizontal_start_point[1]), + ( + horizontal_start_point[0], + -(horizontal_start_point[1] + thickness), + ), + ( + vertical_mid_point[0] + thickness, + -(horizontal_start_point[1] + thickness), + ), + ( + vertical_mid_point[0] + thickness, + horizontal_start_point[1] + thickness, + ), + ( + horizontal_start_point[0], + horizontal_start_point[1] + thickness, + ), + horizontal_start_point, + ] + + # adds any vertical displacement and the connection type to the points + points = [(point[0], point[1] + vertical_displacement, "straight") for point in points] + + wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj) + solid = wire.extrude(until=distance / 2, both=True) + solid = rotate_solid(angles=azimuthal_placement_angles, solid=solid) + + if with_inner_leg: + inner_leg_connection_points = [ + (points[0][0], points[0][1], "straight"), + (points[1][0], points[1][1], "straight"), + (points[4][0], points[4][1], "straight"), + (points[5][0], points[5][1], "straight"), + (points[0][0], points[0][1], "straight"), + ] + inner_wire = create_wire_workplane_from_points( + points=inner_leg_connection_points, plane=plane, origin=origin, obj=obj + ) + inner_solid = inner_wire.extrude(until=distance / 2, both=True) + inner_solid = rotate_solid(angles=azimuthal_placement_angles, solid=inner_solid) + solid = solid.union(inner_solid) + + solid.name = name + solid.color = color + return solid diff --git a/src/paramak/workplanes/u_shaped_dome.py b/src/paramak/workplanes/u_shaped_dome.py new file mode 100644 index 000000000..0de535d55 --- /dev/null +++ b/src/paramak/workplanes/u_shaped_dome.py @@ -0,0 +1,105 @@ +import typing + +from paramak import center_column_shield_cylinder, constant_thickness_dome + +from ..utils import create_wire_workplane_from_points + + +def u_shaped_dome( + radius: float = 310, + reference_point: tuple = ("lower", 0), + dish_height: typing.Tuple[float, float] = 50, + cylinder_height: float = 400, + thickness: float = 16, + rotation_angle: float = 180, + name: str = "u_shaped_dome", + plane="XZ", + upper_or_lower="upper", +): + """A cylindrical u shaped dome with constant thickness. + + Arguments: + radius: the radius from which the centres of the vessel meets the outer + circumference. + reference_point: the x,z coordinates to build te vessel from. Can be + either the 'center' with a value or 'lower' with a + value. For example + dish_height: the height of the lower and upper dish sections. + cylinder_height: the height of the cylindrical section of the vacuum + vessel. + thickness: the radial thickness of the vessel in cm. + """ + + if not isinstance(radius, (float, int)): + raise ValueError(f"radius must be a number. Not {type(radius)}") + if radius <= 0: + msg = "radius must be a positive number above 0. " f"Not {radius}" + raise ValueError(msg) + + if not isinstance(thickness, (float, int)): + msg = f"VacuumVessel.thickness must be a number. Not {type(thickness)}" + raise ValueError(msg) + if thickness <= 0: + msg = f"VacuumVessel.thickness must be a positive number above 0. Not {value}" + raise ValueError(msg) + + # + # - - + # - + # - - - + # - - + # - - + # - | + # | | + # | | + # | | + # c,p | | + # | | + + if reference_point[0] == "center": + center_height = reference_point[1] + lower_chord_center_height = reference_point[1] - 0.5 * cylinder_height + upper_chord_center_height = reference_point[1] + 0.5 * cylinder_height + elif reference_point[0] == "lower": + center_height = reference_point[1] + thickness + dish_height + 0.5 * cylinder_height + lower_chord_center_height = reference_point[1] + thickness + dish_height + upper_chord_center_height = reference_point[1] + thickness + dish_height + cylinder_height + else: + raise ValueError('reference_point should be a tuple where the first value is either "center" or "lower"') + + cylinder_section = center_column_shield_cylinder( + height=cylinder_height, + inner_radius=radius - thickness, + thickness=thickness, + reference_point=("center", center_height), + rotation_angle=rotation_angle, + plane=plane, + ) + + if upper_or_lower == "upper": + dome_section = constant_thickness_dome( + thickness=thickness, + chord_center_height=upper_chord_center_height, + chord_width=(radius - thickness) * 2, + chord_height=dish_height, + upper_or_lower="upper", + rotation_angle=rotation_angle, + plane=plane, + ) + + elif upper_or_lower == "lower": + dome_section = constant_thickness_dome( + thickness=thickness, + chord_center_height=lower_chord_center_height, + chord_width=(radius - thickness) * 2, + chord_height=dish_height, + upper_or_lower="lower", + rotation_angle=rotation_angle, + plane=plane, + ) + else: + raise ValueError(f'upper_or_lower must be either "lower" or "upper" not {upper_or_lower}') + + dome_section.name = name + cylinder_section.name = name + return dome_section, cylinder_section diff --git a/tests/test_assemblies/__init__.py b/tests/test_assemblies/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/test_assemblies/test_spherical_tokamak.py b/tests/test_assemblies/test_spherical_tokamak.py new file mode 100644 index 000000000..bf926ff1d --- /dev/null +++ b/tests/test_assemblies/test_spherical_tokamak.py @@ -0,0 +1,102 @@ +from pathlib import Path + +import pytest +from cad_to_dagmc import CadToDagmc + +import paramak + +from .test_utils import transport_particles_on_h5m_geometry + + +@pytest.mark.parametrize("rotation_angle", [30, 360]) +def test_transport_with_magnets(rotation_angle): + poloidal_field_coils = [] + for case_thickness, height, width, center_point in zip( + [10, 15], + [20, 50], + [20, 50], + [(500, 300), (590, 100)], + ): + poloidal_field_coils.append( + paramak.poloidal_field_coil( + height=height, width=width, center_point=center_point, rotation_angle=rotation_angle + ) + ) + poloidal_field_coils.append( + paramak.poloidal_field_coil_case( + coil_height=height, + coil_width=width, + casing_thickness=case_thickness, + rotation_angle=rotation_angle, + center_point=center_point, + ) + ) + + my_reactor = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ("gap", 10), + ], + elongation=2, + triangularity=0.55, + rotation_angle=rotation_angle, + add_extra_cut_shapes=poloidal_field_coils, + ) + my_reactor.save(f"spherical_tokamak_with_magnets_{rotation_angle}.step") + assert Path(f"spherical_tokamak_with_magnets_{rotation_angle}.step").exists() + + my_model = CadToDagmc() + material_tags = ["mat1"] * 11 # rear wall is being split into 2 parts by the magnet that is cut out + my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) + my_model.export_dagmc_h5m_file(min_mesh_size=2, max_mesh_size=30.0) + + h5m_filename = "dagmc.h5m" + flux = transport_particles_on_h5m_geometry( + h5m_filename=h5m_filename, + material_tags=material_tags, + nuclides=["H1"] * len(material_tags), + cross_sections_xml="tests/cross_sections.xml", + ) + assert flux > 0.0 + + +def test_transport_without_magnets(): + reactor = paramak.spherical_tokamak_from_plasma( + radial_builds=[ + ("gap", 10), + ("layer", 50), + ("layer", 15), + ("gap", 50), + ("plasma", 300), + ("gap", 60), + ("layer", 15), + ("layer", 60), + ("layer", 10), + ("gap", 10), + ], + elongation=2, + triangularity=0.55, + ) + reactor.save("spherical_tokamak.step") + + my_model = CadToDagmc() + material_tags = ["mat1"] * 6 + my_model.add_cadquery_object(cadquery_object=reactor, material_tags=material_tags) + + my_model.export_dagmc_h5m_file(filename="dagmc.h5m", min_mesh_size=10.0, max_mesh_size=100.0) + + flux = transport_particles_on_h5m_geometry( + h5m_filename="dagmc.h5m", + material_tags=material_tags, + nuclides=["H1"] * len(material_tags), + cross_sections_xml="tests/cross_sections.xml", + ) + assert flux > 0.0 diff --git a/tests/test_assemblies/test_utils.py b/tests/test_assemblies/test_utils.py new file mode 100644 index 000000000..d9f691c0e --- /dev/null +++ b/tests/test_assemblies/test_utils.py @@ -0,0 +1,102 @@ +def transport_particles_on_h5m_geometry( + h5m_filename: str, + material_tags: list, + nuclides: list = None, + cross_sections_xml: str = None, +): + """A function for testing the geometry file with particle transport in + DAGMC OpenMC. Requires openmc and either the cross_sections_xml to be + specified. Returns the flux on each volume + + Arg: + h5m_filename: The name of the DAGMC h5m file to test + material_tags: The + nuclides: + cross_sections_xml: + + """ + + import openmc + from openmc.data import NATURAL_ABUNDANCE + + if nuclides is None: + nuclides = list(NATURAL_ABUNDANCE.keys()) + + materials = openmc.Materials() + for i, material_tag in enumerate(set(material_tags)): + # simplified material definitions have been used to keen this example minimal + mat_dag_material_tag = openmc.Material(name=material_tag) + mat_dag_material_tag.add_nuclide(nuclides[i], 1, "ao") + mat_dag_material_tag.set_density("g/cm3", 0.01) + + materials.append(mat_dag_material_tag) + + if cross_sections_xml: + openmc.config["cross_sections"] = cross_sections_xml + + else: + with open("cross_sections.xml", "w") as file: + file.write( + """ + + + + + + """ + ) + + openmc.config["cross_sections"] = "cross_sections.xml" + + dag_univ = openmc.DAGMCUniverse(filename=h5m_filename) + bound_dag_univ = dag_univ.bounded_universe() + geometry = openmc.Geometry(root=bound_dag_univ) + + # initializes a new source object + my_source = openmc.IndependentSource() + + center_of_geometry = ( + (dag_univ.bounding_box[0][0] + dag_univ.bounding_box[1][0]) / 2, + (dag_univ.bounding_box[0][1] + dag_univ.bounding_box[1][1]) / 2, + (dag_univ.bounding_box[0][2] + dag_univ.bounding_box[1][2]) / 2, + ) + # sets the location of the source which is not on a vertex + center_of_geometry_nudged = ( + center_of_geometry[0] + 0.1, + center_of_geometry[1] + 0.1, + center_of_geometry[2] + 0.1, + ) + + my_source.space = openmc.stats.Point(center_of_geometry_nudged) + # sets the direction to isotropic + my_source.angle = openmc.stats.Isotropic() + # sets the energy distribution to 100% 14MeV neutrons + my_source.energy = openmc.stats.Discrete([14e6], [1]) + + # specifies the simulation computational intensity + settings = openmc.Settings() + settings.batches = 10 + settings.particles = 10000 + settings.inactive = 0 + settings.run_mode = "fixed source" + settings.source = my_source + + # adds a tally to record the heat deposited in entire geometry + cell_tally = openmc.Tally(name="flux") + cell_tally.scores = ["flux"] + + # groups the two tallies + tallies = openmc.Tallies([cell_tally]) + + # builds the openmc model + my_model = openmc.Model(materials=materials, geometry=geometry, settings=settings, tallies=tallies) + + # starts the simulation + output_file = my_model.run() + + # loads up the output file from the simulation + statepoint = openmc.StatePoint(output_file) + + my_flux_cell_tally = statepoint.get_tally(name="flux") + + return my_flux_cell_tally.mean.flatten()[0] diff --git a/tests/test_utils.py b/tests/test_utils.py new file mode 100644 index 000000000..29ecfa202 --- /dev/null +++ b/tests/test_utils.py @@ -0,0 +1,1047 @@ +import pytest + +from paramak.utils import ( + ValidationError, + get_gap_after_plasma, + get_plasma_value, + sum_after_gap_following_plasma, + sum_up_to_plasma, + validate_divertor_radial_build, + validate_plasma_radial_build, +) + + +def test_validate_divertor_radial_build_valid(): + radial_build = [("gap", 10), ("lower_divertor", 20)] + assert validate_divertor_radial_build(radial_build) is None + + +def test_validate_divertor_radial_build_invalid_length(): + radial_build = [("gap", 10)] + with pytest.raises(ValidationError, match="should only contain two entries"): + validate_divertor_radial_build(radial_build) + + +def test_validate_divertor_radial_build_invalid_tuple_length(): + radial_build = [("gap", 10, 5), ("lower_divertor", 20)] + with pytest.raises(ValidationError, match="should only contain tuples of length 2"): + validate_divertor_radial_build(radial_build) + + +def test_validate_divertor_radial_build_invalid_second_entry(): + radial_build = [("gap", 10), ("divertor", 20)] + with pytest.raises(ValidationError, match='should be either "lower_divertor" or "upper_divertor"'): + validate_divertor_radial_build(radial_build) + + +def test_validate_divertor_radial_build_invalid_first_entry(): + radial_build = [("layer", 10), ("lower_divertor", 20)] + with pytest.raises(ValidationError, match='should be a "gap"'): + validate_divertor_radial_build(radial_build) + + +def test_validate_divertor_radial_build_non_positive_thickness(): + radial_build = [("gap", -10), ("lower_divertor", 20)] + with pytest.raises(ValidationError, match="should both be positive values"): + validate_divertor_radial_build(radial_build) + + +def test_validate_divertor_radial_build_invalid_thickness_type(): + radial_build = [("gap", "10"), ("lower_divertor", 20)] + with pytest.raises(ValidationError, match="should both be integers or floats"): + validate_divertor_radial_build(radial_build) + + +def test_get_plasma_value(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + assert get_plasma_value(radial_build) == 50 + + +def test_get_plasma_value_not_found(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValueError, match="'plasma' entry not found"): + get_plasma_value(radial_build) + + +def test_valid_case(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + validate_plasma_radial_build(radial_build) # Should not raise an error + + +def test_plasma_not_preceded_by_gap(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="'plasma' entry must be preceded and followed by a 'gap'"): + validate_plasma_radial_build(radial_build) + + +def test_plasma_not_followed_by_gap(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "gap", + 5, + ), + ( + "plasma", + 50, + ), + ( + "layer", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="'plasma' entry must be preceded and followed by a 'gap'"): + validate_plasma_radial_build(radial_build) + + +def test_missing_plasma(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="'plasma' entry not found or found multiple times"): + validate_plasma_radial_build(radial_build) + + +def test_multiple_plasma(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="Multiple 'plasma' entries found"): + validate_plasma_radial_build(radial_build) + + +def test_first_entry_not_string(): + radial_build = [ + ( + 10, + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="Invalid tuple structure at index 0"): + validate_plasma_radial_build(radial_build) + + +def test_second_entry_not_number(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + "50", + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="Invalid tuple structure at index 1"): + validate_plasma_radial_build(radial_build) + + +def test_invalid_string(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "invalid", + 5, + ), # Invalid string + ( + "gap", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="Invalid string 'invalid' at index 2"): + validate_plasma_radial_build(radial_build) + + +def test_plasma_first_entry(): + radial_build = [ + ( + "plasma", + 50, + ), + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="'plasma' entry must have at least one entry before and after it"): + validate_plasma_radial_build(radial_build) + + +def test_plasma_last_entry(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ( + "plasma", + 50, + ), + ] + with pytest.raises(ValidationError, match="'plasma' entry must have at least one entry before and after it"): + validate_plasma_radial_build(radial_build) + + +def test_non_positive_values(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 0, + ), # Non-positive value + ( + "gap", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValidationError, match="Non-positive value '0' at index 2"): + validate_plasma_radial_build(radial_build) + + +def test_sum_up_to_plasma_middle(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + assert sum_up_to_plasma(radial_build) == 115 + + +def test_sum_up_to_plasma_first(): + radial_build = [ + ( + "plasma", + 50, + ), + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + assert sum_up_to_plasma(radial_build) == 0 + + +def test_sum_up_to_plasma_last(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ( + "plasma", + 50, + ), + ] + assert sum_up_to_plasma(radial_build) == 191 + + +def test_sum_up_to_plasma_not_present(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "layer", + 5, + ), + ( + "gap", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + assert sum_up_to_plasma(radial_build) == 191 + + +def test_sum_up_to_plasma_empty(): + radial_build = [] + assert sum_up_to_plasma(radial_build) == 0 + + +def test_sum_up_to_plasma_multiple_entries(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 20, + ), + ( + "layer", + 30, + ), + ( + "gap", + 40, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + assert sum_up_to_plasma(radial_build) == 100 + + +def test_get_gap_after_plasma_not_followed_by_gap(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "gap", + 5, + ), + ( + "plasma", + 50, + ), + ( + "layer", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValueError, match="'plasma' entry is not followed by a 'gap'"): + get_gap_after_plasma(radial_build) + + +def test_get_gap_after_plasma_not_found(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "gap", + 5, + ), + ( + "layer", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValueError, match="'plasma' entry not found"): + get_gap_after_plasma(radial_build) + + +def test_sum_after_gap_following_plasma(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "gap", + 5, + ), + ( + "plasma", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + assert sum_after_gap_following_plasma(radial_build) == 16 + + +def test_sum_after_gap_following_plasma_not_found(): + radial_build = [ + ( + "gap", + 10, + ), + ( + "layer", + 50, + ), + ( + "gap", + 5, + ), + ( + "layer", + 50, + ), + ( + "gap", + 60, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "layer", + 2, + ), + ( + "gap", + 10, + ), + ] + with pytest.raises(ValueError, match="'plasma' entry not found"): + sum_after_gap_following_plasma(radial_build) diff --git a/tests/test_workplanes/__init__.py b/tests/test_workplanes/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/test_workplanes/test_blanket_constant_thickness_arc_h.py b/tests/test_workplanes/test_blanket_constant_thickness_arc_h.py new file mode 100644 index 000000000..d2f66c943 --- /dev/null +++ b/tests/test_workplanes/test_blanket_constant_thickness_arc_h.py @@ -0,0 +1,37 @@ +import math + +import paramak + + +def test_relative_shape_volume(): + """Creates two blankets using the blanket_constant_thickness_arc_h + parametric component and checks that their relative volumes and face areas + are correct.""" + + test_shape_360 = paramak.blanket_constant_thickness_arc_h( + inner_lower_point=(300, -200), + inner_mid_point=(500, 0), + inner_upper_point=(300, 200), + thickness=20, + rotation_angle=360, + ) + + test_shape_180 = paramak.blanket_constant_thickness_arc_h( + inner_lower_point=(300, -200), + inner_mid_point=(500, 0), + inner_upper_point=(300, 200), + thickness=20, + rotation_angle=180, + ) + + assert test_shape_360.val().Volume() > 1000 + assert math.isclose(test_shape_360.val().Volume(), 2 * test_shape_180.val().Volume()) + + assert len(test_shape_360.val().Faces()) == 4 + + areas = [face.Area() for face in test_shape_360.val().Faces()] + assert len(set([round(i) for i in areas])) == 3 + + areas = [face.Area() for face in test_shape_180.val().Faces()] + assert len(areas) == 6 + assert len(set([round(i) for i in areas])) == 4 diff --git a/tests/test_workplanes/test_blanket_from_plasma.py b/tests/test_workplanes/test_blanket_from_plasma.py new file mode 100644 index 000000000..8150fc6d8 --- /dev/null +++ b/tests/test_workplanes/test_blanket_from_plasma.py @@ -0,0 +1,154 @@ +import os +from pathlib import Path + +import pytest +from cadquery import exporters + +import paramak + +plasma = paramak.plasma_simplified( + major_radius=450, minor_radius=150, triangularity=0.55, elongation=2, rotation_angle=160 +) + +test_shape = paramak.blanket_from_plasma( + thickness=150, + start_angle=-90, + stop_angle=240, +) + + +def test_creation_plasma(): + """Checks that a cadquery solid can be created by passing a plasma to + the BlanketFP parametric component.""" + + assert test_shape.vals()[0].Volume() > 1000 + + +def test_faces(): + """creates a blanket using the BlanketFP parametric component and checks + that a solid with the correct number of faces is created""" + + test_shape = paramak.blanket_from_plasma(thickness=150, start_angle=-90, stop_angle=240, rotation_angle=360) + assert len(test_shape.vals()[0].Faces()) == 4 + + test_shape = paramak.blanket_from_plasma(thickness=150, start_angle=-90, stop_angle=240, rotation_angle=180) + assert len(test_shape.vals()[0].Faces()) == 6 + + +def test_creation_variable_thickness_from_tuple(): + """Checks that a cadquery solid can be created using the BlanketFP + parametric component when a tuple of thicknesses is passed as an + argument.""" + + test_shape = paramak.blanket_from_plasma(start_angle=-90, stop_angle=240, thickness=(100, 200)) + + assert test_shape.vals()[0].Volume() > 1000 + + +def test_creation_variable_thickness_from_2_lists(): + """Checks that a cadquery solid can be created using the BlanketFP + parametric component when a list of angles and a list of thicknesses + are passed as an argument.""" + + test_shape = paramak.blanket_from_plasma(start_angle=-90, stop_angle=240, thickness=[(-90, 240), [10, 30]]) + + assert test_shape is not None + + +def test_creation_variable_thickness_function(): + """Checks that a cadquery solid can be created using the BlanketFP + parametric component when a thickness function is passed as an + argument.""" + + def thickness(theta): + return 10 + 0.1 * theta + + test_shape = paramak.blanket_from_plasma(start_angle=-90, stop_angle=240, thickness=thickness) + + assert test_shape.vals()[0].Volume() > 1000 + + +def test_creation_variable_offset_from_tuple(): + """Checks that a cadquery solid can be created using the BlanketFP + parametric component when a tuple of offsets is passed as an + argument.""" + + test_shape = paramak.blanket_from_plasma(thickness=150, start_angle=-90, stop_angle=240, offset_from_plasma=(0, 10)) + + assert test_shape.vals()[0].Volume() > 1000 + + +def test_creation_variable_offset_from_2_lists(): + """Checks that a cadquery solid can be created using the BlanketFP + parametric component when a list of offsets and a list of angles are + passed as an argument.""" + + test_shape = paramak.blanket_from_plasma( + thickness=150, start_angle=90, stop_angle=270, offset_from_plasma=[[270, 100, 90], [0, 5, 10]] + ) + + assert test_shape is not None + + +def test_creation_variable_offset_error(): + """Checks that an error is raised when two lists with different + lengths are passed in offset_from_plasma as an argument.""" + + def test_different_lengths(): + with pytest.raises(ValueError) as excinfo: + paramak.blanket_from_plasma( + thickness=150, start_angle=90, stop_angle=270, offset_from_plasma=[[270, 100, 90], [0, 5, 10, 15]] + ) + assert "maximum recursion" in str(excinfo.value) + + +def test_creation_variable_offset_function(): + """Checks that a cadquery solid can be created using the BlanketFP + parametric component when an offset function is passed.""" + + def offset(theta): + return 10 + 0.1 * theta + + test_shape = paramak.blanket_from_plasma(thickness=150, start_angle=-90, stop_angle=240, offset_from_plasma=offset) + + assert test_shape is not None + assert test_shape.vals()[0].Volume() > 1000 + + +def test_full_cov_stp_export(): + """Creates a blanket using the BlanketFP parametric component with full + coverage and checks that an stp file can be exported using the export_stp + method.""" + + test_shape = paramak.blanket_from_plasma(thickness=150, rotation_angle=180, start_angle=0, stop_angle=360) + + exporters.export(test_shape, "test_blanket_full_cov.step") + assert Path("test_blanket_full_cov.step").exists() + os.system("rm test_blanket_full_cov.step") + + +def test_full_cov_full_rotation(): + """Creates a blanket using the BlanketFP parametric component with full + coverage and full rotation and checks that an stp file can be exported using + the export_stp method.""" + + test_shape = paramak.blanket_from_plasma(thickness=150, rotation_angle=360, start_angle=0, stop_angle=360) + + exporters.export(test_shape, "test_blanket_full_cov_full_rot.step") + assert Path("test_blanket_full_cov_full_rot.step").exists() + os.system("rm test_blanket_full_cov_full_rot.step") + + +def test_overlapping(): + """Creates an overlapping geometry and checks that a warning is raised.""" + + with pytest.warns(UserWarning, match="blanket_from_plasma: Some points with negative R"): + paramak.blanket_from_plasma( + major_radius=100, + minor_radius=100, + triangularity=0.5, + elongation=2, + thickness=200, + stop_angle=360, + start_angle=0, + ) diff --git a/tests/test_workplanes/test_center_column_shield_cylinder.py b/tests/test_workplanes/test_center_column_shield_cylinder.py new file mode 100644 index 000000000..902732643 --- /dev/null +++ b/tests/test_workplanes/test_center_column_shield_cylinder.py @@ -0,0 +1,29 @@ +import cadquery as cq + +import paramak + + +def test_creation(): + test_shape = paramak.center_column_shield_cylinder(height=100, inner_radius=10, thickness=20, rotation_angle=90) + + cq.exporters.export(test_shape, "center_column_shield_cylinder.step") + + +def test_reference_point(): + test_shape = paramak.center_column_shield_cylinder( + height=100, inner_radius=10, thickness=20, rotation_angle=90, reference_point=("center", 0) + ) + assert test_shape.val().BoundingBox().zmin == -50 + assert test_shape.val().BoundingBox().zmax == 50 + + test_shape = paramak.center_column_shield_cylinder( + height=100, inner_radius=10, thickness=20, rotation_angle=90, reference_point=("lower", 200) + ) + assert test_shape.val().BoundingBox().zmin == 200 + assert test_shape.val().BoundingBox().zmax == 300 + + test_shape = paramak.center_column_shield_cylinder( + height=100, inner_radius=10, thickness=20, rotation_angle=90, reference_point=("center", -200) + ) + assert test_shape.val().BoundingBox().zmin == -250 + assert test_shape.val().BoundingBox().zmax == -150 diff --git a/tests/test_workplanes/test_constant_thickness_dome.py b/tests/test_workplanes/test_constant_thickness_dome.py new file mode 100644 index 000000000..682e883b2 --- /dev/null +++ b/tests/test_workplanes/test_constant_thickness_dome.py @@ -0,0 +1,9 @@ +import cadquery as cq + +import paramak +import paramak.workplanes + +a = paramak.constant_thickness_dome(upper_or_lower="upper") +cq.exporters.export(a, "upper.step") +a = paramak.constant_thickness_dome(upper_or_lower="lower") +cq.exporters.export(a, "lower.step") diff --git a/tests/test_workplanes/test_dished_vacuum_vessel.py b/tests/test_workplanes/test_dished_vacuum_vessel.py new file mode 100644 index 000000000..409d1015a --- /dev/null +++ b/tests/test_workplanes/test_dished_vacuum_vessel.py @@ -0,0 +1,30 @@ +import cadquery as cq + +import paramak + + +def test_creation(): + lower_dome_section, cylinder_section, upper_dome_section = paramak.dished_vacuum_vessel() + cq.exporters.export(lower_dome_section, "lower_dome_section.step") + cq.exporters.export(cylinder_section, "cylinder_section.step") + cq.exporters.export(upper_dome_section, "upper_dome_section.step") + + +def test_reference_point(): + lower_dome_section, cylinder_section, upper_dome_section = paramak.dished_vacuum_vessel( + dish_height=(50, 50), cylinder_height=400, thickness=15, reference_point=("center", 0) + ) + assert lower_dome_section.val().BoundingBox().zmin == -265 + assert upper_dome_section.val().BoundingBox().zmax == 265 + + lower_dome_section, cylinder_section, upper_dome_section = paramak.dished_vacuum_vessel( + dish_height=(50, 50), cylinder_height=400, thickness=15, reference_point=("lower", 200) + ) + assert lower_dome_section.val().BoundingBox().zmin == 200 + assert upper_dome_section.val().BoundingBox().zmax == 730 + + lower_dome_section, cylinder_section, upper_dome_section = paramak.dished_vacuum_vessel( + dish_height=(50, 50), cylinder_height=400, thickness=15, reference_point=("center", -200) + ) + assert lower_dome_section.val().BoundingBox().zmin == -465 + assert upper_dome_section.val().BoundingBox().zmax == 65 diff --git a/tests/test_workplanes/test_plasma_simplified.py b/tests/test_workplanes/test_plasma_simplified.py new file mode 100644 index 000000000..7dbea1ce0 --- /dev/null +++ b/tests/test_workplanes/test_plasma_simplified.py @@ -0,0 +1,45 @@ +from pathlib import Path + +import cadquery as cq +import pytest +from cad_to_dagmc import CadToDagmc + +import paramak + +from .utils import transport_particles_on_h5m_geometry + + +@pytest.mark.parametrize("rotation_angle", [30, 90, 180, 360]) +def test_creation_different_angles(rotation_angle): + test_shape = paramak.plasma_simplified(rotation_angle=rotation_angle) + + cq.exporters.export(test_shape, f"plasma_simplified_{rotation_angle}.step") + + assert Path(f"plasma_simplified_{rotation_angle}.step").exists() + + +@pytest.mark.parametrize("rotation_angle", [60, 360]) +def test_transport_different_angles(rotation_angle): + test_shape = paramak.plasma_simplified(rotation_angle=rotation_angle) + + my_model = CadToDagmc() + material_tags = ["mat1"] + my_model.add_cadquery_object( + cadquery_object=test_shape, + material_tags=material_tags, + ) + + h5m_filename = "dagmc.h5m" + my_model.export_dagmc_h5m_file( + filename=h5m_filename, + min_mesh_size=10, + max_mesh_size=200, + ) + + flux = transport_particles_on_h5m_geometry( + h5m_filename=h5m_filename, + material_tags=material_tags, + nuclides=["H1"], + cross_sections_xml="tests/cross_sections.xml", + ) + assert flux > 0.0 diff --git a/tests/test_workplanes/test_poloidal_field_coil.py b/tests/test_workplanes/test_poloidal_field_coil.py new file mode 100644 index 000000000..aee6ac771 --- /dev/null +++ b/tests/test_workplanes/test_poloidal_field_coil.py @@ -0,0 +1,20 @@ +import math + +from cadquery import exporters + +import paramak + + +def test_construction(): + poloidal_field_coil = paramak.poloidal_field_coil(height=20, width=10, center_point=(100, 50), rotation_angle=360) + + exporters.export(poloidal_field_coil, "poloidal_field_coil.step") + + +def test_volume_rotation_angle(): + poloidal_field_coil = paramak.poloidal_field_coil(height=20, width=10, center_point=(100, 50), rotation_angle=360) + half_poloidal_field_coil = paramak.poloidal_field_coil( + height=20, width=10, center_point=(100, 50), rotation_angle=180 + ) + + assert math.isclose(half_poloidal_field_coil.val().Volume(), poloidal_field_coil.val().Volume() / 2) diff --git a/tests/test_workplanes/test_poloidal_field_coil_case.py b/tests/test_workplanes/test_poloidal_field_coil_case.py new file mode 100644 index 000000000..c2c3210d3 --- /dev/null +++ b/tests/test_workplanes/test_poloidal_field_coil_case.py @@ -0,0 +1,24 @@ +import math + +from cadquery import exporters + +import paramak + + +def test_construction(): + poloidal_field_coil_case = paramak.poloidal_field_coil_case( + coil_height=20, coil_width=10, center_point=(100, 50), rotation_angle=360, casing_thickness=10 + ) + + exporters.export(poloidal_field_coil_case, "poloidal_field_coil_case.step") + + +def test_volume_rotation_angle(): + poloidal_field_coil_case = paramak.poloidal_field_coil_case( + coil_height=20, coil_width=10, center_point=(100, 50), casing_thickness=10, rotation_angle=360 + ) + half_poloidal_field_coil_case = paramak.poloidal_field_coil_case( + coil_height=20, coil_width=10, center_point=(100, 50), casing_thickness=10, rotation_angle=180 + ) + + assert math.isclose(half_poloidal_field_coil_case.val().Volume(), poloidal_field_coil_case.val().Volume() / 2) diff --git a/tests/test_workplanes/test_toroidal_field_coil_rectangle.py b/tests/test_workplanes/test_toroidal_field_coil_rectangle.py new file mode 100644 index 000000000..bf4a2958c --- /dev/null +++ b/tests/test_workplanes/test_toroidal_field_coil_rectangle.py @@ -0,0 +1,24 @@ +import math + +from cadquery import exporters + +import paramak + + +def test_construction(): + toroidal_field_coil_rectangle = paramak.toroidal_field_coil_rectangle(azimuthal_placement_angles=[0, 90, 180]) + + exporters.export(toroidal_field_coil_rectangle, "toroidal_field_coil_rectangle.step") + + +def test_volume_rotation_angle(): + toroidal_field_coil_rectangle = paramak.toroidal_field_coil_rectangle( + azimuthal_placement_angles=[0, 90], horizontal_start_point=(50, 200) + ) + half_toroidal_field_coil_rectangle = paramak.toroidal_field_coil_rectangle( + azimuthal_placement_angles=[0], horizontal_start_point=(50, 200) + ) + + assert math.isclose( + half_toroidal_field_coil_rectangle.val().Volume(), toroidal_field_coil_rectangle.val().Volume() / 2 + ) diff --git a/tests/test_workplanes/test_u_shaped_dome.py b/tests/test_workplanes/test_u_shaped_dome.py new file mode 100644 index 000000000..2845b55ed --- /dev/null +++ b/tests/test_workplanes/test_u_shaped_dome.py @@ -0,0 +1,17 @@ +from cadquery import exporters + +import paramak + +# TODO try OCP 7.7.1 and newer cadquery +# lower with rotation_angle less than 360 can fail, OCP 7.7.0 + + +def test_construction(): + dome_section, cylinder_section = paramak.u_shaped_dome(upper_or_lower="lower", rotation_angle=360) + + exporters.export(dome_section, "dome_section.step") + exporters.export(cylinder_section, "cylinder_section.step") + + domes = paramak.u_shaped_dome(upper_or_lower="upper", rotation_angle=360) + domes = paramak.u_shaped_dome(upper_or_lower="upper", rotation_angle=360, reference_point=("lower", 10)) + domes = paramak.u_shaped_dome(upper_or_lower="upper", rotation_angle=360, reference_point=("center", 20)) diff --git a/tests/test_workplanes/utils.py b/tests/test_workplanes/utils.py new file mode 100644 index 000000000..d9f691c0e --- /dev/null +++ b/tests/test_workplanes/utils.py @@ -0,0 +1,102 @@ +def transport_particles_on_h5m_geometry( + h5m_filename: str, + material_tags: list, + nuclides: list = None, + cross_sections_xml: str = None, +): + """A function for testing the geometry file with particle transport in + DAGMC OpenMC. Requires openmc and either the cross_sections_xml to be + specified. Returns the flux on each volume + + Arg: + h5m_filename: The name of the DAGMC h5m file to test + material_tags: The + nuclides: + cross_sections_xml: + + """ + + import openmc + from openmc.data import NATURAL_ABUNDANCE + + if nuclides is None: + nuclides = list(NATURAL_ABUNDANCE.keys()) + + materials = openmc.Materials() + for i, material_tag in enumerate(set(material_tags)): + # simplified material definitions have been used to keen this example minimal + mat_dag_material_tag = openmc.Material(name=material_tag) + mat_dag_material_tag.add_nuclide(nuclides[i], 1, "ao") + mat_dag_material_tag.set_density("g/cm3", 0.01) + + materials.append(mat_dag_material_tag) + + if cross_sections_xml: + openmc.config["cross_sections"] = cross_sections_xml + + else: + with open("cross_sections.xml", "w") as file: + file.write( + """ + + + + + + """ + ) + + openmc.config["cross_sections"] = "cross_sections.xml" + + dag_univ = openmc.DAGMCUniverse(filename=h5m_filename) + bound_dag_univ = dag_univ.bounded_universe() + geometry = openmc.Geometry(root=bound_dag_univ) + + # initializes a new source object + my_source = openmc.IndependentSource() + + center_of_geometry = ( + (dag_univ.bounding_box[0][0] + dag_univ.bounding_box[1][0]) / 2, + (dag_univ.bounding_box[0][1] + dag_univ.bounding_box[1][1]) / 2, + (dag_univ.bounding_box[0][2] + dag_univ.bounding_box[1][2]) / 2, + ) + # sets the location of the source which is not on a vertex + center_of_geometry_nudged = ( + center_of_geometry[0] + 0.1, + center_of_geometry[1] + 0.1, + center_of_geometry[2] + 0.1, + ) + + my_source.space = openmc.stats.Point(center_of_geometry_nudged) + # sets the direction to isotropic + my_source.angle = openmc.stats.Isotropic() + # sets the energy distribution to 100% 14MeV neutrons + my_source.energy = openmc.stats.Discrete([14e6], [1]) + + # specifies the simulation computational intensity + settings = openmc.Settings() + settings.batches = 10 + settings.particles = 10000 + settings.inactive = 0 + settings.run_mode = "fixed source" + settings.source = my_source + + # adds a tally to record the heat deposited in entire geometry + cell_tally = openmc.Tally(name="flux") + cell_tally.scores = ["flux"] + + # groups the two tallies + tallies = openmc.Tallies([cell_tally]) + + # builds the openmc model + my_model = openmc.Model(materials=materials, geometry=geometry, settings=settings, tallies=tallies) + + # starts the simulation + output_file = my_model.run() + + # loads up the output file from the simulation + statepoint = openmc.StatePoint(output_file) + + my_flux_cell_tally = statepoint.get_tally(name="flux") + + return my_flux_cell_tally.mean.flatten()[0]