From 6d6407343908208f5e1755cfc5d38e380dfd397f Mon Sep 17 00:00:00 2001 From: je-cook <81617086+je-cook@users.noreply.github.com> Date: Tue, 14 May 2024 14:30:37 +0100 Subject: [PATCH] =?UTF-8?q?=E2=9C=A8=202D=20Neutronics=20feature=20branch?= =?UTF-8?q?=20(#2656)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * 2D Neutronics module (#2344) * Initial commit to repository * Create basic README.md * Editing README.md * Code styling: Black, isort and whitespace/trailing new lines (#1) * blacken and isort * whitespace and new lines * Rearrage and refactor (#2) * Adding plasma source, DCLL and WCLL. Material defintions and geometries still need tweaking. * Changing from toroidal first wall profile to match points and adding extra source geometry plots. * Adding source code for parametric plasma source * Adjusting point too close to plasma. * Minor format tidy up. Moving function. * Passing geometry_variables to get_fw_points * Update make_materials.py Changing tungsten isotopes and tidying export_materials() * Update make_materials.py W isotopes corrected. * Made most of it OOP, especially the summary, so that the data can be recalled. * Added vscode into gitignore * Finished reformatting the Summary object so now it prints the entire run properly. Not documented yet, and kwargs not passed on yet. * Finished some reformatting. Todo list is at the top. * Some work done on load_fw_points, but incomplete. TODO list updated. * Added some documentations, but found a lot of magic numbers and 'things that don't belong in this module'. I'll deal with them later. * Changed the ordering and naming of the arguments to make_geometry in make_geometry.py, and added some docstring. * Fxied the types of Blankets available using BlanketType in make_materials.py. * Fixed the naming for cells in quick_tbr_heating.py - renamed them to cells_and_cell_lists * Updated the todo list * Cleaned up create_tallies, and in the process, created MaterialsLibrary in make_materials.py so that it becomes modular enough for testing. * Formatted the 'print_df' booleans away into a decorator - can change into printing to bluemira log later. * Added documentations to the rest of the important classes. * Updated todo list * one line change to make pandas print more prettily, before merge commit * moved all relevant files into the neutronics directory. * Moved files from the 2d-neutronics repo to here. * restored comment style fix * moved the neutronics files to the correct directory. * changed the import relative path * relocated files to the correct directory * changed the import path * minor changes before complete overhaul (the latter is for adhering to bluemira code style) * moved some constants to constants.py * Finished one task: placed the relevant constants and conversion units into constants.py * updated to-do list * Further broke down the make_materials.py file into materials_definition.py file * Refactored all of materials. * updated .gitignore * Moved the results formatting functions to a separate file * removed unused file * Improved styling * Broke functions down into appropriate modules * Split into appropriate modules * Added typing information * Updated the import diagram to reflect the import relationships between modules. * Fixed all fixable formatting issues according to flake8 * Fixed all flake8 messages! * installs script * cleaned up script * added complete apt install deps * Updated install script, removed pps_api * Updated to-do list, improved result presentation workflow * .npy requires conversion to BluemiraWire * ➕ Use ninja and install libopenmpi-dev * ➕ Add openmc to conda * ✨ Start using wire as input * 🚨 Some ruff fixes * 🎨 Move ifmain to example * 🚚 Move npy * 🎨 Fix example * 🎨 Get constants from base * 🎨 Fix imports * 🚨 Some ruff fixes * added openmc version to install script, added pps-isotropic package import to quick tbr * removed my path from config * 🎨 Cleanup printing * 🎨 Unit fix * 🎨 Unit fix2 * 🎨 Unit fix3 * 🎨 Use openmc enum * Preparing for unit changes * quick commit before rebase, units half-done * 🙈 Remove noise * 🔥 Remove so file * Removed the not-sensible default constants and placed thme into neutronics.ex.py instead. * Added TODO comment * 🎨 Cell filter cleanup * 🎨 Cleanup tallies * Changed all units * minor fix of import statement * Monkey patched successfully * monkeypatched successfully, and separated plasma geometry variables from tokamak geometry variables. This is the commit BEFORE applying the DataclassUnitConverter pattern. * Minor update to add comments. * Removed the duplicate codes inside params.py. This is the commit AFTER applying the DataclassUnitConverter pattern. * Bug fix - missing 'self,' argument fixed. * Updated todos * 🎨 Change raw_uc percent use and runtime var * 🚧 Make geometry dataclass refactor * 🚧 Slight cleanup of params * 🚨 Ignores * 📝 Docs * 📝 Fixes and docs * 📝 Fix documentation build * Update neutronics.ex.py * Update pyproject.toml * 🐛 Fix volume calculation * 🎨 Cleanup vol calc * 🚨 Fix black * 🎨 Small fixes * 🚚 Move quick tbr to neutronics_axisymmetric * 🎨 Add stochastic volume to output results * Added comments on why only two quantities are excluded in the conversion * pps-isotropic's API changed, making appropriate changes here to maintain usability * Added markdown docs on neutronics.ex.py * 💩 Hack for old version of openmc * 👥 Add contributors --------- Co-authored-by: jamesnha <92076761+jamesnha@users.noreply.github.com> Co-authored-by: je-cook <81617086+je-cook@users.noreply.github.com> Co-authored-by: Oliver Funk * ✨ OpenMC data downloader (#2618) * ✨ OpenMC data downloader * 🎨 Parse known args * 🎨 Ifmain * 🎨 Rename file * 🎨 Write them all then delete all * 🐛 Overindent * 🚨 Cleanup and linting * 🐛 No output * 🐛 Fix neutronics json location * 🐛 Hack to force data into a given directory * 🎨 Spelling * 🏷️ Typing fix * 🚨 Ruff fixes * 🙈 Add to gitignore * Minor error fixes to the outputted tables (#2689) * Fixed unit conversion error in photon heat flux; and turned off Bremsstrahlung heating by default * Undo accidental indent * Kept all statepoint files' units as openmc as cgs, so that the conversion to SI is done purely in python. Fixed misinterpretation of photon flux. * added comments to explain tally score choice. * Minor bug fix to pass code quality checks * 🚨 Ruff fixes * first move to BM materials, lots of refactoring needed but the example runs * 🎨 Basic cleanup * 🎨 Dataclassify * 🐛 Undo Be12Ti density change * 💡 Minor fixes * 📄 Copyright fixes after rebase * ✨ Allow '*' to extract all XS data * ⬆️ Openmc 0.13.3 * temperature docstring. * removed rounding * calculate the DT energy more accurately using plasma physics. * Simplified code to use n_DT_reactions instead. * 1. extracted (major_r, minor_r, elong, triang) into a different class than TokamakGeometry, and named this class PlasmaGeometry. 2. Removed the 'CGS' and 'PPS' from the 3 classes in params.py, and instead appended their baseclasses with 'Base', so that it becomes clearer that these dataclasses store SI units unless the .csg or .plasma_physics_units property is accessed. * (minor bug fix) * Updated the Shafranov shift description. * 🚨 Fix linter * Minor updates to docstrings * 🚨 Ruff & python3.10 fixes * Increase CSG speed without compromising on flexibility (#3079) * :Changed params.py according to Matti's suggestion * Made temporary script * Added json data to be loaded as test case bluemirawires * Added post_init checks. * Created basic script to translate from existing dataclass TokamakGeometry (make_csg) * WIP: making pre-cell. * WIP: making pre-cell. * Bodging serializing issue - fix later. * Deleted failed files. * Test script in attempt to replace the magic numbers/magic procedure to something more generalizable * quick commit * Script to load the divertor. * Added capability to automatically slice up the blanket according to the first wall panels * Updated script for creating neutronics csg * extended the functionality of signed distance to check points as well. * Added function based on signed_distance to check for overlaps. * Add tolerance to account for low sampling rate/ float errors. * Adding module to optimally slice the neutronics model. * Updated module that chooses better surfaces for faster neutronics runs. * Made the volume optimization converge. * Made a better approximation of the reactor. Added tools necessary to make such approximation * Added better documentations. Formatting. * Made a very rigorous structure to create openmc.Cells. * removed unused code * Updated test case * Half-working state * Finally finished the blanket bug fix. * Updated docstrings * Updated equations * Updated comment * Interited from collections.abc.Sequence to avoid having to make def __iter__(). * Cleaned up excess code * minor changes to naming and docstring * Updated the run modes into context managers. * Enforced self.setup() to be ran before self.run() in the contextmanagers. * Added ability to slice divertor wire * Variable renaming * Removed redundant checks * Improved typing (Iterables vs Sequence), and added classes to prepare for tallying. * Added tallying functionality * Developing divertor cells automatic conversion into pre-cells... * Added type hints for return types. * Fixed import * Successfully exported into DivertorPreCellArray. * Updated script correspondingly * split up WireInfo etc. into its own file. * Made make_divertor_pre_cell_array more robust (able to straighten out the cell walls better before passing onto DivertorPreCell). * Typing and docstring improvement * Added function to flip a direction vector if it is pointing the wrong way; Moved plane functionalities to geometry/plane.py * Removed plane creation functionality (migrated to bluemira.geometry.plane) * Minor change to type hinting in make_csg.py * Minor changes to type hints in radial_wall.py * Slight change to variable name * Added plan to make the divertor slicing more robust. * Moved the warning to an appropriate place * Improved the interior_wire_offset method * Fixed minor mistake in direction * Working version: produces pre-cell arrays for both divertor and blanket * Pointed the chord in the correct direction again * Removed unused data loading files. * Improved file structure so each class has separate responsibility. Surfaces that aren't shared among different BlanketCell stacks won't be created in BlanketCellArray. * Minor bug fixes related to types and WireInfo * Moved constants into constants.py * Docstring updates * Simplified boolean expression; made plasma cell. * Deleted unused statement. * updated test script to plot better. * Named file more sensibly. * Fixed some ID control issues * docstring rearrangement * Fixed up the get_all_hollow_merged_cells method so the test script won't have to write something as long to plot the overall-cells. * Fixed the thickness issue - now instead we have absolute thickness. * Got as far as getting some segfaults out of openmc. * Finally working. Tallies need fixing * Fixed floating point issue * Added missing conversion function * Changed thickness. * Modified to allow for a variable number of batches. * Delete some of the obsolete classes * Updated TODO tags in test script. * Updated TODOs * Added missing method back in * Updated documentations slightly to make it more readable. * Updated docstring to conform usage of anti- vs counter- * Ensured backward (neutronics.ex.py) compatibility * Calculated volumes without relying on openmc. * Added function to support volume calculation of the openmc.Cell's .volumes. * Optimized volume calculation (commented them out where they're not needed) * Working version. * Added TODO on reducing the number of negated regions. * ruff fixes * Some updates to how surfaces are get * fixed ONE TODO tag. * Updated zip strict = True almost everywhere, pairwise where applicable. * Cleaned up code slightly at make_csg * We can now control the thickness of the central solenoid from parameter frame directly. * Attempting to use vacuum vessel wire as well. * Further attempt at incorporating the vacuum vessel interior curve into the pre-cells * VV now successfully sandwhiched between the two curvesgit add . * Minor formatting updates * Finished porting over the materials library * Cleaned out all low-lying fruit TODO tags * Made note about convexity. * Neutronics Solver interface (#3231) * 🚧 OpenMC designer * 🚧 Working designer * ✨ Openmc solver * ✨ Fixup from upstream * 🚚 Move stuff * Neutronics old code removal and architectural adjustments (#3232) * 🌐 Bye Z * ♻️ Initial cleanup * ♻️ More cleanup * ♻️ More cleanup * ♻️ More cleanup * 🔥 Destroying Vertices pt1 * 🔥 Destroying Vertices pt2 * ♻️ More cleanup * 🔥 Remove Sequence inheritance * 🎨 Types * 🚨 Remove ruff ignores * ♻️ More cleanup * ♻️ Remove global hangar variable * ♻️ More cleanup * ♻️ More cleanup * 🚨 More cleanup * 📝 Docs cleanup * 📝 Docs cleanup * 🐛 Fix double counting of divertor material * 🐛 Fix 0 volumetric heating * 🎨 Convert SingleNullTokamak to a function * 🏷️ Typing and general cleanup * ♻️ Cleanup dimensions * 🎨 Order functions * ✨ Split of Neutroics pre cell csg and openmc running (#3243) * 🎨 Neutronics PreCell / openmc separation * 🚚 Move stuff about * 🎨 Reduce surface that requires openmc * 🚚 Move materials about * 🎨 Initial material rearrangement * 🌐 More internationalisation * ♻️ Neutronics material movement * 🎨 More materials stuff * ♻️ Cleanup * 🎨 Materials cleanup * Neutronics EUDEMO integration (#3239) * inital go * Saving * hard part done * pass vals to Solver * Updated WireInfoList so it doesn't use horrible mixed types. * expected div wire * Fixed the issue of 'sometimes misordering wires', and used recursion to shorted the break_wire_into_convex_chunks method for better readability with no loss in speed. * 🚧 Tmp get it working * Fixed the .reverse() method, which didn't flip the sign of the tangent. * Improved plotting in plot_2d and plot_surfaces (mainly used for debugging) * 🐛 Fix inner panel point * 🔥 Reduce interface * ♻️ More cleanup * 🎨 More cleanup * 🎨 Remove while loop from break wire into convex chunks * 🎨 Reduce interface * 🎨 More cleanup * 🔥 Remove readme * 🐛 Typing and divertor start point fix * 🎨 Fix extra error and cleanup some reactor mechanics * 🎨 Improve file save location * 🎨 Only pass in ivcshapes not divertor * 🎨 Feed in correct blanket cut point * 🎨 Cleanup runmode * Fixed leaking plasma void. * 🎨 Cleanup inputs * Added the the half_bounding_box function to get the minimum positive r. * 🔥 Remove old data files * ♻️ Run through cleanup * Fixed the final overlap issue at the bottom of the divertor * 🚚 Move some geometry functions to geometry * 🐛 Fix tfcoil size * 🐛 Fix blanket test * 🚨 Linter fix * 🚚 Move some things about * 🚚 Move openmc to codes * 🎨 Isolate eudemo from openmc * 🚚 Move neutronics under radiation_transport * Added volume tag * Added checks to ensure convexity. * Forced the pair of first and last pre-cells in the blanket to mate with the first and last pre-cell pair of divertor. * Fixed the convexity proble; but a different issue arose: floating point precision of the exterior points are all wrong. * Removed debug statement * Added tolerances back in * Floating point precision issue FIXEDgit add make_pre_cell.py openmc/make_csg.py Now we can use tighter tolerances. * 🔥 Remove old examples * integrated rad shield into the neutronics run * move set_volumes to CellStage --------- Co-authored-by: ocean Co-authored-by: james * 🎨 Radiation shield initial material * 📝 Remove neutronics example reference * 📝 Add notice for use and update environment * Improve tallies (#3275) * 🎨 Allow custom filter for tallies * 🐛 Fix tally function for eudemo run * 🐛 Fix remove reindex to avoid nans * Changed tallies to get the required data (need to multiply by appropriate constants in output.py to finish the job) * minor bug fix (missed a variable) * 🐛 Fix filter creation * 🐛 Fix output for temporary output --------- Co-authored-by: ocean * ⬆️ Upgrade openmc --------- Co-authored-by: Ocean Co-authored-by: jamesnha <92076761+jamesnha@users.noreply.github.com> Co-authored-by: Oliver Funk Co-authored-by: Oliver Funk --- .gitignore | 6 + CITATION.cff | 8 + bluemira/codes/openmc/__init__.py | 10 + bluemira/codes/openmc/make_csg.py | 2247 +++++++++++++++++ bluemira/codes/openmc/material.py | 113 + bluemira/codes/openmc/output.py | 339 +++ bluemira/codes/openmc/params.py | 165 ++ bluemira/codes/openmc/solver.py | 519 ++++ bluemira/codes/openmc/sources.py | 77 + bluemira/codes/openmc/tallying.py | 90 + bluemira/codes/wrapper.py | 28 + bluemira/equilibria/physics.py | 2 +- bluemira/geometry/coordinates.py | 52 + bluemira/geometry/plane.py | 64 +- bluemira/geometry/tools.py | 191 +- bluemira/radiation_transport/error.py | 8 + .../neutronics/__init__.py | 8 + .../neutronics/blanket_data.py | 595 +++++ .../neutronics/constants.py | 108 + .../neutronics/geometry.py | 152 ++ .../neutronics/make_pre_cell.py | 767 ++++++ .../neutronics/materials.py | 37 + .../neutronics/neutronics_axisymmetric.py | 287 +++ .../neutronics/radial_wall.py | 320 +++ .../radiation_transport/neutronics/slicing.py | 937 +++++++ .../radiation_transport/neutronics/wires.py | 182 ++ conda/environment.yml | 1 + eudemo/README.md | 4 + eudemo/config/build_config.json | 22 +- eudemo/config/params.json | 66 + eudemo/eudemo/blanket/__init__.py | 30 +- eudemo/eudemo/blanket/designer.py | 46 +- eudemo/eudemo/ivc/__init__.py | 7 + eudemo/eudemo/neutronics/__init__.py | 6 + eudemo/eudemo/neutronics/run.py | 99 + eudemo/eudemo/params.py | 13 + eudemo/eudemo/reactor.py | 13 +- eudemo/eudemo_tests/blanket/test_blanket.py | 27 +- eudemo/eudemo_tests/blanket/test_designer.py | 2 +- .../radiation_calculation_solver_DEMO.py | 4 +- pyproject.toml | 7 +- requirements.txt | 1 + scripts/install-openmc.sh | 45 + scripts/openmc_data_download/__init__.py | 1 + .../multithreaded_download.py | 132 + .../neutronics_data_downloader.py | 248 ++ .../nuclear_data_isotopes.json | 9 + .../nuclear_data_isotopes_example.json | 386 +++ 48 files changed, 8439 insertions(+), 42 deletions(-) create mode 100644 bluemira/codes/openmc/__init__.py create mode 100644 bluemira/codes/openmc/make_csg.py create mode 100644 bluemira/codes/openmc/material.py create mode 100644 bluemira/codes/openmc/output.py create mode 100644 bluemira/codes/openmc/params.py create mode 100644 bluemira/codes/openmc/solver.py create mode 100644 bluemira/codes/openmc/sources.py create mode 100644 bluemira/codes/openmc/tallying.py create mode 100644 bluemira/radiation_transport/neutronics/__init__.py create mode 100644 bluemira/radiation_transport/neutronics/blanket_data.py create mode 100644 bluemira/radiation_transport/neutronics/constants.py create mode 100644 bluemira/radiation_transport/neutronics/geometry.py create mode 100644 bluemira/radiation_transport/neutronics/make_pre_cell.py create mode 100644 bluemira/radiation_transport/neutronics/materials.py create mode 100644 bluemira/radiation_transport/neutronics/neutronics_axisymmetric.py create mode 100644 bluemira/radiation_transport/neutronics/radial_wall.py create mode 100644 bluemira/radiation_transport/neutronics/slicing.py create mode 100644 bluemira/radiation_transport/neutronics/wires.py create mode 100644 eudemo/eudemo/neutronics/__init__.py create mode 100644 eudemo/eudemo/neutronics/run.py create mode 100755 scripts/install-openmc.sh create mode 100644 scripts/openmc_data_download/__init__.py create mode 100644 scripts/openmc_data_download/multithreaded_download.py create mode 100644 scripts/openmc_data_download/neutronics_data_downloader.py create mode 100644 scripts/openmc_data_download/nuclear_data_isotopes.json create mode 100644 scripts/openmc_data_download/nuclear_data_isotopes_example.json diff --git a/.gitignore b/.gitignore index 17973d904c..06f2d6a2e6 100644 --- a/.gitignore +++ b/.gitignore @@ -179,4 +179,10 @@ autoapi generated_data !generated_data/naughty_geometry/README.md tests/bluemira/test_generated_data +cross_section_data plasmod_*.dat + +# openmc outputs +*.svg +*.xml +*.out diff --git a/CITATION.cff b/CITATION.cff index 66c840df10..3b0f7a8531 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -60,6 +60,10 @@ authors: given-names: G. A. affiliation: United Kingdom Atomic Energy Authority, Culham Science Centre, Abingdon, Oxfordshire OX14 3DB, United Kingdom + - family-names: Hagues + given-names: J. + affiliation: United Kingdom Atomic Energy Authority, Culham Science Centre, Abingdon, Oxfordshire OX14 3DB, United Kingdom + - family-names: Humphrey given-names: L. affiliation: United Kingdom Atomic Energy Authority, Culham Science Centre, Abingdon, Oxfordshire OX14 3DB, United Kingdom @@ -88,6 +92,10 @@ authors: given-names: D. affiliation: United Kingdom Atomic Energy Authority, Culham Science Centre, Abingdon, Oxfordshire OX14 3DB, United Kingdom + - family-names: Wong + given-names: O. + affiliation: United Kingdom Atomic Energy Authority, Culham Science Centre, Abingdon, Oxfordshire OX14 3DB, United Kingdom + message: "If you use this software, please cite it using these references." references: - type: article diff --git a/bluemira/codes/openmc/__init__.py b/bluemira/codes/openmc/__init__.py new file mode 100644 index 0000000000..2e89083aed --- /dev/null +++ b/bluemira/codes/openmc/__init__.py @@ -0,0 +1,10 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""OpenMC interface""" + +from bluemira.codes.openmc.solver import OpenMCNeutronicsSolver as Solver + +__all__ = ["Solver"] diff --git a/bluemira/codes/openmc/make_csg.py b/bluemira/codes/openmc/make_csg.py new file mode 100644 index 0000000000..c243d19ff9 --- /dev/null +++ b/bluemira/codes/openmc/make_csg.py @@ -0,0 +1,2247 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +""" +Create csg geometry by converting from bluemira geometry objects made of wires. All units +in this module are in SI (distrance:[m]) unless otherwise specified by the docstring. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from itertools import chain, pairwise +from typing import TYPE_CHECKING + +import matplotlib.pyplot as plt +import numpy as np +import openmc +import openmc.region + +from bluemira.codes.openmc.material import CellType +from bluemira.geometry.constants import D_TOLERANCE +from bluemira.geometry.error import GeometryError +from bluemira.geometry.solid import BluemiraSolid +from bluemira.geometry.tools import ( + is_convex, + make_circle_arc_3P, + make_polygon, + polygon_revolve_signed_volume, + revolve_shape, +) +from bluemira.geometry.wire import BluemiraWire +from bluemira.radiation_transport.neutronics.constants import ( + DTOL_CM, + to_cm, + to_cm3, + to_m, +) +from bluemira.radiation_transport.neutronics.radial_wall import Vert +from bluemira.radiation_transport.neutronics.wires import CircleInfo + +if TYPE_CHECKING: + from collections.abc import Iterable, Iterator, Sequence + + import numpy.typing as npt + + from bluemira.codes.openmc.material import MaterialsLibrary + from bluemira.geometry.coordinates import Coordinates + from bluemira.radiation_transport.neutronics.geometry import ( + DivertorThickness, + TokamakDimensions, + ) + from bluemira.radiation_transport.neutronics.make_pre_cell import ( + DivertorPreCell, + DivertorPreCellArray, + PreCell, + PreCellArray, + ) + from bluemira.radiation_transport.neutronics.neutronics_axisymmetric import ( + NeutronicsReactor, + ) + from bluemira.radiation_transport.neutronics.wires import ( + StraightLineInfo, + WireInfoList, + ) + + +# Found to work by trial and error. I'm sorry. +SHRINK_DISTANCE = 0.0005 # [m] = 0.05cm = 0.5 mm + + +@dataclass +class CellStage: + """Stage of making cells.""" + + blanket: BlanketCellArray + divertor: DivertorCellArray + tf_coils: list[openmc.Cell] + cs_coil: openmc.Cell + plasma: openmc.Cell + radiation_shield: openmc.Cell + ext_void: openmc.Cell + universe: openmc.region.Intersection + + @property + def cells(self): + """Get the list of all cells.""" + return ( + *chain.from_iterable((*self.blanket, *self.divertor)), + *self.tf_coils, + self.cs_coil, + self.plasma, + self.radiation_shield, + self.ext_void, + ) + + def get_all_hollow_merged_cells(self): + """Blanket and divertor cells""" + return [ + *[openmc.Cell(region=stack.get_overall_region()) for stack in self.blanket], + *[openmc.Cell(region=stack.get_overall_region()) for stack in self.divertor], + ] + + def set_volumes(self): + """ + Sets the volume of the voids. Not necessary/ used anywhere yet. + """ + ext_vertices = exterior_vertices(self.blanket, self.divertor) + total_universe_volume = ( + # top - bottom + (self.universe[0].surface.z0 - self.universe[1].surface.z0) + * np.pi + * self.universe[2].surface.r ** 2 # cylinder + ) # cm^3 + + # is this needed? + # self.universe.volume = total_universe_volume + + outer_boundary_volume = to_cm3( + polygon_revolve_signed_volume(ext_vertices[:, ::2]) + ) + ext_void_volume = total_universe_volume - outer_boundary_volume + if self.tf_coils: + for coil in self.tf_coils: + ext_void_volume -= coil.volume + if self.cs_coil: + ext_void_volume -= self.cs_coil.volume + self.ext_void.volume = ext_void_volume + blanket_volumes = sum(cell.volume for cell in chain.from_iterable(self.blanket)) + divertor_volumes = sum( + cell.volume for cell in chain.from_iterable(self.divertor) + ) + self.plasma.volume = outer_boundary_volume - blanket_volumes - divertor_volumes + + +def is_monotonically_increasing(series): + """Check if a series is monotonically increasing""" # or decreasing + return all(np.diff(series) >= 0) # or all(diff<0) + + +def plot_surfaces(surfaces_list: list[openmc.Surface], ax=None): + """ + Plot a list of surfaces in matplotlib. + """ + if not ax: + ax = plt.axes() + # ax.set_aspect(1.0) # don't do this as it makes the plot hard to read. + for i, surface in enumerate(surfaces_list): + plot_surface_at_1000cm(ax, surface, color_num=i) + ax.legend() + ax.set_aspect("equal") + ax.set_ylim([-1000, 1000]) + ax.set_xlim([-1000, 1000]) + return ax + + +def plot_surface_at_1000cm(ax, surface: openmc.Surface, color_num: int): + """ + In the range [-1000, 1000], plot the RZ cross-section of the ZCylinder/ZPlane/ZCone. + """ + if isinstance(surface, openmc.ZCylinder): + ax.plot( + [surface.x0, surface.x0], + [-1000, 1000], + label=f"{surface.id}: {surface.name}", + color=f"C{color_num}", + ) + elif isinstance(surface, openmc.ZPlane): + ax.plot( + [-1000, 1000], + [surface.z0, surface.z0], + label=f"{surface.id}: {surface.name}", + color=f"C{color_num}", + ) + elif isinstance(surface, openmc.ZCone): + intercept = surface.z0 + slope = 1 / np.sqrt(surface.r2) + + def equation_pos(x): + return slope * np.array(x) + intercept + + def equation_neg(x): + return -slope * np.array(x) + intercept + + y_pos, y_neg = equation_pos([-1000, 1000]), equation_neg([-1000, 1000]) + ax.plot( + [-1000, 1000], + y_pos, + label=f"{surface.id}: {surface.name} (upper)", + color=f"C{color_num}", + ) + ax.plot( + [-1000, 1000], + y_neg, + label=f"{surface.id}: {surface.name} (lower)", + linestyle="--", + color=f"C{color_num}", + ) + + +def get_depth_values( + pre_cell: PreCell, blanket_dimensions: TokamakDimensions +) -> npt.NDArray[np.float64]: + """ + Choose the depth values that this pre-cell is suppose to use, according to where it + is physically positioned (hence is classified as inboard or outboard). + + Parameters + ---------- + pre_cell: + :class:`~PreCell` to be classified as either inboard or outboard + blanket_dimensions: + :class:`bluemira.radiation_transport.neutronics.params.TokamakDimensions` + recording the dimensions of the blanket in SI units (unit: [m]). + + Returns + ------- + depth_series: + a series of floats corresponding to the N-1 interfaces between the N layers. + Each float represents how deep into the blanket (i.e. how many [m] into the + first wall we need to drill, from the plasma facing surface) to hit that + interface layer. + """ + if check_inboard_outboard(pre_cell, blanket_dimensions): + return blanket_dimensions.inboard.get_interface_depths() + return blanket_dimensions.outboard.get_interface_depths() + + +def check_inboard_outboard( + pre_cell: PreCell, blanket_dimensions: TokamakDimensions +) -> bool: + """If this pre-cell is an inboard, return True. + Otherwise, this pre-cell belongs to outboard, return False + """ + # reference radius + return ( + pre_cell.vertex[0].mean() < blanket_dimensions.inboard_outboard_transition_radius + ) + + +def torus_from_3points( + point1: npt.NDArray[np.float64], + point2: npt.NDArray[np.float64], + point3: npt.NDArray[np.float64], + surface_id: int | None = None, + name: str = "", +) -> openmc.ZTorus: + """ + Make a circular torus centered on the z-axis using 3 points. + All 3 points should lie on the RZ plane AND the surface of the torus simultaneously. + + Parameters + ---------- + point1, point2, point3: + RZ coordinates of the 3 points on the surface of the torus. + surface_id, name: + See openmc.Surface + """ + point1 = point1[0], 0, point1[-1] + point2 = point2[0], 0, point2[-1] + point3 = point3[0], 0, point3[-1] + circle = make_circle_arc_3P(point1, point2, point3) + cad_circle = circle.shape.OrderedEdges[0].Curve + center, radius = cad_circle.Center, cad_circle.Radius + return torus_from_circle( + [center[0], center[-1]], radius, surface_id=surface_id, name=name + ) + + +def torus_from_circle( + center: Sequence[float], + minor_radius: float, + surface_id: int | None = None, + name: str = "", +) -> openmc.ZTorus: + """ + Make a circular torus centered on the z-axis. + The circle would lie on the RZ plane AND the surface of the torus simultaneously. + + Parameters + ---------- + minor_radius: + Radius of the cross-section circle, which forms the minor radius of the torus. + center: + Center of the cross-section circle, which forms the center of the torus. + surface_id, name: + See openmc.Surface + """ + return z_torus( + [center[0], center[-1]], minor_radius, surface_id=surface_id, name=name + ) + + +def z_torus( + center: npt.ArrayLike, + minor_radius: float, + surface_id: int | None = None, + name: str = "", +) -> openmc.ZTorus: + """ + A circular torus centered on the z-axis. + The center refers to the major radius and it's z coordinate. + + Parameters + ---------- + center: + The center of the torus' RZ plane cross-section + minor_radius: + minor radius of the torus + + """ + major_radius, height, minor_radius = to_cm([center[0], center[-1], minor_radius]) + return openmc.ZTorus( + z0=height, + a=major_radius, + b=minor_radius, + c=minor_radius, + surface_id=surface_id, + name=name, + ) + + +def choose_halfspace( + surface: openmc.Surface, choice_points: npt.NDArray +) -> openmc.Halfspace: + """ + Simply take the centroid point of all of the choice_points, and choose the + corresponding half space + + Parameters + ---------- + surface: + an openmc surface + + """ + pt = np.mean(to_cm(choice_points), axis=0) + value = surface.evaluate(pt) + if value > 0: + return +surface + if value < 0: + return -surface + raise GeometryError("Centroid point is directly on the surface") + + +def choose_plane_cylinders( + surface: openmc.ZPlane | openmc.ZCylinder, choice_points: npt.NDArray +) -> openmc.Halfspace: + """ + choose a side of the Halfspace in the region of ZPlane and ZCylinder. + + Parameters + ---------- + surface + :class:`openmc.surface.Surface` of a openmc.ZPlane or openmc.ZCylinder + choice_points: np.ndarray of shape (N, 3) + a list of points representing the vertices of a convex polygon in RZ plane + + """ + x, y, z = np.array(to_cm(choice_points)).T + values = surface.evaluate([x, y, z]) + threshold = DTOL_CM + if isinstance(surface, openmc.ZCylinder): + threshold = 2 * DTOL_CM * surface.r + DTOL_CM**2 + + if all(values >= -threshold): + return +surface + if all(values <= threshold): + return -surface + + raise GeometryError(f"There are points on both sides of this {type(surface)}!") + + +# simplfying openmc.Intersection by associativity +def flat_intersection(region_list: Iterable[openmc.Region]) -> openmc.Intersection: + """ + Get the flat intersection of an entire list of regions. + e.g. (a (b c)) becomes (a b c) + """ + return openmc.Intersection( + intersection_dictionary(openmc.Intersection(region_list)).values() + ) + + +def intersection_dictionary(region: openmc.Region) -> dict[str, openmc.Region]: + """Get a dictionary of all of the elements that shall be intersected together, + applying the rule of associativity + """ + if isinstance(region, openmc.Halfspace): # termination condition + return {region.side + str(region.surface.id): region} + if isinstance(region, openmc.Intersection): + final_intersection = {} + for _r in region: + final_intersection.update(intersection_dictionary(_r)) + return final_intersection + return {str(region): region} + + +# simplfying openmc.Union by associativity +def flat_union(region_list: Iterable[openmc.Region]) -> openmc.Union: + """ + Get the flat union of an entire list of regions. + e.g. (a | (b | c)) becomes (a | b | c) + """ + return openmc.Union(union_dictionary(openmc.Union(region_list)).values()) + + +# TODO: Raise issue/papercut to check if simplifying the boolean expressions can yield +# speedup or not, and if so, we should attempt to simplify it further. +# E.g. the expression (-1 ((-1107 -1) | -1108)) can be simplified to (-1107 | -1108) -1; +# And don't even get me started on how much things can get simplified when ~ is involved. +# It is possible that boolean expressions get condensed appropriately before getting +# parsed onto openmc. I can't tell either way. + + +def union_dictionary(region: openmc.Region) -> dict[str, openmc.Region]: + """Get a dictionary of all of the elements that shall be unioned together, + applying the rule of associativity + """ + if isinstance(region, openmc.Halfspace): # termination condition + return {region.side + str(region.surface.id): region} + if isinstance(region, openmc.Union): + final_intersection = {} + for _r in region: + final_intersection.update(union_dictionary(_r)) + return final_intersection + return {str(region): region} + + +def round_up_next_openmc_ids(surface_step_size: int = 1000, cell_step_size: int = 100): + """ + Make openmc's surfaces' and cells' next IDs to be incremented to the next + pre-determined interval. + """ + openmc.Surface.next_id = ( + int(max(openmc.Surface.used_ids) / surface_step_size + 1) * surface_step_size + 1 + ) + openmc.Cell.next_id = ( + int(max(openmc.Cell.used_ids) / cell_step_size + 1) * cell_step_size + 1 + ) + + +def exterior_vertices(blanket, divertor) -> npt.NDArray: + """ + Get the 3D coordinates of every point at the outer boundary of the tokamak's + poloidal cross-section. + + Returns + ------- + coordinates + array of shape (N+1+n*M, 3), where N = number of blanket pre-cells, + M = number of divertor pre-cells, n = discretisation_level used when chopping + up the divertor in + :meth:`bluemira.radiation_transport.neutronics.DivertorWireAndExteriorCurve.make_divertor_pre_cell_array` + """ + return np.concatenate([ + blanket.exterior_vertices(), + divertor.exterior_vertices()[::-1], + ]) + + +def interior_vertices(blanket, divertor) -> npt.NDArray: + """ + Get the 3D coordinates of every point at the interior boundary of the tokamak's + poloidal cross-section + + Returns + ------- + coordinates + array of shape ((N+1)+sum(number of interior points of the divertor), 3), + where N = number of blanket pre-cells, M = number of divertor pre-cells. + Runs clockwise, beginning at the inboard blanket-divertor joining point. + """ + return np.concatenate([ + blanket.interior_vertices(), + divertor.interior_vertices()[::-1], + ]) + + +def make_universe_box( + csg, z_min: float, z_max: float, r_max: float, *, control_id: bool = False +): + """Box up the universe in a cylinder (including top and bottom).""" + bottom = csg.find_suitable_z_plane( + z_min, + boundary_type="vacuum", + surface_id=998 if control_id else None, + name="Universe bottom", + ) + top = csg.find_suitable_z_plane( + z_max, + boundary_type="vacuum", + surface_id=999 if control_id else None, + name="Universe top", + ) + universe_cylinder = openmc.ZCylinder( + r=to_cm(r_max), + surface_id=1000 if control_id else None, + boundary_type="vacuum", + name="Max radius of Universe", + ) + return -top & +bottom & -universe_cylinder + + +def make_radiation_shield_box( + csg, + z_min: float, + z_max: float, + r_max: float, + universe: openmc.region.Intersection, + materials: MaterialsLibrary, +): + """Define the radiation shield wall as a hollow of the universe box.""" + bottom_inner = csg.find_suitable_z_plane( + z_min, + name="Radiation shield bottom inner wall", + ) + + top_inner = csg.find_suitable_z_plane( + z_max, + name="Radiation shield top inner wall", + ) + + radial_cylinder_inboard = openmc.ZCylinder( + r=to_cm(r_max), + name="Radiation shield wall radial inboard", + ) + + rad_shield_inner = -top_inner & +bottom_inner & -radial_cylinder_inboard + + return openmc.Cell( + name="Radiation shield wall", + fill=materials.match_material(CellType.RadiationShield), + region=universe & ~rad_shield_inner, + ) + + +def make_coils( + csg, + solenoid_radius: float, + tf_coil_thick: float, + z_min: float, + z_max: float, + materials, +) -> tuple[openmc.Cell, list[openmc.Cell]]: + """ + Make tf coil and the central solenoid. The former wraps around the latter. + + Parameters + ---------- + solenoid_radius: + Central solenoid radius [m] + tf_coil_thick: + Thickness of the tf-coil, wrapped around the central solenoid [m] + z_max: + z-coordinate of the the top z-plane shared by both cylinders + (cs and tf coil) + z_min + z-coordinate of the the bottom z-plane shared by both cylinders + (cs and tf coil) + """ + if tf_coil_thick <= 0 or solenoid_radius <= 0: + raise GeometryError( + "Centrol column TF Coils and solenoid must have positive thicknesses!" + ) + solenoid = openmc.ZCylinder(r=to_cm(solenoid_radius)) + central_tf_coil = openmc.ZCylinder(r=to_cm(tf_coil_thick + solenoid_radius)) + top = csg.find_suitable_z_plane( + z_max, + [z_max - D_TOLERANCE, z_max + D_TOLERANCE], + name="Top of central solenoid", + ) + bottom = csg.find_suitable_z_plane( + z_min, + [z_min - D_TOLERANCE, z_min + D_TOLERANCE], + name="Bottom of central solenoid", + ) + central_solenoid = openmc.Cell( + name="Central solenoid", + fill=materials.match_material(CellType.CSCoil), + region=+bottom & -top & -solenoid, + ) + tf_coils = [ + openmc.Cell( + name="TF coil (sheath around central solenoid)", + fill=materials.match_material(CellType.TFCoil), + region=+bottom & -top & +solenoid & -central_tf_coil, + ) + ] + central_solenoid.volume = (top.z0 - bottom.z0) * np.pi * solenoid.r**2 + tf_coils[0].volume = ( + (top.z0 - bottom.z0) * np.pi * (central_tf_coil.r**2 - solenoid.r**2) + ) + return central_solenoid, tf_coils + + +def make_dividing_surface(csg, component): + """Surface that marks the end of the divertor/blanket's exterior.""" + exterior_pts = component.exterior_vertices()[:, ::2] + return csg.surface_from_2points(exterior_pts[0], exterior_pts[-1]) + + +def blanket_and_divertor_outer_regions( + csg, blanket, divertor, *, control_id: bool = False +) -> openmc.Region: + """ + Get the entire tokamak's poloidal cross-section (everything inside + self.geom.boundary) as an openmc.Region. + """ + dividing_surface = make_dividing_surface(csg, blanket) + blanket_outer = csg.region_from_surface_series( + [*blanket.exterior_surfaces(), dividing_surface], + blanket.exterior_vertices(), + control_id=control_id, + ) + divertor_outer = csg.region_from_surface_series( + [*divertor.exterior_surfaces(), dividing_surface], + divertor.exterior_vertices(), + control_id=control_id, + ) + return blanket_outer, divertor_outer + + +def plasma_void(csg, blanket, divertor, *, control_id: bool = False) -> openmc.Region: + """Get the plasma chamber's poloidal cross-section""" + blanket_interior_pts = blanket.interior_vertices() + dividing_surface = make_dividing_surface(csg, blanket) + if not is_convex(blanket_interior_pts): + raise GeometryError( + f"{blanket} interior (blanket_outline's inner_curve) needs to be convex!" + ) + plasma = csg.region_from_surface_series( + [*blanket.interior_surfaces(), dividing_surface], + blanket_interior_pts, + control_id=control_id, + ) + + divertor_exterior_vertices = divertor.exterior_vertices() + if not is_convex(divertor_exterior_vertices): + raise GeometryError(f"{divertor} exterior needs to be convex!") + exhaust_including_divertor = csg.region_from_surface_series( + [*divertor.exterior_surfaces(), dividing_surface], + divertor_exterior_vertices, + control_id=control_id, + ) + + return flat_union([plasma, exhaust_including_divertor]) & ~divertor.exclusion_zone( + control_id=control_id + ) + + +def make_void_cells( + csg, + universe: openmc.region.Intersection, + blanket: BlanketCellArray, + divertor: DivertorCellArray, + central_solenoid: openmc.Cell, + tf_coils: list[openmc.Cell] | None, + rad_shield: openmc.Cell | None = None, + *, + control_id: bool = False, +): + """Make the plasma chamber and the outside ext_void. This should be called AFTER + the blanket and divertor cells are created. + """ + blanket_silhouette, divertor_silhouette = blanket_and_divertor_outer_regions( + csg, blanket, divertor, control_id=control_id + ) + void_region = universe & ~blanket_silhouette & ~divertor_silhouette + if tf_coils: + void_region &= ~tf_coils[0].region + if central_solenoid: + void_region &= ~central_solenoid.region + if rad_shield: + void_region &= ~rad_shield.region + + return ( + openmc.Cell( + region=plasma_void(csg, blanket, divertor, control_id=control_id), + fill=None, + name="Plasma void", + ), + openmc.Cell( + region=flat_intersection(void_region), + fill=None, + name="Exterior void", + ), + ) + + +def make_cell_arrays( + pre_cell_reactor: NeutronicsReactor, + csg: BluemiraNeutronicsCSG, + materials: MaterialsLibrary, + *, + control_id: bool = False, +) -> CellStage: + """Make pre-cell arrays for the blanket and the divertor. + + Parameters + ---------- + materials: + library containing information about the materials + tokamak_dimensions: + A parameter + :class:`bluemira.radiation_transport.neutronics.params.TokamakDimensions`, + Specifying the dimensions of various layers in the blanket, divertor, and + central solenoid. + control_id: bool + Whether to set the blanket Cells and surface IDs by force or not. + With this set to True, it will be easier to understand where each cell came + from. However, it will lead to warnings and errors if a cell/surface is + generated to use a cell/surface ID that has already been used respectively. + Keep this as False if you're running openmc simulations multiple times in one + session. + """ + # determine universe_box + + z_max, z_min, r_max, r_min = pre_cell_reactor.half_bounding_box + + z_min_adj = z_min - D_TOLERANCE + z_max_adj = z_max + D_TOLERANCE + r_max_adj = r_max + D_TOLERANCE + + rad_shield_wall_tk = pre_cell_reactor.tokamak_dimensions.rad_shield.wall + + # make the universe box, incorporates the radiation shield wall + universe = make_universe_box( + csg, + z_min_adj - rad_shield_wall_tk, + z_max_adj + rad_shield_wall_tk, + r_max_adj + rad_shield_wall_tk, + control_id=control_id, + ) + + blanket = BlanketCellArray.from_pre_cell_array( + pre_cell_reactor.blanket, + materials, + pre_cell_reactor.tokamak_dimensions, + csg, + control_id=control_id, + ) + + # change the cell and surface id register before making the divertor. + # (ids will only count up from here.) + if control_id: + round_up_next_openmc_ids() + + divertor = DivertorCellArray.from_pre_cell_array( + pre_cell_reactor.divertor, + materials, + pre_cell_reactor.tokamak_dimensions.divertor, + csg=csg, + override_start_end_surfaces=(blanket[0].ccw_surface, blanket[-1].cw_surface), + # ID cannot be controlled at this point. + ) + + # make the plasma cell and the exterior void. + if control_id: + round_up_next_openmc_ids() + + cs, tf = make_coils( + csg, + r_min - pre_cell_reactor.tokamak_dimensions.cs_coil.thickness, + pre_cell_reactor.tokamak_dimensions.cs_coil.thickness, + z_min_adj, + z_max_adj, + materials, + ) + # make the radiation shield wall + # which is a hollow of the universe box + rad_shield = make_radiation_shield_box( + csg, + z_min_adj, + z_max_adj, + r_max_adj, + universe, + materials, + ) + plasma, ext_void = make_void_cells( + csg, + universe=universe, + blanket=blanket, + divertor=divertor, + central_solenoid=cs, + tf_coils=tf, + rad_shield=rad_shield, + control_id=control_id, + ) + + cell_stage = CellStage( + blanket=blanket, + divertor=divertor, + tf_coils=tf, + cs_coil=cs, + plasma=plasma, + radiation_shield=rad_shield, + ext_void=ext_void, + universe=universe, + ) + cell_stage.set_volumes() + + return cell_stage + + +class BluemiraNeutronicsCSG: + """Container for CSG planes to enable reuse of planes, very eco friendly""" + + def __init__(self): + # it's called a hangar because it's where the planes are parked ;) + self.hangar = {} + + def surface_from_2points( + self, + point1: npt.NDArray[np.float64], + point2: npt.NDArray[np.float64], + surface_id: int | None = None, + name: str = "", + ) -> openmc.Surface | openmc.model.ZConeOneSided | None: + """ + Create either a cylinder, a cone, or a surface from 2 points using only the + rz coordinates of any two points on it. + + Parameters + ---------- + point1, point2: + any two non-trivial (i.e. cannot be the same) points on the rz cross-section + of the surface, each containing the r and z coordinates + Units: [m] + surface_id, name: + see openmc.Surface + + Returns + ------- + surface: + if the two points provided are redundant: don't return anything, as this is a + single point pretending to be a surface. This will come in handy for handling + the creation of BlanketCells made with 3 surfaces rather than 4. + """ + point1, point2 = to_cm((point1, point2)) + dr, dz = point2 - point1 + if np.isclose(dr, 0, rtol=0, atol=DTOL_CM): + if np.isclose(dz, 0, rtol=0, atol=DTOL_CM): + return None + return openmc.ZCylinder(r=point1[0], surface_id=surface_id, name=name) + if np.isclose(dz, 0, rtol=0, atol=DTOL_CM): + z = point1[-1] + self.hangar[z] = openmc.ZPlane(z0=z, surface_id=surface_id, name=name) + return self.hangar[z] + slope = dz / dr + z_intercept = -slope * point1[0] + point1[-1] + return openmc.ZCone( + z0=z_intercept, r2=slope**-2, surface_id=surface_id, name=name + ) + + def surface_from_straight_line( + self, + straight_line_info: StraightLineInfo, + surface_id: int | None = None, + name: str = "", + ): + """Create a surface to match the straight line info provided.""" + start_end = np.array(straight_line_info[:2])[:, ::2] + return self.surface_from_2points(*start_end, surface_id=surface_id, name=name) + + def surfaces_from_info_list(self, wire_info_list: WireInfoList, name: str = ""): + """ + Create a list of surfaces using a list of wire infos. + + Parameters + ---------- + wire_info_list + List of wires + name + This name will be *reused* across all of the surfaces created in this list. + """ + surface_list = [] + for wire in wire_info_list.info_list: + info = wire.key_points + plane_cone_cylinder = self.surface_from_straight_line(info, name=name) + if isinstance(info, CircleInfo): + torus = torus_from_circle(info.center, info.radius, name=name) + # will need the openmc.Union of these two objects later. + surface_list.append((plane_cone_cylinder, torus)) + else: + surface_list.append((plane_cone_cylinder,)) + return tuple(surface_list) + + def find_suitable_z_plane( + self, + z0: float, + z_range: Iterable[float] | None = None, + surface_id: int | None = None, + name: str = "", + **kwargs, + ): + """ + Find a suitable z from the hangar, or create a new one if no matches are found. + + Parameters + ---------- + z0: + The height of the plane, if we need to create it. Unit: [m] + z_range: + If we a suitable z-plane already exists, then we only accept it if it lies + within this range of z. Unit: [m] + surface_id, name: + See openmc.Surface + """ + if z_range is not None: + z_min, z_max = min(z_range), max(z_range) + for key in self.hangar: + if z_min <= key <= z_max: + return self.hangar[key] # return the first match + self.hangar[z0] = openmc.ZPlane( + z0=to_cm(z0), surface_id=surface_id, name=name, **kwargs + ) + return self.hangar[z0] + + def choose_region_cone( + self, + surface: openmc.ZCone, + choice_points: npt.NDArray, + *, + control_id: bool = False, + ) -> openmc.Region: + """ + choose the region for a ZCone. + When reading this function's code, bear in mind that a Z cone can be separated + into 3 parts: + + A. the upper cone (evaluates to negative), + B. outside of the cone (evaluates to positive), + C. the lower cone (evaluates to negative). + + We have to account for the following cases: + + +------------+---------+------------+ + | upper cone | outside | lower cone | + +============+=========+============+ + | Y | N | N | + +------------+---------+------------+ + | Y | Y | N | + +------------+---------+------------+ + | N | Y | N | + +------------+---------+------------+ + | N | Y | Y | + +------------+---------+------------+ + | N | N | Y | + +------------+---------+------------+ + + All other cases should raise an error. + The tricky part to handle is the floating point precision problem. + It's possible that th every point used to create the cone does not lie on the + cone/ lies on the wrong side of the cone. + Hence the first step is to shrink the choice_points by SHRINK_DISTANCE + towards the centroid. + + Parameters + ---------- + surface: + where all points are expected to be excluded from at least one of its two + one-sided cones. + choice_points: + An array of points that, after choosing the appropriate region, + should all lie in the chosen region. + control_id: + When an ambiguity plane is needed, we ned to create a surface. if + control_id = True, then this would force the surface to be created with + id = 1000 + the id of the cone. This is typically only used so that + we have full control of (and easily understandable records of) + every surfaces' ID. + Thus elsewhere in the code, most other classes/methods turns control_id + on when cell_ids are also provided (proving intention on controlling IDs + of OpenMC objects). + + Returns + ------- + region + openmc.Region, specifically (openmc.Halfspace) or + (openmc.Union of 2 openmc.Halfspaces) + """ + # shrink to avoid floating point number comparison imprecision issues. + # Especially important when the choice point sits exactly on the surface. + centroid = np.mean(choice_points, axis=0) + step_dir = centroid - choice_points + unit_step_dir = (step_dir.T / np.linalg.norm(step_dir, axis=1)).T + shrunken_points = choice_points + SHRINK_DISTANCE * unit_step_dir + + # # Alternative + # shrunken_points = (choice_points + 0.01 * centroid) / 1.01 + # take one step towards the centroid = 0.1 cm + + x, y, z = np.array(to_cm(shrunken_points)).T + values = surface.evaluate([x, y, z]) + middle = values > 0 + if all(middle): # exist outside of cone + return +surface + + z_dist = z - surface.z0 + upper_cone = np.logical_and(~middle, z_dist > 0) + lower_cone = np.logical_and(~middle, z_dist < 0) + + # all of the following cases requires a plane to be made. + if control_id: + surface_id = 1000 + surface.id + else: + surface_id = max(openmc.Surface.used_ids) + surface_id += 1000 + surface_id % 1000 + + if all(upper_cone): + # everything in the upper cone. + # the highest we can cut is at the lowest z. + plane = self.find_suitable_z_plane( + to_m(surface.z0), + to_m([surface.z0 - DTOL_CM, min(z) - DTOL_CM]), + surface_id=surface_id, + name=f"Ambiguity plane for cone {surface.id}", + ) + return -surface & +plane + if all(lower_cone): + # everything in the lower cone + # the lowest we can cut is at the highest z. + plane = self.find_suitable_z_plane( + to_m(surface.z0), + to_m([max(z) + DTOL_CM, surface.z0 + DTOL_CM]), + surface_id=surface_id, + name=f"Ambiguity plane for cone {surface.id}", + ) + return -surface & -plane + if all(np.logical_or(upper_cone, lower_cone)): + raise GeometryError( + "Both cones have vertices lying inside! Cannot compute a contiguous " + "region that works for both. Check if polygon is convex?" + ) + # In this rare case, make its own plane. + plane = self.find_suitable_z_plane( + to_m(surface.z0), + to_m([surface.z0 + DTOL_CM, surface.z0 - DTOL_CM]), + surface_id=surface_id, + name=f"Ambiguity plane for cone {surface.id}", + ) + if all(np.logical_or(upper_cone, middle)): + return +surface | +plane + if all(np.logical_or(lower_cone, middle)): + return +surface | -plane + raise GeometryError("Can't have points in upper cone, lower cone AND outside!") + + def choose_region( + self, + surface: openmc.Surface + | tuple[openmc.Surface] + | tuple[openmc.Surface | openmc.ZTorus], + vertices_array: npt.NDArray, + *, + control_id: bool = False, + ) -> openmc.Region: + """ + Pick the correct region of the surface that includes all of the points in + vertices_array. + + Parameters + ---------- + surface + Either a :class:`openmc.Surface`, or a 1-tuple or 2-tuple of + :class:`openmc.Surface`. If it is a tuple, the first element is always a + :class:`openmc.ZPlane`/:class:`openmc.ZCone`/:class:`openmc.ZCylinder`; + the second element (if present) is always :class:`openmc.ZTorus`. + vertices_array + array of shape (?, 3), that the final region should include. + control_id + Passed as argument onto + :meth:`~bluemira.radiation_transport.neutronics.make_csg.BluemiraNeutronicsCSG.choose_region_cone` + + Returns + ------- + An openmc.Region built from surface provided and includes all of these + """ + if isinstance(surface, openmc.ZPlane | openmc.ZCylinder): + return choose_plane_cylinders(surface, vertices_array) + + if isinstance(surface, openmc.ZCone): + return self.choose_region_cone( + surface, vertices_array, control_id=control_id + ) + + if isinstance(surface, tuple): + chosen_first_region = self.choose_region( + surface[0], vertices_array, control_id=control_id + ) + # cone, or cylinder, or plane, or (cone | ambiguity_surface) + if len(surface) == 1: + return chosen_first_region + # -surface[1] = inside of torus + return flat_union((chosen_first_region, -surface[1])) + raise NotImplementedError( + f"Surface {type(surface)} is not ready for use in the axis-symmetric case!" + ) + + def region_from_surface_series( + self, + series_of_surfaces: Sequence[ + openmc.Surface | tuple[openmc.Surface, openmc.ZTorus | None] | None + ], + vertices_array: npt.NDArray, + *, + control_id: bool = False, + ) -> openmc.Intersection: + """ + Switch between choose_region() and choose_region_from_tuple_of_surfaces() + depending on the type of each element in the series_of_surfaces. + + Parameters + ---------- + series_of_surfaces + Each of them can be a None, a 1-tuple of surface, a 2-tuple of surfaces, or a + surface. For the last 3 options, see + :func:`~bluemira.radiation_transport.neutronics.make_csg.BluemiraNeutronicsCSG.choose_region` + for more. + + vertices_array + array of shape (?, 3), where every single point should be included by, or at + least on the edge of the returned Region. + control_id + Passed as argument onto + :meth:`~bluemira.radiation_transport.neutronics.make_csg.BluemiraNeutronicsCSG.choose_region_cone` + + Returns + ------- + intersection_region: + openmc.Intersection of a list of + [(openmc.Halfspace) or (openmc.Union of openmc.Halfspace)] + """ + return flat_intersection([ + self.choose_region(surface, vertices_array, control_id=control_id) + for surface in series_of_surfaces + if surface is not None + ]) + + +class BlanketCell(openmc.Cell): + """ + A generic blanket cell that forms the base class for the five specialised types of + blanket cells. + + It's a special case of openmc.Cell, in that it has 3 to 4 surfaces + (mandatory surfaces: exterior_surface, ccw_surface, cw_surface; + optional surface: interior_surface), and it is more wieldy because we don't have to + specify the relevant half-space for each surface; instead the corners of the cell + is provided by the user, such that the appropriate regions are chosen. + """ + + def __init__( + self, + exterior_surface: openmc.Surface, + ccw_surface: openmc.Surface, + cw_surface: openmc.Surface, + interior_surface: openmc.Surface | None, + vertices: Coordinates, + csg: BluemiraNeutronicsCSG, + cell_id: int | None = None, + name: str = "", + fill: openmc.Material | None = None, + ): + """ + Create the openmc.Cell from 3 to 4 surfaces and an example point. + + Parameters + ---------- + exterior_surface + Surface on the exterior side of the cell + ccw_surface + Surface on the ccw wall side of the cell + cw_surface + Surface on the cw wall side of the cell + vertices + A list of points. Could be 2D or 3D. + interior_surface + Surface on the interior side of the cell + cell_id + see :class:`openmc.Cell` + name + see :class:`openmc.Cell` + fill + see :class:`openmc.Cell` + """ + self.exterior_surface = exterior_surface + self.ccw_surface = ccw_surface + self.cw_surface = cw_surface + self.interior_surface = interior_surface + self.vertex = vertices + self.csg = csg + + super().__init__( + cell_id=cell_id, + name=name, + fill=fill, + region=self.csg.region_from_surface_series( + [exterior_surface, ccw_surface, cw_surface, interior_surface], + self.vertex.T, # We just assume it is convex + control_id=bool(cell_id), + ), + ) + + self.volume = to_cm3(polygon_revolve_signed_volume(vertices[::2].T)) + if self.volume <= 0: + raise GeometryError("Wrong ordering of vertices!") + + +class BlanketCellStack: + """ + A stack of openmc.Cells, first cell is closest to the interior and last cell is + closest to the exterior. They should all be situated at the same poloidal angle. + """ + + def __init__(self, cell_stack: list[BlanketCell]): + """ + The shared surface between adjacent cells must be the SAME one, i.e. same id and + hash, not just identical. + They must share the SAME counter-clockwise surface and the SAME clockwise surface + (left and right side surfaces of the stack, for the stack pointing straight up). + + Because the bounding_box function is INCORRECT, we can't perform a check on the + bounding box to confirm that the stack is linearly increasing/decreasing in xyz + bounds. + + Series of cells: + + SurfaceCell + FirstWallCell + BreedingZoneCell + ManifoldCell + VacuumVesselCell + """ + self.cell_stack = cell_stack + for int_cell, ext_cell in pairwise(cell_stack): + if int_cell.exterior_surface is not ext_cell.interior_surface: + raise ValueError("Expected a contiguous stack of cells!") + + def __len__(self) -> int: + """Number of cells in stack""" + return len(self.cell_stack) + + def __getitem__(self, index_or_slice) -> list[BlanketCell] | BlanketCell: + """Get cell from stack""" + return self.cell_stack[index_or_slice] + + def __iter__(self) -> Iterator[BlanketCell]: + """Iterator for BlanketCellStack""" + return iter(self.cell_stack) + + def __repr__(self) -> str: + """String representation""" + return ( + super() + .__repr__() + .replace(" at ", f" of {len(self.cell_stack)} BlanketCells at ") + ) + + @staticmethod + def check_cut_point_ordering( + cut_point_series: npt.NDArray[np.float64], + direction_vector: npt.NDArray[np.float64], + location_msg: str = "", + ): + """ + Parameters + ---------- + cut_point_series: + array of shape (M+1, 2) where M = number of cells in the blanket cell stack + (i.e. number of layers in the blanket). Each point has two dimensions + direction_vector: + direction that these points are all supposed to go towards. + """ + direction = direction_vector / np.linalg.norm(direction_vector) + projections = np.dot(np.array(cut_point_series)[:, [0, -1]], direction) + if not is_monotonically_increasing(projections): + raise GeometryError(f"Some surfaces crosses over each other! {location_msg}") + + @property + def interior_surface(self): + """Get interior surface""" + return self.cell_stack[0].interior_surface + + @property + def exterior_surface(self): + """Get exterior surface""" + return self.cell_stack[-1].exterior_surface + + @property + def ccw_surface(self): + """Get counter clockwise surface""" + return self.cell_stack[0].ccw_surface + + @property + def cw_surface(self): + """Get clockwise surface""" + return self.cell_stack[0].cw_surface + + @property + def interfaces(self): + """ + All of the radial surfaces, including the innermost (exposed to plasma) and + outermost (facing vacuum vessel); arranged in that order (from innermost to + outermost). + """ + if not hasattr(self, "_interfaces"): + self._interfaces = [cell.interior_surface for cell in self.cell_stack] + self._interfaces.append(self.cell_stack[-1].exterior_surface) + return self._interfaces + + def get_overall_region( + self, csg: BluemiraNeutronicsCSG, *, control_id: bool = False + ) -> openmc.Region: + """ + Calculate the region covering the entire cell stack. + + Parameters + ---------- + control_id + Passed as argument onto + :meth:`~bluemira.radiation_transport.neutronics.make_csg.BluemiraNeutronicsCSG.region_from_surface_series` + """ + vertices = np.vstack(( + self.cell_stack[0].vertex.T[(1, 2),], + self.cell_stack[-1].vertex.T[(3, 0),], + )) + if not is_convex(vertices): + raise GeometryError(f"{self}'s vertices need to be convex!") + return csg.region_from_surface_series( + [ + self.exterior_surface, + self.ccw_surface, + self.cw_surface, + self.interior_surface, + ], + vertices, + control_id=control_id, + ) + + @classmethod + def from_pre_cell( + cls, + pre_cell: PreCell, + ccw_surface: openmc.Surface, + cw_surface: openmc.Surface, + depth_series: npt.NDArray, + csg: BluemiraNeutronicsCSG, + fill_lib: MaterialsLibrary, + *, + inboard: bool, + blanket_stack_num: int | None = None, + ): + """ + Create a CellStack using a precell and TWO surfaces that sandwiches that precell. + + Parameters + ---------- + pre_cell + An instance of :class:`~PreCell` + ccw_surf + An instance of :class:`openmc.surface.Surface` + cw_surf + An instance of :class:`openmc.surface.Surface` + depth_series + a series of floats corresponding to the N-2 interfaces between the N-1 + layers, whereas the N-th layer is the vacuum vessel (and the pre-cell has + already stored the thickness for that). + Each float represents how deep into the blanket (i.e. how many [cm] into the + first wall we need to drill, from the plasma facing surface) to hit that + interface layer. + fill_lib + :class:`~MaterialsLibrary` so that it separates into .inboard, .outboard, + .divertor, .tf_coil_windings, etc. + inboard + boolean denoting whether this cell is inboard or outboard + blanket_stack_num + An optional number indexing the current stack. Used for labelling. + If None: we will not be controlling the cell and surfaces id. + """ + # check exterior wire is correct + ext_curve_comp = pre_cell.exterior_wire.shape.OrderedEdges + if len(ext_curve_comp) != 1: + raise TypeError("Incorrect type of BluemiraWire parsed in.") + if not ext_curve_comp[0].Curve.TypeId.startswith("Part::GeomLine"): + raise NotImplementedError("Not ready to make curved-line cross-section yet!") + + # 1. Calculate cut points required to make the surface stack, without actually + # creating the surfaces. + # shape (M+1, 2, 2) + wall_cut_pts = np.asarray([ + pre_cell.cell_walls.starts, + *( + pre_cell.get_cell_wall_cut_points_by_thickness(interface_depth) + for interface_depth in depth_series + ), + np.array([pre_cell.vv_point[:, 0], pre_cell.vv_point[:, 1]]), + pre_cell.cell_walls.ends, + ]) + # 1.1 perform sanity check + directions = np.diff(pre_cell.cell_walls, axis=1) # shape (2, 1, 2) + dirs = directions[:, 0, :] + i = "(unspecified)" if blanket_stack_num is None else blanket_stack_num + cls.check_cut_point_ordering( + wall_cut_pts[:, 0], + dirs[0], + location_msg=f"Occuring in cell stack {i}'s CCW wall", + ) + cls.check_cut_point_ordering( + wall_cut_pts[:, 1], + dirs[1], + location_msg=f"Occuring in cell stack {i}'s CW wall", + ) + # 2. Accumulate the corners of each cell. + vertices = [ + np.array([ + [out[1, 0], inn[1, 0], inn[0, 0], out[0, 0]], + np.full(4, 0), + [out[1, 1], inn[1, 1], inn[0, 1], out[0, 1]], + ]) + for inn, out in pairwise(wall_cut_pts) + ] + # shape (M, 2, 2) + projection_ccw = wall_cut_pts[:, 0] @ dirs[0] / np.linalg.norm(dirs[0]) + projection_cw = wall_cut_pts[:, 1] @ dirs[1] / np.linalg.norm(dirs[1]) + layer_mask = np.array( + [ + not (ccw_depth <= DTOL_CM and cw_depth <= DTOL_CM) + for (ccw_depth, cw_depth) in zip( + np.diff(projection_ccw), np.diff(projection_cw), strict=True + ) + ], + dtype=bool, + ) # shape (M,) + + # 3. Choose the ID of the stack's surfaces and cells. + if blanket_stack_num is not None: + # Note: all IDs must be natural number, i.e. integer > 0. + # So we're using an indexing scheme that starts from 1. + # len=M + cell_ids = [10 * blanket_stack_num + v + 1 for v in range(len(vertices))] + # left (ccw_surface) surface had already been created, and since our indexing + # scheme starts from 1, therefore we're using +2 in the following line. + # len=M+1 + surface_ids = [ + 10 * blanket_stack_num + v + 2 for v in range(len(wall_cut_pts)) + ] + else: + cell_ids = [None] * len(vertices) # len=M + surface_ids = [None] * len(wall_cut_pts) # len=M+1 + + # 4. create the actual openmc.Surfaces and Cells. + int_surf = ( + csg.surface_from_2points( + *wall_cut_pts[0], + surface_id=surface_ids[0], + name=f"plasma-facing surface of blanket cell stack {i}", + ) + if pre_cell.interior_wire + else None + ) # account for the case. + + cell_stack = [] + # k = range(0, M - layer_mask == False) + cell_types = np.array( + ( + CellType.BlanketSurface, + CellType.BlanketFirstWall, + CellType.BlanketBreedingZone, + CellType.BlanketManifold, + CellType.VacuumVessel, + ), + dtype=object, + )[layer_mask] + for k, (points, cell_type) in enumerate( + zip(wall_cut_pts[1:][layer_mask], cell_types, strict=False) + ): + j = k + 1 # = range(1, M+1) + if j > 1: + int_surf.name = ( + f"{cell_type}-{cell_type.name} " + f"interface boundary of blanket cell stack {i}" + ) + ext_surf = csg.surface_from_2points( + *points, + surface_id=surface_ids[j], # up to M+1 + ) + cell_stack.append( + BlanketCell( + exterior_surface=ext_surf, + ccw_surface=ccw_surface, + cw_surface=cw_surface, + interior_surface=int_surf, + vertices=vertices[k], # up to M + csg=csg, + cell_id=cell_ids[k], # up to M + name=f"{cell_type.name} of blanket cell stack {i}", + fill=fill_lib.match_material(cell_type, inboard=inboard), + ) + ) + int_surf = ext_surf + int_surf.name = "vacuum-vessel-facing surface" + + return cls(cell_stack) + + +class BlanketCellArray: + """ + An array of BlanketCellStack. Interior and exterior curve are both assumed convex. + + Parameters + ---------- + blanket_cell_array + a list of BlanketCellStack + + """ + + def __init__( + self, blanket_cell_array: list[BlanketCellStack], csg: BluemiraNeutronicsCSG + ): + """ + Create array from a list of BlanketCellStack + """ + self.blanket_cell_array = blanket_cell_array + self.poloidal_surfaces = [self.blanket_cell_array[0].ccw_surface] + self.radial_surfaces = [] + self.csg = csg + for i, this_stack in enumerate(self.blanket_cell_array): + self.poloidal_surfaces.append(this_stack.cw_surface) + self.radial_surfaces.append(this_stack.interfaces) + + # check neighbouring cells share the same lateral surface + if i != len(self.blanket_cell_array) - 1: + next_stack = self.blanket_cell_array[i + 1] + if this_stack.cw_surface is not next_stack.ccw_surface: + raise GeometryError( + f"Neighbouring BlanketCellStack [{i}] and " + f"[{i + 1}] are not aligned!" + ) + + def __len__(self) -> int: + """Number of cell stacks""" + return len(self.blanket_cell_array) + + def __getitem__(self, index_or_slice) -> list[BlanketCellStack] | BlanketCellStack: + """Get cell stack""" + return self.blanket_cell_array[index_or_slice] + + def __iter__(self) -> Iterator[BlanketCellStack]: + """Iterator for BlanketCellArray""" + return iter(self.blanket_cell_array) + + def __repr__(self) -> str: + """String representation""" + return ( + super() + .__repr__() + .replace(" at ", f" of {len(self.blanket_cell_array)} BlanketCellStacks at ") + ) + + def exterior_vertices(self) -> npt.NDArray: + """ + Returns all of the tokamak's poloidal cross-section's outside corners' + coordinates, in 3D. + + Returns + ------- + exterior_vertices: + array of shape (N+1, 3) arranged clockwise (inboard to outboard). + """ + return np.asarray([ + self.blanket_cell_array[0][-1].vertex.T[Vert.exterior_start], + *( + stack[-1].vertex.T[Vert.exterior_end] + for stack in self.blanket_cell_array + ), + ]) + + def interior_vertices(self) -> npt.NDArray: + """ + Returns all of the tokamak's poloidal cross-section's inside corners' + coordinates, in 3D. + + Parameters + ---------- + interior_vertices: + array of shape (N+1, 3) arranged clockwise (inboard to outboard). + """ + return np.asarray([ + self.blanket_cell_array[0][0].vertex.T[Vert.interior_start], + *(stack[0].vertex.T[Vert.interior_end] for stack in self.blanket_cell_array), + ]) + + def interior_surfaces(self) -> list[openmc.Surface]: + """ + Get all of the innermost (plasm-facing) surface. + Runs clockwise. + """ + return [surf_stack[0] for surf_stack in self.radial_surfaces] + + def exterior_surfaces(self) -> list[openmc.Surface]: + """ + Get all of the outermost (vacuum-vessel-facing) surface. + Runs clockwise. + """ + return [surf_stack[-1] for surf_stack in self.radial_surfaces] + + def exclusion_zone(self, *, control_id: bool = False) -> openmc.Region: + """ + Get the exclusion zone AWAY from the plasma. + Usage: plasma_region = openmc.Union(..., ~self.exclusion_zone(), ...) + Assumes that all of the panels (interior surfaces) together forms a convex hull. + + Parameters + ---------- + control_id + Passed as argument onto + :meth:`~bluemira.radiation_transport.neutronics.make_csg.BluemiraNeutronicsCSG.region_from_surface_series`. + """ + exclusion_zone_by_stack = [] + for stack in self.blanket_cell_array: + stack_vertices = np.vstack(( + stack[0].vertex.T[(1, 2),], + stack[-1].vertex.T[(3, 0),], + )) + if not is_convex(stack_vertices): + raise GeometryError(f"{self}'s vertices need to be convex!") + + exclusion_zone_by_stack.append( + self.csg.region_from_surface_series( + [stack.cw_surface, stack.ccw_surface, stack.interior_surface], + stack_vertices, + control_id=control_id, + ) + ) + return openmc.Union(exclusion_zone_by_stack) + + @classmethod + def from_pre_cell_array( + cls, + pre_cell_array: PreCellArray, + materials: MaterialsLibrary, + blanket_dimensions: TokamakDimensions, + csg: BluemiraNeutronicsCSG, + *, + control_id: bool = False, + ) -> BlanketCellArray: + """ + Create a BlanketCellArray from a + :class:`~bluemira.radiation_transport.neutronics.make_pre_cell.PreCellArray`. + This method assumes itself is the first method to be run to create cells in the + :class:`~openmc.Universe.` + + Parameters + ---------- + pre_cell_array + PreCellArray + materials + :class:`~MaterialsLibrary` so that it separates into .inboard, .outboard, + .divertor, .tf_coil_windings, etc. + blanket_dimensions + :class:`bluemira.radiation_transport.neutronics.params.TokamakDimensions` + recording the dimensions of the blanket in SI units (unit: [m]). + control_id + Passed as argument onto + :meth:`~bluemira.radiation_transport.neutronics.make_csg.BluemiraNeutronicsCSG.region_from_surface_series`. + """ + cell_walls = pre_cell_array.cell_walls + # TODO: when contorl_id, we're forced to start at id=0 + + # left wall + ccw_surf = csg.surface_from_2points( + *cell_walls[0], + surface_id=1 if control_id else None, + name="Blanket cell wall 0", + ) + cell_array = [] + for i, (pre_cell, cw_wall) in enumerate( + zip(pre_cell_array.pre_cells, cell_walls[1:], strict=True) + ): + # right wall + cw_surf = csg.surface_from_2points( + *cw_wall, + surface_id=1 + 10 * (i + 1) if control_id else None, + name=f"Blanket cell wall of blanket cell stack {i + 1}", + ) + + stack = BlanketCellStack.from_pre_cell( + pre_cell, + ccw_surf, + cw_surf, + get_depth_values(pre_cell, blanket_dimensions), + csg=csg, + fill_lib=materials, + inboard=check_inboard_outboard(pre_cell, blanket_dimensions), + blanket_stack_num=i if control_id else None, + ) + cell_array.append(stack) + ccw_surf = cw_surf + + return cls(cell_array, csg) + + +class DivertorCell(openmc.Cell): + """ + A generic Divertor cell forming either the (inner target's/outer target's/ + dome's) (surface/ bulk). + """ + + def __init__( + self, + exterior_surfaces: list[tuple[openmc.Surface]], + cw_surface: openmc.Surface, + ccw_surface: openmc.Surface, + interior_surfaces: list[tuple[openmc.Surface]], + exterior_wire: WireInfoList, + interior_wire: WireInfoList, + csg: BluemiraNeutronicsCSG, + subtractive_region: openmc.Region | None = None, + cell_id: int | None = None, + name: str = "", + fill: openmc.Material | None = None, + ): + """Create a cell from exterior_surface""" + self.exterior_surfaces = exterior_surfaces + self.cw_surface = cw_surface + self.ccw_surface = ccw_surface + self.interior_surfaces = interior_surfaces + self.exterior_wire = exterior_wire + self.interior_wire = interior_wire + self.csg = csg + + region = self.csg.region_from_surface_series( + [ + self.cw_surface, + self.ccw_surface, + *self.exterior_surfaces, + *self.interior_surfaces, + ], + self.get_all_vertices(), + control_id=bool(cell_id), + ) + + if subtractive_region: + region &= ~subtractive_region + super().__init__(cell_id=cell_id, name=name, fill=fill, region=region) + self.volume = self.get_volume() + + @property + def outline(self): + """ + Make the outline into a BluemiraWire. This method is created solely for the + purpose of calculating the volume. + + This is slow but it is accurate and works well. + """ + if not hasattr(self, "_outline"): + wire_list = [ + self.interior_wire.restore_to_wire(), + self.exterior_wire.restore_to_wire(), + ] + # make the connecting wires so that BluemiraWire doesn't moan about + # having a discontinuous wires. + i = self.interior_wire.get_3D_coordinates() + e = self.exterior_wire.get_3D_coordinates() + if not np.array_equal(e[-1], i[0]): + wire_list.insert(0, make_polygon([e[-1], i[0]])) + if not np.array_equal(i[-1], e[0]): + wire_list.insert(-1, make_polygon([i[-1], e[0]])) + self._outline = BluemiraWire(wire_list) + return self._outline + + def get_volume(self): + """ + Get the volume using the BluemiraWire of its own outline. + """ + half_solid = BluemiraSolid(revolve_shape(self.outline)) + cm3_volume = to_cm3(half_solid.volume * 2) + if cm3_volume <= 0: + raise GeometryError("Volume (as calculated by FreeCAD) is negative!") + return cm3_volume + + def get_all_vertices(self) -> npt.NDArray: + """ + Get all of the vertices of this cell, which should help us find its convex hull. + """ + return np.concatenate([ + self.exterior_wire.get_3D_coordinates(), + self.interior_wire.get_3D_coordinates(), + ]) + + def exclusion_zone( + self, + *, + away_from_plasma: bool = True, + control_id: bool = False, + additional_test_points: npt.NDArray | None = None, + ) -> openmc.Region: + """ + Get the exclusion zone of a semi-CONVEX cell. + + This can only be validly used: + + If away_from_plasma=True, then the interior side of the cell must be convex. + If away_from_plasma=False, then the exterior_side of the cell must be convex. + + Usage: + + next_cell_region = flat_intersection(..., ~this_cell.exclusion_zone()) + + Parameters + ---------- + control_id + Passed as argument onto + :func:`~bluemira.radiation_transport.neutronics.make_csg.region_from_surface_series` + """ + if away_from_plasma: + vertices_array = self.interior_wire.get_3D_coordinates() + if additional_test_points is not None: + vertices_array = np.concatenate([additional_test_points, vertices_array]) + + if not is_convex(vertices_array): + raise GeometryError( + f"{self} (excluding the surface)'s vertices needs to be convex!" + ) + return self.csg.region_from_surface_series( + [self.cw_surface, self.ccw_surface, *self.interior_surfaces], + vertices_array, + control_id=control_id, + ) + # exclusion zone facing towards the plasma + vertices_array = self.exterior_wire.get_3D_coordinates() + if additional_test_points is not None: + vertices_array = np.concatenate([additional_test_points, vertices_array]) + + if not is_convex(vertices_array): + raise GeometryError( + f"{self} (excluding the vacuum vessel)'s vertices needs to be convex!" + ) + return self.csg.region_from_surface_series( + [self.cw_surface, self.ccw_surface, *self.exterior_surfaces], + vertices_array, + control_id=control_id, + ) + + +class DivertorCellStack: + """ + A CONVEX object! i.e. all its exterior points together should make a convex hull. + A stack of DivertorCells (openmc.Cells), first cell is closest to the interior and + last cell is closest to the exterior. They should all be situated on the same + poloidal angle. + """ + + def __init__( + self, divertor_cell_stack: list[DivertorCell], csg: BluemiraNeutronicsCSG + ): + self.cell_stack = divertor_cell_stack + self.csg = csg + # This check below is invalid because of how we subtract region instead. + # for int_cell, ext_cell in pairwise(self.cell_stack): + # if int_cell.exterior_surfaces is not ext_cell.interior_surfaces: + # raise ValueError("Expected a contiguous stack of cells!") + + @property + def interior_surfaces(self): + """Get interior surfaces""" + return self.cell_stack[0].interior_surfaces + + @property + def exterior_surfaces(self): + """Get exterior surfaces""" + return self.cell_stack[-1].exterior_surfaces + + @property + def ccw_surface(self): + """Get counter clockwise surface""" + return self.cell_stack[-1].ccw_surface + + @property + def cw_surface(self): + """Get clockwise surface""" + return self.cell_stack[-1].cw_surface + + @property + def exterior_wire(self): + """Alias to find the outermost cell's exterior_wire""" + return self.cell_stack[-1].exterior_wire + + @property + def interior_wire(self): + """Alias to find the innermost cell's interior_wire""" + return self.cell_stack[0].interior_wire + + @property + def interfaces(self): + """ + All of the radial surfaces, including the innermost (exposed to plasma) and + outermost (facing the vacuum vessel); arranged in that order (from innermost to + outermost). + """ + if not hasattr(self, "_interfaces"): + self._interfaces = [cell.interior_surfaces for cell in self.cell_stack] + self._interfaces.append(self.cell_stack[-1].exterior_surfaces) + return self._interfaces # list of list of (1- or 2-tuple of) surfaces. + + def __len__(self) -> int: + """Length of DivertorCellStack""" + return len(self.cell_stack) + + def __getitem__(self, index_or_slice) -> list[DivertorCell] | DivertorCell: + """Get item for DivertorCellStack""" + return self.cell_stack[index_or_slice] + + def __iter__(self) -> Iterator[DivertorCell]: + """Iterator for DivertorCellStack""" + return iter(self.cell_stack) + + def __repr__(self) -> str: + """String representation""" + return super().__repr__().replace(" at ", f" of {len(self)} DivertorCells at ") + + def get_all_vertices(self) -> npt.NDArray: + """ + Returns + ------- + vertices_array + shape = (N+M, 3) + """ + return np.concatenate([ + self.interior_wire.get_3D_coordinates(), + self.exterior_wire.get_3D_coordinates(), + ]) + + def get_overall_region(self, *, control_id: bool = False) -> openmc.Region: + """ + Get the region that this cell-stack encompasses. + + Parameters + ---------- + control_id + Passed as argument onto + :func:`~bluemira.radiation_transport.neutronics.make_csg.region_from_surface_series` + """ + all_vertices = self.get_all_vertices() + if not is_convex(all_vertices): + raise GeometryError(f"overall_region of {self} needs to be convex!") + + return self.csg.region_from_surface_series( + [ + self.cw_surface, + self.ccw_surface, + *self.interior_surfaces, + *self.exterior_surfaces, + ], + all_vertices, + control_id=control_id, + ) + + @classmethod + def from_divertor_pre_cell( + cls, + divertor_pre_cell: DivertorPreCell, + cw_surface: openmc.Surface, + ccw_surface: openmc.Surface, + materials: MaterialsLibrary, + csg: BluemiraNeutronicsCSG, + armour_thickness: float = 0, + stack_num: str | int = "", + ) -> DivertorCellStack: + """ + Create a stack from a single pre-cell and two poloidal surfaces sandwiching it. + + Parameters + ---------- + stack_num: + A string or number to identify the cell stack by. + """ + # I apologise that this is still hard to read. + # I'm trying to make cell_stack a 3-element list if armour_thickness>0, + # but a 2-element list if armour_thickness==0. + + outermost_wire = divertor_pre_cell.exterior_wire + outermost_surf = csg.surfaces_from_info_list(outermost_wire) + outer_wire = divertor_pre_cell.vv_wire + outer_surf = csg.surfaces_from_info_list(outer_wire) + if armour_thickness > 0: + inner_wire = divertor_pre_cell.offset_interior_wire(armour_thickness) + inner_surf = csg.surfaces_from_info_list(inner_wire) + innermost_wire = divertor_pre_cell.interior_wire + innermost_surf = csg.surfaces_from_info_list(innermost_wire) + else: + inner_wire = divertor_pre_cell.interior_wire + inner_surf = csg.surfaces_from_info_list(inner_wire) + # make the middle cell + cell_stack = [ + # The middle cell is the only cell guaranteed to be convex. + # Therefore it is the first cell to be made. + DivertorCell( + # surfaces: ext, cw, ccw, int. + outer_surf, + cw_surface, + ccw_surface, + inner_surf, + # wires: ext, int. + outer_wire, + inner_wire, + csg=csg, + name=f"Bulk of divertor in diver cell stack {stack_num}", + fill=materials.match_material(CellType.DivertorBulk), + ), + ] + # make the vacuum vessel cell + cell_stack.append( + DivertorCell( + # surfaces: ext, cw, ccw, int. + outermost_surf, + cw_surface, + ccw_surface, + inner_surf, + # wires: ext, int. + outermost_wire, + outer_wire.reverse(), + csg=csg, + subtractive_region=cell_stack[0].exclusion_zone( + away_from_plasma=False, + additional_test_points=innermost_wire.get_3D_coordinates() + if armour_thickness > 0 + else inner_wire.get_3D_coordinates(), + ), + name="Vacuum Vessel behind the divertor in divertor cell stack" + f"{stack_num}", + fill=materials.match_material(CellType.DivertorFirstWall), + ) + ) + # Unfortunately, this does mean that the vacuum vessel has a larger ID than the + # divertor cassette. + if armour_thickness > 0: + # exterior of bulk becomes the interior of surface cell. + cell_stack.insert( + 0, + DivertorCell( + # surfaces: ext, cw, ccw, int. + # Same ext surfaces as before. + # It'll be handled by subtractive_region later. + outer_surf, + cw_surface, + ccw_surface, + innermost_surf, + # wires: ext, int. + inner_wire.reverse(), + innermost_wire, + csg=csg, + name=f"Divertor armour in divertor cell stack {stack_num}", + # subtract away everything in the first cell. + subtractive_region=cell_stack[0].exclusion_zone( + away_from_plasma=True, + additional_test_points=outermost_wire.get_3D_coordinates(), + ), + fill=materials.match_material(CellType.DivertorSurface), + ), + ) + # again, unfortunately, this does mean that the surface armour cell has the + # largest ID. + return cls(cell_stack, csg) + + +class DivertorCellArray: + """Turn the divertor into a cell array""" + + def __init__(self, cell_array: list[DivertorCellStack]): + """Create array from a list of DivertorCellStack.""" + self.cell_array = cell_array + self.poloidal_surfaces = [self.cell_array[0].cw_surface] + self.radial_surfaces = [] + # check neighbouring cells have the same cell stack. + for i, this_stack in enumerate(self.cell_array): + self.poloidal_surfaces.append(this_stack.ccw_surface) + self.radial_surfaces.append(this_stack.interfaces) + + # check neighbouring cells share the same lateral surface + if i != len(self.cell_array) - 1: + next_stack = self.cell_array[i + 1] + if this_stack.ccw_surface is not next_stack.cw_surface: + raise GeometryError( + f"Neighbouring DivertorCellStack {i} and {i + 1} are expected to" + " share the same poloidal wall." + ) + + def __len__(self) -> int: + """Length of DivertorCellArray""" + return len(self.cell_array) + + def __getitem__(self, index_or_slice) -> list[DivertorCellStack] | DivertorCellStack: + """Get item for DivertorCellArray""" + return self.cell_array[index_or_slice] + + def __iter__(self) -> Iterator[DivertorCellStack]: + """Iterator for DivertorCellArray""" + return iter(self.cell_array) + + def __repr__(self) -> str: + """String representation""" + return ( + super().__repr__().replace(" at ", f" of {len(self)} DivertorCellStacks at") + ) + + def interior_surfaces(self) -> list[openmc.Surface]: + """ + Get all of the innermost (plasm-facing) surface. + Runs clockwise. + """ + return list( + chain.from_iterable([surf_stack[0] for surf_stack in self.radial_surfaces]) + ) + + def exterior_surfaces(self) -> list[openmc.Surface]: + """ + Get all of the outermost (vacuum-vessel-facing) surface. + Runs clockwise. + """ + return list( + chain.from_iterable([ + surf_stack[-1][::-1] for surf_stack in self.radial_surfaces + ]) + ) + + def exterior_vertices(self) -> npt.NDArray: + """ + Returns all of the tokamak's poloidal cross-section's outside corners' + coordinates, in 3D. + + Returns + ------- + exterior_vertices: npt.NDArray of shape (N+1, 3) + Arranged counter-clockwise (inboard to outboard). + """ + return np.concatenate([ + stack.exterior_wire.get_3D_coordinates()[::-1] for stack in self.cell_array + ]) + + def interior_vertices(self) -> npt.NDArray: + """ + Returns all of the tokamak's poloidal cross-section's inside corners' + coordinates, in 3D. + + Parameters + ---------- + interior_vertices: npt.NDArray of shape (N+1, 3) + Arranged counter-clockwise (inboard to outboard). + """ + return np.concatenate([ + stack.interior_wire.get_3D_coordinates() for stack in self.cell_array + ]) + + def exclusion_zone(self, *, control_id: bool = False) -> openmc.Region: + """ + Get the exclusion zone AWAY from the plasma. + Usage: plasma_region = openmc.Union(..., ~self.exclusion_zone(), ...) + Assumes every single cell-stack is made of an interior surface which itself forms + a convex hull. + + Parameters + ---------- + control_id + Passed as argument onto + :func:`~bluemira.radiation_transport.neutronics.make_csg.region_from_surface_series` + """ + return openmc.Union([ + stack[0].exclusion_zone( + away_from_plasma=True, + control_id=control_id, + additional_test_points=stack.exterior_wire.get_3D_coordinates(), + ) + for stack in self.cell_array + ]) + + @classmethod + def from_pre_cell_array( + cls, + pre_cell_array: DivertorPreCellArray, + materials: MaterialsLibrary, + divertor_thickness: DivertorThickness, + csg: BluemiraNeutronicsCSG, + override_start_end_surfaces: tuple[openmc.Surface, openmc.Surface] | None = None, + ) -> DivertorCellArray: + """ + Create the entire divertor from the pre-cell array. + + Parameters + ---------- + pre_cell_array + The array of divertor pre-cells. + materials + container of openmc.Material + divertor_thickness + A parameter + :class:`bluemira.radiation_transport.neutronics.params.DivertorThickness`. + For now it only has one scalar value stating how thick the + divertor armour should be. + override_start_end_surfaces + openmc.Surfaces that would be used as the first cw_surface and last + ccw_surface + """ + stack_list = [] + + def get_final_surface() -> openmc.Surface: + """Generate the final surface on-the-fly so that it gets the correct id.""" + if override_start_end_surfaces: + return override_start_end_surfaces[-1] + return csg.surface_from_straight_line( + pre_cell_array[-1].ccw_wall[-1].key_points + ) + + if override_start_end_surfaces: + cw_surf = override_start_end_surfaces[0] + else: + cw_surf = csg.surface_from_straight_line( + pre_cell_array[0].cw_wall[0].key_points + ) + for i, dpc in enumerate(pre_cell_array): + if i == (len(pre_cell_array) - 1): + ccw_surf = get_final_surface() + else: + ccw_surf = csg.surface_from_straight_line(dpc.ccw_wall[-1].key_points) + stack_list.append( + DivertorCellStack.from_divertor_pre_cell( + dpc, + cw_surf, + ccw_surf, + materials, + csg, + divertor_thickness.surface, + i + 1, + ) + ) + cw_surf = ccw_surf + return cls(stack_list) + + def get_hollow_merged_cells(self) -> list[openmc.Cell]: + """ + Returns a list of cells (unnamed, unspecified-ID) where each corresponds to a + cell-stack. + """ + return [ + openmc.Cell(region=stack.get_overall_region()) for stack in self.cell_array + ] diff --git a/bluemira/codes/openmc/material.py b/bluemira/codes/openmc/material.py new file mode 100644 index 0000000000..58c6d28e77 --- /dev/null +++ b/bluemira/codes/openmc/material.py @@ -0,0 +1,113 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""OpenMC CSG neutronics materials""" + +from __future__ import annotations + +from dataclasses import asdict, dataclass, fields +from enum import Enum, auto +from typing import TYPE_CHECKING + +import openmc + +if TYPE_CHECKING: + from pathlib import Path + + from bluemira.radiation_transport.neutronics.materials import NeutronicsMaterials + + +class CellType(Enum): + """ + The five layers of the blanket as used in the neutronics simulation. + + Variables + --------- + Surface + The surface layer of the first wall. + First wall + Typically made of tungsten or Eurofer + BreedingZone + Where tritium is bred + Manifold + The pipe works and supporting structure + VacuumVessel + The vacuum vessel keeping the plasma from mixing with outside air. + RadiationShield + The radiation shield surrounding the reactor, also called a bio shield. + + """ + + BlanketSurface = auto() + BlanketFirstWall = auto() + BlanketBreedingZone = auto() + BlanketManifold = auto() + VacuumVessel = auto() + DivertorBulk = auto() + DivertorFirstWall = auto() + DivertorSurface = auto() + TFCoil = auto() + CSCoil = auto() + RadiationShield = auto() + + +@dataclass +class MaterialsLibrary: + """A dictionary of materials according to the type of blanket used""" + + inb_vv_mat: openmc.Material + inb_fw_mat: openmc.Material + inb_bz_mat: openmc.Material + inb_mani_mat: openmc.Material + divertor_mat: openmc.Material + div_fw_mat: openmc.Material + outb_fw_mat: openmc.Material + outb_bz_mat: openmc.Material + outb_mani_mat: openmc.Material + outb_vv_mat: openmc.Material + tf_coil_mat: openmc.Material + container_mat: openmc.Material + inb_sf_mat: openmc.Material + outb_sf_mat: openmc.Material + div_sf_mat: openmc.Material + rad_shield: openmc.Material + + @classmethod + def from_neutronics_materials(cls, materials_lib: NeutronicsMaterials): + """Initialise from neutronics materials""" + return cls(**{ + field.name: getattr(materials_lib, field.name).to_openmc_material() + for field in fields(materials_lib) + }) + + def match_material(self, cell_type: CellType, *, inboard: bool = False): + """Choose the appropriate blanket material for the given blanket cell type.""" + match cell_type: + case CellType.BlanketSurface: + return self.inb_sf_mat if inboard else self.outb_sf_mat + case CellType.BlanketFirstWall: + return self.inb_fw_mat if inboard else self.outb_fw_mat + case CellType.BlanketBreedingZone: + return self.inb_bz_mat if inboard else self.outb_bz_mat + case CellType.BlanketManifold: + return self.inb_mani_mat if inboard else self.outb_mani_mat + case CellType.VacuumVessel: + return self.inb_vv_mat if inboard else self.outb_vv_mat + case CellType.CSCoil: + return self.container_mat + case CellType.TFCoil: + return self.tf_coil_mat + case CellType.DivertorBulk: + return self.divertor_mat + case CellType.DivertorFirstWall: + return self.div_fw_mat + case CellType.DivertorSurface: + return self.div_sf_mat + case CellType.RadiationShield: + return self.rad_shield + + def export_to_xml(self, path: str | Path = "materials.xml"): + """Exports material defintions to xml""" + return openmc.Materials(asdict(self).values()).export_to_xml(path) diff --git a/bluemira/codes/openmc/output.py b/bluemira/codes/openmc/output.py new file mode 100644 index 0000000000..11e187e1fc --- /dev/null +++ b/bluemira/codes/openmc/output.py @@ -0,0 +1,339 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Functions to present the results prettily +(Including both printed/logged texts and images) +""" + +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import openmc +from tabulate import tabulate + +from bluemira.base.constants import raw_uc +from bluemira.base.look_and_feel import bluemira_debug +from bluemira.radiation_transport.neutronics.constants import DPACoefficients + + +def get_percent_err(row): + """ + Calculate a percentage error to the required row, + assuming cells had been filled out. + + Parameters + ---------- + row: pd.Series object + It should have the "mean" and "std. dev." + row['mean']: float + row['std. dev.']: float + + Returns + ------- + fractional_error: float + + Usage + ----- + dataframe.apply(get_percent_err), + where dataframe must have one row named "std. dev." and another named "mean". + """ + # if percentage error > 1E7:) + if np.isclose(row["mean"], 0.0, rtol=0.0, atol=row["std. dev."] / 100000): + return np.nan + # else: normal mode of operation: divide std by mean, then multiply by 100. + return row["std. dev."] / row["mean"] * 100.0 + + +@dataclass +class OpenMCResult: + """ + Class that looks opens up the openmc universe from the statepoint file, + so that the dataframes containing the relevant results + can be generated and reformatted by its methods. + """ + + tbr: float + tbr_err: float + heating: dict + neutron_wall_load: dict + """Neutron wall load (eV)""" + + photon_heat_flux: dict + """Photon heat flux""" + + universe: openmc.Universe + src_rate: float + statepoint: openmc.StatePoint + statepoint_file: str + cell_names: dict + cell_vols: dict # [m^3] + mat_names: dict + + @classmethod + def from_run( + cls, + universe: openmc.Universe, + src_rate: float, + statepoint_file: str = "", + ): + """Create results class from run statepoint""" + # Create cell and material name dictionaries to allow easy mapping to dataframe + cell_names = {} + mat_names = {} + for cell_id, _cell in universe.cells.items(): + cell_names[cell_id] = _cell.name + if _cell.fill is not None: # if this cell is not a void + mat_names[_cell.fill.id] = _cell.fill.name + + # Creating cell volume dictionary to allow easy mapping to dataframe + # provided by openmc in cm^3, but we want to save it in m^3 + cell_vols = { + cell_id: raw_uc(universe.cells[cell_id].volume, "cm^3", "m^3") + if isinstance(universe.cells[cell_id].volume, float) + else None + for cell_id in universe.cells + } + # Loads up the output file from the simulation + statepoint = openmc.StatePoint(statepoint_file) + tbr, tbr_err = cls._load_tbr(statepoint) + + return cls( + universe=universe, + src_rate=src_rate, + statepoint_file=statepoint_file, + statepoint=statepoint, + cell_names=cell_names, + cell_vols=cell_vols, + mat_names=mat_names, + tbr=tbr, + tbr_err=tbr_err, + heating=cls._load_heating(statepoint, mat_names, src_rate), + neutron_wall_load=cls._load_neutron_wall_loading( + statepoint, cell_names, cell_vols, src_rate + ), + photon_heat_flux=cls._load_photon_heat_flux( + statepoint, cell_names, cell_vols, src_rate + ), + ) + + @staticmethod + def _load_volume_calculation_from_file( + volume_file_path: Path, cell_names: list[str] + ): + """ + Load the volume file to record as volume information. + + Parameters + ---------- + volume_file_path + + Cell_names + indicative names to print. + """ + if volume_file_path.is_file(): + vol_results = openmc.VolumeCalculation.from_hdf5(volume_file_path) + vols = vol_results.volumes + ids = list(vols.keys()) + cell_volumes = { + "cell": ids, + "cell_names": [cell_names[i] for i in ids], + "Stochastic Volumes": list(raw_uc(list(vols.values()), "cm^3", "m^3")), + } + + else: + bluemira_debug("No volume file found") + vol_results = None + cell_volumes = None + + return vol_results, cell_volumes + + @staticmethod + def _load_dataframe_from_statepoint(statepoint, tally_name: str): + return statepoint.get_tally(name=tally_name).get_pandas_dataframe() + + @staticmethod + def _convert_dict_contents(dataset: dict[str, dict[int, list[str | float]]]): + for k, v in dataset.items(): + vals = list(v.values()) if isinstance(v, dict) else v + dataset[k] = vals if isinstance(vals[0], str) else np.array(vals) + return dataset + + @classmethod + def _load_tbr(cls, statepoint): + """Load the TBR value and uncertainty.""" + tbr_df = cls._load_dataframe_from_statepoint(statepoint, "TBR") + return tbr_df["mean"].sum(), tbr_df["std. dev."].sum() + + @classmethod + def _load_heating(cls, statepoint, mat_names, src_rate): + """Load the heating (sorted by material) dataframe""" + # mean and std. dev. are given in eV per source particle, + # so we don't need to show them to the user. + heating_df = cls._load_dataframe_from_statepoint(statepoint, "Total power") + heating_df["material_name"] = heating_df["material"].map(mat_names) + heating_df["mean(W)"] = raw_uc( + heating_df["mean"].to_numpy() * src_rate, "eV/s", "W" + ) + heating_df["err."] = raw_uc( + heating_df["std. dev."].to_numpy() * src_rate, "eV/s", "W" + ) + heating_df["%err."] = heating_df.apply(get_percent_err, axis=1) + # rearrange dataframe into this desired order + heating_df = heating_df[ + [ + "material", + "material_name", + "nuclide", + "score", + "mean(W)", + "err.", + "%err.", + ] + ] + hdf = heating_df.to_dict() + return cls._convert_dict_contents(hdf) + + @classmethod + def _load_neutron_wall_loading(cls, statepoint, cell_names, cell_vols, src_rate): + """Load the neutron wall load dataframe""" + dfa_coefs = DPACoefficients() # default assumes iron (Fe) is used. + n_wl_df = cls._load_dataframe_from_statepoint( + statepoint, "neutron flux in every cell" + ) + n_wl_df["cell_name"] = n_wl_df["cell"].map(cell_names) + n_wl_df["vol (m^3)"] = n_wl_df["cell"].map(cell_vols) + total_displacements_per_second = ( + n_wl_df["mean"] * dfa_coefs.displacements_per_damage_eV * src_rate + ) # "mean" has units "eV per source particle" + # total number of atomic displacements per second in the cell. + num_atoms_in_cell = n_wl_df["vol (m^3)"] * raw_uc( + dfa_coefs.atoms_per_cc, "1/cm^3", "1/m^3" + ) + n_wl_df["dpa/fpy"] = raw_uc( + total_displacements_per_second.to_numpy() / num_atoms_in_cell.to_numpy(), + "1/s", + "1/year", + ) + + n_wl_df["%err."] = n_wl_df.apply(get_percent_err, axis=1) + # keep only the surface cells: + n_wl_df = n_wl_df.drop( + n_wl_df[~n_wl_df["cell_name"].str.contains("Surface")].index + ) + # DataFrame columns rearrangement + n_wl_df = n_wl_df[ + [ + "cell", + "cell_name", + "particle", + "nuclide", + "score", + "vol (m^3)", + "dpa/fpy", + "%err.", + ] + ] + + return cls._convert_dict_contents(n_wl_df.to_dict()) + + @classmethod + def _load_photon_heat_flux(cls, statepoint, cell_names, cell_vols, src_rate): + """Load the photon heaat flux dataframe""" + p_hf_df = cls._load_dataframe_from_statepoint(statepoint, "photon heating") + p_hf_df["cell_name"] = p_hf_df["cell"].map(cell_names) + + p_hf_df["%err."] = p_hf_df.apply(get_percent_err, axis=1) + p_hf_df["vol (m^3)"] = p_hf_df["cell"].map(cell_vols) + p_hf_df["heating (W)"] = photon_heating = raw_uc( + p_hf_df["mean"].to_numpy() * src_rate, "eV/s", "W" + ) + p_hf_df["heating std.dev."] = photon_heating_stddev = raw_uc( + p_hf_df["std. dev."].to_numpy() * src_rate, "eV/s", "W" + ) + p_hf_df["vol. heating (W/m3)"] = photon_heating / p_hf_df["vol (m^3)"] + p_hf_df["vol. heating std.dev."] = photon_heating_stddev / p_hf_df["vol (m^3)"] + + # Scaling first wall results by factor to surface results + surface_total = p_hf_df.loc[ + p_hf_df["cell_name"].str.contains("Surface"), "heating (W)" + ].sum() + cell_total = p_hf_df.loc[ + ~p_hf_df["cell_name"].str.contains("Surface"), "heating (W)" + ].sum() + + surface_factor = surface_total / cell_total + # in-place modification + p_hf_df["vol. heating (W/m3)"] = np.where( + ~p_hf_df["cell_name"].str.contains( + "Surface" + ), # modify the matching entries, + p_hf_df["heating (W)"] * surface_factor, + p_hf_df["heating (W)"], # otherwise leave it unchanged. + ) + p_hf_df["vol. heating (W/m3)"] = np.where( + ~p_hf_df["cell_name"].str.contains("Surface"), + p_hf_df["heating std.dev."] * surface_factor, + p_hf_df["heating std.dev."], + ) + p_hf_df["vol. heating (W/m3)"] = np.where( + ~p_hf_df["cell_name"].str.contains("Surface"), + p_hf_df["vol. heating (W/m3)"] * surface_factor, + p_hf_df["vol. heating (W/m3)"], + ) + p_hf_df["vol. heating (W/m3)"] = np.where( + ~p_hf_df["cell_name"].str.contains("Surface"), + p_hf_df["vol. heating std.dev."] * surface_factor, + p_hf_df["vol. heating std.dev."], + ) + # DataFrame columns rearrangement + p_hf_df = p_hf_df[ + [ + "cell", + "cell_name", + "particle", + "nuclide", + "score", + "vol (m^3)", + "heating (W)", + "heating std.dev.", + "vol. heating (W/m3)", + "vol. heating std.dev.", + "%err.", + ] + ] + p_hf_dict = p_hf_df.to_dict() + return cls._convert_dict_contents(p_hf_dict) + + def __str__(self): + """String representation""" + ret_str = f"TBR\n{self.tbr:.3f}±{self.tbr_err:.3f}" + for title, data in zip( + ("Heating (W)", "Neutron Wall Load (eV)", "Photon Heat Flux (W m)"), + ( + self.heating, + self.neutron_wall_load, + self.photon_heat_flux, + ), + strict=True, + ): + ret_str += f"\n{title}\n{self._tabulate(data)}" + + return ret_str + + @staticmethod + def _tabulate( + records: dict[str, str | float], + tablefmt: str = "fancy_grid", + floatfmt: str = ".3g", + ) -> str: + return tabulate( + zip(*records.values(), strict=False), + headers=records.keys(), + tablefmt=tablefmt, + showindex=False, + numalign="right", + floatfmt=floatfmt, + ) diff --git a/bluemira/codes/openmc/params.py b/bluemira/codes/openmc/params.py new file mode 100644 index 0000000000..75176c890d --- /dev/null +++ b/bluemira/codes/openmc/params.py @@ -0,0 +1,165 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""dataclasses containing parameters used to set up the openmc model.""" + +from __future__ import annotations + +from dataclasses import dataclass, fields + +from bluemira.base.constants import raw_uc +from bluemira.base.look_and_feel import bluemira_warn +from bluemira.base.parameter_frame import Parameter, ParameterFrame +from bluemira.geometry.error import GeometryError + + +@dataclass +class OpenMCNeutronicsSolverParams(ParameterFrame): + """ + + Parameters + ---------- + major_radius: + Major radius of the machine + aspect_ratio: + aspect ratio of the machine + elongation: + elongation of the plasma + triangularity: + triangularity of the plasma + reactor_power: + total reactor (thermal) power when operating at 100% + peaking_factor: + (max. heat flux on fw)/(avg. heat flux on fw) + temperature: + plasma temperature (assumed to be uniform throughout the plasma) + shaf_shift: + Shafranov shift + shift of the centre of flux surfaces, i.e. + mean(min radius, max radius) of the LCFS, + towards the outboard radial direction. + vertical_shift: + how far (upwards) in the z direction is the centre of the plasma + shifted compared to the geometric center of the poloidal cross-section. + """ + + R_0: Parameter[float] + """Major Radius""" + A: Parameter[float] + """Aspect ratio""" + kappa: Parameter[float] + """Plasma elongation""" + delta: Parameter[float] + """Plasma triangularity""" + reactor_power: Parameter[float] # [W] + peaking_factor: Parameter[float] # [dimensionless] + T_e: Parameter[float] + """Average plasma electron temperature [J]""" + shaf_shift: Parameter[float] + """Shafranov shift""" + vertical_shift: Parameter[float] # [m] + + +@dataclass(frozen=True) +class PlasmaSourceParameters: + """ + Parameters describing the plasma source, + i.e. where the plasma is positioned (and therefore where the power is concentrated), + and what temperature the plasma is at. + + Parameters + ---------- + reactor_power: + total reactor (thermal) power when operating at 100% + peaking_factor: + (max. heat flux on fw)/(avg. heat flux on fw) + temperature: + plasma temperature (assumed to be uniform throughout the plasma) + shaf_shift: + Shafranov shift + shift of the centre of flux surfaces, i.e. + mean(min radius, max radius) of the LCFS, + towards the outboard radial direction. + vertical_shift: + how far (upwards) in the z direction is the centre of the plasma + shifted compared to the geometric center of the poloidal cross-section. + plasma_physics_units: + Plasma_physics_units converted variables + """ + + major_radius: float # [m] + aspect_ratio: float # [dimensionless] + elongation: float # [dimensionless] + triangularity: float # [dimensionless] + reactor_power: float # [W] + peaking_factor: float # [dimensionless] + temperature: float # [K] + shaf_shift: float # [m] + vertical_shift: float # [m] + plasma_physics_units: PlasmaSourceParameters | None = None + + def __post_init__(self): + """Check dimensionless variables are sensible.""" + if self.peaking_factor < 1.0: + raise ValueError( + "Peaking factor (peak heat load/avg. heat load) " + "must be larger than 1, by definition." + ) + if self.aspect_ratio < 1.0: + raise GeometryError( + "By construction, tokamak aspect ratio " "can't be smaller than 1." + ) + if self.elongation < 1.0: + raise GeometryError("Elongation can't be smaller than 1") + if abs(self.triangularity) > 1.0: + # triangularity <0 is known as reversed/ negative triangularity. + bluemira_warn( + "Triangularity with magnitude >1 implies that the difference" + "between the major radius and R_upper is larger than the minor radius." + ) + + @property + def minor_radius(self): + """Calculate minor radius from + aspect_ratio = major_radius/minor_radius + """ + return self.major_radius / self.aspect_ratio + + @classmethod + def from_parameterframe(cls, params: ParameterFrame): + """ + Convert from si units dataclass + :class:`~bluemira.radiation_transport.neutronics.params.PlasmaSourceParameters` + + This gives the illusion that self.cgs.x = scale_factor*self.x + We rely on the 'frozen' nature of this dataclass so these links don't break. + """ + conversion = { + "major_radius": ("m", "cm"), + "reactor_power": ("W", "MW"), + "temperature": ("J", "keV"), + "shaf_shift": ("m", "cm"), + "vertical_shift": ("m", "cm"), + } + mapping = { + "aspect_ratio": "A", + "major_radius": "R_0", + "elongation": "kappa", + "triangularity": "delta", + "temperature": "T_e", + } + param_convert_dict = {} + param_dict = {} + for k in fields(cls): + if k.name == "plasma_physics_units": + continue + val = getattr(params, mapping.get(k.name, k.name)).value + param_dict[k.name] = val + if k.name in conversion: + param_convert_dict[k.name] = raw_uc(val, *conversion[k.name]) + else: + param_convert_dict[k.name] = val + + return cls(**param_dict, plasma_physics_units=cls(**param_convert_dict)) diff --git a/bluemira/codes/openmc/solver.py b/bluemira/codes/openmc/solver.py new file mode 100644 index 0000000000..f91b1e7abc --- /dev/null +++ b/bluemira/codes/openmc/solver.py @@ -0,0 +1,519 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""OpenMC designer""" + +from __future__ import annotations + +from collections.abc import Callable +from contextlib import contextmanager +from dataclasses import dataclass, fields +from enum import auto +from operator import attrgetter +from pathlib import Path +from typing import Literal + +import numpy as np +import openmc + +from bluemira.base.constants import raw_uc +from bluemira.base.look_and_feel import bluemira_debug +from bluemira.base.parameter_frame import ParameterFrame, make_parameter_frame +from bluemira.base.tools import _timing +from bluemira.codes.interface import ( + BaseRunMode, + CodesSetup, + CodesSolver, + CodesTask, + CodesTeardown, +) +from bluemira.codes.openmc.make_csg import ( + BlanketCellArray, + BluemiraNeutronicsCSG, + DivertorCellArray, + make_cell_arrays, +) +from bluemira.codes.openmc.material import MaterialsLibrary +from bluemira.codes.openmc.output import OpenMCResult +from bluemira.codes.openmc.params import ( + OpenMCNeutronicsSolverParams, + PlasmaSourceParameters, +) +from bluemira.codes.openmc.tallying import filter_cells +from bluemira.plasma_physics.reactions import n_DT_reactions + + +class OpenMCRunModes(BaseRunMode): + """OpenMC run modes""" + + RUN = openmc.settings.RunMode.FIXED_SOURCE.value + RUN_AND_PLOT = auto() + PLOT = openmc.settings.RunMode.PLOT.value + VOLUME = openmc.settings.RunMode.VOLUME.value + + +OPENMC_NAME = "OpenMC" + + +@dataclass +class OpenMCSimulationRuntimeParameters: + """Parameters used in the actual simulation + + Parameters + ---------- + particles: + Number of neutrons emitted by the plasma source per batch. + batches: + How many batches to simulate. + photon_transport: + Whether to simulate the transport of photons (i.e. gamma-rays created) or not. + electron_treatment: + The way in which OpenMC handles secondary charged particles. + 'thick-target bremsstrahlung' or 'local energy deposition' + 'thick-target bremsstrahlung' accounts for the energy carried away by + bremsstrahlung photons and deposited elsewhere, whereas 'local energy + deposition' assumes electrons deposit all energies locally. + (the latter is expected to be computationally faster.) + run_mode: + see below for details: + https://docs.openmc.org/en/stable/usersguide/settings.html#run-modes + openmc_write_summary: + whether openmc should write a 'summary.h5' file or not. + cross_section_xml: + Where the xml file for cross-section is stored locally. + """ + + particles: int # number of particles used in the neutronics simulation + cross_section_xml: str | Path + batches: int = 2 + photon_transport: bool = True + # Bremsstrahlung only matters for very thin objects + electron_treatment: Literal["ttb", "led"] = "led" + run_mode: str = OpenMCRunModes.RUN.value + openmc_write_summary: bool = False + parametric_source: bool = True + plot_axis: str = "xz" + plot_pixel_per_metre: int = 100 + + +class Setup(CodesSetup): + """Setup task for OpenMC solver""" + + def __init__( + self, + out_path: str, + codes_name: str, + cross_section_xml: str, + source, + cell_arrays, + pre_cell_model, + materials, + ): + super().__init__(None, codes_name) + + self.out_path = out_path + self.cells = cell_arrays.cells + self.cross_section_xml = cross_section_xml + self.source = source + self.blanket_cell_array = cell_arrays.blanket + self.divertor_cell_array = cell_arrays.divertor + self.pre_cell_model = pre_cell_model + self.materials = materials + self.matlist = attrgetter( + "outb_sf_mat", + "outb_fw_mat", + "outb_bz_mat", + "outb_mani_mat", + "outb_vv_mat", + "divertor_mat", + "div_fw_mat", + "tf_coil_mat", + ) + + @contextmanager + def _base_setup(self, run_mode, *, debug: bool = False): + from openmc.config import config # noqa: PLC0415 + + self.files_created = set() + folder = run_mode.name.lower() + cwd = Path(self.out_path, folder) + cwd.mkdir(parents=True, exist_ok=True) + + config["cross_sections"] = self.cross_section_xml + + self.settings = openmc.Settings( + run_mode=run_mode.value, output={"summary": False} + ) + self.universe = openmc.Universe(cells=self.cells) + self.geometry = openmc.Geometry(self.universe) + self.settings.verbosity = 10 if debug else 6 + try: + yield + finally: + for obj, pth in ( + (self.settings, Path(self.out_path, folder, "settings.xml")), + (self.geometry, Path(self.out_path, folder, "geometry.xml")), + (self.materials, Path(self.out_path, folder, "materials.xml")), + ): + obj.export_to_xml(pth) + self.files_created.add(pth) + + def _set_tallies( + self, + run_mode, + tally_function: TALLY_FUNCTION_TYPE, + blanket_cell_array: BlanketCellArray, + divertor_cell_array: DivertorCellArray, + material_list: list[openmc.Material], + ): + out_path = Path(self.out_path, run_mode.name.lower(), "tallies.xml") + tallies_list = [] + for name, scores, filters in tally_function( + material_list, blanket_cell_array, divertor_cell_array + ): + tally = openmc.Tally(name=name) + tally.scores = [scores] + tally.filters = filters + tallies_list.append(tally) + + openmc.Tallies(tallies_list).export_to_xml(out_path) + self.files_created.add(out_path) + + def run( + self, + run_mode, + runtime_params, + source_params, + tally_function, + *, + debug: bool = False, + ): + """Run stage for setup openmc""" + with self._base_setup(run_mode, debug=debug): + self.settings.particles = runtime_params.particles + self.settings.source = self.source(source_params) + self.settings.batches = int(runtime_params.batches) + self.settings.photon_transport = runtime_params.photon_transport + self.settings.electron_treatment = runtime_params.electron_treatment + + self._set_tallies( + run_mode, + tally_function, + self.blanket_cell_array, + self.divertor_cell_array, + self.matlist(self.materials), + ) + self.files_created.add(f"statepoint.{runtime_params.batches}.h5") + self.files_created.add("tallies.out") + + def plot( + self, + run_mode, + runtime_params, + _source_params, + _tally_function, + *, + debug: bool = False, + ): + """Plot stage for setup openmc""" + with self._base_setup(run_mode, debug=debug): + z_max, _z_min, r_max, _r_min = self.pre_cell_model.bounding_box + plot_width_0 = r_max * 2.1 + plot_width_1 = z_max * 3.1 + plot = openmc.Plot() + plot.basis = runtime_params.plot_axis + plot.pixels = [ + int(plot_width_0 * runtime_params.plot_pixel_per_metre), + int(plot_width_1 * runtime_params.plot_pixel_per_metre), + ] + plot.width = raw_uc([plot_width_0, plot_width_1], "m", "cm") + plot.show_overlaps = True + + plot_pth = Path(self.out_path, run_mode.name.lower(), "plots.xml") + openmc.Plots([plot]).export_to_xml(plot_pth) + self.files_created.add(plot_pth) + + def volume( + self, + run_mode, + runtime_params, + _source_params, + _tally_function, + *, + debug: bool = False, + ): + """Stochastic volume stage for setup openmc""" + z_max, z_min, r_max, r_min = self.pre_cell_model.bounding_box + + min_xyz = (r_min, r_min, z_min) + max_xyz = (r_max, r_max, z_max) + + with self._base_setup(run_mode, debug=debug): + self.settings.volume_calculations = openmc.VolumeCalculation( + self.cells, + runtime_params.particles, + raw_uc(min_xyz, "m", "cm"), + raw_uc(max_xyz, "m", "cm"), + ) + self.files_created.add(Path(self.out_path, run_mode.name.lower(), "summary.h5")) + # single batch + self.files_created.add( + Path(self.out_path, run_mode.name.lower(), "statepoint.1.h5") + ) + + +class Run(CodesTask): + """Run task for OpenMC solver""" + + def __init__(self, out_path: Path, codes_name: str): + super().__init__(None, codes_name) + + self.out_path = out_path + + def _run(self, run_mode, *, debug: bool = False): + """Run openmc""" + folder = run_mode.name.lower() + cwd = Path(self.out_path, folder) + cwd.mkdir(parents=True, exist_ok=True) + _timing( + openmc.run, + "Executed in", + f"Running OpenMC in {folder} mode", + debug_info_str=False, + )( + output=debug, + threads=None, + geometry_debug=False, + restart_file=None, + tracks=False, + cwd=cwd, + openmc_exec="openmc", + mpi_args=None, + event_based=False, + path_input=None, + ) + + def run(self, run_mode, *, debug: bool = False): + """Run stage for run task""" + self._run(run_mode, debug=debug) + + def plot(self, run_mode, *, debug: bool = False): + """Plot stage for run task""" + self._run(run_mode, debug=debug) + + def volume(self, run_mode, *, debug: bool = False): + """Stochastic volume stage for run task""" + self._run(run_mode, debug=debug) + + +class Teardown(CodesTeardown): + """Teardown task for OpenMC solver""" + + def __init__( + self, + cells, + out_path: str, + codes_name: str, + ): + super().__init__(None, codes_name) + + self.out_path = out_path + self.cells = cells + + @staticmethod + def delete_files(files_created): + """Remove files generated during the run (mainly .xml files.)""" + removed_files, failed_to_remove_files = [], [] + for file_name in files_created: + if (f := file_name).exists(): + f.unlink() + removed_files.append(file_name) + else: + failed_to_remove_files.append(file_name) + + if removed_files: + bluemira_debug(f"Removed files {removed_files}") + if failed_to_remove_files: + bluemira_debug( + f"Attempted to remove files {failed_to_remove_files} but " + "they don't exists." + ) + + def run( + self, + universe, + files_created, + source_params, + statepoint_file, + *, + delete_files: bool = False, + ): + """Run stage for Teardown task""" + result = OpenMCResult.from_run( + universe, + n_DT_reactions(source_params.plasma_physics_units.reactor_power), + statepoint_file, + ) + if delete_files: + self.delete_files(files_created) + return result + + def plot( + self, + _universe, + files_created, + *_args, + delete_files: bool = False, + ): + """Plot stage for Teardown task""" + if delete_files: + self.delete_files(files_created) + + def volume( + self, + _universe, + files_created, + _source_params, + _statepoint_file, + *, + delete_files: bool = False, + ) -> dict[int, float]: + """Stochastic volume stage for teardown task""" + if delete_files: + self.delete_files(files_created) + return { + cell.id: raw_uc( + np.nan if cell.volume is None else cell.volume, "cm^3", "m^3" + ) + for cell in self.cells + } + + +TALLY_FUNCTION_TYPE = Callable[ + [list[openmc.Material], BlanketCellArray, DivertorCellArray], + tuple[ + str, + str, + list[openmc.CellFilter | openmc.MaterialFilter | openmc.ParticleFilter], + ], +] + + +class OpenMCNeutronicsSolver(CodesSolver): + """OpenMC 2D neutronics solver""" + + name: str = OPENMC_NAME + param_cls: type[OpenMCNeutronicsSolverParams] = OpenMCNeutronicsSolverParams + params: OpenMCNeutronicsSolverParams + run_mode_cls: type[OpenMCRunModes] = OpenMCRunModes + setup_cls: type[Setup] = Setup + run_cls: type[Run] = Run + teardown_cls: type[Teardown] = Teardown + + def __init__( + self, + params: dict | ParameterFrame, + build_config: dict, + neutronics_pre_cell_model, + source: Callable[[PlasmaSourceParameters], openmc.source.SourceBase], + tally_function: TALLY_FUNCTION_TYPE | None = None, + ): + self.params = make_parameter_frame(params, self.param_cls) + self.build_config = build_config + + self.out_path = self.build_config.get("neutronics_output_path", Path.cwd()) + + self.source = source + + self.pre_cell_model = neutronics_pre_cell_model + self.materials = MaterialsLibrary.from_neutronics_materials( + self.pre_cell_model.material_library + ) + + self.cell_arrays = make_cell_arrays( + self.pre_cell_model, BluemiraNeutronicsCSG(), self.materials, control_id=True + ) + + self.tally_function = filter_cells if tally_function is None else tally_function + + @property + def source(self) -> Callable[[PlasmaSourceParameters], openmc.Source]: + """Source term for OpenMC""" + return self._source + + @source.setter + def source(self, value: Callable[[PlasmaSourceParameters], openmc.Source]): + self._source = value + + @property + def tally_function(self) -> TALLY_FUNCTION_TYPE: + """Function used to set up tallies""" + return self._tally_function + + @tally_function.setter + def tally_function(self, value: TALLY_FUNCTION_TYPE): + self._tally_function = value + + def execute(self, *, debug=False) -> OpenMCResult | dict[int, float]: + """Execute the setup, run, and teardown tasks, in order.""" + run_mode = self.build_config.get("run_mode", self.run_mode_cls.RUN) + if isinstance(run_mode, str): + run_mode = self.run_mode_cls.from_string(run_mode) + + source_params = PlasmaSourceParameters.from_parameterframe(self.params) + runtime_params = OpenMCSimulationRuntimeParameters(**{ + k.name: self.build_config[k.name] + for k in fields(OpenMCSimulationRuntimeParameters) + if k.name in self.build_config + }) + if run_mode is OpenMCRunModes.RUN_AND_PLOT: + for run_mode in (OpenMCRunModes.PLOT, OpenMCRunModes.RUN): + result = self._single_run( + run_mode, source_params, runtime_params, debug=debug + ) + return result + return self._single_run(run_mode, source_params, runtime_params, debug=debug) + + def _single_run( + self, + run_mode: OpenMCRunModes, + source_params: PlasmaSourceParameters, + runtime_params: OpenMCSimulationRuntimeParameters, + *, + debug=False, + ) -> OpenMCResult | dict[int, float]: + self._setup = self.setup_cls( + self.out_path, + self.name, + str(self.build_config["cross_section_xml"]), + self.source, + self.cell_arrays, + self.pre_cell_model, + self.materials, + ) + self._run = self.run_cls(self.out_path, self.name) + self._teardown = self.teardown_cls( + self.cell_arrays.cells, self.out_path, self.name + ) + + result = None + if setup := self._get_execution_method(self._setup, run_mode): + result = setup( + run_mode, runtime_params, source_params, self.tally_function, debug=debug + ) + if run := self._get_execution_method(self._run, run_mode): + result = run(run_mode, debug=debug) + if teardown := self._get_execution_method(self._teardown, run_mode): + result = teardown( + self._setup.universe, + self._setup.files_created, + source_params, + Path( + self.out_path, + run_mode.name.lower(), + f"statepoint.{runtime_params.batches}.h5", + ), + ) + return result diff --git a/bluemira/codes/openmc/sources.py b/bluemira/codes/openmc/sources.py new file mode 100644 index 0000000000..17d72dd457 --- /dev/null +++ b/bluemira/codes/openmc/sources.py @@ -0,0 +1,77 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Neutronics sources""" + +import numpy as np +import openmc + +from bluemira.base.constants import raw_uc +from bluemira.codes.openmc.params import PlasmaSourceParameters +from bluemira.radiation_transport.error import SourceError +from bluemira.radiation_transport.neutronics.constants import dt_neutron_energy + +try: + from pps_isotropic.source import create_parametric_plasma_source + + PPS_ISO_INSTALLED = True +except ImportError: + PPS_ISO_INSTALLED = False + + +def make_pps_source(source_parameters: PlasmaSourceParameters) -> openmc.Source: + """Make a plasma source""" + if not PPS_ISO_INSTALLED: + raise SourceError("pps_isotropic installation not found") + return create_parametric_plasma_source( + # tokamak geometry + major_r=source_parameters.plasma_physics_units.major_radius, + minor_r=source_parameters.plasma_physics_units.minor_radius, + elongation=source_parameters.plasma_physics_units.elongation, + triangularity=source_parameters.plasma_physics_units.triangularity, + # plasma geometry + peaking_factor=source_parameters.plasma_physics_units.peaking_factor, + temperature=source_parameters.plasma_physics_units.temperature, + radial_shift=source_parameters.plasma_physics_units.shaf_shift, + vertical_shift=source_parameters.plasma_physics_units.vertical_shift, + # plasma type + mode="DT", + ) + + +def make_ring_source(source_parameters: PlasmaSourceParameters) -> openmc.Source: + """Create the ring source""" + return create_ring_source( + source_parameters.major_radius, source_parameters.shaf_shift + ) + + +def create_ring_source(major_r_cm: float, shaf_shift_cm: float) -> openmc.Source: + """ + Creating simple line ring source lying on the Z=0 plane, + at r = major radius + shafranov shift, + producing 14.1 MeV neutrons with no variation in energy. + A more accurate source will slightly affect the wall loadings and dpa profiles. + + Parameters + ---------- + major_r_cm: + major radius [cm] + shaf_shift_cm: + shafranov shift [cm] + """ + ring_source = openmc.Source() + source_radii_cm = openmc.stats.Discrete([major_r_cm + shaf_shift_cm], [1]) + source_z_values = openmc.stats.Discrete([0], [1]) + source_angles = openmc.stats.Uniform(a=0.0, b=2 * np.pi) + ring_source.space = openmc.stats.CylindricalIndependent( + r=source_radii_cm, phi=source_angles, z=source_z_values, origin=(0.0, 0.0, 0.0) + ) + ring_source.angle = openmc.stats.Isotropic() + ring_source.energy = openmc.stats.Discrete( + [raw_uc(dt_neutron_energy, "J", "eV")], [1] + ) + + return ring_source diff --git a/bluemira/codes/openmc/tallying.py b/bluemira/codes/openmc/tallying.py new file mode 100644 index 0000000000..b58ec7131a --- /dev/null +++ b/bluemira/codes/openmc/tallying.py @@ -0,0 +1,90 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Functions for creating the openmc tallies.""" + +from itertools import chain + +import openmc + +from bluemira.codes.openmc.make_csg import BlanketCellArray, DivertorCellArray + + +def filter_cells( + material_list, + blanket_cell_array: BlanketCellArray, + divertor_cell_array: DivertorCellArray, +): + """ + Create scores and the filter for the scores. Give them names. + + Returns + ------- + TBR + Achieved by (n,Xt) reaction, which counts the number of tritium-producing + nuclear reactions per neutron emitted at the source. + + We used the (n,Xt) score because the Lithium produces a maximum of 1 Tritium per + reaction, so there won't be any concerns about uncer-counting the TBR. + + Powers + Measures the nuclear heating in various locations and materials, and interpret + this as power. "damage-energy" is given by eV per source neutron. + Multiply by neutron source rate, and then divide by (number of atoms and + threshold displacement energy) to get the DPA. + + Fluence + Measures # of neutrons streaming through. + "flux" is given in # per source particle, so multiply by # of source neutrons to + get the total fluence over the simulation. + Divide by area to get fluence in unit: cm^-2. + + """ + blanket_cells = [*chain.from_iterable(blanket_cell_array)] + div_cells = [*chain.from_iterable(divertor_cell_array)] + cells = blanket_cells + div_cells + fw_surf_cells = [ + *(stack[0] for stack in blanket_cell_array), + *(stack[1] for stack in blanket_cell_array), + ] + vv_cells = [ + *(stack[-1] for stack in blanket_cell_array), + *(stack[-1] for stack in divertor_cell_array), + ] + bz_cells = [stack[2] for stack in blanket_cell_array] + + # Cell filters + # blanket_cell_filter = openmc.CellFilter(blanket_cells) + div_cell_filter = openmc.CellFilter(div_cells) + cell_filter = openmc.CellFilter(cells) + fw_surf_filter = openmc.CellFilter(fw_surf_cells) + vv_filter = openmc.CellFilter(vv_cells) + bz_filter = openmc.CellFilter(bz_cells) + + # material filters + mat_filter = openmc.MaterialFilter(material_list[:-1]) + eurofer_filter = openmc.MaterialFilter([material_list[-1]]) + neutron_filter = openmc.ParticleFilter(["neutron"]) + photon_filter = openmc.ParticleFilter(["photon"]) + + # name, scores, filters + return ( + ("TBR", "(n,Xt)", []), # theoretical maximum TBR only, obviously. + # Powers + ("Total power", "heating", [mat_filter]), + ("divertor power", "heating", [div_cell_filter]), + ("vacuum vessel power", "heating", [vv_filter]), + ("breeding blanket power", "heating", [bz_filter]), + # Fluence + ("neutron flux in every cell", "flux", [cell_filter, neutron_filter]), + ("photon heating", "heating", [fw_surf_filter, photon_filter]), + # ("neutron flux in 2d mesh", "flux", [cyl_mesh_filter, neutron_filter]), + # TF winding pack does not exits yet, so this will have to wait + # DPA + ("eurofer damage", "damage-energy", [cell_filter, eurofer_filter]), + # used to get the EUROFER OBMP + ("divertor damage", "damage-energy", [div_cell_filter, mat_filter]), + ("vacuum vessel damage", "damage-energy", [vv_filter]), + ) diff --git a/bluemira/codes/wrapper.py b/bluemira/codes/wrapper.py index 1c9ca37ae1..72d3150205 100644 --- a/bluemira/codes/wrapper.py +++ b/bluemira/codes/wrapper.py @@ -105,3 +105,31 @@ def transport_code_solver( """ transp = get_code_interface(module) return transp.Solver(params, build_config) + + +def neutronics_code_solver( + params: ParameterFrame, + build_config: BuildConfig, + neutronics_model, + source, + tally_function=None, + module: str = "OPENMC", +) -> CodesSolver: + """ + Neutronics solver + + Parameters + ---------- + params: + ParameterFrame for neutronics code + build_config: + build configuration dictionary + module: + Module to use + + Returns + ------- + The solver object to be run + """ + neutron = get_code_interface(module) + return neutron.Solver(params, build_config, neutronics_model, source, tally_function) diff --git a/bluemira/equilibria/physics.py b/bluemira/equilibria/physics.py index cbb44b34d8..75fd7014d9 100644 --- a/bluemira/equilibria/physics.py +++ b/bluemira/equilibria/physics.py @@ -217,7 +217,7 @@ def calc_energy(eq: Equilibrium) -> float: return volume_integral(Bp**2 * mask, eq.x, eq.dx, eq.dz) / (2 * MU_0) -def calc_Li(eq: Equilibrium) -> float: # noqa: N802 +def calc_Li(eq: Equilibrium) -> float: """ Calculates the internal inductance of the plasma [H] diff --git a/bluemira/geometry/coordinates.py b/bluemira/geometry/coordinates.py index 86e47d3a34..2bf8a9b4fd 100644 --- a/bluemira/geometry/coordinates.py +++ b/bluemira/geometry/coordinates.py @@ -940,6 +940,43 @@ def vector_intersect( return point +def get_bisection_line( + p1: npt.NDArray[float], + p2: npt.NDArray[float], + p3: npt.NDArray[float], + p4: npt.NDArray[float], +) -> tuple[npt.NDArray[float], npt.NDArray[float]]: + """ + Find the bisection line between two lines. + + Parameters + ---------- + p1: + The first point on the first vector (shape: (2,)). + p2: + The second point on the first vector (shape: (2,)). + p3: + The first point on the second vector (shape: (2,)). + p4: + The second point on the second vector (shape: (2,)). + + Returns + ------- + origin: + A point on that bisection line. (shape: (2,)) + direction: + A normal vector that the bisection line points in (shape: (2,)) + """ + origin = vector_intersect(p1, p2, p3, p4) + da = p2 - p1 + db = p4 - p3 + normed_da = da / np.linalg.norm(da) + normed_db = db / np.linalg.norm(db) + dc = normed_da + normed_db + direction = dc / np.linalg.norm(dc) + return origin, direction + + # ============================================================================= # Coordinate class and parsers # ============================================================================= @@ -1722,3 +1759,18 @@ def join_intersect( args.append(coords1.argmin([x, 0, z])) return list(set(args)) return None + + +def choose_direction( + vector: npt.NDArray[float], + lower_pt: npt.NDArray[float], + higher_pt: npt.NDArray[float], +): + """ + Flip the vector to the correct side (multiply by +1 or -1) so that + when lower_pt is projected onto the vector, it has a smaller value than + when higher_pt is projected onto the vector. + """ + if (vector @ lower_pt) > (vector @ higher_pt): + return -vector + return vector diff --git a/bluemira/geometry/plane.py b/bluemira/geometry/plane.py index de21f41638..050104c1e5 100644 --- a/bluemira/geometry/plane.py +++ b/bluemira/geometry/plane.py @@ -15,14 +15,16 @@ import numpy as np +import bluemira.codes._freecadapi as cadapi +from bluemira.geometry.constants import VERY_BIG +from bluemira.geometry.face import BluemiraFace + if TYPE_CHECKING: from collections.abc import Iterable - from bluemira.geometry.placement import BluemiraPlacement + from numpy import typing as npt -import bluemira.codes._freecadapi as cadapi -from bluemira.geometry.constants import VERY_BIG -from bluemira.geometry.face import BluemiraFace + from bluemira.geometry.placement import BluemiraPlacement __all__ = ["BluemiraPlane"] @@ -164,3 +166,57 @@ def to_placement(self) -> BluemiraPlacement: from bluemira.geometry.placement import BluemiraPlacement # noqa: PLC0415 return BluemiraPlacement._create(cadapi.placement_from_plane(self._shape)) + + +def xz_plane_from_2_points( + point1: npt.NDArray[float], point2: npt.NDArray[float] +) -> BluemiraPlane: + """ + Make a plane that is perpendicular to the RZ plane using only 2 points. + + Parameters + ---------- + point1, point2: + npt.NDArray of shape (3,) + """ + # Draw an extra leg in the Y-axis direction to make an L shape. + point1 = np.array([point1[0], 0, point1[-1]]) + point2 = np.array([point2[0], 0, point2[-1]]) + point3 = point1 + np.array([0, 1, 0]) + return BluemiraPlane.from_3_points(point1, point2, point3) + + +def x_plane(x: float): + """Make a vertical plane (perpendicular to X) at a specified value of x.""" + # Simply draw an L shape in the YZ plane + return BluemiraPlane.from_3_points([x, 0, 0], [x, 0, 1], [x, 1, 0]) + + +def z_plane(z: float): + """Make a horizontal plane (perpendicular to Z) at a specified value of z.""" + # Simply draw an L shape in the XY plane + return BluemiraPlane.from_3_points([0, 0, z], [0, 1, z], [1, 0, z]) + + +def calculate_plane_dir( + start_point, end_point +) -> tuple[BluemiraPlane, npt.NDArray[float]]: + """ + Calculate the cutting plane and the direction of the cut from 2 points. + Both points must lie on the RZ plane. + + Parameters + ---------- + start_point, end_point: + 3D arrays of single points (shape = (3,)) + + Returns + ------- + plane: BluemiraPlane + A plane in 3D. + cut_direction: npt.NDArray[float] + A 3D vector (shape = (3,)) + """ + plane = xz_plane_from_2_points(start_point, end_point) + cut_direction = end_point - start_point + return plane, cut_direction diff --git a/bluemira/geometry/tools.py b/bluemira/geometry/tools.py index d4aa0f6d4f..3fb6ccc3f6 100644 --- a/bluemira/geometry/tools.py +++ b/bluemira/geometry/tools.py @@ -20,6 +20,7 @@ import numba as nb import numpy as np +from numpy import typing as npt from scipy.spatial import ConvexHull from bluemira.base.constants import EPS @@ -208,6 +209,13 @@ def _make_vertex(point: Iterable[float]) -> cadapi.apiVertex: ------- Vertex at the point """ + if isinstance(point, Coordinates): + if np.shape(point) != (3, 1): + raise GeometryError( + "Can only cast the 3D coordinates of a single point" + "into a cadapi vertex!" + ) + point = point.points[0] if len(point) != 3: # noqa: PLR2004 raise GeometryError("Points must be of dimension 3.") @@ -759,6 +767,89 @@ def convex_hull_wires_2d( return make_polygon(hull_coords, closed=True) +# # ============================================================================= +# # Volume function +# # ============================================================================= +def polygon_revolve_signed_volume(polygon: npt.ArrayLike) -> float: + """ + Revolve a polygon along the z axis, and return the volume. + + A polygon placed in the RHS of the z-axis in the xz plane would have positive volume + if it runs clockwise, and negative volume if it runs counter-clockwise. + + Similarly a polygon placed on the LHS of the z-axis in the xz plane would have + negative volume if it runs clockwise, positive volume if it runs counter-clockwise. + + Parameters + ---------- + polygon: + Stores the x-z coordinate pairs of the four coordinates. + + Notes + ----- + Consider one edge of the polygon, which has two vertices, $p$ and $c$. + TODO: insert graphics + + When revolved around the z-axis, this trapezium forms a the frustum of a cone. + The expression for the volume of this frustrum needs to be modified to avoid + ZeroDivisionError, thus it is recast into the following (also the simplest) form: + :math:`V = \\frac{\\pi}{3} (p_z - c_z) (p_x^2 + p_x c_x + c_x^2)`. + + Adding together the signed volume of all edges, the excess negative volume from one + side would cancel out the excess positive volume from the other, such that + abs(signed volume)= the volume of the polygon after being revolved around the z-axis. + """ + polygon = np.array(polygon) + if np.ndim(polygon) != 2 or np.shape(polygon)[1] != 2: # noqa: PLR2004 + raise ValueError("This function takes in an np.ndarray of shape (N, 2).") + previous_points, current_points = polygon, np.roll(polygon, -1, axis=0) + px, pz = previous_points[:, 0], previous_points[:, -1] + cx, cz = current_points[:, 0], current_points[:, -1] + volume_3_over_pi = (pz - cz) * (px**2 + px * cx + cx**2) + return np.pi / 3 * sum(volume_3_over_pi) + + +def partial_diff_of_volume( + three_vertices: Sequence[Sequence[float]], + normalised_direction_vector: Iterable[float], +) -> float: + """ + Gives the relationship between how the the solid volume varies with the position of + one of its verticies. More precisely, it gives gives the the partial derivative of + the volume of the solid revolved out of a polygon when one vertex of that polygon + is moved in the direction specified by normalised_direction_vector. + + Parameters + ---------- + three_vertices: + Contain (x, z) coordinates of the polygon. It extracts only the vertex being + moved, and the two vertices around it. three_vertices[0] and three_vertices[2] + are anchor vertices that cannot be adjusted. shape (3, 2) + normalised_direction_vector: + Direction that the point is allowed to move in. shape = (2,) + + Notes + ----- + Let there be 3 points, :math:`q`, :math:`r`, and :math:`s`, forming two edges of a + polygon. When r is moved, the polygon's revolved solid volume changes. + After a hefty amount of derivation, everything cancels out to give the expression + .. math:: + + \\frac{dV}{d r_z} = q_z q_x - r_z q_x + 2 q_z r_x - 2 s_z r_x + r_z s_x - s_z s_x + \\frac{dV}{d r_x} = (q_x + r_x + s_x) (s_x - q_x) + + + The dot product between the direction of motion and the vector :math:`\\frac{dV}{dr}` + gives the required scalar derivative showing "how much does the volume change when + r is moved in a certain direction by one unit length". + """ + (qx, qz), (rx, rz), (sx, sz) = three_vertices + x_component = qz * qx - rz * qx + 2 * qz * rx - 2 * sz * rx + rz * sx - sz * sx + z_component = (qx + rx + sx) * (sx - qx) + xz_derivatives = np.array([x_component, z_component]).T + return np.pi / 3 * np.dot(normalised_direction_vector, xz_derivatives) + + # # ============================================================================= # # Shape operation # # ============================================================================= @@ -1059,6 +1150,38 @@ def slice_shape( return None +def get_wire_plane_intersect( + convex_bm_wire: BluemiraWire, plane: BluemiraPlane, cut_direction: npt.NDArray[float] +) -> npt.NDArray[float]: + """ + Cut a wire using a plane. + + Parameters + ---------- + convex_bm_wire: + The wire that we're interested in cutting. + plane: + Plane that is cutting the wire. + cut_direction: + np.ndarray with shape==(3,) + + Returns + ------- + intersection point: + np.ndarray with shape==(3,) + """ + intersection_points = slice_shape(convex_bm_wire, plane) + if len(intersection_points) > 1: + if len(intersection_points) > 2: # noqa: PLR2004 + bluemira_warn( + "convex_bm_wire expected to be a convex hull, but isn't.\n" + "Proceeding by choosing the final intersection point..." + ) + final_intersection = np.argmax(np.dot(intersection_points, cut_direction)) + return intersection_points[final_intersection] + return intersection_points[0] + + def circular_pattern( shape: BluemiraGeo, origin: tuple[float, float, float] = (0, 0, 0), @@ -1124,6 +1247,29 @@ def mirror_shape( return convert(cadapi.mirror_shape(shape.shape, base, direction), label=label) +def is_convex(points: npt.NDArray): + """ + Check that the the list of xz points are strictly convex, i.e. + Not even collinear points are allowed. + + However, repeated points are allowed, as the repeated point would be ignored; + and points are allowed to be entered in 3D (xyz), but the y component would be + ignored as well. + + Parameters + ---------- + points + A list of points that we want to check the convexity for. Shape = (n, 2/3) + + Returns + ------- + boolean + """ + if np.shape(points)[1] == 3: # noqa: PLR2004 + points = np.array(points)[:, ::2] # squash 3D points down to 2D + return len(np.unique(points, axis=0)) == ConvexHull(points).nsimplex + + # # ============================================================================= # # Save functions # # ============================================================================= @@ -1297,30 +1443,35 @@ def signed_distance_2D_polygon( return d -def signed_distance(wire_1: BluemiraWire, wire_2: BluemiraWire) -> float: +def signed_distance( + origin: BluemiraWire | Coordinates, target: BluemiraWire | Coordinates +) -> float: """ Single-valued signed "distance" function between two wires. Will return negative - values if wire_1 does not touch or intersect wire_2, 0 if there is one intersection, + values if origin does not touch or intersect target, 0 if there is one intersection, and a positive estimate of the intersection length if there are overlaps. Parameters ---------- - wire_1: - Subject wire - wire_2: - Target wire + origin: + a 0D/1D set of points + target: + a 0D/1D set of points Returns ------- - Signed distance from wire_1 to wire_2 + Closest distance between origin and target Notes ----- This is not a pure implementation of a distance function, as for overlapping wires a metric of the quantity of overlap is returned (a positive value). This nevertheless enables the use of such a function as a constraint in gradient-based optimisers. + + This function has been extended to allow the target wire to be a point + (:class:`~bluemira.geometry.coordinates.Coordinates`) as well """ - d, vectors = distance_to(wire_1, wire_2) + d, vectors = distance_to(origin, target) # Intersections are exactly 0.0 if d == 0.0: if len(vectors) <= 1: @@ -1345,6 +1496,30 @@ def signed_distance(wire_1: BluemiraWire, wire_2: BluemiraWire) -> float: return -d +def raise_error_if_overlap( + origin: BluemiraWire | Coordinates, + target: BluemiraWire | Coordinates, + origin_name: str = "", + target_name: str = "", +): + """ + Raise an error if two wires/points intersects overlaps. + """ + check_overlaps = signed_distance(origin, target) + if check_overlaps < -D_TOLERANCE: + return + if not origin_name: + origin_name = "origin " + origin.__class__.__name__ + if not target_name: + target_name = "target " + target.__class__.__name__ + if -D_TOLERANCE <= check_overlaps <= 0: + # Sometimes intersecting lines can still appears to separate (negative), + # but only by just a little. So a small negative number is included in the check. + raise GeometryError(f"{origin_name} likely intersects {target_name} !") + if check_overlaps > 0: + raise GeometryError(f"{origin_name} and {target_name} partially/fully overlaps!") + + # ====================================================================================== # Boolean operations # ====================================================================================== diff --git a/bluemira/radiation_transport/error.py b/bluemira/radiation_transport/error.py index 449708f47b..f96fe5ff42 100644 --- a/bluemira/radiation_transport/error.py +++ b/bluemira/radiation_transport/error.py @@ -21,3 +21,11 @@ class AdvectionTransportError(RadiationTransportError): """ Error class for advective transport solver. """ + + +class NeutronicsError(RadiationTransportError): + """Error for neutronics""" + + +class SourceError(NeutronicsError): + """Error for sources""" diff --git a/bluemira/radiation_transport/neutronics/__init__.py b/bluemira/radiation_transport/neutronics/__init__.py new file mode 100644 index 0000000000..5783dd2ea4 --- /dev/null +++ b/bluemira/radiation_transport/neutronics/__init__.py @@ -0,0 +1,8 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +""" +Neutronics simulation. +""" diff --git a/bluemira/radiation_transport/neutronics/blanket_data.py b/bluemira/radiation_transport/neutronics/blanket_data.py new file mode 100644 index 0000000000..123c78d870 --- /dev/null +++ b/bluemira/radiation_transport/neutronics/blanket_data.py @@ -0,0 +1,595 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later + +"""Create specific materials from known blanket data.""" + +from copy import deepcopy +from dataclasses import dataclass +from enum import Enum, auto + +from bluemira.base.constants import raw_uc +from bluemira.materials.material import MassFractionMaterial +from bluemira.materials.mixtures import HomogenisedMixture, MixtureFraction +from bluemira.radiation_transport.neutronics.materials import NeutronicsMaterials + +# Elements +he_cool_mat = MassFractionMaterial( + name="helium", nuclides={"He4": 1.0}, density=0.008867, percent_type="ao" +) + +tungsten_mat = MassFractionMaterial( + name="tungsten", + nuclides={ + "W182": 0.266, + "W183": 0.143, + "W184": 0.307, + "W186": 0.284, + }, + percent_type="ao", + density=19.3, + density_unit="g/cm3", +) + +water_mat = MassFractionMaterial( + name="water", + nuclides={"H1": 2, "O16": 1}, + percent_type="ao", + density=0.866, + density_unit="g/cm3", +) + +al2o3_mat = MassFractionMaterial( + name="Aluminium Oxide", + nuclides={"Al27": 2, "O16": 3}, + percent_type="ao", + density=3.95, + density_unit="g/cm3", +) + + +# alloys +eurofer_mat = MassFractionMaterial( + name="eurofer", + elements={"Fe": 0.9006, "Cr": 0.0886}, + nuclides={ + "W182": 0.0108 * 0.266, + "W183": 0.0108 * 0.143, + "W184": 0.0108 * 0.307, + "W186": 0.0108 * 0.284, + }, + percent_type="wo", + density=7.78, + density_unit="g/cm3", +) + + +Be12Ti_mat = MassFractionMaterial( + name="Be12Ti", + elements={"Be": 12, "Ti": 1}, + percent_type="ao", + density=2.25, + density_unit="g/cm3", +) + + +# Lithium-containing materials +def make_PbLi_mat(li_enrich_ao) -> MassFractionMaterial: + """Make PbLi according to the enrichment fraction inputted.""" + return MassFractionMaterial( + name="PbLi", + elements={"Pb": 0.83, "Li": 0.17}, + percent_type="ao", + enrichment=li_enrich_ao, + enrichment_target="Li6", + enrichment_type="ao", + density=9.4, + density_unit="g/cm3", + ) + + +def make_Li4SiO4_mat(li_enrich_ao) -> MassFractionMaterial: + """Making enriched Li4SiO4 from elements with enrichment of Li6 enrichment""" + return MassFractionMaterial( + name="lithium_orthosilicate", + elements={"Li": 4}, + nuclides={"Si28": 1, "O16": 4}, + percent_type="ao", + enrichment=li_enrich_ao, + enrichment_target="Li6", + enrichment_type="ao", + density=2.247 + 0.078 * (100.0 - li_enrich_ao) / 100.0, + density_unit="g/cm3", + ) + + +def make_Li2TiO3_mat(li_enrich_ao) -> MassFractionMaterial: + """Make Li2TiO3 according to the enrichment fraction inputted.""" + return MassFractionMaterial( + name="lithium_titanate", + elements={"Li": 2, "Ti": 1}, + nuclides={"O16": 3}, + percent_type="ao", + enrichment=li_enrich_ao, + enrichment_target="Li6", + enrichment_type="ao", + density=3.28 + 0.06 * (100.0 - li_enrich_ao) / 100.0, + density_unit="g/cm3", + ) + + +# mixture of existing materials +lined_euro_mat = HomogenisedMixture( + name="Eurofer with Al2O3 lining", + materials=[ + MixtureFraction(eurofer_mat, 2.0 / 2.4), + MixtureFraction(al2o3_mat, 0.4 / 2.4), + ], + percent_type="vo", +) + + +# Lithium-containing material that is also a mixture of existing materials +def make_KALOS_ACB_mat(li_enrich_ao) -> HomogenisedMixture: + """Ref: Current status and future perspectives of EU ceramic breeder development""" + return HomogenisedMixture( + name="kalos_acb", # optional name of homogeneous material + materials=[ # molar combination adjusted to atom fractions + MixtureFraction( + make_Li4SiO4_mat(li_enrich_ao), 9 * 0.65 / (9 * 0.65 + 6 * 0.35) + ), + MixtureFraction( + make_Li2TiO3_mat(li_enrich_ao), 6 * 0.35 / (9 * 0.65 + 6 * 0.35) + ), + ], + percent_type="ao", + packing_fraction=0.642, + ) # combination fraction type is by atom fraction + # todo: check if this packing fraction is correct (as set above) + # KALOS_ACB_mat.set_density("g/cm3", 2.52 * 0.642) # applying packing fraction + + +def duplicate_mat_as( + mat_to_clone: MassFractionMaterial | HomogenisedMixture, + new_name: str, + new_id: int | None = None, +) -> MassFractionMaterial | HomogenisedMixture: + """Clones and renames an OpenMC material""" + new_mat = deepcopy(mat_to_clone) + new_mat.material_id = new_id + new_mat.name = new_name + + return new_mat + + +class BlanketType(Enum): + """Types of allowed blankets, named by their acronyms.""" + + DCLL = auto() + HCPB = auto() + WCLL = auto() + + @classmethod + def _missing_(cls, value: str): + try: + return cls[value.upper()] + except KeyError: + raise ValueError( + f"Invalid query: {value}. Choose from: {(*cls._member_names_,)}" + ) from None + + +@dataclass +class ReactorBaseMaterials: + """Minimum set of materials that can create a tokamak. + The rest can be populated by duplication using a priori knowledge, + e.g. inboard material = outboard material etc. + """ + + inb_vv_mat: HomogenisedMixture + inb_fw_mat: HomogenisedMixture + inb_bz_mat: HomogenisedMixture + inb_mani_mat: HomogenisedMixture + divertor_mat: HomogenisedMixture + div_fw_mat: HomogenisedMixture + + +@dataclass +class BreederTypeParameters: + """Dataclass to hold information about the breeder blanket material + and design choices. + """ + + enrichment_fraction_Li6: float + blanket_type: BlanketType + + +def _make_dcll_mats(li_enrich_ao: float) -> ReactorBaseMaterials: + """Creates openmc material definitions for a dcll blanket. + + Parameters + ---------- + --_------- + li_enrich_ao: float + Enrichment of Li-6 as a percentage + to be parsed as argument to openmc.Material.add_element + + Notes + ----- + Divertor definition from Neutronic analyses of the preliminary + design of a DCLL blanket for the EUROfusion DEMO power, 24 March 2016 + Using Eurofer instead of SS316LN + """ + inb_vv_mat = HomogenisedMixture( + name="inb_vacuum_vessel", + material_id=104, + materials=[ + MixtureFraction(eurofer_mat, 0.8), + MixtureFraction(water_mat, 0.2), + ], + percent_type="vo", + ) + + # Making first wall + inb_fw_mat = HomogenisedMixture( + name="inb_first_wall", + material_id=101, + materials=[ + MixtureFraction(tungsten_mat, 2.0 / 27.0), + MixtureFraction(eurofer_mat, 1.5 / 27.0), + MixtureFraction(he_cool_mat, 12.0 / 27.0), + MixtureFraction(lined_euro_mat, 11.5 / 27.0), + ], + percent_type="vo", + ) + + # Making blanket + inb_bz_mat = HomogenisedMixture( + name="inb_breeder_zone", + material_id=102, + materials=[ + MixtureFraction(lined_euro_mat, 0.0605 + 0.9395 * 0.05), + MixtureFraction(make_PbLi_mat(li_enrich_ao), 0.9395 * 0.95), + ], + percent_type="vo", + ) + + return ReactorBaseMaterials( + inb_vv_mat=inb_vv_mat, + inb_fw_mat=inb_fw_mat, + inb_bz_mat=inb_bz_mat, + inb_mani_mat=HomogenisedMixture( + name="inb_manifold", + material_id=103, + materials=[ + MixtureFraction(eurofer_mat, 0.573), + MixtureFraction(inb_bz_mat, 0.426), + ], # 1% void + percent_type="vo", + ), + divertor_mat=duplicate_mat_as(inb_vv_mat, "divertor", 301), + div_fw_mat=duplicate_mat_as(inb_fw_mat, "div_first_wall", 302), + ) + + +def _make_hcpb_mats(li_enrich_ao: float) -> ReactorBaseMaterials: + """Creates openmc material definitions for an hcpb blanket. + + Parameters + ---------- + li_enrich_ao: + Enrichment of Li-6 as a percentage + to be parsed as argument to openmc.Material.add_element + + Notes + ----- + HCPB Design Report, 26/07/2019 + WPBB-DEL-BB-1.2.1-T005-D001 + """ + inb_vv_mat = HomogenisedMixture( + name="inb_vacuum_vessel", # optional name of homogeneous material + material_id=104, + materials=[ + MixtureFraction(eurofer_mat, 0.6), + MixtureFraction(water_mat, 0.4), + ], + percent_type="vo", + ) + + # Making blanket + structural_fraction_vo = 0.128 + multiplier_fraction_vo = 0.493 # 0.647 + breeder_fraction_vo = 0.103 # 0.163 + helium_fraction_vo = 0.276 # 0.062 + + return ReactorBaseMaterials( + inb_vv_mat=inb_vv_mat, + inb_fw_mat=HomogenisedMixture( + name="inb_first_wall", # optional name of homogeneous material + material_id=101, + materials=[ + MixtureFraction(tungsten_mat, 2.0 / 27.0), + MixtureFraction(eurofer_mat, 25.0 * 0.573 / 27.0), + MixtureFraction(he_cool_mat, 25.0 * 0.427 / 27.0), + ], + percent_type="vo", + ), + inb_bz_mat=HomogenisedMixture( + name="inb_breeder_zone", + material_id=102, + materials=[ + MixtureFraction(eurofer_mat, structural_fraction_vo), + MixtureFraction(Be12Ti_mat, multiplier_fraction_vo), + MixtureFraction(make_KALOS_ACB_mat(li_enrich_ao), breeder_fraction_vo), + MixtureFraction(he_cool_mat, helium_fraction_vo), + ], + percent_type="vo", + ), + inb_mani_mat=HomogenisedMixture( + name="inb_manifold", + material_id=103, + materials=[ + MixtureFraction(eurofer_mat, 0.4724), + MixtureFraction(make_KALOS_ACB_mat(li_enrich_ao), 0.0241), + MixtureFraction(he_cool_mat, 0.5035), + ], + percent_type="vo", + ), + divertor_mat=duplicate_mat_as(inb_vv_mat, "divertor", 301), + div_fw_mat=HomogenisedMixture( + name="div_first_wall", + material_id=302, + materials=[ + MixtureFraction(tungsten_mat, 16.0 / 25.0), + MixtureFraction(water_mat, 4.5 / 25.0), + MixtureFraction(eurofer_mat, 4.5 / 25.0), + ], + percent_type="vo", + ), + ) + + +def _make_wcll_mats(li_enrich_ao: float) -> ReactorBaseMaterials: + """Creates openmc material definitions for a wcll blanket + + Parameters + ---------- + li_enrich_ao: + Enrichment of Li-6 as a percentage + to be parsed as argument to openmc.Material.add_element + + Notes + ----- + Ref. D. Nevo and M. Oron-Carl, WCLL Design Report 2018, Eurofusion, + WPBB-DEL-BB-3.2.1-T005-D001, June 2019. + """ + PbLi_mat = make_PbLi_mat(li_enrich_ao) + + # Divertor definition from Neutronic analyses of the preliminary + # design of a DCLL blanket for the EUROfusion DEMO power, 24 March 2016 + # Using Eurofer instead of SS316LN + inb_fw_mat = HomogenisedMixture( + name="inb_first_wall", + material_id=101, + materials=[ + MixtureFraction(tungsten_mat, 0.0766), + MixtureFraction(water_mat, 0.1321), + MixtureFraction(eurofer_mat, 0.7913), + ], + percent_type="vo", + ) + + return ReactorBaseMaterials( + inb_vv_mat=HomogenisedMixture( + name="inb_vacuum_vessel", + material_id=104, + materials=[ + MixtureFraction(eurofer_mat, 0.6), + MixtureFraction(water_mat, 0.4), + ], + percent_type="vo", + ), + inb_fw_mat=inb_fw_mat, + inb_bz_mat=HomogenisedMixture( + name="inb_breeder_zone", + material_id=102, + materials=[ + MixtureFraction(tungsten_mat, 0.0004), + MixtureFraction(PbLi_mat, 0.8238), + MixtureFraction(water_mat, 0.0176), + MixtureFraction(eurofer_mat, 0.1582), + ], + percent_type="vo", + ), + inb_mani_mat=HomogenisedMixture( + name="inb_manifold", + material_id=103, + materials=[ + MixtureFraction(PbLi_mat, 0.2129), + MixtureFraction(water_mat, 0.2514), + MixtureFraction(eurofer_mat, 0.5357), + ], + percent_type="vo", + ), + divertor_mat=duplicate_mat_as(eurofer_mat, "divertor", 301), + div_fw_mat=duplicate_mat_as(inb_fw_mat, "div_first_wall", 302), + ) + + +@dataclass(frozen=True) +class TokamakGeometry: + """The thickness measurements for all of the generic components of the tokamak. + + Parameters + ---------- + inb_fw_thick: + inboard first wall thickness [m] + inb_bz_thick: + inboard breeding zone thickness [m] + inb_mnfld_thick: + inboard manifold thickness [m] + inb_vv_thick: + inboard vacuum vessel thickness [m] + tf_thick: + toroidal field coil thickness [m] + outb_fw_thick: + outboard first wall thickness [m] + outb_bz_thick: + outboard breeding zone thickness [m] + outb_mnfld_thick: + outboard manifold thickness [m] + outb_vv_thick: + outboard vacuum vessel thickness [m] + inb_gap: + inboard gap [m] + """ + + inb_fw_thick: float + inb_bz_thick: float + inb_mnfld_thick: float + inb_vv_thick: float + tf_thick: float + outb_fw_thick: float + outb_bz_thick: float + outb_mnfld_thick: float + outb_vv_thick: float + inb_gap: float + + +def get_preset_physical_properties( + blanket_type: str | BlanketType, +) -> tuple[BreederTypeParameters, TokamakGeometry]: + """ + Works as a switch-case for choosing the tokamak geometry + and blankets for a given blanket type. + The allowed list of blanket types are specified in BlanketType. + Currently, the blanket types with pre-populated data in this function are: + {'wcll', 'dcll', 'hcpb'} + """ + blanket_type = BlanketType(blanket_type) + + breeder_materials = BreederTypeParameters( + blanket_type=blanket_type, + enrichment_fraction_Li6=0.60, + ) + + # Geometry variables + + # Break down from here. + # Paper inboard build --- + # Nuclear analyses of solid breeder blanket options for DEMO: + # Status,challenges and outlook, + # Pereslavtsev, 2019 + # + # 0.400, # TF Coil inner + # 0.200, # gap from figures + # 0.060, # VV steel wall + # 0.480, # VV + # 0.060, # VV steel wall + # 0.020, # gap from figures + # 0.350, # Back Supporting Structure + # 0.060, # Back Wall and Gas Collectors Back wall = 3.0 + # 0.350, # breeder zone + # 0.022 # fw and armour + + shared_building_geometry = { # that are identical in all three types of reactors. + "inb_gap": 0.2, # [m] + "inb_vv_thick": 0.6, # [m] + "tf_thick": 0.4, # [m] + "outb_vv_thick": 0.6, # [m] + } + if blanket_type is BlanketType.WCLL: + tokamak_geometry = TokamakGeometry( + **shared_building_geometry, + inb_fw_thick=0.027, # [m] + inb_bz_thick=0.378, # [m] + inb_mnfld_thick=0.435, # [m] + outb_fw_thick=0.027, # [m] + outb_bz_thick=0.538, # [m] + outb_mnfld_thick=0.429, # [m] + ) + elif blanket_type is BlanketType.DCLL: + tokamak_geometry = TokamakGeometry( + **shared_building_geometry, + inb_fw_thick=0.022, # [m] + inb_bz_thick=0.300, # [m] + inb_mnfld_thick=0.178, # [m] + outb_fw_thick=0.022, # [m] + outb_bz_thick=0.640, # [m] + outb_mnfld_thick=0.248, # [m] + ) + elif blanket_type is BlanketType.HCPB: + # HCPB Design Report, 26/07/2019 + tokamak_geometry = TokamakGeometry( + **shared_building_geometry, + inb_fw_thick=0.027, # [m] + inb_bz_thick=0.460, # [m] + inb_mnfld_thick=0.560, # [m] + outb_fw_thick=0.027, # [m] + outb_bz_thick=0.460, # [m] + outb_mnfld_thick=0.560, # [m] + ) + return breeder_materials, tokamak_geometry + + +def create_materials( + breeder_materials: BreederTypeParameters, +) -> NeutronicsMaterials: + """ + Parameters + ---------- + breeder_materials: + dataclass containing attributes: 'blanket_type', 'enrichment_fraction_Li6' + """ + return _create_from_blanket_type( + breeder_materials.blanket_type, + raw_uc(breeder_materials.enrichment_fraction_Li6, "", "%"), + ) + + +def _create_from_blanket_type( + blanket_type: BlanketType, li_enrich_ao: float +) -> NeutronicsMaterials: + """Create Materials Library by specifying just the blanket type + + Parameters + ---------- + blanket_type: + the blanket type + li_enrich_ao: + PERCENTAGE enrichment of Li6 (float between 0 - 100) + """ + if blanket_type is BlanketType.DCLL: + base_materials = _make_dcll_mats(li_enrich_ao) + elif blanket_type is BlanketType.HCPB: + base_materials = _make_hcpb_mats(li_enrich_ao) + elif blanket_type is BlanketType.WCLL: + base_materials = _make_wcll_mats(li_enrich_ao) + return NeutronicsMaterials( + inb_vv_mat=base_materials.inb_vv_mat, + inb_fw_mat=base_materials.inb_fw_mat, + inb_bz_mat=base_materials.inb_bz_mat, + inb_mani_mat=base_materials.inb_mani_mat, + divertor_mat=base_materials.divertor_mat, + div_fw_mat=base_materials.div_fw_mat, + outb_fw_mat=duplicate_mat_as(base_materials.inb_fw_mat, "outb_first_wall", 201), + outb_bz_mat=duplicate_mat_as( + base_materials.inb_bz_mat, "outb_breeder_zone", 202 + ), + outb_mani_mat=duplicate_mat_as( + base_materials.inb_mani_mat, "outb_manifold", 203 + ), + outb_vv_mat=duplicate_mat_as( + base_materials.inb_vv_mat, "outb_vacuum_vessel", 204 + ), + tf_coil_mat=duplicate_mat_as(eurofer_mat, "tf_coil", 401), + container_mat=duplicate_mat_as(base_materials.inb_vv_mat, "container", 501), + # surfaces + inb_sf_mat=duplicate_mat_as(eurofer_mat, "inb_sf", 601), + outb_sf_mat=duplicate_mat_as(eurofer_mat, "outb_sf", 602), + div_sf_mat=duplicate_mat_as(eurofer_mat, "div_sf", 603), + # TODO get shield material + rad_shield=duplicate_mat_as(eurofer_mat, "div_sf", 604), + ) diff --git a/bluemira/radiation_transport/neutronics/constants.py b/bluemira/radiation_transport/neutronics/constants.py new file mode 100644 index 0000000000..70ef1339ad --- /dev/null +++ b/bluemira/radiation_transport/neutronics/constants.py @@ -0,0 +1,108 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""constants used for the neutronics module""" + +from periodictable import elements + +from bluemira.base.constants import ( + ELECTRON_MOLAR_MASS, + HE_MOLAR_MASS, + NEUTRON_MOLAR_MASS, + N_AVOGADRO, + raw_uc, +) +from bluemira.geometry.constants import D_TOLERANCE +from bluemira.plasma_physics.reactions import E_DT_fusion + +DTOL_CM = raw_uc(D_TOLERANCE, "m", "cm") + + +def to_cm(m): + """Converter for m to cm""" + return raw_uc(m, "m", "cm") + + +def to_m(cm): + """Converter for cm to m""" + return raw_uc(cm, "cm", "m") + + +def to_cm3(m3): + """Converter for m3 to cm3""" + return raw_uc(m3, "m^3", "cm^3") + + +# Amount of energy released in a single dt fusion reaction, in MeV. +energy_per_dt = raw_uc(E_DT_fusion(), "eV", "J") + +# Amount of energy carried away by the neutron, which is about 4/5 of that. +ALHPA_MOLAR_MASS = HE_MOLAR_MASS - ELECTRON_MOLAR_MASS + +# ignoring the binding energy of the electron, too minute. +dt_neutron_energy = energy_per_dt * ( + ALHPA_MOLAR_MASS / (ALHPA_MOLAR_MASS + NEUTRON_MOLAR_MASS) +) # [J] + +# Energy required to displace an Fe atom in Fe. See docstring of DPACoefficients +dpa_Fe_threshold_eV = 40 # Source cites 40 eV. + +# how many degrees misalignment tolerated while merging almost-parallel wires into one. +TOLERANCE_DEGREES = 6.0 + +# Default value to discretise the BluemiraWire. +# Set to 10 to preserve speed without too much loss in precision. +DISCRETISATION_LEVEL = 10 + + +# The following material science constants are in cgs. +Fe_molar_mass_g = elements.isotope("Fe").mass +Fe_density_g_cc = elements.isotope("Fe").density + + +class DPACoefficients: + """ + Get the coefficients required + + To convert the number of damage into the number of displacements. + number of atoms in region = avogadro * density * volume / molecular mass + number of atoms in 1 cc = avogadro * density / molecular mass + dpa_per_second_of_operation = src_rate * displacements / atoms + dpa_fpy = dpa_per_second_of_operation / S_TO_YEAR + + Notes + ----- + Shengli Chena, David Bernard + On the calculation of atomic displacements using damage energy + Results in Physics 16 (2020) 102835 + https://doi.org/10.1016/j.rinp.2019.102835 + """ + + def __init__( + self, + density_g_cc: float = Fe_density_g_cc, + molar_mass_g: float = Fe_molar_mass_g, + dpa_threshold_eV: float = dpa_Fe_threshold_eV, + ): + """ + Parameters + ---------- + density_g_cc: float [g/cm^2] + density of the wall material, + where the damage (in DPA) would be calculated later. + molar_mass_g: float [g/mole] + molar mass of the wall material, + where the damage (in DPA) would be calculated later. + dpa_threshold_eV: float [eV/count] + the average amount of energy dispersed + by displacing one atom in the wall material's lattice. + + Attributes/values + ----------------- + atoms_per_cc: number density, given in cgs. + displacements_per_damage_eV: + """ + self.atoms_per_cc = N_AVOGADRO * density_g_cc / molar_mass_g + self.displacements_per_damage_eV = 0.8 / (2 * dpa_threshold_eV) diff --git a/bluemira/radiation_transport/neutronics/geometry.py b/bluemira/radiation_transport/neutronics/geometry.py new file mode 100644 index 0000000000..641a6a041c --- /dev/null +++ b/bluemira/radiation_transport/neutronics/geometry.py @@ -0,0 +1,152 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""dataclasses containing parameters used to set up the openmc model.""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING + +import numpy as np + +if TYPE_CHECKING: + from bluemira.radiation_transport.neutronics.neutronics_axisymmetric import ( + NeutronicsReactorParameterFrame, + ) + + +@dataclass +class BlanketThickness: + """ + Give the depth of the interfaces between blanket layers. + + Parameters + ---------- + surface + Thickness of the surface layer of the blanket. Can be zero. + Only used for tallying purpose, i.e. not a physical component. + first_wall + Thickness of the first wall. + breeding_zone + Thickness of the breedng zone. Could be zero if the breeding zone is absent. + + Note + ---- + Thickness of the vacuum vessel is not required because we we assume it fills up the + remaining space between the manifold's end and the outer_boundary. + """ + + surface: float + first_wall: float + breeding_zone: float + manifold: float + + def get_interface_depths(self): + """Return the depth of the interface layers""" + return np.cumsum([ + self.surface, + self.first_wall, + self.breeding_zone, + ]) + + +@dataclass +class DivertorThickness: + """ + Divertor dimensions. + For now it only has 1 value: the surface layer thickness. + + Parameters + ---------- + surface + The surface layer of the divertor, which we expect to be made of a different + material (e.g. Tungsten or alloy of Tungsten) from the bulk support & cooling + structures of the divertor. + """ + + surface: float + + +@dataclass +class ToroidalFieldCoilDimension: + """ + Gives the toroidal field coil diameters. Working with the simplest assumption, we + assume that the tf coil is circular for now. + + Parameters + ---------- + inner_diameter + (i.e. inner diameter of the windings.) + outer_diameter + Outer diameter of the windings. + """ + + thickness: float + inner_radius: float + + +@dataclass +class RadiationShieldThickness: + """ + Radiation shield dimensions. + For now it only has 1 value: the wall layer thickness. + + Parameters + ---------- + wall + The wall thickness of the radiation shield. + """ + + wall: float + + +@dataclass +class TokamakDimensions: + """ + The dimensions of the simplest axis-symmetric case of the tokamak. + + Parameters + ---------- + inboard + thicknesses of the inboard blanket + outboard + thicknesses of the outboard blanket + divertor + thicknesses of the divertor components + cs_coil + diameters of the toroidal field coil in the + """ + + inboard: BlanketThickness + inboard_outboard_transition_radius: float + outboard: BlanketThickness + divertor: DivertorThickness + cs_coil: ToroidalFieldCoilDimension + rad_shield: RadiationShieldThickness + + @classmethod + def from_parameterframe( + cls, params: NeutronicsReactorParameterFrame, r_inner_cut: float + ): + """Setup tokamak dimensions""" + return cls( + BlanketThickness( + params.fw_blanket_surface_tk.value, + params.inboard_fw_tk.value, + params.inboard_breeding_tk.value, + params.blk_ib_manifold.value, + ), + r_inner_cut, + BlanketThickness( + params.fw_blanket_surface_tk.value, + params.outboard_fw_tk.value, + params.outboard_breeding_tk.value, + params.blk_ob_manifold.value, + ), + DivertorThickness(params.fw_divertor_surface_tk.value), + ToroidalFieldCoilDimension(params.tk_tf_inboard.value, params.r_tf_in.value), + RadiationShieldThickness(params.tk_rs.value), + ) diff --git a/bluemira/radiation_transport/neutronics/make_pre_cell.py b/bluemira/radiation_transport/neutronics/make_pre_cell.py new file mode 100644 index 0000000000..c8ce001a9d --- /dev/null +++ b/bluemira/radiation_transport/neutronics/make_pre_cell.py @@ -0,0 +1,767 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Make pre-cells using bluemira wires.""" + +from __future__ import annotations + +from itertools import pairwise + +import numpy as np +from numpy import typing as npt + +from bluemira.base.constants import EPS +from bluemira.display import plot_2d, show_cad +from bluemira.geometry.coordinates import ( + Coordinates, + choose_direction, + get_bisection_line, +) +from bluemira.geometry.error import GeometryError +from bluemira.geometry.solid import BluemiraSolid +from bluemira.geometry.tools import ( + is_convex, + make_polygon, + raise_error_if_overlap, + revolve_shape, +) +from bluemira.geometry.wire import BluemiraWire +from bluemira.radiation_transport.neutronics.radial_wall import ( + CellWalls, + Vert, +) +from bluemira.radiation_transport.neutronics.wires import ( + CircleInfo, + StraightLineInfo, + WireInfo, + WireInfoList, +) + +CCW_90 = np.array([[0, 0, -1], [0, 1, 0], [1, 0, 0]]) +CW_90 = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) + + +def ratio_of_distances( + point_of_interest: npt.NDArray[np.float64], + anchor1: npt.NDArray[np.float64], + normal1: npt.NDArray[np.float64], + anchor2: npt.NDArray[np.float64], + normal2: npt.NDArray[np.float64], +) -> npt.NDArray[np.float64]: + """ + Find how close a point is to line 1 and line 2, and then express that ratio as a + tuple of floats that sums up to unit. + Each line is given by the user by specifying any point on that line, and a direction + NORMAL to that line. The point_of_interest must lie on the positive side of the line. + + Parameters + ---------- + point_of_interest: + point to which we want to calculate the ratio of distances. + anchor1, anchor2: + Any point on line 1 and line 2 respectively. + normal1, normal2: + The positive distance direction of line 1 and line 2 respectively. + + Returns + ------- + dist_to_1, dist_to_2: + ratio of distances. Sum of these two numbers should yield unity (1.0). + """ + dist_to_1 = (point_of_interest - anchor1) @ normal1 + dist_to_2 = (point_of_interest - anchor2) @ normal2 + if dist_to_1 < -EPS or dist_to_2 < -EPS: + raise GeometryError( + "Expecting point_of_interest to lie on the positive side of both lines!" + ) + total_dist = dist_to_1 + dist_to_2 + return np.array([dist_to_1, dist_to_2]) / total_dist + + +def find_equidistant_point( + point1: npt.NDArray[np.float64], point2: npt.NDArray[np.float64], distance: float +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """Find the two (or 0) points on a 2D plane that are equidistant to each other. + + Parameters + ---------- + point1, point2: + 2D points, each with shape (2,) + distance: + the distance that both points must obey by. + + Returns + ------- + intersection1, intersection2: + The two intersection points of circle1 and circle2. + """ + mid_point = (point1 + point2) / 2 + sep = point2 - point1 + half_sep = np.linalg.norm(sep) / 2 # scalar + if half_sep > distance: + raise GeometryError("The two points are separated by > 2 * distance!") + orth_length = np.sqrt(distance**2 - half_sep**2) + orth_dir = np.array([-sep[1], sep[0]]) + orth_dir /= np.linalg.norm(orth_dir) + return mid_point + (orth_dir * orth_length), mid_point - (orth_dir * orth_length) + + +def pick_higher_point( + point1: npt.NDArray[np.float64], + point2: npt.NDArray[np.float64], + vector: npt.NDArray[np.float64], +) -> npt.NDArray[np.float64]: + """Pick the point that, when projected onto `vector`, gives a higher value.""" + if (vector @ point1) > (vector @ point2): + return point1 + return point2 + + +def calculate_new_circle( + old_circle_info: CircleInfo, new_points: npt.NDArray +) -> CircleInfo: + """ + Calculate how far does the new circle get shifted. + + Parameters + ---------- + old_circle_info: + an object accessed by WireInfoList[i].key_points. + info on circle where the start_point and end_point are each of shape (3,). + new_points: + array of shape (2, 3) + + Returns + ------- + new_circle_info: + An instance of CircleInfo representing the new (scaled) arc of circle. + """ + new_chord_vector = np.diff(new_points, axis=0) + old_chord_vector = np.diff(old_circle_info[:2], axis=0) + scale_factor = np.linalg.norm(new_chord_vector) / np.linalg.norm(old_chord_vector) + new_radius = old_circle_info.radius * scale_factor + possible_centers = find_equidistant_point(*new_points[:, ::2], new_radius) + center1, center2 = np.insert(possible_centers, 1, 0, axis=1) + + old_chord_mid_point = np.mean(old_circle_info[:2], axis=0) + old_radius_vector = np.array(old_circle_info.center) - old_chord_mid_point + # chord should stay on the same side of the center after transformation. + new_center = pick_higher_point(center1, center2, old_radius_vector) + + return CircleInfo(*new_points, new_center, new_radius) + + +class PreCell: + """ + A pre-cell is the BluemiraWire outlining the reactor cross-section + BEFORE they have been simplified into straight-lines. + Unlike a Cell, its outline may be constructed from curved lines. + """ + + def __init__( + self, + interior_wire: BluemiraWire | Coordinates, + vv_wire: BluemiraWire, + exterior_wire: BluemiraWire, + ): + """ + Parameters + ---------- + interior_wire: + Either a wire representing the interior-boundary (i.e. plasma-facing side) + of a blanket's pre-cell, running in the counter-clockwise direction when + viewing the right hand side poloidal cross-section, + i.e. downwards if inboard, upwards if outboard. + or a single Coordinates point, representing a point on the interior-boundary + vv_wire: + A wire representing the interface between the vacuum vessel and the blanket. + exterior_wire: + A wire representing the exterior-boundary of a + blanket's pre-cell, running in the clockwise direction when viewing the + right hand side poloidal cross-section, + i.e. upwards if inboard, downwards if outboard. + """ + self.interior_wire = interior_wire + self.vv_wire = vv_wire + self.exterior_wire = exterior_wire + raise_error_if_overlap( + self.interior_wire, self.vv_wire, "interior wire", "vacuum vessel wire" + ) + raise_error_if_overlap( + self.vv_wire, self.exterior_wire, "vacuum vessel wire", "exterior wire" + ) + raise_error_if_overlap( + self.interior_wire, self.exterior_wire, "exterior wire", "interior wire" + ) + ext_start, ext_end = exterior_wire.start_point(), exterior_wire.end_point() + vv_start, vv_end = vv_wire.start_point(), vv_wire.end_point() + if isinstance(interior_wire, Coordinates): + int_start = int_end = interior_wire + self._inner_curve = make_polygon( + np.array([ext_end, interior_wire, ext_start]).T, closed=False + ) + else: + int_start, int_end = interior_wire.start_point(), interior_wire.end_point() + self._out2in = make_polygon(np.array([ext_end, int_start]).T, closed=False) + self._in2out = make_polygon(np.array([int_end, ext_start]).T, closed=False) + self._inner_curve = BluemiraWire([ + self._out2in, + self.interior_wire, + self._in2out, + ]) + raise_error_if_overlap( + self._out2in, + self._in2out, + "cell-start cutting plane", + "cell-end cutting plane", + ) + + self.vertex = Coordinates([ext_end, int_start, int_end, ext_start]).xz + self.vv_point = Coordinates([vv_start, vv_end]).xz + # The wire outlining the PreCell + self.outline = BluemiraWire([self.exterior_wire, self._inner_curve]) + + @property + def half_solid(self) -> BluemiraSolid: + """ + Create the 180° revolved shape on demand only. + Revolved 180° instead of 360° for easier viewing + """ + if not hasattr(self, "_half_solid"): + self._half_solid = BluemiraSolid(revolve_shape(self.outline)) + return self._half_solid + + @property + def blanket_outline(self) -> BluemiraWire: + """ + Create the outline of the blanket, i.e. the part excluding the vacuum vessel. + """ + if not hasattr(self, "_blanket_outline"): + vv_start, vv_end = self.vv_wire.start_point(), self.vv_wire.end_point() + if isinstance(self.interior_wire, Coordinates): + in_start = in_end = self.interior_wire + inner_curve = make_polygon( + np.array([vv_end, self.interior_wire, vv_start]).T, closed=False + ) + else: + in_start = self.interior_wire.start_point() + in_end = self.interior_wire.end_point() + inner_curve = BluemiraWire([ + make_polygon(np.array([vv_end, in_start]).T, closed=False), + self.interior_wire, + make_polygon(np.array([in_end, vv_start]).T, closed=False), + ]) + self._blanket_outline = BluemiraWire([self.vv_wire, inner_curve]) + return self._blanket_outline + + @property + def blanket_half_solid(self) -> BluemiraSolid: + """Get the volume of the blanket""" + if not hasattr(self, "_blanket_half_solid"): + self._blanket_half_solid = BluemiraSolid(revolve_shape(self.blanket_outline)) + return self._blanket_half_solid + + def plot_2d(self, *args, **kwargs): + """Plot the outline in 2D""" + return plot_2d([self.outline, self.vv_wire], *args, **kwargs) + + def show_cad(self, *args, **kwargs): + """Plot the outline in 3D""" + return show_cad(self.half_solid, *args, **kwargs) + + @property + def cell_walls(self) -> CellWalls: + """ + The side (clockwise side and counter-clockwise) walls of this cell. + Only create it when called, because some instances of PreCell will never use it. + + it is of type + :class:`~bluemira.radiation_transport.neutronics.radial_wall.CellWalls`. + """ + if not hasattr(self, "_cell_walls"): + self._cell_walls = CellWalls( + np.asarray(( + self.vertex.T[(Vert.interior_end, Vert.exterior_start),], + self.vertex.T[(Vert.interior_start, Vert.exterior_end),], + )) + ) + + return self._cell_walls + + @property + def normal_to_interior(self) -> npt.NDArray: + """ + The vector pointing from the interior_wire direction towards the exterior_wire, + specifically, it's perpendicular to the interior_wire. + Also only created when called, because it's not needed + """ + if not hasattr(self, "_normal_to_interior"): + if isinstance(self.interior_wire, Coordinates): + self._normal_to_interior = get_bisection_line( + *self.cell_walls.reshape([4, 2]) + )[1] + else: + interior_vector = self.cell_walls.starts[0] - self.cell_walls.starts[1] + normal = np.array([interior_vector[-1], -interior_vector[0]]) + self._normal_to_interior = normal / np.linalg.norm(normal) + + return self._normal_to_interior + + def get_cell_wall_cut_points_by_fraction(self, fraction: float) -> npt.NDArray: + """ + Find the cut points on the cell's side walls by multiplying the original lengths + by a fraction. When fraction=0, this returns the interior_start and interior_end. + + Parameters + ---------- + fraction: + A scalar value + + Returns + ------- + new end points + The position of the pre-cell wall end points at the required fraction, array + of shape (2, 2) [[cw_wall x, cw_wall z], [ccw_wall x, ccw_wall z]]. + """ + return self.cell_walls.calculate_new_end_points( + self.cell_walls.lengths * fraction + ) + + def get_cell_wall_cut_points_by_thickness(self, thickness: float): + """ + Offset a line parallel to the interior_wire towards the exterior direction. + Then, find where this line intersect the cell's side walls. + + Parameters + ---------- + thickness: + A scalar value + + Returns + ------- + new end points + The position of the pre-cell wall end points at the required thickness, array + of shape (2, 2) [[cw_wall x, cw_wall z], [ccw_wall x, ccw_wall z]]. + """ + return self.cell_walls.calculate_new_end_points( + thickness / (self.cell_walls.directions @ self.normal_to_interior) + ) + + +class PreCellArray: + """ + A list of pre-cells materials + + Parameters + ---------- + list_of_pre_cells: + An adjacent list of pre-cells + """ + + def __init__(self, list_of_pre_cells: list[PreCell]): + """The list of pre-cells must be ajacent to each other.""" + self.pre_cells = list(list_of_pre_cells) + for this_cell, next_cell in pairwise(self.pre_cells): + # perform check that they are actually adjacent + this_wall = ( + this_cell.vertex[:, Vert.exterior_end], + this_cell.vertex[:, Vert.interior_start], + ) + next_wall = ( + next_cell.vertex[:, Vert.exterior_start], + next_cell.vertex[:, Vert.interior_end], + ) + if not np.array_equal(this_wall, next_wall): + raise GeometryError( + "Adjacent pre-cells are expected to have matching" + f"corners; but instead we have {this_wall}!={next_wall}." + ) + self.cell_walls = CellWalls.from_pre_cell_array(self) + + if not is_convex(self.exterior_vertices()[:, ::2]): + raise GeometryError(f"{self} must have convex exterior wires!") + + @property + def volumes(self) -> tuple[float]: + """Create the iterable of volumes on demand.""" + if not hasattr(self, "_volumes"): + # Immutable property, hence wrapped in tuple. + volume_list = [] + for pre_cell in self.pre_cells: + blanket_volume = pre_cell.blanket_half_solid.volume * 2 + vv_volume = pre_cell.half_solid.volume * 2 - blanket_volume + volume_list.append((blanket_volume, vv_volume)) + self._volumes = tuple(volume_list) + return self._volumes + + def straighten_exterior(self, *, preserve_volume: bool = False) -> PreCellArray: + """ + Turn the exterior curves of each cell into a straight edge. + This is done at the PreCellArray level instead of the PreCell level to allow + volume preservation, see Parameters below for more details. + + Parameters + ---------- + preserve_volume: + Whether to preserve the volume of each cell during the transformation from + pre-cell with curved-edge to pre-cell with straight edges. + If True, increase the length of the cut lines appropriately to compensate for + the volume loss due to the straight line approximation. + """ + exterior_walls_copy = self.cell_walls.copy() + interior_walls_copy = CellWalls.from_pre_cell_array_vv(self) + if preserve_volume: + blanket_volumes, vv_volumes = np.array(self.volumes).T + total_volumes = blanket_volumes + vv_volumes + interior_walls_copy.optimise_to_match_individual_volumes(blanket_volumes) + exterior_walls_copy.optimise_to_match_individual_volumes(total_volumes) + new_pre_cells = [] + for i in range(len(self.pre_cells)): + j = i + 1 + straightened_vv_interior = make_polygon( + [ + [interior_walls_copy[i][1, 0], interior_walls_copy[j][1, 0]], + [0, 0], + [interior_walls_copy[i][1, 1], interior_walls_copy[j][1, 1]], + ], # fill it back up to 3D to make the polygon + label="Straight edge approximation of the interior of the vacuum vessel " + f"at pre-cell {i}", + closed=False, + ) + straightened_exterior = make_polygon( + [ + [exterior_walls_copy[i][1, 0], exterior_walls_copy[j][1, 0]], + [0, 0], + [exterior_walls_copy[i][1, 1], exterior_walls_copy[j][1, 1]], + ], # fill it back up to 3D to make the polygon + label=f"straight edge approximation of the exterior of pre-cell {i}", + closed=False, + ) + new_pre_cells.append( + PreCell( + self[i].interior_wire, + straightened_vv_interior, + straightened_exterior, + ) + ) + return PreCellArray(new_pre_cells) + + def plot_2d(self, *args, **kwargs): + """Plot pre cells in 2d""" + return plot_2d( + [ + *(pc.outline for pc in self.pre_cells), + *(pc.vv_wire for pc in self.pre_cells), + ], + *args, + **kwargs, + ) + + def show_cad(self, *args, **kwargs): + """Show pre cell CAD""" + return show_cad([pc.half_solid for pc in self.pre_cells], *args, **kwargs) + + def exterior_vertices(self) -> npt.NDArray: + """ + Returns all of the vertices on the exterior side of the pre-cell array. + + Returns + ------- + exterior_vertices: + array of shape (N+1, 3) arranged clockwise (inboard to outboard). + """ + return np.insert(self.cell_walls[:, 1], 1, 0, axis=-1) + + def interior_vertices(self) -> npt.NDArray: + """ + Returns all of the vertices on the interior side of the pre-cell array. + + Returns + ------- + interior_vertices: + array of shape (N+1, 3) arranged clockwise (inboard to outboard). + """ + return np.insert(self.cell_walls[:, 0], 1, 0, axis=-1) + + def __len__(self) -> int: + """Number of pre cells""" + return len(self.pre_cells) + + def __getitem__(self, index_or_slice) -> list[PreCell] | PreCell: + """Get pre cell""" + return self.pre_cells[index_or_slice] + + def __setitem__(self, index_or_slice, new_pre_cell: PreCell | PreCellArray): + """Set an element to be a new Precell, or a slice to be a new PreCellarray.""" + if isinstance(new_pre_cell, PreCell): + self.pre_cells[index_or_slice] = new_pre_cell + else: + self.pre_cells[index_or_slice] = new_pre_cell.pre_cells + + def __add__(self, other_array) -> PreCellArray: + """Adding two list together to create a new one.""" + if isinstance(other_array, PreCellArray): + return PreCellArray(self.pre_cells + other_array.pre_cells) + raise TypeError( + f"Addition not implemented between PreCellArray and {type(other_array)}" + ) + + def __repr__(self) -> str: + """String representation""" + return super().__repr__().replace(" at ", f" of {len(self)} PreCells at ") + + def copy(self): + """ + NOT a deepcopy, each element of the new_copy.pre_cells list points to the same + items as the self.pre_cells + + Returns + ------- + new_copy + copy of itself + """ + return PreCellArray(self.pre_cells.copy()) + + +class DivertorPreCell: + """ + An intermediate class between the bluemira wire and the final csg product. + A divertor pre-cell is the equivalent of a blanket's pre-cell, but for the divertor. + """ + + def __init__( + self, + interior_wire: WireInfoList, + vv_wire: WireInfoList, + exterior_wire: WireInfoList, + ): + """ + Parameters + ---------- + interior_wire: + WireInfoList of a wire on the interior side of the cell running + counter-clockwise + vv_wire: + WireInfoList of the external surface of the vacuum vessel + exterior_wire: + WireInfoList of a wire on the exterior side of the cell running clockwise + """ + self.interior_wire = interior_wire + self.vv_wire = vv_wire + self.exterior_wire = exterior_wire + # cw_wall and ccw_wall are of type WireInfoLists!! + self.cw_wall = WireInfoList([ + WireInfo.from_2P( + self.exterior_wire.end_point, self.interior_wire.start_point + ) + ]) + self.ccw_wall = WireInfoList([ + WireInfo.from_2P( + self.interior_wire.end_point, self.exterior_wire.start_point + ) + ]) + self.vertex = Coordinates( + np.asarray([ + self.cw_wall.start_point, + self.cw_wall.end_point, + self.ccw_wall.start_point, + self.ccw_wall.end_point, + ]) + ) + + def plot_2d(self, *args, **kwargs): + """Plot 2d precell""" + return plot_2d(self.outline, *args, **kwargs) + + def show_cad(self, *args, **kwargs): + """Show precell CAD""" + return show_cad(self.half_solid, *args, **kwargs) + + @property + def outline(self) -> BluemiraWire: + """ + We don't need the volume value, so we're only going to generate the outline + when the user wants to plot it. + """ + if not hasattr(self, "_outline"): + self._outline = BluemiraWire([ + self.cw_wall.restore_to_wire(), + self.interior_wire.restore_to_wire(), + self.ccw_wall.restore_to_wire(), + self.exterior_wire.restore_to_wire(), + ]) + return self._outline + + @property + def half_solid(self) -> BluemiraSolid: + """ + Create the 180° revolved shape on demand only. + Revolved 180° instead of 360° for easier viewing + """ + if not hasattr(self, "_half_solid"): + self._half_solid = BluemiraSolid(revolve_shape(self.outline)) + return self._half_solid + + def offset_interior_wire(self, thickness: float) -> WireInfoList: + """ + Offset the interior wire towards the exterior_wire. + The true problem of expanding/shrinking a wire is a much more difficult one, so + I've only opted for a simpler (but incorrect) approach of pushing the wire to a + desired direction determined by how close it is to the wall. + + TODO: Re-write this method, as it currently do weird things to circles. + + New method should do the following: + + 1. Find a point to be our new origin. (Name that point "p") + This is likely to be a point near exterior_vertices. + 2. Scale down the interior_wire by x%. + + We should give more thought of how to derive/search for the optimal vector + (p[0], p[-1], x), such that all lines are displaced by + + This proposed method has several new benefits: + + 1. The circles will be scaled correctly (center moves towards the new origin + by x%, radius scaled down by x%). + 2. All tangents are preserved, so no need to change them. + """ + int_wire_pts = np.array([ + *(w.key_points[0] for w in self.interior_wire), + self.interior_wire[-1].key_points[1], + ]) + + # assumed normalised + cw_dir = choose_direction( + self.cw_wall[0].tangents[0], self.cw_wall.end_point, self.cw_wall.start_point + ) + ccw_dir = choose_direction( + self.ccw_wall[0].tangents[1], + self.ccw_wall.start_point, + self.ccw_wall.end_point, + ) + + cw_norm = CCW_90 @ cw_dir + cw_anchor = self.cw_wall.start_point + ccw_norm = CW_90 @ ccw_dir + ccw_anchor = self.ccw_wall.end_point + + def unit(new_dir: npt.NDArray[np.float64]): + return new_dir / np.linalg.norm(new_dir) + + shifted_pts = [ + pt + + unit( + np.array([cw_dir, ccw_dir]).T + @ ratio_of_distances(pt, cw_anchor, cw_norm, ccw_anchor, ccw_norm)[::-1] + ) + * thickness + for pt in int_wire_pts + ] + + info_list = [] + for i, (new_start, new_end) in enumerate(pairwise(shifted_pts)): + old_kp = self.interior_wire[i].key_points + if isinstance(old_kp, StraightLineInfo): + info_list.append(WireInfo.from_2P(new_start, new_end)) + else: + info_list.append( + WireInfo( + calculate_new_circle(old_kp, np.array([new_start, new_end])), + # copy the old tangents, which are not going to be correct, + # but we won't be using them anyways. + self.interior_wire[i].tangents, + ) + ) + return WireInfoList(info_list) + + +class DivertorPreCellArray: + """An array of Divertor pre-cells""" + + def __init__(self, list_of_div_pc: list[DivertorPreCell]): + self.pre_cells = list(list_of_div_pc) + # Perform check that they are adjacent + for prev_cell, curr_cell in pairwise(self.pre_cells): + if not np.array_equal( + prev_cell.vertex.xyz[:, Vert.exterior_start], + curr_cell.vertex.xyz[:, Vert.exterior_end], + ): + raise GeometryError("Expect neighbouring cells to share corners!") + + if not is_convex(self.exterior_vertices()[:, ::2]): + raise GeometryError(f"{self} must have convex exterior vertices!") + + def exterior_vertices(self) -> npt.NDArray: + """ + Returns all of the tokamak's poloidal cross-section's outside corners' + coordinates, in 3D. + + Returns + ------- + exterior_vertices: + aray of shape (N+1, 3) arranged counter-clockwise (inboard to outboard). + """ + # Because cells run counter-clockwise but the exterior_wire themselves runs + # clockwise, we have to invert the wire during extraction to make it run + # without double-backing onto itself. + return np.concatenate([ + stack.exterior_wire.get_3D_coordinates()[::-1] for stack in self.pre_cells + ]) + + def interior_vertices(self) -> npt.NDArray: + """ + Returns all of the tokamak's poloidal cross-section's inside corners' + coordinates, in 3D. + + Parameters + ---------- + interior_vertices: + aray of shape (N+1, 3) arranged counter-clockwise (inboard to outboard). + """ + return np.concatenate([ + stack.interior_wire.get_3D_coordinates() for stack in self.pre_cells + ]) + + def __len__(self) -> int: + """Number of pre cells""" + return len(self.pre_cells) + + def __getitem__(self, index_or_slice) -> list[DivertorPreCell] | DivertorPreCell: + """Get pre cell""" + return self.pre_cells[index_or_slice] + + def __repr__(self) -> str: + """String representation""" + return ( + super().__repr__().replace(" at ", f" of {len(self)} DivertorPreCells at ") + ) + + def plot_2d(self, *args, **kwargs): + """Plot precell array cad in 2d""" + return plot_2d( + [ + *(dpc.outline for dpc in self.pre_cells), + *(dpc.vv_wire.restore_to_wire() for dpc in self.pre_cells), + ], + *args, + **kwargs, + ) + + def show_cad(self, *args, **kwargs): + """Show precell array CAD""" + return show_cad([dpc.half_solid for dpc in self.pre_cells], *args, **kwargs) + + def copy(self): + """ + NOT a deepcopy, each element of the new_copy.pre_cells list points to the same + items as the self.pre_cells + + Returns + ------- + new_copy + copy of itself + """ + return DivertorPreCellArray(self.pre_cells.copy()) diff --git a/bluemira/radiation_transport/neutronics/materials.py b/bluemira/radiation_transport/neutronics/materials.py new file mode 100644 index 0000000000..13583bc0b5 --- /dev/null +++ b/bluemira/radiation_transport/neutronics/materials.py @@ -0,0 +1,37 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Create the material sets for each type of reactor.""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from bluemira.materials.material import MassFractionMaterial + from bluemira.materials.mixtures import HomogenisedMixture + + +@dataclass +class NeutronicsMaterials: + """A dictionary of materials according to the type of blanket used""" + + inb_vv_mat: HomogenisedMixture | MassFractionMaterial + inb_fw_mat: HomogenisedMixture | MassFractionMaterial + inb_bz_mat: HomogenisedMixture | MassFractionMaterial + inb_mani_mat: HomogenisedMixture | MassFractionMaterial + divertor_mat: HomogenisedMixture | MassFractionMaterial + div_fw_mat: HomogenisedMixture | MassFractionMaterial + outb_fw_mat: HomogenisedMixture | MassFractionMaterial + outb_bz_mat: HomogenisedMixture | MassFractionMaterial + outb_mani_mat: HomogenisedMixture | MassFractionMaterial + outb_vv_mat: HomogenisedMixture | MassFractionMaterial + tf_coil_mat: HomogenisedMixture | MassFractionMaterial + container_mat: HomogenisedMixture | MassFractionMaterial + inb_sf_mat: HomogenisedMixture | MassFractionMaterial + outb_sf_mat: HomogenisedMixture | MassFractionMaterial + div_sf_mat: HomogenisedMixture | MassFractionMaterial + rad_shield: HomogenisedMixture | MassFractionMaterial diff --git a/bluemira/radiation_transport/neutronics/neutronics_axisymmetric.py b/bluemira/radiation_transport/neutronics/neutronics_axisymmetric.py new file mode 100644 index 0000000000..232892891d --- /dev/null +++ b/bluemira/radiation_transport/neutronics/neutronics_axisymmetric.py @@ -0,0 +1,287 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +""" +Axis-symmetric CSG CAD models for neutronics. +""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from dataclasses import dataclass +from typing import TYPE_CHECKING + +import numpy as np + +from bluemira.base.look_and_feel import bluemira_print +from bluemira.base.parameter_frame import Parameter, ParameterFrame, make_parameter_frame +from bluemira.geometry.plane import calculate_plane_dir +from bluemira.geometry.tools import get_wire_plane_intersect, make_polygon +from bluemira.radiation_transport.neutronics.make_pre_cell import PreCell +from bluemira.radiation_transport.neutronics.slicing import ( + DivertorWireAndExteriorCurve, + PanelsAndExteriorCurve, +) + +if TYPE_CHECKING: + from numpy import typing as npt + + from bluemira.base.reactor import ComponentManager + from bluemira.geometry.wire import BluemiraWire + from bluemira.radiation_transport.neutronics.geometry import TokamakDimensions + from bluemira.radiation_transport.neutronics.make_pre_cell import ( + DivertorPreCellArray, + PreCellArray, + ) + from bluemira.radiation_transport.neutronics.materials import NeutronicsMaterials + + +@dataclass +class ReactorGeometry: + """ + Data storage stage + + Parameters + ---------- + divertor_wire: + The plasma-facing side of the divertor. + panel_break_points: + The start and end points for each first-wall panel + (for N panels, the shape is (N+1, 2)). + boundary: + interface between the inside of the vacuum vessel and the outside of the blanket + vacuum_vessel_wire: + The outer-boundary of the vacuum vessel + """ + + divertor_wire: BluemiraWire + panel_break_points: npt.NDArray + boundary: BluemiraWire + vacuum_vessel_wire: BluemiraWire + + +@dataclass +class CuttingStage: + """Stage of making cuts to the exterior curve/ outer boundary.""" + + blanket: PanelsAndExteriorCurve + divertor: DivertorWireAndExteriorCurve + + +class PreCellStage: + """Stage of making pre-cells""" + + def __init__(self, blanket: PreCellArray, divertor: DivertorPreCellArray): + """Check convexity after initialization""" + self.blanket = blanket.copy() + self.divertor = divertor.copy() + # 1. stretch first blanket cell in PreCellArray to reach div_start_wire + div_start_wire = self.divertor[0].cw_wall.restore_to_wire() + # pull everything down to: div_start_wire. + # Alternatively, choose div_start_wire=self.divertor[0].outline + old_vv_wire = self.blanket[0].vv_wire + ext_pt, i_high, i_low = np.insert(self.blanket[0].vertex, 1, 0, axis=0).T[:3] + i_end = get_wire_plane_intersect( + div_start_wire, *calculate_plane_dir(i_high, i_low) + ) + # v_end = get_wire_plane_intersect( + # div_start_wire, + # *calculate_plane_dir(old_vv_wire.end_point().xyz.flatten(), + # old_vv_wire.start_point().xyz.flatten()) + # ) + in_wire = make_polygon(np.array([i_high, i_end]).T, closed=False) + vv_wire = make_polygon( + np.array([ + self.divertor[0].vv_wire.end_point, + old_vv_wire.end_point().xyz.flatten(), + ]).T, + closed=False, + ) + ex_wire = make_polygon( + np.array([self.divertor[0].vertex.T[0], ext_pt]).T, closed=False + ) + new_start_cell = PreCell(in_wire, vv_wire, ex_wire) + self.blanket[0] = new_start_cell + + # 2. stretch first blanket cell in PreCellArray to reach div_end_wire + div_end_wire = self.divertor[-1].ccw_wall.restore_to_wire() + old_vv_wire = self.blanket[-1].vv_wire + i_low, i_high, ext_pt = np.insert(self.blanket[-1].vertex, 1, 0, axis=0).T[-3:] + i_start = get_wire_plane_intersect( + div_end_wire, *calculate_plane_dir(i_high, i_low) + ) + # v_end = get_wire_plane_intersect( + # div_end_wire, + # *calculate_plane_dir(old_vv_wire.start_point().xyz.flatten(), + # old_vv_wire.end_point().xyz.flatten()) + # ) + in_wire = make_polygon(np.array([i_start, i_high]).T, closed=False) + vv_wire = make_polygon( + np.array([ + old_vv_wire.start_point().xyz.flatten(), + self.divertor[-1].vv_wire.start_point, + ]).T, + closed=False, + ) + ex_wire = make_polygon( + np.array([ext_pt, self.divertor[-1].vertex.T[-1]]).T, closed=False + ) + new_end_cell = PreCell(in_wire, vv_wire, ex_wire) + self.blanket[-1] = new_end_cell + + # re-initialize so that the cell_walls are re-calculated + self.blanket = self.blanket.copy() + + def external_coordinates(self) -> npt.NDArray: + """ + Get the outermost coordinates of the tokamak cross-section from pre-cell array + and divertor pre-cell array. + Runs clockwise, beginning at the inboard blanket-divertor joint. + """ + return np.concatenate([ + self.blanket.exterior_vertices(), + self.divertor.exterior_vertices()[::-1], + ]) + + def bounding_box(self) -> tuple[float, ...]: + """Get bounding box of pre cell stage""" + all_ext_vertices = self.external_coordinates() + z_min = all_ext_vertices[:, -1].min() + z_max = all_ext_vertices[:, -1].max() + r_max = max(abs(all_ext_vertices[:, 0])) + return z_max, z_min, r_max, -r_max + + def half_bounding_box(self) -> tuple[float, ...]: + """ + Get bounding box of the 2D poloidal cross-section of the right-hand half of the + reactor. + """ + all_ext_vertices = self.external_coordinates() + z_min = all_ext_vertices[:, -1].min() + z_max = all_ext_vertices[:, -1].max() + r_max = max(abs(all_ext_vertices[:, 0])) + r_min = min(abs(all_ext_vertices[:, 0])) + return z_max, z_min, r_max, r_min + + +@dataclass +class NeutronicsReactorParameterFrame(ParameterFrame): + """Neutronics reactor parameters""" + + inboard_fw_tk: Parameter[float] + inboard_breeding_tk: Parameter[float] + outboard_fw_tk: Parameter[float] + outboard_breeding_tk: Parameter[float] + r_tf_in: Parameter[float] + tk_tf_inboard: Parameter[float] + fw_divertor_surface_tk: Parameter[float] + fw_blanket_surface_tk: Parameter[float] + blk_ib_manifold: Parameter[float] + tk_rs: Parameter[float] + blk_ob_manifold: Parameter[float] + + +class NeutronicsReactor(ABC): + """Pre csg cell reactor""" + + param_cls = NeutronicsReactorParameterFrame + + def __init__( + self, + params: dict | ParameterFrame, + divertor: ComponentManager, + blanket: ComponentManager, + vacuum_vessel: ComponentManager, + materials_library: NeutronicsMaterials, + *, + snap_to_horizontal_angle: float = 45, + blanket_discretisation: int = 10, + divertor_discretisation: int = 5, + ): + bluemira_print("Creating axis-symmetric neutronics model") + + self.params = make_parameter_frame(params, self.param_cls) + self.material_library = materials_library + ( + self.tokamak_dimensions, + divertor_wire, + panel_points, + blanket_wire, + vacuum_vessel_wire, + ) = self._get_wires_from_components(divertor, blanket, vacuum_vessel) + + self.geom = ReactorGeometry( + divertor_wire, panel_points, blanket_wire, vacuum_vessel_wire + ) + + self._pre_cell_stage = self._create_pre_cell_stage( + blanket_discretisation, divertor_discretisation, snap_to_horizontal_angle + ) + + def _create_pre_cell_stage( + self, blanket_discretisation, divertor_discretisation, snap_to_horizontal_angle + ): + cutting = CuttingStage( + blanket=PanelsAndExteriorCurve( + self.geom.panel_break_points, + self.geom.boundary, + self.geom.vacuum_vessel_wire, + ), + divertor=DivertorWireAndExteriorCurve( + self.geom.divertor_wire, self.geom.boundary, self.geom.vacuum_vessel_wire + ), + ) + divertor = cutting.divertor.make_divertor_pre_cell_array( + discretisation_level=divertor_discretisation + ) + first, last = divertor.exterior_vertices()[(0, -1),] + + blanket = cutting.blanket.make_quadrilateral_pre_cell_array( + discretisation_level=blanket_discretisation, + starting_cut=first[::2], + ending_cut=last[::2], + snap_to_horizontal_angle=snap_to_horizontal_angle, + ) + + return PreCellStage( + blanket=blanket.straighten_exterior(preserve_volume=True), divertor=divertor + ) + + @property + def bounding_box(self) -> tuple[float, ...]: + """Bounding box of Neutronics reactor""" + return self._pre_cell_stage.bounding_box() + + @property + def half_bounding_box(self) -> tuple[float, ...]: + """Bounding box of the right-hand half of the 2D poloidal cross-section""" + return self._pre_cell_stage.half_bounding_box() + + @property + def blanket(self): + """Blanket pre cell""" + return self._pre_cell_stage.blanket + + @property + def divertor(self): + """Divertor pre cell""" + return self._pre_cell_stage.divertor + + def plot_2d(self, *args, **kwargs): + """Plot neutronics reactor 2d profile""" + show = kwargs.pop("show", True) + ax = kwargs.pop("ax", None) + ax = self.blanket.plot_2d(*args, ax=ax, show=False, **kwargs) + return self.divertor.plot_2d(*args, ax=ax, show=show, **kwargs) + + @abstractmethod + def _get_wires_from_components( + self, + divertor: ComponentManager, + blanket: ComponentManager, + vacuum_vessel: ComponentManager, + ) -> tuple[ + TokamakDimensions, BluemiraWire, npt.NDArray, BluemiraWire, BluemiraWire + ]: ... diff --git a/bluemira/radiation_transport/neutronics/radial_wall.py b/bluemira/radiation_transport/neutronics/radial_wall.py new file mode 100644 index 0000000000..bc21a599e4 --- /dev/null +++ b/bluemira/radiation_transport/neutronics/radial_wall.py @@ -0,0 +1,320 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Defining (and changing) the radial (side) walls of PreCell in PreCellArrays.""" + +from __future__ import annotations + +from enum import IntEnum +from typing import TYPE_CHECKING + +import numpy as np +from numpy import typing as npt +from scipy.optimize import newton + +from bluemira.base.look_and_feel import bluemira_debug, bluemira_warn +from bluemira.geometry.constants import EPS_FREECAD +from bluemira.geometry.error import GeometryError +from bluemira.geometry.tools import partial_diff_of_volume, polygon_revolve_signed_volume + +if TYPE_CHECKING: + from collections.abc import Iterable + + from bluemira.radiation_transport.neutronics.make_pre_cell import PreCellArray + + +class Vert(IntEnum): + """Vertices index for cells""" + + exterior_end = 0 + interior_start = 1 + interior_end = 2 + exterior_start = 3 + + +class CellWalls: + """ + A list of start- and end-location vectors of all of the walls dividing neighbouring + pre-cells. + """ + + def __init__(self, cell_walls: npt.ArrayLike): + """ + Parameters + ---------- + cell_walls + array of shape (N+1, 2, 2) + where N = number of cells, so the axis=0 describes the N+1 walls on either + side of each cell, the axis=1 describes the the start and end points, and + axis=2 describes the r and z coordinates. + """ + self.cell_walls = np.array(cell_walls) + if np.shape(self.cell_walls)[1:] != (2, 2): + raise ValueError( + "Expected N values of start and end xz coordinates, i.e. " + f"shape = (N+1, 2, 2); got {np.shape(self.cell_walls)}." + ) + self._starts = self.cell_walls[:, 0] # shape = (N+1, 2) + self._init_ends = self.cell_walls[:, 1] # shape = (N+1, 2) + vector = self._init_ends - self._starts # shape = (N+1, 2) + self.original_lengths = np.linalg.norm(vector, axis=-1) + self.directions = (vector.T / self.original_lengths).T + self.num_cells: int = len(self) - 1 + self.check_volumes_and_lengths() + + def __len__(self) -> int: + """Number of cell wall panels""" + return len(self.cell_walls) + + def __getitem__(self, index_or_slice) -> npt.NDArray | float: + """Get cell wall panel""" + return self.cell_walls[index_or_slice] + + def __setitem__(self, index_or_slice, new_coordinates: npt.NDArray | float): + """ + self[:, :, :] = ... can completely reset some coordinates. + However, a full-reset should be avoided because we don't want to mess with the + start rz coordinates. + """ + self.cell_walls[index_or_slice] = new_coordinates + + def __repr__(self) -> str: + """String representation""" + return ( + super().__repr__().replace(" at ", f" of {len(self.cell_walls)} walls at ") + ) + + def copy(self) -> CellWalls: + """Copy cell wall""" + return CellWalls(self.cell_walls.copy()) + + @classmethod + def from_pre_cell_array(cls, pre_cell_array: PreCellArray) -> CellWalls: + """Use the corner vertices in an array of pre-cells to make a CellWalls.""" + # cut each coordinates down from having shape (3, 1) down to (2,) + return cls([ + *( + (c.vertex[:, Vert.interior_end], c.vertex[:, Vert.exterior_start]) + for c in pre_cell_array + ), + ( + pre_cell_array[-1].vertex[:, Vert.interior_start], + pre_cell_array[-1].vertex[:, Vert.exterior_end], + ), + ]) + + @classmethod + def from_pre_cell_array_vv(cls, pre_cell_array: PreCellArray) -> CellWalls: + """ + Use the corner vertices and the vacuum vessel vertices of the pre-cell array to + make a CellWall. + """ + return cls([ + *( + (c.vertex[:, Vert.interior_end], c.vv_point[:, 0]) + for c in pre_cell_array + ), + ( + pre_cell_array[-1].vertex[:, Vert.interior_start], + pre_cell_array[-1].vv_point[:, 1], + ), + ]) + + @property + def starts(self) -> npt.NDArray: + """The start points of each cell wall. shape = (N+1, 2)""" + return self._starts # the start points never change value. + + @property + def ends(self) -> npt.NDArray: + """The end point changes value depending on the user-set length.""" + return self.cell_walls[:, 1] # shape = (N+1, 2) + + def calculate_new_end_points( + self, lengths: float | npt.NDArray[np.float64] + ) -> npt.NDArray: + """ + Get the end points of each cell wall if they were changed to have the specified + lengths. This is different to set_length in that the new end points are returned, + and the object itself is not modified by the test length(s). + + Parameters + ---------- + lengths + np.ndarray of shape = (N+1,) where N+1 == number of cell walls. + It can also be a scalar, which would then be broadcasted into (N+1,). + + Returns + ------- + new end points + an array of the same shape as self.starts + """ + return self.starts + (self.directions.T * lengths).T + + def get_length(self, i: int) -> float: + """ + Get the length of the i-th cell-wall + + Returns + ------- + length: float + """ + _end_i, _start_i = self.cell_walls[i] + return np.linalg.norm(_end_i - _start_i) + + def set_length(self, i, new_length): + """Set the length of the i-th cell-wall""" + self.cell_walls[i, 1] = self.cell_walls[i, 0] + self.directions[i] * new_length + + @property + def lengths(self): + """Current lengths of the cell walls.""" + return np.linalg.norm(self.ends - self.starts, axis=-1) # shape = (N+1) + + def get_volume(self, i): + """ + Get the volume of the i-ith cell + + Returns + ------- + volume: float + """ + return polygon_revolve_signed_volume( + np.concatenate([self.cell_walls[i], self.cell_walls[i + 1][::-1]]) + ) + + @property + def volumes(self): + """ + Current volumes of the (simplified) cells created by joining straight lines + between neighbouring cell walls. + """ + # shape = (N+1,) + return np.asarray([self.get_volume(i) for i in range(self.num_cells)]) + + def check_volumes_and_lengths(self): + """ + Ensure all cells have positive volumes, to minimise the risk of self-intersecting + lines and negative lengths + """ + if not (all(self.volumes > 0) and all(self.lengths > 0)): + raise GeometryError("At least 1 cell has non-positive volumes!") + + def volume_of_cells_neighbouring(self, i, test_length): + """ + Get the volume of cell[i-1] and cell[i] when cell_wall[i] is set to the + test_length. + + Returns + ------- + volume: float + """ + start_i, dir_i = self.starts[i], self.directions[i] + new_end = start_i + dir_i * test_length + prev_wall, next_wall = self.cell_walls[i - 1 : i + 2 : 2] + return polygon_revolve_signed_volume([ + prev_wall[0], + prev_wall[1], + new_end, + start_i, + ]) + polygon_revolve_signed_volume([ + start_i, + new_end, + next_wall[1], + next_wall[0], + ]) + + def volume_derivative_of_cells_neighbouring(self, i, test_length): + """ + Measure the derivative on the volume of cell[i-1] and cell[i] w.r.t. to + length of cell_wall[i], at cell_wall[i] length = test_length. + + Returns + ------- + dV/dl: float + """ + start_i, dir_i = self.starts[i], self.directions[i] + new_end = start_i + dir_i * test_length + prev_end, next_end = self.ends[i - 1 : i + 2 : 2] + return partial_diff_of_volume( + [prev_end, new_end, start_i], dir_i + ) + partial_diff_of_volume([start_i, new_end, next_end], dir_i) + + def optimise_to_match_individual_volumes( + self, volume_list: Iterable[float], *, max_iter=1000 + ): + """ + Allow the lengths of each wall to increase, so that the overall volumes are + preserved as much as possible. Assuming the entire exterior curve is convex, + then our linear approximation is only going to under-approximate. Therefore + to achieve better approximation, we only need to increase the lengths. + """ + if self.num_cells <= 1: + return + + target_volumes = np.array(list(volume_list)) + if self.num_cells == 2: # noqa: PLR2004 + # only one single step is required for the optimisation + def volume_excess(new_length): + return self.volume_of_cells_neighbouring(1, new_length) - sum( + target_volumes + ) + + def derivative(new_length): + return self.volume_derivative_of_cells_neighbouring(1, new_length) + + self.set_length(1, newton(volume_excess, self.get_length(1), derivative)) + self.check_volumes_and_lengths() + return + + # if more than 3 walls (more than 2 cells) + step_direction, i, i_min, i_max = 1, 1, 1, self.num_cells - 1 + + num_passes_counter = -1 + forward_pass_result = np.zeros(self.num_cells + 1) + + while num_passes_counter < max_iter: + # do not allow length to decrease beyond their original value. + + def excess_volume(test_length, i=i): + return self.volume_of_cells_neighbouring(i, test_length) - sum( + target_volumes[i - 1 : i + 1] + ) + + def dV_dl(test_length, i=i): # noqa: N802 + return self.volume_derivative_of_cells_neighbouring(i, test_length) + + self.set_length( + i, + newton(excess_volume, self.get_length(i), dV_dl) + if excess_volume(self.original_lengths[i]) < 0 + else self.original_lengths[i], + ) + if i == i_min: + # hitting the left end: bounce to the right + step_direction = 1 + num_passes_counter += 1 + backward_pass_result = self.lengths.copy() + # termination condition + if np.allclose( + backward_pass_result, forward_pass_result, rtol=0, atol=EPS_FREECAD + ): + bluemira_debug( + "Cell volume-matching optimisation successful." + "Terminating iterative cell wall length adjustment after " + f"{num_passes_counter} passes." + ) + self.check_volumes_and_lengths() + return + elif i == i_max: + # hitting the right end: bounce to the left + step_direction = -1 + num_passes_counter += 1 + forward_pass_result = self.lengths.copy() + i += step_direction + bluemira_warn( + "Optimisation failed: Did not converge within" + f"{num_passes_counter} iterations!" + ) diff --git a/bluemira/radiation_transport/neutronics/slicing.py b/bluemira/radiation_transport/neutronics/slicing.py new file mode 100644 index 0000000000..0348fd289b --- /dev/null +++ b/bluemira/radiation_transport/neutronics/slicing.py @@ -0,0 +1,937 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Oversees the conversion from bluemira wires into pre-cells, then into csg.""" + +from __future__ import annotations + +from itertools import chain, pairwise, starmap +from typing import TYPE_CHECKING + +import numpy as np + +from bluemira.base.look_and_feel import bluemira_warn +from bluemira.codes import _freecadapi as cadapi +from bluemira.geometry.constants import EPS_FREECAD +from bluemira.geometry.coordinates import choose_direction, get_bisection_line +from bluemira.geometry.error import GeometryError +from bluemira.geometry.plane import ( + BluemiraPlane, + calculate_plane_dir, + x_plane, + xz_plane_from_2_points, + z_plane, +) +from bluemira.geometry.tools import get_wire_plane_intersect, make_polygon +from bluemira.radiation_transport.neutronics.constants import ( + DISCRETISATION_LEVEL, + TOLERANCE_DEGREES, +) +from bluemira.radiation_transport.neutronics.make_pre_cell import ( + DivertorPreCell, + DivertorPreCellArray, + PreCell, + PreCellArray, +) +from bluemira.radiation_transport.neutronics.wires import ( + CircleInfo, + StraightLineInfo, + WireInfo, + WireInfoList, +) + +if TYPE_CHECKING: + from collections.abc import Iterator, Sequence + + import numpy.typing as npt + + from bluemira.geometry.wire import BluemiraWire + + +def cut_curve( + cut_points: list[npt.NDArray], + wire: BluemiraWire, + discretisation_level: int, + *, + reverse: bool = False, +) -> Iterator[npt.NDArray[np.float64]]: + """ + Generator to cut a closed BluemiraWire to size. + Current implementation is to yield a parameter range for every sequential cut-point + pair. Subject to change of implementation in the future when issue 3038 is fixed. + + Parameters + ---------- + cut_points: + A list of at least 3 points, presumed to be lying on the wire, that we want to + cut the wire at. + wire: + The wire that we want to cut up. + discretisation_level: + We yield a list of points of len==discretisation_level, such that we can build a + wire made of (discretisation_level-1) straight lines to approximate each segment. + reversed: + Whether we want the neutron spectrum to go in the increasing or the decreasing + direction. + + Yields + ------ + params_range: + For each segment, yield the values of the parametric variable t such that + `[wire.value_at(t) for t in params_range]` + is a list of points (sufficient to built up a series of (discretisation_level-1) + straight lines that approximate that segment) + """ + cut_params = [wire.parameter_at(cp) for cp in cut_points] + # for reference: cut_params has range: + # [0, 1] for an open wire, [0, 1) for a closed wire. + + # determine whether we're going up or down in parameter (t) space. + finite_difference = np.diff(cut_params) + if len(finite_difference) <= 2: # noqa: PLR2004 + raise GeometryError( + "Too few points! I.e. discretisation_level parameter too low. " + "Can't determine the cut direction!" + ) + if (finite_difference <= 0).sum() <= 1: + # strictly monotonically increasing except for 1 wrap-around point + increasing = True + elif (finite_difference >= 0).sum() <= 1: + # strictly monotonically decreasing except for 1 wrap-around point + increasing = False + else: # no discernible pattern in the increase/decrease + raise GeometryError("Points are too disordered!") + + # generator function + for alpha, beta in pairwise(cut_params): + used_alpha = alpha + if increasing: + if alpha > beta: # alpha is expected to be smaller than beta + used_alpha -= 1.0 + elif alpha < beta: # alpha is expected to be larger than beta. + used_alpha += 1.0 + param_range = np.linspace(used_alpha, beta, discretisation_level) % 1.0 + yield param_range[::-1] if reverse else param_range + + +def check_and_breakdown_wire(wire: BluemiraWire) -> WireInfoList: + """ + Raise GeometryError if the BluemiraWire has an unexpected data storage structure. + Then, get only the key information (start/end points and tangent) of each segment of + the wire. + """ + wire_container = [] + + def add_line( + edge: cadapi.apiEdge, + wire: BluemiraWire, + start_vec: cadapi.apiVector | npt.NDArray, + end_vec: cadapi.apiVector | npt.NDArray, + ): + """Function to record a line""" + return WireInfo( + StraightLineInfo(np.array(start_vec), np.array(end_vec)), + [edge.tangentAt(edge.FirstParameter), edge.tangentAt(edge.LastParameter)], + wire, + ) + + def add_circle( + edge: cadapi.apiEdge, + wire: BluemiraWire, + start_vec: cadapi.apiVector | npt.NDArray, + end_vec: cadapi.apiVector | npt.NDArray, + ): + """Function to record the arc of a circle.""" + return WireInfo( + CircleInfo( + np.array(start_vec), + np.array(end_vec), + np.array(edge.Curve.Center), + edge.Curve.Radius, + ), + [edge.tangentAt(edge.FirstParameter), edge.tangentAt(edge.LastParameter)], + wire, + ) + + for w_edge in wire.edges: + if len(w_edge.boundary) != 1 or len(w_edge.boundary[0].OrderedEdges) != 1: + raise GeometryError("Expected each boundary to contain only 1 curve!") + edge = w_edge.boundary[0].OrderedEdges[0] + + # Create aliases for easier referring to variables. + # The following line may become `edge.start_point(), edge.end_point()` + # when PR # 3095 is merged + current_start, current_end = edge.firstVertex().Point, edge.lastVertex().Point + + # Get the info about this segment of wire + if isinstance(edge.Curve, cadapi.Part.BSplineCurve | cadapi.Part.BezierCurve): + wire_container.extend( + check_and_breakdown_wire( + make_polygon(w_edge.discretise(DISCRETISATION_LEVEL), closed=False) + ).info_list + ) + continue + + if isinstance(edge.Curve, cadapi.Part.Line | cadapi.Part.LineSegment): + wire_info = add_line(edge, w_edge, current_start, current_end) + + elif isinstance(edge.Curve, cadapi.Part.ArcOfCircle | cadapi.Part.Circle): + wire_info = add_circle(edge, w_edge, current_start, current_end) + elif isinstance(edge.Curve, cadapi.Part.ArcOfEllipse | cadapi.Part.Ellipse): + raise NotImplementedError("Conversion for ellipses are not available yet.") + # TODO: implement this feature + else: + raise NotImplementedError(f"Conversion for {edge.Curve} not available yet.") + + # horrible hack to work around issue #3037, which is still an open issue: + # wire segments are sometimes reversed for no apparent reason. + if len(wire_container) == 0: + wire_container.append(wire_info) + continue + distance_to_start = np.linalg.norm( + wire_container[-1].key_points.end_point - wire_info.key_points.start_point + ) + distance_to_end = np.linalg.norm( + wire_container[-1].key_points.end_point - wire_info.key_points.end_point + ) + if distance_to_end < distance_to_start: + wire_info = wire_info.reverse() + wire_container.append(wire_info) + + return WireInfoList(wire_container) + + +def turned_morethan_180( + xyz_vector1: Sequence[float], xyz_vector2: Sequence[float], direction_sign: int +) -> bool: + """ + Checked if one needs to rotate vector 1 by more than 180° in the specified direction + by more than 180° to align with vector 2. + + Parameters + ---------- + xyz_vector1, xyz_vector2: Sequence[float] + xyz array of where the vector points. + direction_sign: signed integer + +1: evaluate rotation required in the counter-clockwise direction. + -1: evaluate rotation required in the clockwise direction. + """ + if xyz_vector1[1] != 0 or xyz_vector2[1] != 0: + raise GeometryError("Tangent vector points out of plane!") + angle1 = np.arctan2(xyz_vector1[2], xyz_vector1[0]) + angle2 = np.arctan2(xyz_vector2[2], xyz_vector2[0]) + if direction_sign == 1: + if angle2 < angle1: + angle2 += 2 * np.pi + return (angle2 - angle1) >= np.pi + if angle2 > angle1: + angle2 -= 2 * np.pi + return (angle1 - angle2) >= np.pi + + +def deviate_less_than( + xyz_vector1: Sequence[float], xyz_vector2: Sequence[float], threshold_degrees: float +) -> bool: + """Check if two vector's angles less than a certain threshold angle (in degrees).""" + angle1 = np.arctan2(xyz_vector1[2], xyz_vector1[0]) + angle2 = np.arctan2(xyz_vector2[2], xyz_vector2[0]) + return np.rad2deg(abs(angle2 - angle1)) < threshold_degrees + + +def straight_lines_deviate_less_than( + info1: WireInfo, info2: WireInfo, threshold_degrees: float +) -> bool: + """ + Check that both lines are straight lines, then check if deviation is less than + threshold_degrees or not + """ + if not ( + isinstance(info1.key_points, StraightLineInfo) + and isinstance(info2.key_points, StraightLineInfo) + ): + return False + return deviate_less_than(info1.tangents[1], info2.tangents[0], threshold_degrees) + + +def break_wire_into_convex_chunks( + wire: BluemiraWire, curvature_sign: int = -1 +) -> list[WireInfoList]: + """ + Break a wire up into several convex wires. + Merge if they are almost collinear. + + Parameters + ---------- + wire: + curvature_sign: + if it's -1: we allow each convex chunk's wire to turn right only. + if it's 1: we allow each convex chunk's wire to turn left only. + + Returns + ------- + convex_chunks: + a list of WireInfos + """ + wire_segments = list(check_and_breakdown_wire(wire)) + # initialising the first chunk + convex_chunks, this_chunk = [], [] + chunk_start_tangent = wire_segments[0].tangents[0] + + def add_to_chunk(this_seg: WireInfo) -> None: + """Add the info and wire into the current chunk and current wire.""" + if this_chunk and straight_lines_deviate_less_than( + this_chunk[-1], this_seg, TOLERANCE_DEGREES + ): + # modify the previous line directly + this_chunk[-1].key_points = StraightLineInfo( + this_chunk[-1].key_points.start_point, this_seg.key_points.end_point + ) + return + this_chunk.append(WireInfo(this_seg.key_points, this_seg.tangents)) + + for no, w_s in enumerate(wire_segments[1:]): + add_to_chunk(wire_segments[no]) + prev_end_tangent = this_chunk[-1].tangents[-1] + next_start_tangent = w_s.tangents[0] + if deviate_less_than( + this_chunk[-1].tangents[1], w_s.tangents[0], TOLERANCE_DEGREES + ): + continue + interior_curve_turned_over_180 = turned_morethan_180( + chunk_start_tangent, next_start_tangent, curvature_sign + ) + if turned_morethan_180(prev_end_tangent, next_start_tangent, curvature_sign): + convex_chunks.append(WireInfoList(this_chunk.copy())) + this_chunk.clear() + chunk_start_tangent = w_s.tangents[0] + elif interior_curve_turned_over_180: + # curled in on itself too much. + bluemira_warn( + "Divertor wire geometry possibly too extreme for program " + "to handle. Check pre-cell visually by using the .plot_2d() methods " + "on the relevant DivertorPreCell and DivertorPreCellArray." + ) + + add_to_chunk(wire_segments[-1]) + convex_chunks.append(WireInfoList(this_chunk.copy())) + + return convex_chunks + + +class PanelsAndExteriorCurve: + """ + A collection of three objects, the first wall panels, the vacuum vessel interior + curve, and the exterior curve. + + Parameters + ---------- + panel_break_points: + numpy array of shape==(N, 2), showing the RZ coordinates of joining points + between adjacent first wall panels, running in the clockwise direction, + from the inboard divertor to the top, then down to the outboard divertor. + vv_interior: + The BluemiraWire representing the inside surface of the vacuum vessel. + vv_exterior: + The BluemiraWire representing the outside surface of the vacuum vessel. + """ + + def __init__( + self, + panel_break_points: npt.NDArray[np.float64], + vv_interior: BluemiraWire, + vv_exterior: BluemiraWire, + ): + """Instantiate from a break point curve and the vv_interior. + + Parameters + ---------- + panel_break_points: + A series of 2D coordinate (of shape = (N+1, 3)) representing the N panels + of the blanket. It (also) runs clockwise (inboard side to outboard side), + same as vv_interior + vv_interior: + A BluemiraWire that runs clockwise, showing the vacuum vessel on the RHHP + cross-section of the tokamak. + """ + self.vv_interior = vv_interior + self.vv_exterior = vv_exterior + # shape = (N+1, 3) + self.interior_panels = panel_break_points + if len(self.interior_panels[0]) != 3 or np.ndim(self.interior_panels) != 2: # noqa: PLR2004 + raise ValueError( + "Expected an input np.ndarray of breakpoints of shape = " + f"(N+1, 2). Instead received shape = {np.shape(panel_break_points)}." + ) + + def get_bisection_line( + self, index: int + ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """Calculate the bisection line that separates two panels at breakpoint[i]. + + Parameters + ---------- + index: + the n-th breakpoint that we want the bisection line to project from. + (N.B. There are N+1 breakpoints for N panels.) + Thus this is the last point of the n-th panel. + """ + p1, p2, p3 = self.interior_panels[index - 1 : index + 2, ::2] + origin_2d, direction_2d = get_bisection_line(p1, p2, p3, p2) + line_origin = np.insert(origin_2d, 1, 0, axis=-1) + line_direction = np.insert(direction_2d, 1, 0, axis=-1) + return line_origin, line_direction + + def calculate_cut_points( + self, + starting_cut: npt.NDArray[np.float64] | None, + ending_cut: npt.NDArray[np.float64] | None, + snap_to_horizontal_angle: float, + ) -> tuple[list[npt.NDArray[np.float64]], ...]: + """ + Cut the curves according to some specified criteria: + In general, a line would be drawn from each panel break point outwards towards + the curves. This would form the cut line. + + The space between two neighbouring cut line makes a PreCell. + + Usually these cut's angle is the bisection of the normal angles of its + neighbouring panels, but in special rules applies when snap_to_horizontal_angle + is set to a value >0. + + Since there is only one neighbouring panel at the start break point and end break + point, starting_cut and end_cut must be specified. If not, they are chosen to be + horizontal cuts. + + Parameters + ---------- + snap_to_horizontal_angle: + If the cutting plane is less than x degrees (°) away from horizontal, + then snap it to horizontal. + allowed range: [0, 90] + starting_cut, ending_cut: + Each is an ndarray with shape (2,), denoting the destination of the cut + lines. + + The program cannot deduce what angle to cut the curves at these + locations without extra user input. Therefore the user is able to use these + options to specify the cut line's destination point for the first and final + points respectively. + + For the first cut line, the cut line would start from interior_panels[0] + and reach starting_cut. + For the final cut line, the cut line would start from interior_panels[-1] + and reach ending_cut. + Both arguments have shape (2,) if given, representing the RZ coordinates. + + + Returns + ------- + vv_cut_points: + cut points where the vacuum vessel interior wire is split + exterior_cut_points: + cut points where the vacuum vessel exterior wire is split + """ + vv_cut_points, exterior_cut_points = [], [] + + def add_cut_points(cutting_plane: BluemiraPlane, cut_direction: npt.NDArray): + """ + Find where the cutting_plane intersect the self.vv_interior and + self.vv_exterior; these points will eventually be used to cut up + self.vv_interior and self.vv_exterior. + + N.B. These cut points must be sequentially added, i.e. they should + follow the curves in the clockwise direction. + """ + vv_cut_points.append( + get_wire_plane_intersect(self.vv_interior, cutting_plane, cut_direction) + ) + exterior_cut_points.append( + get_wire_plane_intersect(self.vv_exterior, cutting_plane, cut_direction) + ) + + threshold_angle = np.deg2rad(snap_to_horizontal_angle) + + # initial cut point + add_cut_points( + *calculate_plane_dir( + self.interior_panels[0], [starting_cut[0], 0, starting_cut[-1]] + ) + ) + + for i in range(1, len(self.interior_panels) - 1): + origin, c_dir = self.get_bisection_line(i) + + if c_dir[0] == 0: + plane = x_plane(self.interior_panels[i][0]) # vertical cut plane + elif abs(np.arctan(c_dir[-1] / c_dir[0])) < threshold_angle: + plane = z_plane(self.interior_panels[i][-1]) # horizontal cut plane + else: + plane = xz_plane_from_2_points(origin, origin + c_dir) + add_cut_points(plane, c_dir) + + # final cut point + add_cut_points( + *calculate_plane_dir( + self.interior_panels[-1], [ending_cut[0], 0, ending_cut[-1]] + ) + ) + + return vv_cut_points, exterior_cut_points + + def execute_curve_cut( + self, + discretisation_level: int, + starting_cut: npt.NDArray[np.float64] | None = None, + ending_cut: npt.NDArray[np.float64] | None = None, + snap_to_horizontal_angle: float = 30.0, + ) -> tuple[list[BluemiraWire], ...]: + """ + Cut the exterior curve into a series and return them. + This is the slowest part of the entire csg-creation process, because of the + discretisation. + + Parameters + ---------- + snap_to_horizontal_angle, starting_cut, ending_cut: + See + :meth:`~bluemira.radiation_transport.neutronics.slicing.PanelsAndExteriorCurve.calculate_cut_points` + discretisation_level: + how many points to use to approximate the curve. + + Returns + ------- + vv_curve_segments: + segments of the vacuum vessel interior curve forming each pre-cell's interior + curve. + exterior_curve_segments: + segments of the vacuum vessel exterior curve forming each pre-cell's exterior + curve. + """ + # TODO: remove discretisation level when issue #3038 is fixed. The raw wire can + # be used without discretisation then. + vv_cut_points, exterior_cut_points = self.calculate_cut_points( + starting_cut, ending_cut, snap_to_horizontal_angle + ) + + vv_curve_segments = [ + self.approximate_curve( + self.vv_interior, + t_range, + label=f"vacuum vessel interior curve {i + 1}", + ) + for i, t_range in enumerate( + cut_curve( + vv_cut_points, self.vv_interior, discretisation_level, reverse=False + ) + ) + ] + + exterior_curve_segments = [ + self.approximate_curve( + self.vv_exterior, + t_range, + label=f"vacuum vessel exterior curve {i + 1}", + ) + for i, t_range in enumerate( + cut_curve( + exterior_cut_points, + self.vv_exterior, + discretisation_level, + reverse=False, + ) + ) + ] + + return vv_curve_segments, exterior_curve_segments + + @staticmethod + def approximate_curve( + curve: BluemiraWire, param_range: npt.NDArray[np.float64], label="" + ) -> BluemiraWire: + """ + Given params_range (an iterable of floats of len n between 0 to 1), create n-1 + straight lines such that the n-th point corresponds to + `curve.value_at(params_range[n])`. + + This implementation shall be replaced when issue #3038 gets resolved. + """ + return make_polygon( + np.array([curve.value_at(t) for t in param_range]).T, + label=label, + closed=False, # we expect it to form an open curve. + ) + + def make_quadrilateral_pre_cell_array( + self, + discretisation_level: int = DISCRETISATION_LEVEL, + starting_cut: npt.NDArray[np.float64] | None = None, + ending_cut: npt.NDArray[np.float64] | None = None, + snap_to_horizontal_angle: float = 30.0, + ) -> PreCellArray: + """ + Cut the exterior curve up, so that it would act as the exterior side of the + quadrilateral pre-cell. Then, the panel would act as the interior side of the + pre-cell. The two remaining sides are the counter-clockwise and clockwise walls. + + The vacuum vessel interior surface is sandwiched between them. + """ + # make horizontal cuts if the starting and ending cuts aren't provided. + if starting_cut is None: + start_r = 0 # starting cut's final destination's r value + starting_cut = np.array([start_r, self.interior_panels[0][-1]]) + if ending_cut is None: + # ending cut's final destination's r value + end_r = self.interior_panels[-1, 0] * 2 + ending_cut = np.array([end_r, self.interior_panels[-1][-1]]) + + return PreCellArray([ + PreCell( + make_polygon(self.interior_panels[i : i + 2][::-1].T, closed=False), + vv_segment, + exterior_segment, + ) + for i, (vv_segment, exterior_segment) in enumerate( + zip( + *self.execute_curve_cut( + discretisation_level, + starting_cut, + ending_cut, + snap_to_horizontal_angle, + ), + strict=True, + ) + ) + ]) + + +class DivertorWireAndExteriorCurve: + """ + A class to store a wire with a vacuum vessel interior curve and an exterior curve. + """ + + def __init__( + self, + divertor_wire: BluemiraWire, + vv_interior: BluemiraWire, + vv_exterior: BluemiraWire, + ): + """ + Instantiate from a BluemiraWire of the divertor and a BluemiraWire of the + exterior (outside of vacuum vessem) curve. Also save the BluemiraWire that goes + between these two (inside of vacuum vessel). + + Parameters + ---------- + divertor_wire + A BluemiraWire that runs from the inboard side to the outboard side of the + tokamak, representing the plasma facing surface of the divertor. + vv_interior + A BluemiraWire that runs clockwise, showing the inside of the vacuum vessel + on the RHHP cross-section of the tokamak. + vv_exterior + A BluemiraWire that runs clockwise, showing the outside of the vacuum vessel + on the RHHP cross-section of the tokamak. + """ + self.convex_segments = break_wire_into_convex_chunks(divertor_wire) + self.key_points = np.array([ # shape = (N+1, 3) + *(seg.key_points[0] for seg in chain.from_iterable(self.convex_segments)), + self.convex_segments[-1][-1].key_points[1], + ]) + self.tangents = [ # shape = (N, 2, 3) + seg.tangents for seg in chain.from_iterable(self.convex_segments) + ] + self.center_point = np.array([ # center of the divertor + (self.key_points[0, 0] + self.key_points[-1, 0]) / 2, # mean x + self.key_points[:, -1].max(), # highest z + ]) + self.vv_interior = vv_interior + self.vv_exterior = vv_exterior + + def get_bisection_line( + self, prev_idx: int + ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """ + Get the the line bisecting the x and the y, represented as origin and direction. + + Parameters + ---------- + prev_idx + two convex chunks are involved when making a bisection line. This index + refers to the smaller of these two chunks' indices. + + Returns + ------- + line_origin: + a point on this line + line_direction: + a normal vector pointing in the direction of this line. + """ + # Get the vectors as arrays of shape = (2,) + anchor1 = np.array(self.convex_segments[prev_idx][-1].key_points[1])[::2] + anchor2 = np.array(self.convex_segments[prev_idx + 1][0].key_points[0])[::2] + # force the directions to point outwards. + direct1 = choose_direction( + np.array(self.convex_segments[prev_idx][-1].tangents[1])[::2], + self.center_point, + anchor1, + ) + direct2 = choose_direction( + np.array(self.convex_segments[prev_idx + 1][0].tangents[0])[::2], + self.center_point, + anchor2, + ) + origin_2d, direction_2d = get_bisection_line( + anchor1 - direct1, anchor1, anchor2 - direct2, anchor2 + ) + return np.insert(origin_2d, 1, 0, axis=-1), np.insert( + direction_2d, 1, 0, axis=-1 + ) + + def calculate_cut_points( + self, + starting_cut: npt.NDArray[np.float64] | None, + ending_cut: npt.NDArray[np.float64] | None, + ) -> tuple[list[npt.NDArray[np.float64]], ...]: + """ + Cut the curves up into N segments to match the N convex chunks of the + divertor. + The space between two neighbouring cut line makes a DivertorPreCell. + + Parameters + ---------- + starting_cut, ending_cut: + Each is an ndarray with shape (2,), denoting the destination of the cut + lines. + + The program cannot deduce what angle to cut the curves at these + locations without extra user input. Therefore the user is able to use these + options to specify the cut line's destination point for the first and final + points respectively. + + For the first cut line, the cut line would start from interior_panels[0] and + reach starting_cut. + For the final cut line, the cut line would start from interior_panels[-1] and + reach ending_cut. + Both arguments have shape (2,) if given, representing the RZ coordinates. + + + Returns + ------- + self.vv_cut_points + cut points where the vacuum vessel interior wire is split + self.exterior_cut_points + cut points where the vacuum vessel exterior wire is split + """ + vv_cut_points, exterior_cut_points = [], [] + + def add_cut_points(cutting_plane: BluemiraPlane, cut_direction: npt.NDArray): + """ + Find where the cutting_plane intersect the self.vv_interior and + self.vv_exterior; these points will eventually be used to cut up + self.vv_interior and self.vv_exterior. + + N.B. These cut points must be sequentially added, i.e. they should + follow the curves in the clockwise direction. + + While identical to :meth:`~PanelsAndExteriorCurve.add_cut_points`, + this can't be refactored away because they're specific to the class. + """ + vv_cut_points.append( + get_wire_plane_intersect(self.vv_interior, cutting_plane, cut_direction) + ) + exterior_cut_points.append( + get_wire_plane_intersect(self.vv_exterior, cutting_plane, cut_direction) + ) + + # initial cut point + add_cut_points( + *calculate_plane_dir( + np.array(self.convex_segments[0][0].key_points[0]), + [starting_cut[0], 0, starting_cut[-1]], + ) + ) + + for i in range(len(self.convex_segments) - 1): + origin, c_dir = self.get_bisection_line(i) + plane = xz_plane_from_2_points(origin, origin + c_dir) + add_cut_points(plane, c_dir) + + # final cut point + add_cut_points( + *calculate_plane_dir( + np.array(self.convex_segments[-1][-1].key_points[1]), + [ending_cut[0], 0, ending_cut[-1]], + ) + ) + + return vv_cut_points, exterior_cut_points + + def execute_curve_cut( + self, + discretisation_level: int, + starting_cut: npt.NDArray[np.float64] | None, + ending_cut: npt.NDArray[np.float64] | None, + ) -> tuple[list[WireInfoList], ...]: + """ + Cut the vacuum vessel curves into a series return these segments. + + (Obsolete) + ---------- + Use the following table to decide whether the + segments run clockwise or counter-clockwise. + | | INCREASING = True | INCREASING = False | + ------------------------------------------------------------ + | reverse = False | clockwise | counter-clockwise | + | reverse = False | counter-clockwise | clockwise | + + This is the slowest part of the entire csg-creation process, because of the + discretisation. + + Parameters + ---------- + starting_cut, ending_cut: + See + :meth:`~bluemira.radiation_transport.neutronics.slicing.DivertorWireAndExteriorCurve.calculate_cut_points` + discretisation_level: + how many points to use to approximate the curve. + + Returns + ------- + vv_curve_segments + + list of WireInfoList describing the each pre-cell's vacuum vessel interior + curve + + exterior_curve_segments + + list of WireInfoList describing the each pre-cell's vacuum vessel exterior + curve + """ + # TODO: remove discretisation level when issue #3038 is fixed. The raw wire can + # be used without discretisation then. + vv_cut_points, exterior_cut_points = self.calculate_cut_points( + starting_cut, ending_cut + ) + + vv_curve_segments = [ + self.approximate_curve(self.vv_interior, t_range) + for t_range in cut_curve( + vv_cut_points, self.vv_interior, discretisation_level, reverse=True + ) + ] + + exterior_curve_segments = [ + self.approximate_curve(self.vv_exterior, t_range) + for t_range in cut_curve( + exterior_cut_points, self.vv_exterior, discretisation_level, reverse=True + ) + ] + + return vv_curve_segments, exterior_curve_segments + + @staticmethod + def approximate_curve( + curve: BluemiraWire, param_range: npt.NDArray[np.float64] + ) -> WireInfoList: + """ + Given params_range (an iterable of floats of len n between 0 to 1), create n-1 + straight lines such that the n-th point corresponds to + `curve.value_at(params_range[n])`. + + This implementation shall be updated/replaced when issue #3038 gets resolved. + """ + return WireInfoList( + starmap(WireInfo.from_2P, pairwise([curve.value_at(t) for t in param_range])) + ) + + def make_divertor_pre_cell_array( + self, + starting_cut: npt.NDArray[np.float64] | None = None, + ending_cut: npt.NDArray[np.float64] | None = None, + discretisation_level: int = DISCRETISATION_LEVEL, + ): + """ + Cut the exterior curve up, so that would act as the exterior side of the + quadrilateral pre-cell. Then, the panel would act as the interior side of the + pre-cell. The two remaining sides are the counter-clockwise and clockwise walls. + The vv_interior curve would be sandwiched between these two. + + Ordering: The cells should run from inoboard to outboard. + + Parameters + ---------- + starting_cut: + shape = (2,) + ending_cut: + shape = (2,) + discretisation_level: + integer: how many points to use to approximate each curve segment. + """ + # deduced starting_cut and ending_cut from the divertor itself. + if not starting_cut: + first_point = np.array(self.convex_segments[0][0].key_points[0])[::2] + tangent = np.array(self.convex_segments[0][0].tangents[0])[::2] + tangent = choose_direction(tangent, self.center_point, first_point) + starting_cut = first_point + tangent + if not ending_cut: + last_point = np.array(self.convex_segments[-1][-1].key_points[1])[::2] + tangent = np.array(self.convex_segments[-1][-1].tangents[1])[::2] + tangent = choose_direction(tangent, self.center_point, last_point) + ending_cut = last_point + tangent + + pre_cell_list = [] + for i, (vv_segment, exterior_segment) in enumerate( + zip( + *self.execute_curve_cut(discretisation_level, starting_cut, ending_cut), + strict=True, + ) + ): + interior_wire = self.convex_segments[i] + # make a new line joining the start and end. + # merge lines if collinear + cw_line = ( + WireInfoList([interior_wire.pop(0)]) + if np.allclose( + exterior_segment.end_point, + interior_wire.start_point, + rtol=0, + atol=EPS_FREECAD, + ) + else WireInfoList([ + WireInfo.from_2P( + exterior_segment.end_point, interior_wire.start_point + ) + ]) + ) + + while straight_lines_deviate_less_than(cw_line[-1], interior_wire[0], 0.5): + cw_line.end_point = interior_wire.start_point + interior_wire.pop(0) + + # merge lines if collinear + ccw_line = ( + WireInfoList([interior_wire.pop(-1)]) + if np.allclose( + interior_wire.end_point, + exterior_segment.start_point, + rtol=0, + atol=EPS_FREECAD, + ) + else WireInfoList([ + WireInfo.from_2P( + interior_wire.end_point, exterior_segment.start_point + ) + ]) + ) + + while straight_lines_deviate_less_than(interior_wire[-1], ccw_line[0], 0.5): + ccw_line.start_point = interior_wire.end_point + interior_wire.pop(-1) + + pre_cell_list.append( + DivertorPreCell(interior_wire, vv_segment, exterior_segment) + ) + + return DivertorPreCellArray(pre_cell_list) diff --git a/bluemira/radiation_transport/neutronics/wires.py b/bluemira/radiation_transport/neutronics/wires.py new file mode 100644 index 0000000000..390a5f12f8 --- /dev/null +++ b/bluemira/radiation_transport/neutronics/wires.py @@ -0,0 +1,182 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +""" +Info about straight line wires and circles. Made to be simpler to modify than a whole +BluemiraWire. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from itertools import pairwise +from typing import TYPE_CHECKING, NamedTuple + +import numpy as np +from numpy import typing as npt + +from bluemira.geometry.error import GeometryError +from bluemira.geometry.tools import make_circle_arc_3P, make_polygon +from bluemira.geometry.wire import BluemiraWire + +if TYPE_CHECKING: + from collections.abc import Iterable, Sequence + + +class StraightLineInfo(NamedTuple): + """Key information about a straight line""" + + start_point: Iterable[float] # 3D coordinates + end_point: Iterable[float] # 3D coordinates + + def reverse(self) -> StraightLineInfo: + """Flip the wire""" + return StraightLineInfo(self.end_point, self.start_point) + + +class CircleInfo(NamedTuple): + """Arc of a Circle, LESS THAN 180°""" + + start_point: Iterable[float] # 3D coordinates + end_point: Iterable[float] # 3D coordinates + center: Iterable[float] # 3D coordinates + radius: float # scalar + + def reverse(self) -> CircleInfo: + """Flip the wire""" + return CircleInfo(self.end_point, self.start_point, self.center, self.radius) + + +@dataclass +class WireInfo: + """ + A tuple to store: + 1. the key points about this wire (and what kind of wire this is) + 2. The tangent to that wire at the start and end + 3. A copy of the wire itself + """ + + # TODO: Perhaps implement more classes so it also work with splines? Or justify why + # we don't need to. Or merge this invention into an existing issue? + + key_points: StraightLineInfo | CircleInfo # 2 points of xyz/ CircleInfo + tangents: Sequence[Iterable[float]] # 2 normalised directional vectors xyz + wire: BluemiraWire | None = None + + def reverse(self) -> WireInfo: + """Flip the wire""" + return type(self)( + self.key_points.reverse(), [-t for t in self.tangents[::-1]], None + ) + + @classmethod + def from_2P( # noqa: N802 + cls, start_point: npt.NDArray[np.float64], end_point: npt.NDArray[np.float64] + ) -> WireInfo: + """ + Create the WireInfo for a straight line (i.e. one where the key_points is of + instance StraightLineInfo) using only two points. + """ + direction = np.array(end_point) - np.array(start_point) + normed_dir = np.array(direction) / np.linalg.norm(direction) + return cls( + StraightLineInfo(np.array(start_point), np.array(end_point)), + [normed_dir, normed_dir], + ) + + +class WireInfoList: + """A class to store info about a series of wires""" + + def __init__(self, info_list: Iterable[WireInfo]): + self.info_list = list(info_list) + for i, (prev_wire, curr_wire) in enumerate(pairwise(self.info_list)): + if not np.array_equal(prev_wire.key_points[1], curr_wire.key_points[0]): + raise GeometryError(f"wire {i + 1} must start where the wire {i} stops.") + + def __len__(self) -> int: + """Number of wire infos""" + return len(self.info_list) + + def __getitem__(self, index_or_slice) -> list[WireInfo] | WireInfo: + """Get a WireInfo""" + return self.info_list[index_or_slice] + + def __repr__(self) -> str: + """String representation""" + return super().__repr__().replace(" at ", f" of {len(self)} WireInfo at ") + + def pop(self, index): + """Pop one element""" + return self.info_list.pop(index) + + def get_3D_coordinates(self) -> npt.NDArray: + """ + Get of the entire wire. + """ + # assume continuity, which is already enforced during initialisation, so we + # should be fine. + # shape = (N+1, 3) + return np.array( + [ + self.info_list[0].key_points[0], + *(seg.key_points[1] for seg in self.info_list), + ], + dtype=float, + ) + + @property + def start_point(self): + """The start_point for the entire series of wires""" + return self.info_list[0].key_points[0] + + @start_point.setter + def start_point(self, new_start_point: npt.NDArray[np.float64]): + """ + Set the start_point to somewhere new. Note this doesn't change the tangents. + """ + old_kp = self.info_list[0].key_points + self.info_list[0].key_points = type(old_kp)(new_start_point, *old_kp[1:]) + + @property + def end_point(self): + """The end_point for the entire series of wires""" + return self.info_list[-1].key_points[1] + + @end_point.setter + def end_point(self, new_end_point): + """Set the end_point to somewhere new. Note this doesn't change the tangents.""" + old_kp = self.info_list[-1].key_points + self.info_list[0].key_points = type(old_kp)( + old_kp[0], new_end_point, *old_kp[2:] + ) + + def reverse(self) -> WireInfoList: + """Flip this list of wires""" + return WireInfoList([info.reverse() for info in self.info_list[::-1]]) + + def restore_to_wire(self) -> BluemiraWire: + """Re-create a bluemira wire from a series of WireInfo.""" + wire_list = [] + for info in self: + start_end = np.array(info.key_points[:2]) + if info.wire: + # quick way to get the wire back without doing any computation is by + # directly copying. + wire_list.append(info.wire) + continue + if isinstance(info.key_points, StraightLineInfo): + info.wire = make_polygon(start_end.T, closed=False) + wire_list.append(info.wire) + else: + # given two points on the circumference, only makes the SHORTER of the + # two possible arcs of the circle. + chord_intersect = start_end.mean(axis=0) + direction = chord_intersect - info.key_points.center + normed_dir = direction / np.linalg.norm(direction) + middle = info.key_points.center + info.key_points.radius * normed_dir + info.wire = make_circle_arc_3P(start_end[0], middle, start_end[1]) + wire_list.append(info.wire) + return BluemiraWire(wire_list) diff --git a/conda/environment.yml b/conda/environment.yml index 49d5e7ae6a..d5c6066c7e 100644 --- a/conda/environment.yml +++ b/conda/environment.yml @@ -24,6 +24,7 @@ dependencies: - freecad=0.21.2 - h5py - c-blosc2=2.14.3 + - openmc>=0.14.0 - graphviz - pip - pip: diff --git a/eudemo/README.md b/eudemo/README.md index a9d08b5e89..6098e87fe1 100644 --- a/eudemo/README.md +++ b/eudemo/README.md @@ -23,3 +23,7 @@ relative to the `eudemo` directory. In future this will be moved to a separate repository. It should be used as a template for how we expect other reactor repositories to be structured. + +## Neutronics + +To use the axis-symmetric neutronics run you will need to download the required neutronics cross section data into the folder `eudemo/config/cross_section_data` or modify the build_config.json accordingly. Secondly you will need to provide a source, the default used in the `make_pps_source` function is not currently open source but should be available shortly. diff --git a/eudemo/config/build_config.json b/eudemo/config/build_config.json index 54cd5cd8b5..84ac36f19f 100644 --- a/eudemo/config/build_config.json +++ b/eudemo/config/build_config.json @@ -1,7 +1,7 @@ { "params": "$path:params.json", "Radial build": { - "run_mode": "run", + "run_mode": "read", "read_dir": "$path_expand:./", "run_dir": "$path_expand:./", "plot": true @@ -60,7 +60,7 @@ } }, "Free boundary equilibrium": { - "run_mode": "run", + "run_mode": "read", "plot": true, "save": true, "fixed_eq_file_path": "$path_expand:./fixed_boundary_eqdsk.json", @@ -105,8 +105,22 @@ }, "Divertor silhouette": {} }, + "Neutronics": { + "cross_section_xml": "$path_expand:./cross_section_data/cross_sections.xml", + "particles": 16800, + "batches": 3, + "photon_transport": true, + "electron_treatment": "ttb", + "run_mode": "run_and_plot", + "openmc_write_summary": false, + "parametric_source": true, + "blanket_type": "HCPB", + "plot_axis": "xz", + "plot_pixel_per_metre": 100, + "neutronics_output_path": "$path_expand:./neutronics" + }, "TF coils": { - "run_mode": "run", + "run_mode": "read", "file_path": "$path_expand:./TFCoilDesign.json", "plot": true, "param_class": "TripleArc", @@ -130,7 +144,7 @@ } }, "PF coils": { - "run_mode": "run", + "run_mode": "read", "file_path": "$path_expand:./PFCoilDesign.json", "verbose": false, "plot": true, diff --git a/eudemo/config/params.json b/eudemo/config/params.json index 490549cc2f..895c511457 100644 --- a/eudemo/config/params.json +++ b/eudemo/config/params.json @@ -1045,5 +1045,71 @@ "unit": "A/W/m^2", "source": "Input", "long_name": "Electron cyclotron resonce heating current drive efficiency" + }, + "reactor_power": { + "value": 1998, + "unit": "MW", + "source": "Input", + "long_name": "neutronics reactor power" + }, + "peaking_factor": { + "value": 1.508, + "unit": "", + "source": "Input", + "long_name": "neutronics peaking factor" + }, + "vertical_shift": { + "value": 0, + "unit": "m", + "source": "Input", + "long_name": "neutronics plasma vertical shift" + }, + "inboard_fw_tk": { + "value": 0, + "unit": "m", + "source": "Input", + "long_name": "First wall inboard thickness" + }, + "outboard_fw_tk": { + "value": 0, + "unit": "m", + "source": "Input", + "long_name": "First wall outboard thickness" + }, + "fw_blanket_surface_tk": { + "value": 0.01, + "unit": "m", + "source": "Input", + "long_name": "First wall blanket surface thickness" + }, + "fw_divertor_surface_tk": { + "value": 0.1, + "unit": "m", + "source": "Input", + "long_name": "First wall divertor surface thickness" + }, + "inboard_breeding_tk": { + "value": 0, + "unit": "m", + "source": "Input", + "long_name": "Blanket breeding zone inboard thickness" + }, + "outboard_breeding_tk": { + "value": 0, + "unit": "m", + "source": "Input", + "long_name": "Blanket breeding zone outboard thickness" + }, + "blk_ib_manifold": { + "value": 0.02, + "unit": "m", + "source": "Input", + "long_name": "Blanket manifold inboard thickness" + }, + "blk_ob_manifold": { + "value": 0.2, + "unit": "m", + "source": "Input", + "long_name": "Blanket manifold outboard thickness" } } diff --git a/eudemo/eudemo/blanket/__init__.py b/eudemo/eudemo/blanket/__init__.py index 9adc53b84c..8e0ba976e8 100644 --- a/eudemo/eudemo/blanket/__init__.py +++ b/eudemo/eudemo/blanket/__init__.py @@ -5,8 +5,11 @@ # SPDX-License-Identifier: LGPL-2.1-or-later """Builders, designers, and components for an EUDEMO blanket.""" +from bluemira.base.components import Component from bluemira.base.reactor import ComponentManager +from bluemira.geometry.coordinates import Coordinates from bluemira.geometry.face import BluemiraFace +from bluemira.geometry.wire import BluemiraWire from eudemo.blanket.builder import BlanketBuilder from eudemo.blanket.designer import BlanketDesigner @@ -16,14 +19,33 @@ class Blanket(ComponentManager): """Wrapper around a Blanket component tree.""" - def inboard_xz_silhouette(self) -> BluemiraFace: - """The poloidal plane silhouette of the inboard blanket segment.""" + def __init__( + self, component_tree: Component, panel_points: Coordinates, r_inner_cut: float + ): + self.r_inner_cut = r_inner_cut + self._panel_points = panel_points + super().__init__(component_tree) + + def panel_points(self) -> Coordinates: + """The panel points of the blanket.""" + return self._panel_points + + def inboard_xz_face(self) -> BluemiraFace: + """The poloidal plane face of the inboard blanket segment.""" return ( self.component().get_component("xz").get_component(BlanketBuilder.IBS).shape ) - def outboard_xz_silhouette(self) -> BluemiraFace: - """The poloidal plane silhouette of the outboard blanket segment.""" + def outboard_xz_face(self) -> BluemiraFace: + """The poloidal plane face of the outboard blanket segment.""" return ( self.component().get_component("xz").get_component(BlanketBuilder.OBS).shape ) + + def inboard_xz_boundary(self) -> BluemiraWire: + """The toroidal plane silhouette of the inboard blanket segment.""" + return self.inboard_xz_face().boundary[0] + + def outboard_xz_boundary(self) -> BluemiraWire: + """The poloidal plane silhouette of the outboard blanket segment.""" + return self.outboard_xz_face().boundary[0] diff --git a/eudemo/eudemo/blanket/designer.py b/eudemo/eudemo/blanket/designer.py index 45c6ada664..a52616d65b 100644 --- a/eudemo/eudemo/blanket/designer.py +++ b/eudemo/eudemo/blanket/designer.py @@ -9,12 +9,14 @@ from typing import TypeVar import numpy as np +from scipy.spatial.distance import euclidean from bluemira.base.designer import Designer from bluemira.base.error import BuilderError from bluemira.base.look_and_feel import bluemira_warn from bluemira.base.parameter_frame import Parameter, ParameterFrame from bluemira.geometry.constants import VERY_BIG +from bluemira.geometry.coordinates import Coordinates from bluemira.geometry.face import BluemiraFace from bluemira.geometry.tools import boolean_cut, make_polygon from bluemira.geometry.wire import BluemiraWire @@ -50,7 +52,7 @@ class BlanketSegments: outboard_boundary: BluemiraWire -class BlanketDesigner(Designer[tuple[BluemiraFace, BluemiraFace]]): +class BlanketDesigner(Designer[tuple[BluemiraFace, BluemiraFace, Coordinates]]): """ Designer for an EUDEMO-style blanket. @@ -103,7 +105,33 @@ def __init__( ) self.cut_angle = cut_angle - def run(self) -> tuple[BluemiraFace, BluemiraFace]: + @staticmethod + def _remove_gap_and_merge( + ib_panels: Coordinates, ob_panels: Coordinates + ) -> Coordinates: + """Merge two sets of coordinates.""" + sm_dist = np.inf + closest_points = None + for ib_pt in ib_panels.points: + dists = [euclidean(ib_pt, ob_pt) for ob_pt in ob_panels.points] + sm_dist_idx = np.argmin(dists) + d = dists[sm_dist_idx] + if d < sm_dist: + sm_dist = d + ob_pt = ob_panels.points[sm_dist_idx] + closest_points = (ib_pt, ob_pt) + + new_coords = [ + p + for p in ib_panels.points + ob_panels.points + if not ( + np.array_equal(p, closest_points[0]) + or np.array_equal(p, closest_points[1]) + ) + ] + return Coordinates(new_coords) + + def run(self) -> tuple[BluemiraFace, BluemiraFace, Coordinates]: """Run the blanket design problem.""" segments = self.segment_blanket() # Inboard @@ -114,7 +142,19 @@ def run(self) -> tuple[BluemiraFace, BluemiraFace]: ob_panels = self.panel_boundary(segments.outboard_boundary) ob_panels_face = BluemiraFace(ob_panels) cut_ob = boolean_cut(segments.outboard, [ob_panels_face])[0] - return cut_ib, cut_ob + + # to get the panel points, we must remove the gap points + # introduced in the ib <-> ob slice + ib_verts = ib_panels.vertexes.points + ob_verts = ob_panels.vertexes.points + ib_verts[0] = cut_ib.vertexes.points[cut_ib.vertexes.argmin(ib_verts[0])] + ob_verts[-1] = cut_ob.vertexes.points[cut_ob.vertexes.argmin(ob_verts[-1])] + + panel_points = self._remove_gap_and_merge( + Coordinates(ib_verts), Coordinates(ob_verts) + ) + + return cut_ib, cut_ob, panel_points def segment_blanket(self) -> BlanketSegments: """ diff --git a/eudemo/eudemo/ivc/__init__.py b/eudemo/eudemo/ivc/__init__.py index 12eac2b476..633b8aaa62 100644 --- a/eudemo/eudemo/ivc/__init__.py +++ b/eudemo/eudemo/ivc/__init__.py @@ -14,6 +14,7 @@ from bluemira.base.designer import run_designer from bluemira.equilibria import find_OX_points +from bluemira.geometry.tools import boolean_cut from eudemo.ivc.divertor_silhouette import DivertorSilhouetteDesigner from eudemo.ivc.ivc_boundary import IVCBoundaryDesigner from eudemo.ivc.plasma_face import PlasmaFaceDesigner @@ -32,6 +33,7 @@ class IVCShapes: """Collection of geometries used to design/build in-vessel components.""" blanket_face: BluemiraFace + div_internal_boundary: BluemiraWire divertor_face: BluemiraFace outer_boundary: BluemiraWire inner_boundary: BluemiraWire @@ -81,8 +83,13 @@ def design_ivc( # We want the boundary wire and face to start and end at the same # place, so we cut the wire again here. wall_boundary = cut_wall_below_x_point(wall_boundary, plasma_face.bounding_box.z_min) + + div_boundary = divertor_face.boundary[0] + div_internal_boundary = boolean_cut(div_boundary, ivc_boundary)[0] + return IVCShapes( blanket_face=plasma_face, + div_internal_boundary=div_internal_boundary, divertor_face=divertor_face, outer_boundary=ivc_boundary, inner_boundary=wall_boundary, diff --git a/eudemo/eudemo/neutronics/__init__.py b/eudemo/eudemo/neutronics/__init__.py new file mode 100644 index 0000000000..8855268e41 --- /dev/null +++ b/eudemo/eudemo/neutronics/__init__.py @@ -0,0 +1,6 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""EUDEMO neutronics""" diff --git a/eudemo/eudemo/neutronics/run.py b/eudemo/eudemo/neutronics/run.py new file mode 100644 index 0000000000..19d5f499bb --- /dev/null +++ b/eudemo/eudemo/neutronics/run.py @@ -0,0 +1,99 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Test script to make the CSG branch work.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from bluemira.codes.wrapper import neutronics_code_solver +from bluemira.radiation_transport.error import NeutronicsError +from bluemira.radiation_transport.neutronics.blanket_data import ( + create_materials, + get_preset_physical_properties, +) +from bluemira.radiation_transport.neutronics.geometry import TokamakDimensions +from bluemira.radiation_transport.neutronics.neutronics_axisymmetric import ( + NeutronicsReactor, + NeutronicsReactorParameterFrame, +) + +if TYPE_CHECKING: + from collections.abc import Callable + + import numpy.typing as npt + import openmc.source + + from bluemira.base.parameter_frame import ParameterFrame + from bluemira.base.reactor import ComponentManager + from bluemira.codes.openmc.params import PlasmaSourceParameters + from bluemira.geometry.wire import BluemiraWire + from eudemo.blanket import Blanket + from eudemo.ivc import IVCShapes + from eudemo.vacuum_vessel import VacuumVessel + + +class EUDEMONeutronicsCSGReactor(NeutronicsReactor): + """EUDEMO Axis-symmetric neutronics model""" + + def _get_wires_from_components( + self, + ivc_shapes: IVCShapes, + blanket: Blanket, + vacuum_vessel: VacuumVessel, + ) -> tuple[TokamakDimensions, BluemiraWire, npt.NDArray, BluemiraWire, BluemiraWire]: + return ( + TokamakDimensions.from_parameterframe(self.params, blanket.r_inner_cut), + ivc_shapes.div_internal_boundary, + blanket.panel_points().T, + ivc_shapes.outer_boundary, + vacuum_vessel.xz_boundary(), + ) + + +def run_neutronics( + params: dict | ParameterFrame, + build_config: dict, + blanket: ComponentManager, + vacuum_vessel: ComponentManager, + ivc_shapes: IVCShapes, + source: Callable[[PlasmaSourceParameters], openmc.source.SourceBase] | None = None, + tally_function=None, +): + """Runs the neutronics model""" + # TODO get these materials from the componentmanager or something similar + breeder_materials, tokamak_geometry = get_preset_physical_properties( + build_config.pop("blanket_type") + ) + material_library = create_materials(breeder_materials) + + csg_params = NeutronicsReactorParameterFrame.from_config_params(params) + csg_params.update_from_dict({ + "inboard_fw_tk": {"value": tokamak_geometry.inb_fw_thick, "unit": "m"}, + "inboard_breeding_tk": {"value": tokamak_geometry.inb_bz_thick, "unit": "m"}, + "outboard_fw_tk": {"value": tokamak_geometry.outb_fw_thick, "unit": "m"}, + "outboard_breeding_tk": {"value": tokamak_geometry.outb_bz_thick, "unit": "m"}, + }) + neutronics_csg = EUDEMONeutronicsCSGReactor( + csg_params, ivc_shapes, blanket, vacuum_vessel, material_library + ) + if source is None: + try: + from bluemira.codes.openmc.sources import make_pps_source # noqa: PLC0415 + except ImportError: + raise NeutronicsError("Cannot import neutronics source") from None + + solver = neutronics_code_solver( + params, + build_config, + neutronics_csg, + source=source or make_pps_source, + tally_function=tally_function, + ) + + res = solver.execute() + + return neutronics_csg, res diff --git a/eudemo/eudemo/params.py b/eudemo/eudemo/params.py index 435f7f4338..cd0f6eadbf 100644 --- a/eudemo/eudemo/params.py +++ b/eudemo/eudemo/params.py @@ -225,3 +225,16 @@ class EUDEMOReactorParams(ParameterFrame): # First wall panelling fw_a_max: Parameter[float] fw_dL_min: Parameter[float] # noqa: N815 + + # CSG neutronics + reactor_power: Parameter[float] + peaking_factor: Parameter[float] + vertical_shift: Parameter[float] + inboard_fw_tk: Parameter[float] + outboard_fw_tk: Parameter[float] + fw_blanket_surface_tk: Parameter[float] + fw_divertor_surface_tk: Parameter[float] + inboard_breeding_tk: Parameter[float] + outboard_breeding_tk: Parameter[float] + blk_ib_manifold: Parameter[float] + blk_ob_manifold: Parameter[float] diff --git a/eudemo/eudemo/reactor.py b/eudemo/eudemo/reactor.py index 6ffaf9b17b..0bf1910bfe 100644 --- a/eudemo/eudemo/reactor.py +++ b/eudemo/eudemo/reactor.py @@ -78,6 +78,7 @@ ) from eudemo.maintenance.upper_port import UpperPortKOZDesigner from eudemo.model_managers import EquilibriumManager +from eudemo.neutronics.run import run_neutronics from eudemo.params import EUDEMOReactorParams from eudemo.pf_coils import PFCoil, PFCoilsDesigner, build_pf_coils_component from eudemo.power_cycle import SteadyStatePowerCycleSolver @@ -207,9 +208,9 @@ def build_blanket( designer = BlanketDesigner( params, blanket_boundary, blanket_face, r_inner_cut, cut_angle ) - ib_silhouette, ob_silhouette = designer.execute() + ib_silhouette, ob_silhouette, panel_points = designer.execute() builder = BlanketBuilder(params, build_config, ib_silhouette, ob_silhouette) - return Blanket(builder.build()) + return Blanket(builder.build(), panel_points, r_inner_cut) def build_tf_coils(params, build_config, separatrix, vvts_cross_section) -> TFCoil: @@ -479,6 +480,14 @@ def build_radiation_plugs(params, build_config, cr_ports, radiation_xz_boundary) cut_angle, ) + csg_model, neutronics_output = run_neutronics( + reactor_config.params_for("Neutronics"), + reactor_config.config_for("Neutronics"), + blanket=reactor.blanket, + vacuum_vessel=reactor.vacuum_vessel, + ivc_shapes=ivc_shapes, + ) + vv_thermal_shield = build_vacuum_vessel_thermal_shield( reactor_config.params_for("Thermal shield"), reactor_config.config_for("Thermal shield", "VVTS"), diff --git a/eudemo/eudemo_tests/blanket/test_blanket.py b/eudemo/eudemo_tests/blanket/test_blanket.py index 4958ec0ac5..263eaecb73 100644 --- a/eudemo/eudemo_tests/blanket/test_blanket.py +++ b/eudemo/eudemo_tests/blanket/test_blanket.py @@ -4,6 +4,7 @@ # # SPDX-License-Identifier: LGPL-2.1-or-later +from bluemira.geometry.coordinates import Coordinates from eudemo.blanket import Blanket, BlanketBuilder from eudemo_tests.blanket.tools import make_simple_blanket @@ -22,22 +23,32 @@ def make_blanket_component(): ib_silhouette=segments.inboard, ob_silhouette=segments.outboard, ) - return segments, builder.build() + return params, segments, builder.build() class TestBlanket: def test_inboard_xz_silhouette_face_from_BlanketBuilder_component_tree(self): - segments, component = make_blanket_component() + params, segments, component = make_blanket_component() + io_cut = segments.inboard.bounding_box.x_max + panel_points = Coordinates( + segments.inboard_boundary.vertexes.points + + segments.outboard_boundary.vertexes.points + ) - blanket = Blanket(component_tree=component) - ib_face = blanket.inboard_xz_silhouette() + blanket = Blanket(component, panel_points, io_cut) + ib_face = blanket.inboard_xz_face() assert ib_face is segments.inboard def test_outboard_xz_silhouette_face_from_BlanketBuilder_component_tree(self): - segments, component = make_blanket_component() - - blanket = Blanket(component_tree=component) - ob_face = blanket.outboard_xz_silhouette() + params, segments, component = make_blanket_component() + io_cut = segments.inboard.bounding_box.x_max + panel_points = Coordinates( + segments.inboard_boundary.vertexes.points + + segments.outboard_boundary.vertexes.points + ) + + blanket = Blanket(component, panel_points, io_cut) + ob_face = blanket.outboard_xz_face() assert ob_face is segments.outboard diff --git a/eudemo/eudemo_tests/blanket/test_designer.py b/eudemo/eudemo_tests/blanket/test_designer.py index f3dcac0838..b5477982e1 100644 --- a/eudemo/eudemo_tests/blanket/test_designer.py +++ b/eudemo/eudemo_tests/blanket/test_designer.py @@ -136,7 +136,7 @@ def test_segments_cut_using_panels(self, paneller_mock): with mock.patch.object(designer, "segment_blanket") as sb_mock: sb_mock.return_value = blanket - ib, ob = designer.run() + ib, ob, panels = designer.run() # These areas were (painstakingly) worked out by hand panel_trapezium_area = 9 / 2 * (4 * np.sqrt(2) - 5) diff --git a/examples/radiation_transport/radiation_calculation_solver_DEMO.py b/examples/radiation_transport/radiation_calculation_solver_DEMO.py index 519f48ee5c..5a7c755b3e 100644 --- a/examples/radiation_transport/radiation_calculation_solver_DEMO.py +++ b/examples/radiation_transport/radiation_calculation_solver_DEMO.py @@ -195,7 +195,7 @@ if core_filter_in(point): rad_sol_grid[j, i] = interpolated_field_values( x_sol[i], z_sol[j], f_core - ) + ).item() else: rad_sol_grid[j, i] = ( rad_sol_grid[j, i] @@ -203,7 +203,7 @@ * (pfr_down_filter(point) * 1.0) * (pfr_up_filter(point) * 1.0) * (core_filter_out(point) * 1.0) - ) + ).item() func = grid_interpolator(x_sol, z_sol, rad_sol_grid) # Calculate radiation of FW points diff --git a/pyproject.toml b/pyproject.toml index 8c66b535ef..b24ebc9f87 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -79,11 +79,7 @@ pinned = [ "matplotlib==3.8.3", "scipy==1.10.1", ] -openmoc = ["OpenMOC @git+https://github.com/mit-crpg/OpenMOC.git@7940c0b"] -openmc = [ - "OpenMC @git+https://github.com/openmc-dev/openmc.git", - "parametric-plasma-source @git+https://github.com/open-radiation-sources/parametric-plasma-source.git", -] +openmc = ["openmc>=0.14.0", "openmc_data"] polyscope = ["polyscope"] radiation = ["cherab"] @@ -311,6 +307,7 @@ ignore-names = [ "*_F", "*eV", "I_not_dI", + "*Li*" ] [tool.ruff.lint.per-file-ignores] diff --git a/requirements.txt b/requirements.txt index feb50bf7e1..88df75f4b9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -21,6 +21,7 @@ numba==0.59.1 numba-scipy @ git+https://github.com/numba/numba-scipy@1e2f244 numexpr==2.9.0 numpy==1.26.4 +openmc_data==0.2.2 packaging==24.0 pandas==2.2.1 periodictable==1.7.0 diff --git a/scripts/install-openmc.sh b/scripts/install-openmc.sh new file mode 100755 index 0000000000..16285df745 --- /dev/null +++ b/scripts/install-openmc.sh @@ -0,0 +1,45 @@ +#!/bin/bash + +if [ "$1" ] + then + VERSION="$1" +else + VERSION="v0.14.0" +fi + +read -p "Are you in the correct python environment? (y/n) " answer +case ${answer:0:1} in + y|Y ) + ;; + * ) + exit;; +esac + +echo +echo Installing... +echo + +clean_up() { + test -d "$tmp_dir" && rm -rf "$tmp_dir" +} + +tmp_dir=$( mktemp -d -t install-openmc.XXX) +trap "clean_up $tmp_dir" EXIT + +set -euxo pipefail + +cd $tmp_dir + +sudo apt install g++ cmake libhdf5-dev libpng-dev libopenmpi-dev ninja-build -y + +git clone --recurse-submodules https://github.com/openmc-dev/openmc.git +cd openmc +git checkout $VERSION +mkdir build && cd build +cmake -GNinja -DCMAKE_INSTALL_PREFIX=$HOME/.local -DOPENMC_USE_MPI=ON .. +ninja install + +cd .. +pip install . + +echo Finished diff --git a/scripts/openmc_data_download/__init__.py b/scripts/openmc_data_download/__init__.py new file mode 100644 index 0000000000..33290173bb --- /dev/null +++ b/scripts/openmc_data_download/__init__.py @@ -0,0 +1 @@ +"""Neutronics Data Downloader""" diff --git a/scripts/openmc_data_download/multithreaded_download.py b/scripts/openmc_data_download/multithreaded_download.py new file mode 100644 index 0000000000..16aa91c767 --- /dev/null +++ b/scripts/openmc_data_download/multithreaded_download.py @@ -0,0 +1,132 @@ +"""Asynchronous multithreaded file downloader""" + +import asyncio +import concurrent.futures +import functools +import hashlib +import os +import signal +from collections.abc import Callable +from pathlib import Path + +import requests +from rich.progress import Progress + + +async def get_size(url: str, timeout: int = 10) -> int: + """Get size of file""" + response = requests.head( # noqa: ASYNC100 + url, allow_redirects=True, timeout=timeout + ) + return int(response.headers["Content-Length"]) + + +def download_range( + url: str, + start: int, + end: int, + output: Path, + progress: Callable[[], None], + timeout: int = 10, +): + """Download worker""" + if not (output.is_file() and output.stat().st_size == end + 1 - start): + headers = {"User-Agent": "Mozilla/5.0", "Range": f"bytes={start}-{end}"} + response = requests.get(url, headers=headers, timeout=timeout) + with open(output, "wb") as f: + for part in response.iter_content(1024): + f.write(part) + progress() + + +async def _download( + run: Callable, + url: str, + output: Path, + *, + chunk_size: int = 10000000, + timeout: int = 10, +): + """Downloader event loop""" + if not url.startswith(("http:", "https:")): + raise ValueError("Must be an http or https url") + + file_size = await get_size(url) + + if output.is_file() and output.stat().st_size == file_size: + return + + chunks = range(0, file_size, chunk_size) + with Progress() as progress: + updater = functools.partial( + progress.update, + progress.add_task(f"Downloading {output.parts[-1]}", total=len(chunks)), + advance=1, + ) + tasks = [ + run( + download_range, + url, + start, + start + chunk_size - 1, + Path(f"{output}.part{i}"), + updater, + timeout, + ) + for i, start in enumerate(chunks) + ] + + await asyncio.wait(tasks) + + with open(output, "wb") as o: # noqa: ASYNC101 + for i in range(len(chunks)): + chunk_path = Path(f"{output}.part{i}") + o.write(Path(chunk_path).read_bytes()) + + # Write all before deleting + for i in range(len(chunks)): + Path(f"{output}.part{i}").unlink() + + +def downloader( + url: str, + checksum: int | None = None, + as_browser: bool = False, # noqa: ARG001, FBT001, FBT002 + output_path: os.PathLike | None = None, + output_filename: os.PathLike | None = None, + *, + max_workers: int = 5, + timeout: int = 10, +): + """Asynchronous multithreaded downloader""" + if output_filename is None: + output_filename = Path(url.rsplit("/", 1)[-1]) + + if output_path is None: + local_path = Path(output_filename) + else: + Path(output_path).mkdir(parents=True, exist_ok=True) + local_path = Path(output_path, output_filename) + + # Reenable Keyboard Interrupt + signal.signal(signal.SIGINT, signal.SIG_DFL) + + executor = concurrent.futures.ThreadPoolExecutor(max_workers=max_workers) + loop = asyncio.new_event_loop() + run = functools.partial(loop.run_in_executor, executor) + asyncio.set_event_loop(loop) + + try: + loop.run_until_complete(_download(run, url, local_path, timeout=timeout)) + finally: + loop.close() + + if ( + checksum is not None + and hashlib.md5(open(local_path, "rb").read()).hexdigest() # noqa: S324, SIM115 + != checksum + ): + raise OSError( + f"MD5 checksum for {local_path} does not match." + "Please ensure you checksum is correct." + ) diff --git a/scripts/openmc_data_download/neutronics_data_downloader.py b/scripts/openmc_data_download/neutronics_data_downloader.py new file mode 100644 index 0000000000..a6da6ed250 --- /dev/null +++ b/scripts/openmc_data_download/neutronics_data_downloader.py @@ -0,0 +1,248 @@ +"""Neutronics data downloader main script""" + +import argparse +import fnmatch +import functools +import json +import sys +import tarfile +import zipfile +from collections.abc import Callable, Iterable +from os import chdir +from pathlib import Path +from shutil import rmtree +from types import ModuleType +from unittest.mock import patch +from xml.etree import cElementTree # noqa: S405 + +from multithreaded_download import downloader +from rich.progress import track + +from bluemira.base.look_and_feel import ( + bluemira_print, + bluemira_print_flush, + bluemira_warn, +) + + +def state_download_size(download_size: int, uncompressed_size: int, units: str): + """Patch warning to logging""" + bluemira_warn( + f"This script will download up to {download_size} {units} " + "of data. Extracting and processing the data may require as much " + f"as {uncompressed_size} {units} of additional free disk space." + ) + + +def extractor( + compressed_files: list[Path], + extraction_dir: Path, + del_compressed_file: bool, # noqa: FBT001 +): + """Customisable extractor""" + Path.mkdir(extraction_dir, parents=True, exist_ok=True) + + if not isinstance(compressed_files, Iterable): + compressed_files = [compressed_files] + + for f in compressed_files: + suffix = "".join(f.suffixes[-2:]) + if suffix.endswith(".zip"): + with zipfile.ZipFile(f, "r") as zip_handler: + bluemira_print("Getting file list") + file_list = {z.filename: z for z in zip_handler.infolist()} + for m in track( + filter_members(f.parts[-1], file_list), description="Extracting" + ): + file = Path(extraction_dir, m.filename) + if not file.is_file() or file.stat().st_size != m.file_size: + zip_handler.extract(m, path=extraction_dir) + + elif suffix in {".tar.gz", ".tgz", ".tar.bz2", ".tar.xz", ".xz"}: + with tarfile.open(f, "r") as tgz: + bluemira_print("Getting file list") + file_list = {m.get_info()["name"]: m for m in tgz.getmembers()} + for m in track( + filter_members(f.parts[-1], file_list), description="Extracting" + ): + m_inf = m.get_info() + file = Path(extraction_dir, m_inf["name"]) + if not file.is_file() or file.stat().st_size != m_inf["size"]: + tgz.extract(m, path=extraction_dir) + else: + raise ValueError( + f"File type not currently supported by extraction function {f!s}" + ) + + if del_compressed_file: + rmtree(compressed_files, ignore_errors=True) + + +def _filter_members( + file: str, + isotope_file: str, + filename: str, + members: dict[str, tarfile.TarInfo | zipfile.ZipInfo], +) -> list[tarfile.TarInfo] | list[zipfile.ZipInfo]: + """Filter archive contents to only extract wanted files""" + import openmc_data.convert.convert_tendl as tendl # noqa: PLC0415 + from openmc_data import all_release_details as ard # noqa: PLC0415 + + with open(Path(Path(file).parent, isotope_file)) as fh: + isotope_data = json.load(fh) + if filename in ard["tendl"][tendl.args.release]["neutron"]["compressed_files"]: + return _filter( + filename, members, isotope_data["tendl"]["neutron"], lambda m: f"*/{m}" + ) + if filename == ard["endf"]["b7.1"]["neutron"]["compressed_files"][0]: + return _filter( + filename, + members, + isotope_data["endf"]["neutron"], + lambda m: f"*/{m[:-3]}_{m[-3:]}*.ace", + ) + if filename == ard["endf"]["b7.1"]["neutron"]["compressed_files"][1]: + return _filter(filename, members, ["bebeo", "obeo"], lambda m: f"{m}.acer") + if filename in ard["endf"]["b7.1"]["photon"]["compressed_files"]: + return _filter( + filename, members, isotope_data["endf"]["photon"], lambda m: f"*{m}*.endf" + ) + raise ValueError("Unknown archive") + + +def _filter( + filename: str, + members: dict[str, tarfile.TarInfo | zipfile.ZipInfo], + datakeys: list[str], + filt: Callable, +) -> list[tarfile.TarInfo] | list[zipfile.ZipInfo]: + """Filter archive members""" + filtered_members = [] + mem_keys = members.keys() + if datakeys == "*": + return list(members.values()) + for m in datakeys: + if file := fnmatch.filter(mem_keys, filt(m)): + filtered_members.append(members[file[0]]) + else: + bluemira_warn(f"Cant find {m} in {filename}") + return filtered_members + + +def _convertion_progress(*string, **_kwargs): + """Patch print function for logging""" + bluemira_print_flush("".join([str(s) for s in string])) + + +def combine_xml( + lib_names: tuple[str, ...], + thermal_files: tuple[str, ...], + thermal_prefix: Path, +): + """Combine xml files""" + bluemira_print("Removing uneeded files") + for i in thermal_files: + Path(thermal_prefix / f"{i}.h5").unlink() + for file in ["bebeo", "obeo"]: + Path("endf-b7.1-ace", f"{file}.acer").unlink() + + bluemira_print("Combining cross section xml files") + xml_handle = [ + cElementTree.parse(Path(name, "cross_sections.xml")) # noqa: S313 + for name in lib_names + ] + for name, xml in zip(lib_names, xml_handle, strict=False): + data = xml.getroot() + remove_list = [] + for elem in data.iter(): + if elem.tag == "library": + if elem.attrib["materials"] in thermal_files: + remove_list.append(elem) + else: + elem.attrib["path"] = str(Path(name, elem.attrib["path"])) + for elem in remove_list: + data.remove(elem) + + data = xml_handle[0].getroot() + for xml in xml_handle[1:]: + data.extend(list(xml.getroot().iter())[1:]) + + xml_handle[0].write("cross_sections.xml") + + +def download_data( + download: Callable, + libs: tuple[ModuleType, ...], + lib_names: tuple[str, ...], +): + """Download neutronics data""" + for name, lib in zip(lib_names, libs, strict=False): + bluemira_print(f"Downloading {name} cross section data") + lib.state_download_size = state_download_size + lib.download = download + lib.extract = extractor + lib.args.destination = Path(name) + with patch("builtins.print", new=_convertion_progress): + lib.main() + print() + + +def parse_args(args: list[str] | None = None) -> argparse.Namespace: + """Parse arguments""" + parser = argparse.ArgumentParser("Bluemira Neutronics data downloader") + + parser.add_argument("-l", "--location", default=Path.cwd() / "bluemira_openmc_data") + parser.add_argument("--download-threads", type=int, default=5) + parser.add_argument("--isotope-file", type=str, default="nuclear_data_isotopes.json") + p, unknown = parser.parse_known_args(args) + p.location = Path(p.location) + sys.argv = [sys.argv[0], *unknown] + return p + + +class ChgDir: + """Change directory context manager""" + + def __init__(self, path: Path): + self.path = path + self.origin = Path().absolute() + + def __enter__(self): + """Change directory""" + chdir(self.path) + + def __exit__(self, typ, exc, tb): + """Revert directory change""" + chdir(self.origin) + + +def main(p): + """Main function""" + root_folder = p.location + root_folder.mkdir(parents=True, exist_ok=True) + + download = functools.partial(downloader, max_workers=p.download_threads) + + # Imported after parsing arguments because argparse is called on import here... + import openmc_data.convert.convert_endf as endf # noqa: PLC0415 + import openmc_data.convert.convert_tendl as tendl # noqa: PLC0415 + + libs = (tendl, endf) + lib_names = tuple(lib.__name__.split("_")[-1] for lib in libs) + + with ChgDir(root_folder): + download_data(download, libs, lib_names) + + # convert_endf crashes if you dont have these available... + thermal_files = ("c_O_in_BeO", "c_Be_in_BeO") + thermal_prefix = endf.args.destination / "neutron" + + combine_xml(lib_names, thermal_files, thermal_prefix) + + +if __name__ == "__main__": + p = parse_args(None) + filter_members = functools.partial( + _filter_members, str(Path(__file__).resolve()), p.isotope_file + ) + main(p) diff --git a/scripts/openmc_data_download/nuclear_data_isotopes.json b/scripts/openmc_data_download/nuclear_data_isotopes.json new file mode 100644 index 0000000000..4a4a3bc8d2 --- /dev/null +++ b/scripts/openmc_data_download/nuclear_data_isotopes.json @@ -0,0 +1,9 @@ +{ + "tendl": { + "neutron": "*" + }, + "endf": { + "neutron": "*", + "photon": "*" + } +} diff --git a/scripts/openmc_data_download/nuclear_data_isotopes_example.json b/scripts/openmc_data_download/nuclear_data_isotopes_example.json new file mode 100644 index 0000000000..52677ebda7 --- /dev/null +++ b/scripts/openmc_data_download/nuclear_data_isotopes_example.json @@ -0,0 +1,386 @@ +{ + "tendl": { + "neutron": [ + "C012", + "C013", + "Ne020", + "Ne021", + "Ne022", + "O018", + "Os184", + "Os186", + "Os187", + "Os188", + "Os189", + "Os190", + "Os192", + "Pt190", + "Pt192", + "Pt194", + "Pt195", + "Pt196", + "Pt198", + "Yb168", + "Yb170", + "Yb171", + "Yb172", + "Yb173", + "Yb174", + "Yb176" + ] + }, + "endf": { + "neutron": [ + "Ag107", + "Ag109", + "Al027", + "Ar036", + "Ar038", + "Ar040", + "As075", + "Au197", + "B010", + "B011", + "Ba130", + "Ba132", + "Ba134", + "Ba135", + "Ba136", + "Ba137", + "Ba138", + "Be009", + "Bi209", + "Br079", + "Br081", + "C000", + "Ca040", + "Ca042", + "Ca043", + "Ca044", + "Ca046", + "Ca048", + "Cd106", + "Cd108", + "Cd110", + "Cd111", + "Cd112", + "Cd113", + "Cd114", + "Cd116", + "Ce136", + "Ce138", + "Ce140", + "Ce142", + "Cl035", + "Cl037", + "Co059", + "Cr050", + "Cr052", + "Cr053", + "Cr054", + "Cs133", + "Cu063", + "Cu065", + "Dy156", + "Dy158", + "Dy160", + "Dy161", + "Dy162", + "Dy163", + "Dy164", + "Er162", + "Er164", + "Er166", + "Er167", + "Er168", + "Er170", + "Eu151", + "Eu153", + "F019", + "Fe054", + "Fe056", + "Fe057", + "Fe058", + "Ga069", + "Ga071", + "Gd152", + "Gd154", + "Gd155", + "Gd156", + "Gd157", + "Gd158", + "Gd160", + "Ge070", + "Ge072", + "Ge073", + "Ge074", + "Ge076", + "H001", + "H002", + "He003", + "He004", + "Hf174", + "Hf176", + "Hf177", + "Hf178", + "Hf179", + "Hf180", + "Hg196", + "Hg198", + "Hg199", + "Hg200", + "Hg201", + "Hg202", + "Hg204", + "Ho165", + "I127", + "In113", + "In115", + "Ir191", + "Ir193", + "K039", + "K040", + "K041", + "Kr078", + "Kr080", + "Kr082", + "Kr083", + "Kr084", + "Kr086", + "La138", + "La139", + "Li006", + "Li007", + "Lu175", + "Lu176", + "Mg024", + "Mg025", + "Mg026", + "Mn055", + "Mo100", + "Mo092", + "Mo094", + "Mo095", + "Mo096", + "Mo097", + "Mo098", + "N014", + "N015", + "Na023", + "Nb093", + "Nd142", + "Nd143", + "Nd144", + "Nd145", + "Nd146", + "Nd148", + "Nd150", + "Ni058", + "Ni060", + "Ni061", + "Ni062", + "Ni064", + "O016", + "O017", + "P031", + "Pa231", + "Pb204", + "Pb206", + "Pb207", + "Pb208", + "Pd102", + "Pd104", + "Pd105", + "Pd106", + "Pd108", + "Pd110", + "Pr141", + "Rb085", + "Rb087", + "Re185", + "Re187", + "Rh103", + "Ru100", + "Ru101", + "Ru102", + "Ru104", + "Ru096", + "Ru098", + "Ru099", + "S032", + "S033", + "S034", + "S036", + "Sb121", + "Sb123", + "Sc045", + "Se074", + "Se076", + "Se077", + "Se078", + "Se080", + "Se082", + "Si028", + "Si029", + "Si030", + "Sm144", + "Sm147", + "Sm148", + "Sm149", + "Sm150", + "Sm152", + "Sm154", + "Sn112", + "Sn114", + "Sn115", + "Sn116", + "Sn117", + "Sn118", + "Sn119", + "Sn120", + "Sn122", + "Sn124", + "Sr084", + "Sr086", + "Sr087", + "Sr088", + "Ta180", + "Ta181", + "Tb159", + "Te120", + "Te122", + "Te123", + "Te124", + "Te125", + "Te126", + "Te128", + "Te130", + "Th230", + "Th232", + "Ti046", + "Ti047", + "Ti048", + "Ti049", + "Ti050", + "Tl203", + "Tl205", + "Tm169", + "U234", + "U235", + "U238", + "V050", + "V051", + "W180", + "W182", + "W183", + "W184", + "W186", + "Xe124", + "Xe126", + "Xe128", + "Xe129", + "Xe130", + "Xe131", + "Xe132", + "Xe134", + "Xe136", + "Y089", + "Zn064", + "Zn066", + "Zn067", + "Zn068", + "Zn070", + "Zr090", + "Zr091", + "Zr092", + "Zr094", + "Zr096" + ], + "photon": [ + "Ag", + "Al", + "Ar", + "As", + "Au", + "B", + "Ba", + "Be", + "Bi", + "Br", + "C", + "Ca", + "Cd", + "Ce", + "Cl", + "Co", + "Cr", + "Cs", + "Cu", + "Dy", + "Er", + "Eu", + "F", + "Fe", + "Ga", + "Gd", + "Ge", + "H", + "He", + "Hf", + "Hg", + "Ho", + "I", + "In", + "Ir", + "K", + "Kr", + "La", + "Li", + "Lu", + "Mg", + "Mn", + "Mo", + "N", + "Na", + "Nb", + "Nd", + "Ne", + "Ni", + "O", + "Os", + "P", + "Pa", + "Pb", + "Pd", + "Pr", + "Pt", + "Rb", + "Re", + "Rh", + "Ru", + "S", + "Sb", + "Sc", + "Se", + "Si", + "Sm", + "Sn", + "Sr", + "Ta", + "Tb", + "Te", + "Th", + "Ti", + "Tl", + "Tm", + "U", + "V", + "W", + "Xe", + "Y", + "Yb", + "Zn", + "Zr" + ] + } +}