Coordinates

Most of the fields in this package are not most naturally written down in cartesian coordinates. A tokamak equilibrium wants cylindrical or toroidal coordinates, a Solov'ev equilibrium wants cylindrical coordinates normalised to the major radius, and only the simplest fields are genuinely cartesian. ElectromagneticFields.jl therefore lets every equilibrium carry its own chart, and generates all evaluation routines in the coordinates of that chart.

The price of that freedom is that a vector no longer has just three numbers. In a curvilinear chart the coordinate basis is neither orthogonal nor normalised, so there are several inequivalent ways to write the components of one and the same vector. This page describes the three that the package uses, and how they show up in the generated code.

Charts

A chart is a map from curvilinear coordinates $\xi$ to the cartesian coordinates $x$ of the ambient space,

\[F : \xi = (\xi^1, \xi^2, \xi^3) \mapsto x = (x^1, x^2, x^3) ,\]

together with its inverse. An equilibrium declares its chart by overloading the traits x¹, x², x³ for $F$ and ξ¹, ξ², ξ³ for $F^{-1}$, plus the metric coefficients g₁₁ … g₃₃, the volume element J, the orientation, and optionally the bounds minx¹ … maxx³ of the coordinate domain. Everything not overloaded falls back to a cartesian default: identity metric, infinite range, positive orientation.

Five charts are in use:

Note that SolovevSymmetric is a cartesian equilibrium despite its name, and that the Solov'ev chart coordinates are dimensionless, which is where the factors of $R_0^2$ in its metric come from.

The examples on this page all use the cylindrical tokamak, which is a good representative: it is curvilinear, its metric is non-trivial, and it is left-handed.

using ElectromagneticFields
using LinearAlgebra

field = FieldFunctions(AxisymmetricTokamakCylindricalEquilibrium(6.2, 5.3, 2.0))

t = 0.0
ξ = [6.5, 0.5, 0.25]

The chart itself is generated as to_cartesian and from_cartesian, each returning the three components as an SVector:

to_cartesian(field, t, ξ)

3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 6.297930741119191
 1.6081257351543992
 0.5

roundtrip = from_cartesian(field, t, to_cartesian(field, t, ξ)) ≈ ξ
roundtrip

true

The bounds of the coordinate domain come as rangemin and rangemax. For this chart only the toroidal angle is bounded:

rangemin(field, t, ξ), rangemax(field, t, ξ)

([-Inf, -Inf, 0.0], [Inf, Inf, 6.283185307179586])

Tangent Map and Metric

Differentiating the chart gives the tangent map $DF$, and the metric is the pullback of the cartesian one along it:

\[{DF^i}_j = \frac{\partial x^i}{\partial \xi^j} , \qquad g_{ij} = \sum_k {DF^k}_i \, {DF^k}_j , \qquad g^{ij} = (g^{-1})_{ij} .\]

Both are generated as matrix-valued functions: DF and DF̄ for the tangent map and its inverse, g♭ and g♯ for the metric and its inverse. The metric follows the musical isomorphisms — ♭ lowers an index, ♯ raises one — while the inverse of the tangent map carries an overbar, its components being the derivatives of the inverse chart,

\[{\bar{DF}^i}_j = \frac{\partial \xi^i}{\partial x^j} .\]

DF(field, t, ξ)

3×3 StaticArraysCore.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 0.968912  0.0  -1.60813
 0.247404  0.0   6.29793
 0.0       1.0   0.0

The three identities relating them hold pointwise:

checks = (DF̄(field, t, ξ) ≈ inv(DF(field, t, ξ)),
          DF(field, t, ξ)' * DF(field, t, ξ) ≈ g♭(field, t, ξ),
          g♯(field, t, ξ) ≈ inv(g♭(field, t, ξ)))
checks

(true, true, true)

Derivatives of the metric are generated as well, Dg♭ and Dg♯ for the first and DDg♭ and DDg♯ for the second, which is what a geometric integrator needs to assemble Christoffel symbols. They are rank-three and rank-four SArrays, so Dg♭(field, t, ξ)[i,j,k] is $\partial_k g_{ij}$.

Covariant, Contravariant and Physical Components

A vector field and a one-form are different objects, but in the presence of a metric they carry the same information, and the package generates both representations of every field it computes. The metric raises and lowers the indices:

\[v_i = g_{ij} \, v^j , \qquad v^i = g^{ij} \, v_j .\]

Components with a lower index are called covariant, components with an upper index contravariant. In the code the two conversions are one-liners, covariant_to_contravariant and contravariant_to_covariant, and the naming follows the musical isomorphisms that perform them: B♭ is covariant, B♯ contravariant.

