Bringing 3D GVEC equilibria to STRUPHY

STRUcture-Preserving HYbrid code (STRUPHY) is a code base for kinetic-fluid plasma simulations developed at the Max Planck Institute for Plasma Physics in Garching. In this package we add to STRUPHY several general Python interfaces for reading magneto-hydrodynamic (MHD) equilibrium data for Tokamaks and Stellarators. In particular, interfaces to the equilibrium code Galerkin Variational Equilibrium Code (GVEC) (and in the future, to the eqdsk file format) are provided. From the data, a flux-aligned coordinate system is created, and the MHD quantities are transformed into differential forms for the use in STRUPHY.

These new features are tested by a variety of benchmark simulations, comparing STRUPHY to well-established hybrid codes, such as HMGC and MEGA. The focus of the benchmarks is on resonant excitation of Alfvén eigenmodes by energetic particles in Tokamak equilibria with Shafranov shift.

Documentation: GitLab Pages (Not yet accessible)

How to use

First and foremost, one need to setup the package. It can be accomplished as simple as a pip install -e ./. See Section "Build and install" below for detailed instructions.

This code needs to read in a GVEC MHD equilibrium in JSON format. Therefore, the .dat GVEC output data has to be converted into .json using a utility as demonstrated below. This only needs to be done once, per GVEC equilibrium.

from  gvec_to_python.reader.gvec_reader import GVEC_Reader

# Spline coefficients are not provided by GVEC by default.
read_filepath = 'GVEC/testcases/circ_tok/'
read_filename = 'CIRC_TOK_State_0000_00010000.dat'
save_filepath = 'GVEC/testcases/circ_tok/'
save_filename = 'CIRC_TOK_State_0000_00010000.json'
reader = GVEC_Reader(read_filepath, read_filename, save_filepath, save_filename, with_spl_coef=False)

# Special case with spline coefficients.
read_filepath = 'GVEC/testcases/ellipstell/'
read_filename = 'GVEC_ellipStell_profile_update_State_0000_00010000.dat'
save_filepath = 'GVEC/testcases/ellipstell/'
save_filename = 'GVEC_ellipStell_profile_update_State_0000_00010000.json'
reader = GVEC_Reader(read_filepath, read_filename, save_filepath, save_filename, with_spl_coef=True)

The main features can then be used simply as follows:

import numpy as np
from gvec_to_python import GVEC, Form, Variable

gvec_data_filepath = 'GVEC/testcases/ellipstell/'
gvec_data_filename = 'GVEC_ellipStell_profile_update_State_0000_00010000.json'
gvec = GVEC(gvec_data_filepath, gvec_data_filename)

# The functions can accept various types of inputs:
# - Individual point in logical coordinates.
(s, u, v) = (0.1, 0.2, 0.3)
# - 1D arrays of a grid.
s_range = np.linspace(0, 1, 11)
u_range = np.linspace(0, 1, 21)
v_range = np.linspace(0, 1, 31)
# - A sparse 3D meshgrid.
s_sparse, u_sparse, v_sparse = np.meshgrid(s_range, u_range, v_range, indexing='ij', sparse=True)
# - A dense 3D meshgrid.
s_dense,  u_dense,  v_dense  = np.meshgrid(s_range, u_range, v_range, indexing='ij', sparse=False)

# Function parameters can either be floats gvec.F(s, u, v), 
# or 1D arrays gvec.F(s_range, u_range, v_range), 
# or 3D arrays from sparse meshgrid gvec.F(s_sparse, u_sparse, v_sparse), 
# or dense meshgrid gvec.F(s_dense, u_dense, v_dense).
(x, y, z) = gvec.F(s, u, v)
jacobian  = gvec.DF(s, u, v)
pressure  = gvec.get_variable(s, u, v, variable=Variable.PRESSURE, form=Form.ZERO)
pressure  = gvec.P(s, u, v) # Alternatively, a shorthand is available.
p_3       = gvec.P_3(s, u, v) # 3-form pressure.

