Geometry (mapped domains)

In Struphy, a computational domain is defined by a mapping from logical coordinates to physical space. Concretely, we use a map

\[F:(0,1)^3 \rightarrow \Omega,\qquad \boldsymbol{\eta}=(\eta_1,\eta_2,\eta_3) \mapsto \mathbf{x}=(x,y,z),\]

where the logical unit cube is easy to mesh, and the mapping bends this mesh into the physical geometry of interest.

This is why mapped domains are so useful: we can keep structured discretizations in logical space while representing curved geometries (cylinders, tori, tokamaks, stellarators) in physical space.

A practical challenge is the polar singularity (the magnetic axis). In Struphy, two strategies are common:

1. Use polar splines in the de Rham setup (more expensive).

2. Cut out a small hole around the axis (cheap, but particles cannot live in that excluded region).

In this tutorial we focus on the second strategy and build intuition from simple to more physics-driven mappings.

HollowCylinder

We start with HollowCylinder, the simplest curvilinear example: a cylindrical shell obtained from the logical cube by an analytical mapping.

from struphy import domains

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

A 3D rendering helps to connect the mapping parameters to geometry: inner radius, outer radius, and length become visually obvious.

domain.show_3d()

2026-06-16 08:07:34.978 (   1.485s) [    7FBCBDDA4B80]vtkXOpenGLRenderWindow.:1460  WARN| bad X server connection. DISPLAY=
2026-06-16 08:07:34.979 (   1.485s) [    7FBCBDDA4B80]vtkOpenGLRenderWindow.c:645   WARN| Failed to load EGL! Please install the EGL library from your distribution's package manager.
2026-06-16 08:07:36.988 (   3.495s) [    7FBCBDDA4B80]vtkOpenGLRenderWindow.c:645   WARN| Failed to load EGL! Please install the EGL library from your distribution's package manager.

Let us inspect the default parameters of HollowCylinder. Read them as: a1 = inner radius, a2 = outer radius, Lz = axial length, and poc controls angular periodicity (“piece of cake”).

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

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

poc lets us model only a sector of the full cylinder and reconstruct periodicity through the mapping. This is useful when symmetry allows reducing the simulated angular extent.

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

Some key attributes to interpret are:

• kind_map: internal mapping family identifier (20-29 are cylinder/torus analytical mappings).

• pole: True when the mapping includes a polar singularity on axis.

• periodic_eta3: tells whether the logical direction \(\eta_3\) is periodic in this mapping.

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

20
False
False

It is also useful to inspect callable methods of a domain object: this shows what geometric operations (mapping, Jacobian, metric-related helpers, visualization, etc.) are available for downstream workflows.

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

The domain object is itself callable: evaluating domain(...) applies the mapping \(F(\eta_1,\eta_2,\eta_3)=(x,y,z)\) from logical to physical coordinates.

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.

Now we vary a1 to control how much of the axis neighborhood is removed. Increasing a1 enlarges the excluded core and avoids singular behavior at the center.

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

If we set a1=0, the hole disappears and the axis is included. In that case, domain.pole becomes True, signaling the polar singularity discussed in the introduction.

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

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

20
True
False

HollowTorus

HollowTorus extends the same idea to toroidal geometry: the logical cube is mapped to a torus with optional inner hole around the magnetic axis.

Compared to HollowCylinder, important new parameters are:

• R0: major radius (distance from torus center to tube center),

• sfl: whether to use straight-field-line poloidal parametrization,

• tor_period: built-in toroidal periodicity (e.g. 3 means one third of a full torus in the mapped geometry).

Let us create HollowTorus with default parameters:

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

For this mapping, periodic_eta3 is True because the toroidal angle is encoded in \(\eta_3\) and treated periodically by construction.

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

22
False
True

Let us inspect the defaults. In particular, compare a1/a2 (minor radii), R0 (major radius), and periodicity-related parameters to understand the torus shape and angular parametrization choices.

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

Now we modify three knobs at once:

• a1 to change the hole around the axis,

• sfl=True to switch poloidal coordinates to straight-field-line style,

• tor_period=1 to represent the full torus period in the map.

This helps build intuition for how coordinate choices alter the geometry and periodicity.

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

Tokamak

For realistic tokamak geometry, Struphy uses mappings tied to an axisymmetric MHD equilibrium. Instead of prescribing a simple analytic torus, the domain is built from a poloidal flux function \(\psi\) and field-line tracing.

So conceptually, the workflow is: equilibrium model \(\rightarrow\) flux surfaces \(\rightarrow\) mapped computational domain.

Tokamak is the domain class implementing this idea.

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 0x7fbbfadac760>, no rescaling performed in EQDSK output.
  warnings.warn(

domain.show_3d()

2026-06-16 08:07:44.244 (  10.751s) [    7FBCBDDA4B80]vtkOpenGLRenderWindow.c:645   WARN| Failed to load EGL! Please install the EGL library from your distribution's package manager.

Let us inspect the default tokamak mapping parameters. Notice how they now include options linked to equilibrium/flux-surface construction rather than only simple geometric radii.

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

A Tokamak domain is always coupled to an AxisymmMHDequilibrium, which provides the flux function \(\psi\) used to construct the map.

By default, this is AdhocTorus, a convenient analytic equilibrium for testing. We can also switch to EQDSKequilibrium, which represents a more realistic equilibrium source.

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 0x7fbcb0e8b6d0>, no rescaling performed in EQDSK output.
  warnings.warn(

Finally, we shrink the inner excluded region by adjusting psi_shifts, which changes how close the mapped domain gets to the magnetic axis.

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

Stellarator mappings

Tokamak mappings above are axisymmetric. Stellarators, in contrast, are intrinsically 3D and typically require equilibrium data generated by dedicated external codes.

Struphy provides interfaces to read equilibrium geometry from:

• GVEC equilibrium code

• DESC equilibrium code

GVEC interface

The corresponding Struphy domain class is 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 0x7fbbfaac9330>, no rescaling performed in GVEC output.
  warnings.warn(

Inspecting defaults is useful here as well: they describe both geometry controls and how the external equilibrium information is interpreted by the mapping.

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)

Next, we modify equilibrium-side parameters to open a hole around the magnetic axis (rmin) and use the whole stellarator (use_nfp=False). These settings are passed through 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 0x7fbbec9bc8b0>, no rescaling performed in GVEC output.
  warnings.warn(

DESC interface

The analogous interface for DESC is DESCunit. As with GVEC, equilibrium choices are specified through DESCequilibrium, then consumed by the domain mapping:

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

domain.show()