Guiding-center ORbit Integration with Local Linearization Approach
GORILLA computes guiding-center orbits for charged particles of given mass, charge and energy in toroidal fusion devices with three-dimensional field geometry. This guiding-center orbit code is designed to be used in scientific plasma physics simulations in the field of magnetic confinement fusion.
The guiding-center orbits are traced via a quasi-geometric integration method described in Ref. [1]. There, high order interpolation of electromagnetic fields in space is replaced by a special linear interpolation, leading to locally linear Hamiltonian equations of motion with piecewise constant coefficients. The underlying formulation treats the motion in the piecewise linear fields exactly. This further leads to conservation of total energy, magnetic moment and phase space volume. Furthermore, the approach reduces computational effort and noise sensitivity. Guiding-center orbits are computed without taking collisions into account.
Due to its formulation in general curvilinear coordinates, GORILLA is not limited by the field topology. That means that the computation domain of GORILLA covers both the closed field line region (i.e. the plasma core) and the open field line region (i.e. the scrape-off layer).
For various simulations in magnetic confinement fusion, direct modeling of guiding-center particle orbits is utilized, e.g. global kinetic computations of quasi-steady plasma parameters or fast alpha particle loss estimation for stellarator optimization. In such complex simulations a simple interface for the guiding-center orbit integration part is needed. Namely, the initial condition in five-dimensional phase space is provided (i.e. guiding-center position, parallel and perpendicular velocity) and the main interest is in the condition after a prescribed time step while the integration process itself is irrelevant. Such a pure “orbit time step routine” acting as an interface with a plasma physics simulation is provided (orbit_timestep_gorilla
).
However, the integration process itself can be of high interest as well, thus, a program allowing the detailed analysis of guiding-center orbits, the time evolution of their respective invariants of motion and Poincaré plots is at disposal as well (gorilla_plot
).
Both applications are realized for demonstration in the program (test_gorilla_main
).
The code is free to use and modify under the MIT License and links to Runge-Kutta-Fehlberg routines in
SRC/contrib/rkf45.f90
from https://people.sc.fsu.edu/~jburkardt/f_src/rkf45/rkf45.html under the GNU LGPL License.
The magnetic field can be provided by magnetohydrodynamics (MHD) equilibria with nested magnetic flux surfaces either in 2D (e.g. EFIT) or in 3D (e.g. VMEC). Supported equilibria are in the g-file or NetCDF format, respectively.
For both equilibria formats, test files for the limited purpose of computing guiding-center orbits are provided.
The g-file test equilibrium (g_file_for_test
) was provided by the ASDEX Upgrade Team for testing purposes and corresponds to the axisymmetric tokamak field configuration of ASDEX Upgrade (shot 26884 at 4300 ms) described in Ref. [3].
The VMEC NetCDF test equlibrium (netcdf_file_for_test.nc
) was provided by Michael Drevlak for testing purposes and corresponds to the stellarator field configuration described in Ref. [4], namely, a quasi-isodynamic reactor-scale device with five toroidal field periods and a major radius of 25 m.
The g-file test equilibrium in WEST geometry (g_file_for_test_WEST
) was provided by the SOLEDGE3X-EIRENE Team for testing purposes and corresponds to the axisymmetric tokamak field configuration of WEST for the code SOLEDGE3X-EIRENE (shot 54903 at 8 s) described in Ref. [5].
A detailed description of the working principle of GORILLA can be found in DOCUMENTATION/GORILLA_DOC.pdf
.
The following supplemental material is available in DOCUMENTATION/SUPPLEMENTAL_MATERIAL
:
- arXiv preprint of Ref. [1]
- Master's thesis of M. Eder (preliminary work with some detailed explanations referenced in
GORILLA_DOC.pdf
) - Master's thesis of L. Bauer (preliminary work with some detailed explanations referenced in
GORILLA_DOC.pdf
) - Bachelor's thesis of D. Forstenlechner (unit testing with pFUnit)
GORILLA can be built with make
or cmake
.
- GNU Fortan
- NetCDF
- LAPACK/BLAS
To install requirements on Ubuntu Linux use
sudo apt install wget unzip git gfortran make cmake liblapack-dev libnetcdff-dev
To install requirements on macOS, install Homebrew, then
brew install wget unzip git gcc make cmake netcdf netcdf-fortran libomp
To install pFUnit, follow the instructions on the linked github project page.
To install lcov on Ubuntu Linux use
sudo apt install lcov
To install lcov on macOS use
brew install lcov
N. Flocke, “Algorithm 954: An Accurate and Efficient Cubic and Quartic Equation Solver for Physical Applications” https://doi.org/10.1145/2699468
- Download supplemental material
954.zip
from above webpage.
wget -O 954.zip "https://dl.acm.org/action/downloadSupplement?doi=10.1145%2F2699468&file=954.zip&download=true"
- Copy
954/F90/Src/Polynomial234RootSolvers.f90
toGORILLA/SRC/contrib/
and overwrite existing file. (Existing file with identical name is a placeholder which is necessary for compilation.)
unzip 954.zip
cp 954/F90/Src/Polynomial234RootSolvers.f90 SRC/contrib/
- GORILLA can be run without this external library. The computation of guiding-center orbits is then limited to the numerical Runge-Kutta option of GORILLA.
cd /path/to/GORILLA
make
This will produce test_gorilla_main.x
required to run the code. To specify the location of
NetCDF includes and libraries, one has to set the NCINC
and NCLIB
variable during make
.
To build GORILLA with cmake
, use build.sh
.
cd /path/to/GORILLA
./build.sh
This will produce test_gorilla_main.x
in the folder BUILD/SRC/ required to run the code.
To get additional tests and code coverage, build GORILLA with build_coverage.sh
. This requires the additional tools and a correct set PFUNIT_DIR. See the pFUnit github project page for additional information.
GORILLA currently runs on a single node with OpenMP shared memory parallelization with one particle per thread and background fields residing in main memory.
The main executable is test_gorilla_main.x
.
As an input it takes ....
... the following input files which can be found in the folder INPUT/
tetra_grid.inp
(Input file for settings of the tetrahedronal grid used in GORLLA)gorilla.inp
(Input file for settings of GORILLA)gorilla_plot.inp
(Input file for the program for the analysis of guiding-center orbits)field_divB0.inp
(Input file for loading g-file equilibria - Do not change this file.)preload_for_SYNCH.inp
(Input file for splining magnetic field data of g-file equilibria - Do not change this file.)
... and the MHD equilibrium files which can be found in the folder MHD_EQUILIBRIA/
netcdf_file_for_test.nc
: VMEC NetCDF equlibrium (File name can be specified intetra_grid.inp
.)g_file_for_test
org_file_for_test_WEST
: g-file equilibrium (File name can be specified intetra_grid.inp
.)
For compability with WEST geometry of SOLEDGE3X-EIRENE, additional input files describing the original 2D mesh are needed. Those can be found in MHD_EQUILIBRIA/MESH_SOLEDGE3X_EIRENE
-
knots_for_test.dat
: coordinates ($R$ ,$Z$ ) of the vertices making up the 2D grid (File name can be specified intetra_grid.inp
.) -
triangles_for_test.dat
: association of above mentioned vertices to triangles (triples of vertices) covering the 2D plane (File name can be specified intetra_grid.inp
.)
To produce these files (including the g-file equilibrium) oneself from files provided by SOLEDGE3X-EIRENE, a set of prepocessing MATLAB scripts are at disposal in REPROCESSING/SOLEDGE3X_EIRENE/MESH
and REPROCESSING/SOLEDGE3X_EIRENE/MHD_EQUILIBRIUM
respectively.
A tutorial for running GORILLA and plotting Poincaré cuts, full guiding-center orbits and the appropriate time evolution of invariants of motion is realized redundantly in both MATLAB and Python.
- MATLAB Live Script with the name
plotting_tutorial.mlx
is at disposal inMATLAB
as a step-by-step tutorial for all plotting features of GORILLA.
- Jupyter Notebook with the name
plotting_tutorial.ipynb
is at disposal inPYTHON
as a step-by-step tutorial for all plotting features of GORILLA.
For the Jupyter Notebook as well as the Python scripts used for the examples, some Python packages are needed, including the fortran namelist package f90nml. To install the necessary packages use for example pip.
python -m pip install --upgrade pip
pip install f90nml numpy matplotlib
However, GORILLA itself can be run without these packages. They are only used to incorporate GORILLA properly into Python scripts.
Seven examples for plotting Poincaré cuts, full guiding-center orbits (in plasma core or edge regions) and the appropriate time evolution of invariants of motion can be found in EXAMPLES/example_1
- EXAMPLES/example_7
. There, the necessary soft links are already created and the input files are given, and runs are started with
./test_gorilla_main.x #if the build was done with make
or
./test_gorilla_main_cmake.x #if the build was done with cmake
To avoid hyperthreading issues, it is beneficial to limit the number of threads to
the number of actual CPU cores via the environment variable $OMP_NUM_THREADS
.
Detailed descriptions of the respective input files can be found in INPUT
.
After appropriate compilation of GORILLA, the code can be executed in all of these 7 example folders, respectively.
For the visualization of the output of these seven examples, appropriate plotting methods for Python 3 are at disposal at PYTHON/plot_example_1.py
- PYTHON/plot_example_7.py
.
- Compute a collisionless guiding-center orbit with GORILLA for a trapped Deuterium particle.
- Use a field-aligned grid for a non-axisymmetric VMEC MHD equilibrium.
- Use the GORILLA polynomial option with order K = 4.
- Create a figure with the Poincaré sections (
$v_\parallel = 0$ ) in cylindrical and symmetry flux coordinates. - Compute the normalized parallel adiabatic invariant as a function of banana bounces.
- Compute a collisionless guiding-center orbit with GORILLA for a passing Deuterium particle.
- Use a field-aligned grid for a non-axisymmetric VMEC MHD equilibrium.
- Use the GORILLA Runge-Kutta 4 option.
- Create a figure with the Poincaré plots (
$\varphi = 0$ ) in cylindrical and symmetry flux coordinates. - Compute the normalized total energy as a function of toroidal mappings.
- Compute a collisionless guiding-center orbit with GORILLA for a passing Deuterium particle.
- Use a field-aligned grid for an axisymmetric tokamak equilibrium (g-file)
- Use the GORILLA polynomial option with order K = 3.
- Create a figure with the Poincaré plots (
$\varphi$ = 0) in cylindrical and symmetry flux coordinates. - Compute the normalized toroidal angular momentum as a function of toroidal mappings.
- Compute collisionless guiding-center orbits with GORILLA for two passing and one trapped Deuterium particle.
- Use a field-aligned grid for an axisymmetric tokamak equilibrium (g-file)
- Create a figure with the Poincaré plots (
$\varphi = 0$ ) in cylindrical and symmetry flux coordinates.
- Compute collisionless guiding-center orbits with GORILLA for a trapped Deuterium particle.
- Use a field-aligned grid for an axisymmetric tokamak equilibrium (g-file).
- Plot the plasma boundary, the guiding-center orbits, and the resulting Poincare plot (
$\varphi = 0$ ).
- Compute collisionless guiding-center orbits with GORILLA for passing Deuterium particles in scrape-off layer (WEST geometry).
- Construct a 3D extension of the SOLEDGE3X-EIRENE 2D-mesh for an axisymmetric tokamak equilibrium (g-file).
- Plot the 2D projection of the guiding-center orbits on the original SOLEDGE3X-EIRENE grid.
- Compute collisionless guiding-center orbit with GORILLA for a trapped Deuterium particle with adaptive scheme.
- Use a field-aligned grid for a non-axisymmetric VMEC MHD equilibrium.
- Create a figure with the Poincaré plots (
$\varphi$ = 0) in cylindrical and symmetry flux coordinates. - Compute the normalized parallel adiabatic invariant as a function of banana bounces.
- Plot fluctuation and evolution of energy over the bounces.
- Using a potential that is a scaled poloidal flux (ASDEX/SOLEDGE3X-mesh), calculate the corresponding electric field via central differences.
- Perform two runs, one with and one without the inclusion of the additional terms in the strong electric field Lagrangian using ionised Tungsten as particle species.
- Plot the poincare cross-section of both runs, their respective fluctuations of total energy and toroidal momentum. as well as the electric field and ExB-drift.
- The definition of total energy and toroidal momentum is different in the two runs. Each definition should yield conserved quantities in their respective case.
A detailed explanation of all examples (1-7) including the generation of the appropriate input files (including the example folders in EXAMPLES/MATLAB_RUN
and EXAMPLES/PYTHON_RUN
) and plotting of the results with MATLAB and Python can be found in the folders MATLAB
and PYTHON
, respectively.
Here, the results of GORILLA with different polynominal orders K=2,3,4 and Runge-Kutta 4 are compared in case of examples 1-3. For examples 5-6 orbits for both trapped and passing particles are calculated. For example 7 an additional, in-depth comparison between adaptive and non-adaptive scheme is performed. The last example, example 8, is currently only available via the MATLAB script example_8.m
. A redundant version in PYTHON, as well as a corresponding example-folder with an appropriate script for plotting will follow in the future.
Tests are implemented with pFUnit. The generation of coverage files are done with the compiler option --coverage
and the evaluation of the files is implemented with lcov. To see the coverage report, build with build_coverage.sh
and open index.html
in the folder /BUILD/COVERAGE or take a look in the github workflows Ubuntu
or Mac
.
The tests performed can be found in the folder SRC/TESTS
. Descriptions of the tests are given in the individiual source files alongside some background information on the module/part of the program covered by the tests. The checks focus on low level function testing, providing routines with trial/default inputs and comparing the output/behaviour against expected results. Further tests are planned to be implemented in future releases. Incomplete tests not being able to meet the outlined policy due to restrictions in the current program structure are disabled, but still included in SRC/TESTS
to be implemented at a later date.
If you have questions regarding any aspect of the software then please get in touch with the main developer Michael Eder via email - [email protected]. Alternatively you can create an issue on the Issue Tracker.
There may be bugs. If you think you've caught one, please report it on the Issue Tracker. This is also the place to propose new ideas for features or ask questions about the design of GORILLA.
We welcome help in improving and extending GORILLA. This is managed through Github pull requests; for external contributions we prefer the "fork and pull" workflow while core developers use branches in the main repository:
- First open an Issue to discuss the proposed contribution.
- Make your own project fork and implement the changes there.
- Open a pull request to merge the changes into the main project. A more detailed discussion can take place there before the changes are accepted.
- Michael Eder1
- Christopher G. Albert1
- Lukas M. P. Bauer1
- Georg S. Graßler1
- Daniel Forstenlechner1
- Sergei V. Kasilov1,2,3
- Winfried Kernbichler1
- Markus Meisterhofer1
- Michael Scheidt1
-
Fusion@OEAW, Institut für Theoretische Physik - Computational Physics, Technische Universität Graz, Petersgasse 16, 8010 Graz, Austria
-
Institute of Plasma Physics, National Science Center, “Kharkov Institute of Physics and Technology”, Akademicheskaya str. 1, 61108 Kharkov, Ukraine
-
Department of Applied Physics and Plasma Physics, V. N. Karazin Kharkov National University, Svobody sq. 4, 61022 Kharkov, Ukraine
The authors would like to thank Michael Drevlak, the ASDEX Upgrade Team and the SOLEDGE3X-EIRENE Team for respectively providing the stellarator field configuration, the ASDEX Upgrade MHD equilibrium for shot 26884 at 4300 ms and the WEST MHD equilibrium for shot 54903 at 8 s. Further thanks to Martin Heyn, Philipp Ulbl, Rico Buchholz, Patrick Lainer, and Markus Richter for useful discussions. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
When using this code for scientific publications, please cite both references [1] and [2]:
[1] M. Eder, C.G. Albert, L.M.P. Bauer, S.V. Kasilov and W. Kernbichler “Quasi-geometric integration of guiding-center orbits in piecewise linear toroidal fields” Physics of Plasmas 27, 122508 (2020) https://doi.org/10.1063/5.0022117 Preprint: https://arxiv.org/abs/2007.08151
[2] M. Eder, C.G. Albert, L.M.P. Bauer, Georg S. Graßler, S.V. Kasilov, W. Kernbichler, M. Meisterhofer and M.Scheidt “GORILLA: Guiding-center ORbit Integration with Local Linearization Approach” submitted to Journal of Open Source Software Preprint: https://github.com/openjournals/joss-papers/blob/joss.03116/joss.03116/10.21105.joss.03116.pdf
[3] M. F. Heyn, I. B. Ivanov, S. V. Kasilov, W. Kernbichler, P. Leitner, and V. Nemov, “Quasilinear modelling of RMP interaction with a tokamak plasma: application to ASDEX Upgrade ELM mitigation experiments” Nucl. Fusion 54, 064005 (2014). https://doi.org/10.1088/0029-5515/54/6/064005
[4] M. Drevlak, C. D. Beidler, J. Geiger, P. Helander, and Y. Turkin “Quasi-Isodynamic Configuration with Improved Confinement” 41st EPS Conference on Plasma Physics ECA (2014), Vol. 38F, p. P1.070. http:https://ocs.ciemat.es/EPS2014PAP/pdf/P1.070.pdf
[5] H. Bufferand et al “Progress in edge plasma turbulence modelling—hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry” Nucl. Fusion 61, 116052 (2021) https://doi.org/10.1088/1741-4326/ac2873