Harmonic fields in a hollow torus
This example is one of the tests in MRX. It computes harmonic fields in a hollow toroidal shape and checks that the computed values for the poloidal and toroidal fluxes agree with Ampère’s law and approximate analytical expressions in the thin-torus limit.
Parameters
We begin by defining a number of parameters as well as the mapping function \(\Phi\) that maps the logical domain \([0, 1]^3\) to a hollow toroidal shape in \(\mathbb{R}^3\).
Ip = 1.95
It = 2.46
Is = jnp.array([Ip, It])
n = 6
p = 3
q = p + 2
ɛ = 0.1
π = jnp.pi
μ0 = 1.0
@jax.jit
def F(x):
"""Hollow toroid."""
r, θ, ζ = x
R = 1 + ɛ * (r + 1)/2 * jnp.cos(2 * π * θ)
return jnp.array([R * jnp.cos(2 * π * ζ),
-R * jnp.sin(2 * π * ζ),
ɛ * (r + 1)/2 * jnp.sin(2 * π * θ)])
de Rham sequence
Next, we set up the finite element spaces. Using the convenience method Seq.assemble_all() is only marginally more expensive than assembling only the needed operators.
ns = (n, n, n)
ps = (p, p, p)
types = ("clamped", "periodic", "periodic")
Seq = DeRhamSequence(ns, ps, q, types, F, polar=False, dirichlet=True)
Seq.evaluate_1d()
Seq.assemble_all()
Harmonic fields
There are two non-trivial harmonic fields in a hollow torus: One is linked to the current \(I_p\) flowing through the cavity/tunnel in the torus toroidally, the other is linked to the current \(I_t\) flowing poloidally around the torus and through the “donut” hole in the middle. This is our first sanity check (generalized hermitian eigenvalue problems are not yet implemented in jax.numpy.linalg, so we use scipy.linalg.eigh here):
evs, evecs = sp.linalg.eigh(Seq.M2 @ Seq.dd2, Seq.M2)
assert jnp.sum(evs < 1e-11) == 2 # two harmonic fields
assert jnp.min(evs) > -1e-11 # no negative eigenvalues
Note that Seq.M2 @ Seq.dd2 is the bilinear form corresponding to the weak formulation of the Laplace-de Rham operator on 2-forms. Seq.dd2 is the operator that corresponds to it and acts on the DoFs of 2-forms. In other words, when \(u_h = \sum_i \mathtt{u}_i \Lambda_i^2 \in V_2\):
The harmonic fields are the corresponding Seq.M2-orthonormal eigenvectors:
h1_dof = evecs[:, 0]
h2_dof = evecs[:, 1]
h1 = jax.jit(DiscreteFunction(h1_dof, Seq.Λ2, Seq.E2))
h2 = jax.jit(DiscreteFunction(h2_dof, Seq.Λ2, Seq.E2))
Ampère’s law
Next, we compute contour integrals of the harmonic fields along two non-contractible loops in the torus: one poloidal and one toroidal. These integrals give us the fluxes associated with each harmonic field. According to Ampère’s law, we should have
We compute the line integrals in the logical domain: \(\hat c_1\) and \(\hat c_2\) have simple expressions:
# contour around the enclosed tunnel poloidally:
def c1(θ): return jnp.array([1e-6, θ, 0])
# contour around the center tunnel toroidally:
def c2(ζ): return jnp.array([1 - 1e-6, 0.5, ζ])
The small offsets 1e-6 ensure that the contours lie in the interior of the physical domain. If \(h_1\) and \(h_2\) were one-forms, we could directly evaluate them along the curves in the logical domain since \(\int E \cdot \mathrm dl = \int \hat E \cdot \mathrm d \hat l\) for one-forms under pushforward \(E = \Phi_*^1 \hat E\).
However, they are two-forms, so instead it holds that
In code, the operation \((h, c) \mapsto \oint_c h \cdot \mathrm dl\) is given by:
def h_dl(twoform, curve):
def oneform(x):
DF = jax.jacfwd(F)(x)
return DF.T @ DF @ twoform(x) / jnp.linalg.det(DF)
def integrand(χ):
x = curve(χ)
dx = jax.jacfwd(curve)(χ).reshape(-1)
v = oneform(x).reshape(-1)
return jnp.dot(v, dx)
return integrand
# Integrate h1 along contours using trapezoidal rule
n_q = 256
_χ = jnp.linspace(0, 1, n_q, endpoint=False)
_w = jnp.ones(n_q) * (1/n_q)
Next, we build the matrix \(\mathbb P\) with entires \(\mathbb P_{ij} = \oint_{c_j} h_i \cdot \mathrm dl_j\). We can now solve for the magnetic field using Ampère’s law. It holds that
The ordering of \(h_1\) and \(h_2\) is of course arbitrary, it is this calculation that connects them to \(I_p\) and \(I_t\). As expected, the matrix \(\mathbb P\) only has two entries that are non-zero and these two differ by a multiplicative factor \(\varepsilon\).
P = jnp.array([
[jax.vmap(h_dl(h, c))(_χ) @ _w for c in (c1, c2)]
for h in (h1, h2)
])
b = jnp.linalg.solve(P.T, μ0 * Is)
b_dofs = b[0] * h1_dof + b[1] * h2_dof
Final checks
Finally, we can run a few assertions: The field should be divergence-free and curl-free.
curl_b_dofs = Seq.weak_curl @ b_dofs
div_b_dofs = Seq.strong_div @ b_dofs
assert (curl_b_dofs @ Seq.M1 @ curl_b_dofs)**0.5 < 1e-10
assert (div_b_dofs @ Seq.M3 @ div_b_dofs)**0.5 < 1e-10
Furthermore, we know that \(h_t = \mu_0 I_t \mathbf e_\phi / (2 \pi R)\) - this is simply the vacuum field of a solid torid with the correct magnitude to match \(I_t\). For the poloidal field, \(h_p \approx \mu_0 I_p \mathbf e_\theta / (2 \pi d)\), where \(d\) is the distance to the centerline of the enclosed tunnel assuming \(\varepsilon \ll 1\).
At \(\hat x = (r, θ, ζ) = (0.5, 0, 0)\), we have \(B_r = B_x\), \(B_θ = B_z\), and \(B_ζ = -B_y\) in Cartesian coordinates. Thus, we can check the values of the computed magnetic field \(B\) at this point: The error in \(B_e\) should be essentially zero since neither \(h_1\) nor \(h_2\) have a radial component. The error in \(B_ζ\) is determined by the resolution, while that in \(B_θ\) is dominated by the thin-torus approximation.
def B_expected(x):
r, θ, ζ = x
d = ɛ * (r + 1) / 2
R = 1 + d * jnp.cos(2 * π * θ)
sζ, cζ = jnp.sin(2 * π * ζ), jnp.cos(2 * π * ζ)
sθ, cθ = jnp.sin(2 * π * θ), jnp.cos(2 * π * θ)
B_ζ = μ0 * It / (2 * π * R)
B_θ = μ0 * Ip / (2 * π * d)
Bx = -B_ζ * sζ - B_θ * sθ * cζ
By = -B_ζ * cζ + B_θ * sθ * sζ
Bz = B_θ * cθ
return jnp.array([Bx, By, Bz])
B_computed = jax.jit(Pushforward(DiscreteFunction(
b_dofs, Seq.Λ2, Seq.E2), Seq.F, 2))
y = jnp.array([0.5, 0.0, 0.0])
B_diff = B_computed(y) - B_expected(y)
assert jnp.abs(B_diff[0]) < 1e-12
assert jnp.abs(B_diff/B_expected(y))[1] < (1/n)**p
assert jnp.abs(B_diff/B_expected(y))[2] < ɛ