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()

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)

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)

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)

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__)

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)

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)

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)

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()

[ ]:

domain.show_3d()

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)

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()

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 interface#

The corresponding Struphy domain class is GVECunit. Let us create an instance with default parameters:

[ ]:

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

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()

DESC interface#

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

[ ]:

import contextlib
import io

# DESC is verbose while loading: silence its output, but let errors through.
with contextlib.redirect_stdout(io.StringIO()), contextlib.redirect_stderr(io.StringIO()):
    desc_equil = equils.DESCequilibrium(use_nfp=False)
    domain = domains.DESCunit(desc_equil)

[ ]:

domain.show()