6 - Mapped domains with polar singularity

This tutorial provides access to Struphy domains which can be defined from analytical formulas or through third-party software (VMEC, GVEC, DESC, etc.).

Of particular interest are maps that have a polar singularity (“magnetic axis”). Such singularities are challenging numerically; in Struphy two solutions are possible:

1. Use polar splines when setting up the de Rham sequence (=more expensive)

2. Cut out a small hole of the domain around the singularity (=cheap, but problematic for particles).

We shall focus on the second possibility in this notebook.

HollowCylinder

Let us create the domain HollowCylinder with default parameters:

from struphy import domains

domain = domains.HollowCylinder()
domain.show()

A 3D view of the domain can be displayed with

domain.show_3d()

2026-05-08 07:43:45.231 (   1.781s) [    7FF821BAEB80]vtkXOpenGLRenderWindow.:1460  WARN| bad X server connection. DISPLAY=
2026-05-08 07:43:45.231 (   1.782s) [    7FF821BAEB80]vtkOpenGLRenderWindow.c:645   WARN| Failed to load EGL! Please install the EGL library from your distribution's package manager.
2026-05-08 07:43:46.460 (   3.010s) [    7FF821BAEB80]vtkOpenGLRenderWindow.c:645   WARN| Failed to load EGL! Please install the EGL library from your distribution's package manager.

The default parameters of HollowCylinder are:

for key, val in domain.params.items():
    print(key, "=", val)

a1 = 0.2
a2 = 1.0
Lz = 4.0
poc = 1

We can use the piece-of-cake parameter poc to use only a part of the disk and complete it periodically:

domain_poc = domains.HollowCylinder(poc=3)
domain_poc.show()

Some relevant domain attributes are:

print(domain.kind_map)
print(domain.pole)
print(domain.periodic_eta3)

20
False
False

The domain methods are also quite important:

for attr in dir(domain):
    if callable(getattr(domain, attr)) and "__" not in attr and attr[0] != "_":
        print(attr)

create_geometry_mesh
export_geometry
from_dict
get_params_numpy
jacobian
jacobian_det
jacobian_inv
metric
metric_inv
prepare_arg
prepare_eval_pts
pull
push
show
show_3d
to_dict
transform

Aside from these, the domain object itself is callable:

help(domain.__call__)

Help on method __call__ in module struphy.geometry.base:

__call__(*etas, change_out_order=False, squeeze_out=False, remove_outside=True, identity_map=False) method of struphy.geometry.domains.HollowCylinder instance
    Evaluates the mapping :math:`F : (0, 1)^3 \to \mathbb R^3,\, \boldsymbol \eta \mapsto \mathbf x`.

    Logical coordinates outside of :math:`(0, 1)^3` are evaluated to -1.
    The type of evaluation depends on the shape of the input ``etas``.

    Parameters
    ----------
    *etas : array-like | tuple
        Logical coordinates at which to evaluate. Two cases are possible:
            1. 2d numpy array, where coordinates are taken from eta1 = etas[:, 0], eta2 = etas[:, 1], etc. (like markers).
            2. list/tuple (eta1, eta2, ...), where eta1, eta2, ... can be float or array-like of various shapes.

    change_out_order : bool
        If True, the axis corresponding to x, y, z coordinates in the output array is the last one, otherwise the first one.

    squeeze_out : bool
        Whether to remove singleton dimensions in output array.

    remove_outside : bool
        If True, logical coordinates outside of (0, 1)^3 are NOT evaluated to -1 and are removed in the output array.

    identity_map : bool
        If True, not the mapping F, but the identity map (0, 1)^3 --> (0, 1)^3 is evaluated

    Returns
    -------
    out : ndarray | float
        The Cartesian coordinates corresponding to the given logical ones.

Let us change the size of the hole around the pole:

domain = domains.HollowCylinder(a1=0.05)
domain.show()

Note that if we set the inner radius a1 to zero, the attribute domain.pole becomes True:

domain = domains.HollowCylinder(a1=0.0)
domain.show()

print(domain.kind_map)
print(domain.pole)
print(domain.periodic_eta3)

20
True
False

HollowTorus

