Usage

Working with ElectromagneticFields.jl has two steps. First an equilibrium is constructed, which is a small immutable object holding nothing but the parameters of the field. Then the code for evaluating that field is generated from it into a FieldFunctions object, and that is what the actual computation uses.

Constructing an Equilibrium

Every field is a type with sensible defaults, so the shortest way to get an equilibrium is

using ElectromagneticFields

equ = AxisymmetricTokamakCylindricalEquilibrium()

Axisymmetric Tokamak Equilibrium in (R,Z,ϕ) Coordinates with
  R₀ = 1.0
  B₀ = 1.0
  q₀ = 2.0

All parameters can be passed explicitly, here the major radius $R_0$, the magnetic field strength $B_0$ at the magnetic axis, and the safety factor $q_0$:

equ = AxisymmetricTokamakCylindricalEquilibrium(6.2, 5.3, 2.0)

Axisymmetric Tokamak Equilibrium in (R,Z,ϕ) Coordinates with
  R₀ = 6.2
  B₀ = 5.3
  q₀ = 2.0

Some fields ship named configurations as well, e.g. the Solov'ev equilibrium comes with parameter sets for ITER, NSTX and a field reversed configuration:

SolovevEquilibriumITER()

SolovevEquilibrium Equilibrium with
  R₀ = 6.2
  B₀ = 5.3
  ϵ  = 0.32
  κ  = 1.7
  δ  = 0.33
  α  = -0.155

Generating the Evaluation Code

An equilibrium object on its own cannot be evaluated. The evaluation routines are traced symbolically from its chart, metric and vector potential, differentiated where necessary, and compiled. FieldFunctions does all of that and returns the result as a value:

field = FieldFunctions(equ)

This is an ordinary object. It can be built inside a function, stored in a Dict, or passed to a vector field that depends on the field — typically through a GeometricEquations problem's parameters, from which the equation function evaluates the components it needs.

Every quantity is reached through a generic that takes the field first. They all take the time as their second argument, followed by the coordinates of the evaluation point, which for this equilibrium are $(R, Z, \phi)$:

t = 0.0
x = [6.5, 0.5, 0.0]

B(field, t, x)

5.060970858985749

Both a vector and a splatted call are defined, so B(field, t, x) and B(field, t, x...) are the same thing:

same = B(field, t, x...) == B(field, t, x)
same

true

A perturbation is combined with the equilibrium before any code is generated, so a perturbed field is a FieldFunctions like any other:

perturbed = FieldFunctions(ThetaPinchEquilibrium(), EzCosZPerturbation())
φ(perturbed, t, [0.5, 0.5, 0.25])

0.15915494309189535

What Gets Generated

Rather more than just the magnetic field, and each quantity comes as a whole tensor rather than one function per component. A vector-valued quantity returns an SVector{3}, a Jacobian or the metric an SMatrix{3,3}, and the higher derivatives an SArray of the matching rank, so a single component is an ordinary index:

B♭(field, t, x)[1]

-0.20384615384615384

The naming follows the musical isomorphisms of differential geometry: ♭ lowers an index and gives the covariant components, ♯ raises one and gives the contravariant components, and ♮ gives the components in the physical (orthonormal) frame. They are typed \flat, \sharp and \natural. The bare letter is the magnitude, so B is $|B|$. See Coordinates for what the three representations are and Fields for how the quantities below are derived from the vector potential.

Four things are data rather than functions. The parameters of the equilibrium, which the generated code has baked in as literals:

parameters(field)

(R₀ = 6.2, B₀ = 5.3, q₀ = 2.0)

the equilibrium's own coordinate helpers, whose names vary from field to field:

keys(coordinates(field))

(:X, :Y, :Z, :R, :r, :θ, :ϕ, :r²)

which coordinates of the chart are periodic, periodic(field), and the handedness of the chart. J is the volume element and DF the Jacobian matrix; the two are related through orientation, and this equilibrium uses a left-handed $(R, Z, \phi)$ chart, so the determinant of DF is -J:

using LinearAlgebra

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

(-1, true)

functions(field) returns all of the generated functions at once, as a NamedTuple keyed by the names in the tables above.

Evaluating on a Grid

Since the generated functions are ordinary compiled Julia functions, sampling a field is just a comprehension:

R₀ = parameters(field).R₀

Rgrid = LinRange(0.5 * R₀, 1.5 * R₀, 100)
Zgrid = LinRange(-0.5 * R₀, 0.5 * R₀, 120)

Bfield = [B(field, t, Rgrid[i], Zgrid[j], 0.0)
          for i in eachindex(Rgrid), j in eachindex(Zgrid)]

size(Bfield)

(100, 120)

Note the index order: Bfield[i,j] holds the value at (Rgrid[i], Zgrid[j]), which is also the convention Makie's contour expects, so such an array can be passed straight to a plotting call. See Plotting for how to turn this into a figure.