ElectromagneticFields.jl – Modules

code(equ, pert = ZeroPerturbation(); export_parameters = true, escape = false, output = 0, cse = true)

Generate code for evaluating analytic equilibria: an Expr defining the field functions of equ, which @code splices into the calling module and load_equilibrium evaluates into one.

• export_parameters also emits the equilibrium's scalar parameters as constants.
• escape escapes the generated names, which a macro splicing this into another module needs.
• output is 0, 1 or 2: silent, one line per function, or the generated body as well.
• cse names each repeated subexpression once instead of emitting it in full every time it occurs. This is on by default and is value-preserving to the last bit — subexpressions are named, never rewritten — but the generated code is easier to read against a paper with it off. It matters: SymEngine shares nothing, so db₁dx₁ of the ITER Solov'ev equilibrium with an X-point contains 108 separate evaluations of log(ξ₁) without it, and is 6.8 times slower.

The @code macros take these as key = value arguments, after the equilibrium's parameters: ThetaPinch.@code(B₀, cse = false), or just ThetaPinch.@code cse = false for the defaults.

Every generated function takes (t, ξ₁, ξ₂, ξ₃) and (t, ξ), except orientation(), which takes no arguments: it is the sign of the chart's handedness, constant in both t and ξ. See orientation for what it means and why J alone is not enough.

code_arguments(args)

Split the arguments a @code macro was called with into the equilibrium's parameters and the options destined for code, returned as (parameters, options).

Macros cannot take keyword arguments, so an option is written as key = value among the arguments — @code cse = false or @code(R₀, B₀, q₀, cse = false). A ; key = value tail is accepted too.

Returns the l-th component of the Levi-Civita connection

Returns the i-th component of the one-form corresponding to the vector v

Returns the i-th component of the physical coordinate representation of the one-form α

Returns the i-th component of the vector corresponding to the one-form α

Returns the i-th component of the physical coordinate representation of the one-form α

Returns the m-th component of the cross-product between the vectors v and w

eliminate_common_subexpressions(body::Expr; prefix = "_cse")

Rewrite a generated function body so that every repeated subexpression is evaluated once and read from a local afterwards. Returns a begin ... end block whose value is that of body, and whose locals are named prefix * n.

Value-preserving to the last bit: subexpressions are named, not rewritten.

Generate functions for evaluating analytic equilibria.

Returns the m-th component of the one-form corresponding to the two-form β

load_equilibrium(equ, pert = ZeroPerturbation(); target_module = Main, output = 0, cse = true)
load_equilibrium(f, equ, pert = ZeroPerturbation(); target_module = Main, output = 0, cse = true)

Evaluate functions for evaluating analytic equilibria: generate the field functions of equ with code and evaluate them into target_module, which is returned.

output and cse are passed on to code — in particular cse = false emits every repeated subexpression in full, which is slower but easier to read against a paper.

The first form is only usable from top level. The methods it evaluates are defined in a world age newer than the frame that called it, so that frame cannot call them: from inside a function they are "too new", which Julia reports as a MethodError on a name that plainly exists. At top level — the REPL, a notebook cell, the body of a script — the world age advances between statements, so the generated functions are there for everything that follows.

The second form is for everywhere else. It passes target_module to f through Base.invokelatest, so the body runs in the world age the generated methods live in:

load_equilibrium(equ; target_module = MyModule) do mod
    mod.orientation()
    mod.DF(t, ξ)
end

Both the name lookup and the call have to happen inside f — on Julia 1.12 and newer a binding created by the evaluation is as invisible to the calling frame as a method is. f's value is returned. For a caller that cannot be restructured around a callback, Base.invokelatest(mod.B, t, ξ) does the same thing one call at a time.

The @code macros are not affected either way: they splice the same definitions in at expansion time, so they are already there when the surrounding code is compiled.

Returns the length of the vector v

Normalises the vector v in the metric g

Normalises the vector v in the metric g

orientation(equ)

Sign of the chart's orientation: +1 if (ξ¹, ξ², ξ³) is right-handed, -1 if left-handed.

J(x, equ) is the volume element √|g| = |det DF|, which is what most of the machinery here wants and what functions["J"] exports. The Hodge star and the cross product, however, are orientation-dependent and need the signed determinant det DF = orientation(equ) * J(x, equ).

Several of the charts in this package are left-handed — (R, Z, ϕ) and (r, θ, ϕ) both have det DF < 0, since the right-handed orderings would be (R, ϕ, Z) and (r, ϕ, θ) — so this is not a corner case. Getting it wrong reverses B without any other visible symptom: the same vector potential yields a magnetic field antiparallel to the one the cartesian chart gives at the same physical point.

Orientation of the charts in this package

Every equilibrium built on CartesianEquilibrium uses the identity map (x, y, z), so DF = I, J = 1 and the chart is right-handed. The four families below define their own coordinates, and all four are left-handed — in each case because the toroidal angle ϕ sits in the third slot where the right-handed ordering would put the second poloidal coordinate.

Concretely: AxisymmetricTokamakCartesian, SolovevSymmetric, ThetaPinch, Dipole, ABC, Singular, SymmetricQuadratic, QuadraticPotentials and the three Penning traps are right-handed; AxisymmetricTokamakCylindrical, AxisymmetricTokamakToroidal, AxisymmetricTokamakToroidalRegularization and every Solovev* equilibrium other than SolovevSymmetric (which is a cartesian chart despite the name) are left-handed.

code emits a zero-argument orientation() into the generated module — the one generated function that takes no arguments, since the sign depends on neither t nor ξ — so a loaded equilibrium can recover det DF = orientation() * J(t, ξ) without the equilibrium object.

test_analytic.jl asserts det(DF) ≈ orientation() * J for each of the charts above, reading the generated function, and separately that it equals the trait it was generated from. A new chart that declares the wrong sign therefore fails immediately rather than silently flipping its own B, and so does a generator that stops tracking the trait.

Returns the Christoffel symbol Γⱼₖˡ