Caveat: The mapping function F() returns in the shape (3, #s, #u, #v) so that we can do (x, y, z) = gvec.F(s, u, v), but the Jacobian function DF() returns in the shape (#s, #u, #v, 3, 3) in order to comply with numpy's convention.

Under the hood, each step of the workflow is demonstrated in various Jupyter Notebooks at the repository root. Examples include:

• parsing GVEC output into .json format;
• transforming GVEC basis (X^1, X^2, \zeta) to Cartesian (x, y, z);
• a SymPy derivation of the Jacobian of the transform mapping;
• converting the output data to ParaView's .vtk format.

and so on.

The test files are also helpful in understanding low-level usage, without the clutter of experiment features in the Notebooks.

For now, only torus mapping is implemented (see suv_to_xyz.py), where (x, y, z) = (X^1 cos(\zeta), - X^1 sin(\zeta), X^2). If one wishes to extend to other forms of mapping, such as a periodic cylinder, where (x, y, z) = (X^1, - \frac{L_c}{2 \pi} \zeta, X^2), this is where the modification should be done.

Directory oragnization

The code is arranged in the following manner:

• ./GVEC
  Gvec outputs in various file formats.
• ./literature
  As titled.
• ./logs
  An autocreated directory that stores the execution log. Can be deleted after use.
• ./src/gvec_to_python/base
  B-spline (radial) \times Fourier (angular) basis of GVEC.
  Depends on hylife.
• ./src/gvec_to_python/hmap
  Mapping from logical coordinates (s, u, v), to GVEC basis (X^1, X^2, \xi), to Cartesian coordinates (x, y, z).
  Depends on base.
• ./src/gvec_to_python/hylife
  Code "borrowed" from STRUPHY code base, primarily for constructing 1D B-splines basis for GVEC's coordinate system.
• ./src/gvec_to_python/reader
  Parse GVEC output into .json, and read that .json.
• ./src/gvec_to_python/util
  Well, utils.
• ./src/gvec_to_python/writer
  Writing GVEC output to ParaView .vtk file.
• ./src/gvec_to_python/GVEC_functions.py
  Main class GVEC which combines all functionalities above.
• ./test
  Test functions used in PyTest.
  Run these from the project root.
• ./*.ipynb
  Notebooks used in development, but also demonstrate the workflow of using this code.
• ./.coveragerc
  Configuration for test coverage report.
• ./.gitlab-ci.yml
  GitLab CI (Continuous Integration).
• ./build_and_test.sh
  Emulate part of the GitLab CI process locally.
• ./pyproject.toml, ./setup.cfg
  Python package build settings.

Build and install

The simplest way to use this code in your own project is to git clone the repo, then run pip install -e . from the project root.

To build the project into an installable package, one needs to structure the package inside a src folder, which is already done in this repository. pyproject.toml and setup.cfg are kept at the top-level parent:

[ARBITRARY_NAME]/  
├── LICENSE
├── pyproject.toml  
├── README.md
├── setup.cfg  
├── src/  
│   └── gvec_to_python/  
│       ├── __init__.py
│       ├── ......
└── tests/

From here, one can build the project:
python -m pip install -U build
python -m build

The output will be like this:

[ARBITRARY_NAME]/  
└── dist/
    ├── gvec_to_python-0.1.0-py3-none-any.whl
    └── gvec_to_python-0.1.0.tar.gz

One can then pip install this wheel file directly without uploading the package to PyPI.
i.e. pip install -U dist/gvec_to_python-0.1.0-py3-none-any.whl --force-reinstall
The latest wheel file can be downloaded from the Artifacts on GitLab's Pipeline.
Afterwards, this code can be imported and used like others.

Alternatively, one can clone this source from git and do a dev install with pip install -e ..

GVEC code information

• When browsing the GVEC repository, always look at the branch develop.
• File readin: in src/readstate/readstate.f90 with all container variables defined in src/readstate/readstate_vars.f90
• File output: in src/restart/restart.f90 subroutine WriteStateToASCII
• Preparation of visualization data is found in src/functionals/mhd3d/mhd3d_visu.f90.
  • EXCELLENT EXAMPLE to reference, for the evaluation of the field and other quantities on a (s,\theta,\zeta) grid.
  • Generated example of paraview visualization files can be found in in GVEC/testcases/ellipstell/*.vtu.
• An example for hmap can be found in src/functionals/mhd3d/hmap_RZ.f90 (for cylinder hmap, length of cylinder would be needed).
  • FUNCTION hmap_RZ_eval(sf, q_in)RESULT(x_out) converts coordinate q_in = (R,Z,\zeta) = (X^1,X^2,\zeta) into x_out = (x,y,z).
    • wp stands for working precision: single or double.
  • FUNCTION hmap_RZ_eval_dxdq(sf, q_in, q_vec) RESULT(dxdq_qvec) evaluates total derivatives at q_in = (R,Z,\zeta) = (X^1,X^2,\zeta).
  • FUNCTION hmap_RZ_eval_Jh(sf, q_in) RESULT(Jh) evaluates the Jacobian J_h.
    • Further derivatives of this function, with hmap_RZ_eval_Jh_dq1 and hmap_RZ_eval_Jh_dq2, are not relevant for our purpose.
  • Note FUNCTION hmap_RZ_eval_gij(sf, qL_in, q_G, qR_in) RESULT(g_ab) in case we will need it.
• See also src/gvec_to_hopr/gvec_to_hopr.f90 or src/gvec_to_gene/gvec_to_gene.f90 for examples on evaluating profiles and coordinates.
  • Note SUBROUTINE gvec_to_hopr_SFL() for examples on how to evaluate derivatives and quantities such as \phi,\chi,\iota.
• Fourier bases can be found in src/base/fbase.f90.
  • Similarly, the radial basis can be found in src/base/sbase.f90.
  • The coefficients are stored as a single array in the .dat output file. In the case that the _sincos_ option is specified, the number of modes for _sin_ and _cos_ coefficients is not a known from the output file.
    • Thus, after read-in of the coefficients along with their m,n mode number, the m_max,n_max must be inferred and:
    • modes_sin=(n_max)+m_max*(2*n_max+1) , modes_cos=(n_max+1)-mn_excl+m_max*(2*n_max+1)
    • with mn_excl=0/1 if the cosine m=n=0mode is included/excluded.

GVEC output file

• ASCII format
  • Example file of a simple 3D stellarator case, run with GVEC for 10000 iterations:
    • GVEC/testcases/ellipstell/GVEC_ellipStell_profile_update_State_0000_00010000.dat where relevant variables are marked with (*).
• Sections are separated with a headerline starting with ##
• List of data entries given in List_of_data_entries.md
• Content:
  • 1D grid points in s-direction, can be non-uniform (npoints=nElems+1) (=knots of the spline, without their repetition at the boundary, since they are clamped).
  • Global information (which nfp,hmap).
  • 3D solution variables (X^1,X^2,\lambda(s,\theta,\zeta)).
    • Fourier coefficients x B-spline coefficients.
    • ALL variables use the same radial grid.
    • EACH variable has its own spline basis in (s). (defined by polynomial degree (continuity must be degree-1!), knots from radial 1D grid) and its own Fourier basis in (theta,zeta) (number of modes in m_max,n_max , cos or sin or cos&sin ,NFP...)
  • 1D profiles interpolated at the Greville points of the spline basis of X1 (currently the spline basis of X1 is used, but that might have to be decoupled in the future, since the profiles might need at different radial resolution/degree ...)

Numba vs Pyccel

Runtime comparison of pytest as of commit 48457d88:

To reduce code complexity, support for pyccel is dropped in favour of numba, even though it is slightly faster.