Let us create the domain HollowTorus with default parameters:

domain = domains.HollowTorus()
domain.show()

Note that the attribute periodic_eta3 is True for this mapping:

print(domain.kind_map)
print(domain.pole)
print(domain.periodic_eta3)

22
False
True

The default parameters are:

for key, val in domain.params.items():
    print(key, "=", val)

a1 = 0.1
a2 = 1.0
R0 = 3.0
sfl = False
pol_period = 1
tor_period = 3

Let us change the size of the hole around the pole, the poloidal angle parametrization (sfl) and the tor_period:

domain = domains.HollowTorus(a1=0.05, sfl=True, tor_period=1)
domain.show()

Tokamak

Tokamak is the class for mappings for Tokamak MHD equilibria constructed via field-line tracing of a poloidal flux function \(\psi\).

Let us create a Tokamak with default parameters:

domain = domains.Tokamak()
domain.show()

/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/struphy/fields_background/equils.py:1766: UserWarning: self.units =<struphy.physics.physics.Units object at 0x7ff75a9e8550>, no rescaling performed in EQDSK output.
  warnings.warn(

domain.show_3d()

2026-05-08 07:43:53.269 (   9.819s) [    7FF821BAEB80]vtkOpenGLRenderWindow.c:645   WARN| Failed to load EGL! Please install the EGL library from your distribution's package manager.

The default parameters are:

for key, val in domain.params.items():
    if "cx" not in key and "cy" not in key:
        print(key, "=", val)

equilibrium = EQDSKequilibrium
    rel_path:      True
    file:          None
    data_type:     0
    degree_for_psi:(3, 3)
    psi_resolution:(25.0, 6.25)
    degree_for_flux:3
    flux_resolution:50.0
    n1:            2.0
    n2:            1.0
    na:            0.2
    base_units:    None
num_elements = (8, 32)
degree = (2, 3)
psi_power = 0.75
psi_shifts = (0.01, 2.0)
xi_param = equal_angle
r0 = 0.3
num_elements_pre = (64, 256)
p_pre = (3, 3)
tor_period = 1

The Tokamak domain is always related to an AxisymmMHDequilibrium, which provides the flux function \(\psi\). In the default parameters this is AdhocTorus. Instead, we could also look at the default EQDSKequilibrium:

from struphy import equils

mhd_eq = equils.EQDSKequilibrium()

domain = domains.Tokamak(equilibrium=mhd_eq)
domain.show()

/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/struphy/fields_background/equils.py:1766: UserWarning: self.units =<struphy.physics.physics.Units object at 0x7ff7801d0070>, no rescaling performed in EQDSK output.
  warnings.warn(

Let us shrink the hole:

domain = domains.Tokamak(equilibrium=mhd_eq, psi_shifts=[0.2, 2])
domain.show()

Stellarator mappings

Struphy can read data produced by

• GVEC equilibrium code

• DESC equilibrium code

GVEC interface

The interface is the class GVECunit.

Let us create an instance with default parameters:

domain = domains.GVECunit()
domain.show()

/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/struphy/fields_background/equils.py:2136: UserWarning: self.units =<struphy.physics.physics.Units object at 0x7ff75501b220>, no rescaling performed in GVEC output.
  warnings.warn(

The default parameters are:

for key, val in domain.params.items():
    if "cx" not in key and "cy" not in key and "cz" not in key:
        print(key, "=", val)

Let us put a domain hole around the magnetic axis and use the whole Stellarator (use_nfp=False). The parameters must be passed through the GVECequilibrium:

gvec_equil = equils.GVECequilibrium(rmin=0.1, use_nfp=False)
domain = domains.GVECunit(gvec_equil)
domain.show()

/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/struphy/fields_background/equils.py:2136: UserWarning: self.units =<struphy.physics.physics.Units object at 0x7ff74bf09ab0>, no rescaling performed in GVEC output.
  warnings.warn(

DESC interface

The interface is the class DESCunit. The parameters must be passed through the DESCequilibrium:

%%capture
desc_equil = equils.DESCequilibrium(use_nfp=False)
domain = domains.DESCunit(desc_equil)

domain.show()