Neither of them is what a measurement returns. Covariant and contravariant components are taken with respect to the coordinate basis and its dual, and in a curvilinear chart those basis vectors are neither unit length nor mutually orthogonal — the covariant toroidal component of $B$ below is larger than $|B|$ by nearly a factor of $R$, and carries different units than the other two components. What is usually quoted in the literature, and what this package calls the physical components, are the components in an orthonormal frame:

\[v_{(i)} = {DF^i}_j \, v^j = {\bar{DF}^j}_i \, v_j .\]

That is, physical components are obtained by pushing the vector forward to the ambient cartesian frame. They all carry the same units, and their euclidean norm is the length of the vector. They are written with the musical natural, B♮, the sign that cancels a flat or a sharp.

Some texts define physical components instead by normalising the coordinate basis vectors, $v_{\langle i \rangle} = \sqrt{g_{ii}} \, v^i$. For an orthogonal chart that is also an orthonormal frame, so it gives the same magnitude, but its components differ from the ones here by the rotation relating the local frame to the cartesian one.

The unit vector along the magnetic field is available in all three representations, as b♭ (covariant), b♯ (contravariant) and b♮ (physical):

[b♭(field, t, ξ) b♯(field, t, ξ) b♮(field, t, ξ)]

3×3 StaticArraysCore.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.0402781  -0.0402781  -0.286157
  0.0241668   0.0241668   0.957878
  6.49283     0.153676    0.0241668

The three columns are the same vector. Lowering, raising and pushing forward take one into another:

conversions = (g♭(field, t, ξ) * b♯(field, t, ξ) ≈ b♭(field, t, ξ),
               g♯(field, t, ξ) * b♭(field, t, ξ) ≈ b♯(field, t, ξ),
               b♮(field, t, ξ) ≈ DF(field, t, ξ) * b♯(field, t, ξ),
               b♮(field, t, ξ) ≈ DF̄(field, t, ξ)' * b♭(field, t, ξ))
conversions

(true, true, true, true)

Only two of the three give the length of the vector directly. Contracting a covariant with a contravariant vector is a metric-free operation and needs no $g$; the same contraction of the physical components is an ordinary euclidean dot product. Contracting covariant with covariant, on the other hand, is meaningless.

b♯(field, t, ξ)' * b♭(field, t, ξ), norm(b♮(field, t, ξ))

(1.0000000000000002, 1.0)

Summarising the notation:

Volume Element and Orientation

J is the volume element $\sqrt{|g|} = |\det DF|$. It is not computed from the metric but declared by the chart, so that the generated expressions stay in the closed form the chart author intended rather than passing through a symbolic square root.

Being an absolute value, J says nothing about handedness, and handedness matters: the Hodge star that turns the magnetic two-form into a one-form, and the cross product used to build the perpendicular frame, both need the signed determinant

\[\det DF = \mathrm{orientation} \cdot J .\]

Four of the five charts above are left-handed, because the toroidal angle sits in the third slot where the right-handed ordering would put the second poloidal coordinate — $(R, Z, \phi)$ rather than $(R, \phi, Z)$. This is not a corner case, and getting it wrong has no visible symptom other than a magnetic field pointing the wrong way. Each equilibrium therefore declares its handedness explicitly, and the field carries it as orientation(field) — a stored value rather than a generated function, since the sign depends on neither time nor position. See orientation for the full discussion.

orientation(field), J(field, t, ξ), det(DF(field, t, ξ))

(-1, 6.5, -6.5)

signed = det(DF(field, t, ξ)) ≈ orientation(field) * J(field, t, ξ)
signed

true

Adding a Chart

A new chart is defined by overloading, for the new equilibrium type,

• x¹, x², x³ — the chart map to cartesian coordinates,
• ξ¹, ξ², ξ³ — its inverse,
• g₁₁ … g₃₃ — the metric coefficients, in the chart's own coordinates,
• J — the volume element,
• orientation — +1 or -1, and
• optionally minx¹ … maxx³ — the bounds of the coordinate domain.

Everything else — the inverse metric, the tangent map and its inverse, all derivatives, and the whole tower of field representations described in Fields — is derived symbolically from these by the code generator. A perturbation added to an equilibrium must live on the same chart; the generator asserts that the two agree on J and g.

The consistency of a new chart is checked by the test suite, which asserts $J = \sqrt{\det(DF^T DF)}$ and $\det DF = \mathrm{orientation} \cdot J$ for every equilibrium, so a chart whose metric, volume element and handedness do not fit together fails immediately.