LAGRANGIAN RELAXATION SCHEMES FOR CALCULATING FORCE-FREE MAGNETIC FIELDS, AND THEIR LIMITATIONS

In a numerical relaxation experiment using braided initial fields (Wilmot-Smith et al. 2009) we came across an inconsistency of the resulting numerical force-free state, which is best explained with the help of the following example. Consider a magnetic field obtained from the homogenous field by a simple twisting deformation as shown in Figure 1(a). Obviously, an ideal relaxation toward a force-free state must end again in a homogenous state (J = 0). During this process the Lagrangian relaxation leads to a deformation of the initial computational mesh which exactly cancels the initial deformation applied to the homogenous field. This is a well defined setup in which we know exactly the initial and final states. We now employ the implicit (alternating direction implicit—ADI) relaxation scheme detailed by Craig & Sneyd (1986) to relax our twisted field to a force-free equilibrium. The magnetic field is line tied on all boundaries (${\bf B}\cdot {\hat{\bf n}}=0$ on x and y boundaries). The J × B force as calculated by the numerical scheme decreases monotonically to an arbitrarily small value (e.g., 10⁻⁶–10⁻⁸), giving the appearance that the scheme converges (in an iterative sense) to a force-free equilibrium (to any desired accuracy, down to machine precision). However, when plotting α, the force-free proportionality factor, along a field line it shows variations which are by orders of magnitude higher than would be expected from |J × B| < 10⁻⁸. It is this inconsistency that we investigate in what follows. As we will discuss the convergence of the numerical scheme in what follows, it is worth emphasizing here the distinction between iterative convergence (at fixed resolution N) and real convergence, i.e., convergence toward a "correct" solution as the resolution N becomes sufficiently large.

Zoom In Zoom Out Reset image size

In order to investigate the source of the inconsistency described in the previous section, we consider the test problem outlined there. Specifically, we begin with an initially uniform magnetic field $({\bf B}=b_0{\bf {\hat{z}}})$, and superimpose two regions of toroidal field, centered on the z-axis at ±z₀, with exactly the same functional form, but of opposite signs:

$$\begin{equation} \mathbf {B} = b_{0} \hat{\bf {z}}\,\, +\, \sum _{i=1}^{2} \frac{2 b_{0} \phi _i}{\pi a_r} \exp\!\left(\!- \frac{x^2+y^2}{a_r^{2}} - \frac{\left(z-L_i\right)^{2}}{a_z^{2}}\! \right)\!\left(- y {\bf {\hat{x}}}+ x {\bf {\hat{y}}}\right), \end{equation} \tag{ 3 }$$

with b₀ = 1, $a_r=\sqrt{2}$, a_z = 2 and ϕ₁ = π, ϕ₂ = −π, L₁ = −4, L₂ = 4. We refer to this field in the following as T2. The field T2 (see Figure 1(a)) is constructed such that the two regions of twisted field, which are of opposite sign, should exactly cancel one another under an ideal relaxation, approaching the uniform field (with J = 0) as the equilibrium. Note that |ϕ| is the maximum turning angle of field lines around the z-axis.

One of the great advantages of an ideal Lagrangian relaxation is that it is possible to extract the paths of the magnetic field lines of the final state if one knows them in the initial state, simply by interpolating over the mesh displacement. Calculating the field lines in this way, no error is accumulated as occurs when integrating numerically along B. Given knowledge of the field line paths, one can test the quality of the force-free approximation by plotting α along field lines. For a force-free field

$$\begin{equation} \nabla \times {\bf B} = \alpha {\bf B}, \end{equation} \tag{ 4 }$$

and α should be constant along field lines since taking the divergence of the above yields

$$\begin{equation} {\bf B}\cdot \nabla \alpha =0. \end{equation} \tag{ 5 }$$

We begin by defining the variable α*, motivated by Equation (4), as

$$\begin{equation} \alpha ^* = \frac{ J_{\Vert } }{|{\bf B}|}. \end{equation} \tag{ 6 }$$

We find that for the magnetic field calculated by the relaxation scheme, the value of α* changes dramatically along field lines. Of course the relaxation gives a magnetic field for which J × B is not identically zero. So for a given value of J × B, what is the maximum possible variation in α* along a field line?

Consider

$$\begin{eqnarray} &&\nabla \cdot {\bf J} = \nabla \cdot (J_{\Vert } {\hat{\bf B}}) + \nabla \cdot {\bf J}_{\perp } = \delta,\nonumber \end{eqnarray}$$

say, where δ is representative of the error in calculating J. Using Equation (6) to replace J_|| gives

$$\begin{equation} {\bf B} \cdot \nabla \alpha ^* = - \nabla \cdot {\bf J}_{\perp } + \delta. \end{equation} \tag{ 7 }$$

or

$$\begin{equation} \frac{d\alpha ^*}{dl} = - \frac{\nabla \cdot {\bf J}_{\perp }}{|{\bf B}|} + \frac{\delta }{|{\bf {B}}|}. \end{equation} \tag{ 8 }$$

where l is a parameter along a magnetic field line with units of length. Now suppose that |J × B|/|B|² < within our domain. This implies that |J_⊥| <  |B|, so that

$$\begin{eqnarray} &&|\nabla \cdot {\bf J}_{\perp }| < \frac{\epsilon \,|{\bf B}|}{ d},\nonumber \end{eqnarray}$$

where d is the length scale of variations perpendicular to the magnetic field. Then from Equation (8)

$$\begin{equation} \left| \frac{d \alpha ^*}{dl}\right| < \frac{\epsilon }{d} + \frac{|\delta |}{|{\bf {B}}|}. \end{equation} \tag{ 9 }$$

Returning to our relaxation results, we have, for example, = 10⁻⁶, with |B| ≈ 1, $d \approx \sqrt{2}$. However, we find that |dα*/dl|_max ≈ 0.02. The discrepancy between this figure and the value of must come from the final term in Equation (9). This has been checked by interpolating the data onto a rectangular mesh and approximating ∇ · J using standard finite differences. We find ∇ · J ∼ O(10⁻²), and it therefore appears that the residual currents parallel to B are not relaxed because ∇ · J ≠ 0. As demonstrated below, this error does however decrease as the resolution is increased (see Tables 1–4).

It turns out that the appearance of the errors is related to the way in which J is calculated within the scheme, via a combination of first and second derivatives of the deformation matrix. These derivatives are calculated via finite differences in the numerical scheme, and it is here that these discretization errors arise. This is demonstrated below.

To ascertain the source of the errors, we take our initial state T2 and instead of performing the relaxation procedure, we artificially apply a deformation to the mesh which we can write down as a closed form expression, and moreover for which we can obtain the derivatives of the mesh displacement, and thus the resultant B and J fields, as closed form expressions. Motivated by the results of the relaxation, we impose a similar rotational distortion of the mesh which acts to "untwist" the field, via the transformation

$$\begin{equation} (x,\,y,\,z) \longrightarrow (x\,\cos \,\theta -y\,\sin \,\theta, y\,\cos \,\theta +x\,\sin \,\theta, z), \end{equation} \tag{ 10 }$$

where

$$\begin{equation} \theta = \psi \exp \left(-\frac{x^2+y^2}{4} -\frac{z^2}{16} \right), \end{equation} \tag{ 11 }$$

ψ constant. We now apply this transformation to an initially rectangular mesh on which B is given by T2, and compare the numerical and exact values for each entry in the mesh deformation Jacobian, and each component of B and J. Results are shown for three different values of the parameter ψ in Table 1.

It is clear that there are large errors in the current calculated by the numerical scheme. While errors in B and in each individual term in the mesh distortion Jacobian are much smaller, it turns out that the combination in which they are multiplied, summed, and divided to calculate J incurs large errors. The calculation has been meticulously checked such that we are certain that the error appears not due to mathematical or coding error, but rather due to an accumulation of numerical truncation errors in the process of calculating J (via Equation (2.10) in Craig & Sneyd 1986). Note that the errors increase as the mesh distortion (ψ) increases, and decrease with resolution (N).

Clearly, the accuracy of the force-free approximation is impaired by the accuracy of the curl operation performed in the numerical scheme. One way to increase the accuracy of spatial derivatives could be to use higher-order finite-difference expressions. Existing versions of the scheme use conventional second-order centered differences involving two nearest-neighbor values. If we expect smooth solutions (without grid-scale features) then increasing to fourth-order finite-difference expressions (using four nearest neighbors) is expected to increase the accuracy.

Recall that in the frictional Lagrangian method fluid displacements are determined from an equation of the form

$$\begin{eqnarray*} &&\frac{{\partial x_i}}{\partial t} = A_{i \alpha \beta \gamma } \, \, x_{\alpha, \beta \gamma } + C_i, \end{eqnarray*}$$

where A and C are prescribed tensor and vector functions and summation over repeated (Greek) indices is assumed. Note that the partial differentiation with respect to the background Cartesian coordinates (X₁, X₂, X₃) is indicated using the comma notation (i.e., ∂²f/∂β∂γ = f_,βγ).

The important point for us is that the method involves two spatial derivatives of the Lagrangian variables x_α. This reflects the fact that the Lorentz force is computed using first- and second-order derivatives of the Lagrangian mesh. Now "diagonal derivatives" such as x_i,jj are relatively easy to compute using finite differences: they involve the point itself and two/four nearest neighbors depending on whether the scheme is second or fourth order. It is these derivatives that are handled implicitly (via tridiagonal and pentadiagonal implementations of the ADI method) to provide, formally at least, the unconditional stability of the numerical scheme. Note that the computation of the mixed derivatives can be more complicated: terms such as x_i,jk involve 16 terms in the fourth-order scheme, as opposed to just four when the method is second order. However, irrespective of the formal accuracy, perhaps the main drawback of the scheme is that, unlike ∇ · B, the numerical evaluation of ∇ · J is not guaranteed to vanish identically. Thus, the accuracy of the relaxed solution can be compromised by the presence of rogue currents, especially in regions where the mesh is highly distorted and truncation errors are significant. The examples presented below show explicitly that this can restrict the convergence of the solution with resolution N.

Here, we present an algorithm for calculating the curl of a vector field (say J = ∇ × B) on a nonuniform mesh. This algorithm gives promising results, as discussed below, and is based on rewriting the curl operation, via Stokes' theorem, as a line integral:

$$\begin{equation} \oint _C {\bf B}\cdot {\bf dr} = \int _{U(C)} {\bf {J}}\cdot \hat{{\bf n}} dS = \int _{U(C)} {\bf J}_n dS, \end{equation} \tag{ 12 }$$

where the surface U(C), with unit normal ${\hat{{\bf n}}}$, has the closed curve C as its boundary. The idea is similar to that of Hyman & Shashkov (1997) who have applied such so-called mimetic numerical methods to solving Maxwell's equations (Hyman & Shashkov 1999). Our algorithm differs in some ways from theirs.

Equation (12) can be discretized as follows. Suppose that we want to calculate J at the mesh point X_i,j,k. There are three mesh surfaces that intersect at this point. The first is the ith mesh level in the first index direction. Consider the circuit in this surface shown in Figure 2, and let the nearest-neighbor points to X_i,j,k be denoted x_I, x_II, x_III, x_IV. Then, defining dx₁ = x_II − x_I, dx₂ = x_III − x_II, etc. and B₁ = (B(x_I) + B(x_II))/2, B₂ = (B(x_II) + B(x_III))/2, etc., we can approximate the left-hand side of Equation (12) by

$$\begin{equation} \mathcal {I} = {\bf {B}}_1\cdot {\bf {d}}{\bf {x}}_1 + {\bf {B}}_2\cdot {\bf {d}}{\bf {x}}_2 + {\bf {B}}_3\cdot {\bf {d}}{\bf {x}}_3 + {\bf {B}}_4\cdot {\bf {d}}{\bf {x}}_4. \end{equation} \tag{ 13 }$$

Zoom In Zoom Out Reset image size

Furthermore, the area of the enclosed quadrilateral is

$$\begin{equation} A= \case{1}{4}\, |{\bf dx}_1 \times {\bf dx}_2| + \case{1}{4}\, |{\bf dx}_2 \times {\bf dx}_3| + \case{1}{4}\, |{\bf dx}_3 \times {\bf dx}_4| + \case{1}{4}\, |{\bf dx}_4 \times {\bf dx}_1|. \end{equation} \tag{ 14 }$$

We can now define the direction perpendicular to this mesh surface as

$$\begin{eqnarray} {\bf {n}}^{(1)} &=& \frac{1}{4}\bigg(\frac{{\bf dx}_1 \times {\bf dx}_2}{|{\bf dx}_1 \times {\bf dx}_2|} + \frac{{\bf dx}_2 \times {\bf dx}_3}{|{\bf dx}_2 \times {\bf dx}_3|} + \frac{{\bf dx}_3 \times {\bf dx}_4}{|{\bf dx}_3 \times {\bf dx}_4|} \nonumber\\ &&+ \frac{{\bf dx}_4 \times {\bf dx}_1}{|{\bf dx}_4 \times {\bf dx}_1|} \bigg), \end{eqnarray} \tag{ 15 }$$

and (comparing Equations (12)–(14)) the current in this direction as

$$\begin{equation} J_n^{(1)} = {\mathcal {I}}/{A}. \end{equation} \tag{ 16 }$$

In the same way, we can define J⁽²⁾_n and J⁽³⁾_n perpendicular to the other two mesh surfaces, with normal vectors n⁽²⁾ and n⁽³⁾, that pass through X_i,j,k.

Now, the J^(p)_n are projections of the current we require (denoted J_s) in such a way that

$$\begin{equation} {\bf {n}}^{(p)}\cdot {\bf {J}}_s = J_n^{(p)},\qquad p=1,\,2,\,3. \end{equation} \tag{ 17 }$$

Denoting by $\mathcal {N}$ the matrix whose rows are the row vectors n⁽¹⁾, n⁽²⁾, and n⁽³⁾, and by J_n the vector with components J⁽¹⁾_n,  J⁽²⁾_n,  J⁽³⁾_n we can re-write Equation (17) as $\mathcal {N} {\bf J}_s = {\bf J}_n$. Finally, we obtain the current at point X_i,j,k via

$$\begin{equation} {\bf {J}}_s = \mathcal {N}^{-1} {\bf {J}}_n. \end{equation} \tag{ 18 }$$

Note that $\mathcal {N}$ is always invertible assuming that the n^(p) are linearly independent, i.e., as long as the grid cells have nonzero volume. It is straightforward to verify that this procedure reduces to the standard second-order centered difference expression in each direction for a rectangular (undeformed) mesh.

Since this method makes use of Stokes' theorem, which is "topological" in the sense that is does not depend on a deformation of the loop we are integrating over, we expect the method to be more robust and accurate for deformed grids (see also Hyman & Shashkov 1997). The algorithm has not yet been implemented in the full numerical scheme, as it requires a complete rewriting of the implicit time stepping routine, and a simple explicit implementation turns out to be prohibitively slow to run at reasonable resolution. However, the algorithm is tested and used as a diagnostic in what follows.

First, to demonstrate that in principle the original (second-order) relaxation scheme is sound, we perform the relaxation "experiment" in the following way. We begin with a uniform field $({\bf B}=b_0 {\bf {\hat{z}}})$ on a uniform mesh, and then deform this mesh in such a way that the resulting magnetic field is T2 (with ϕ_i = ±π). This configuration is then relaxed, to a level where the Lorentz force calculated by the numerical scheme is reduced to J × B = 10⁻⁶. The equilibrium field corresponding to T2 is the uniform field ${\bf B}=b_0{\bf {\hat{z}}}$, and in this case (due to the frozen-in condition) the straight field should correspond to the rectangular mesh. (The choice |ϕ| = π is motivated by the braiding example discussed in Wilmot-Smith et al. 2009 since this is the minimum level of twist which yields a magnetic field whose field lines are truly braided in that case.)

To diagnose the success of the relaxation, referring back to Equation (9), we calculate

$$\begin{equation} \epsilon ^* = d \, \frac{\Delta \alpha ^*}{\Delta l}, \end{equation} \tag{ 19 }$$

taking $d=\sqrt{2}$ (radius of twist regions) and Δα* and Δl as the maximum change in α* over a given length, which occurs along the central field line—the z-axis—by symmetry. This expression puts a true value on the "quality" * of the force-free approximation: for a given variation in α*, * provides a lower bound for the maximum value of |J ×  B|/|B|² within our domain, for a current free of errors (i.e., setting δ = 0 in Equation (9)). We obtain a value of * = 9.5 × 10⁻⁷ for resolution N = 61, demonstrating that discretization errors in the scheme are very small when the relaxed state has an approximately rectangular mesh.

We now investigate the promise of the two extensions to the scheme described in the previous section. In our test case (T2) the topology of the field is simple, and the equilibrium field known, permitting the above approach (i.e., relaxation toward a uniform mesh). However, to investigate magnetic fields with non-trivial topology we must approach the problem in a different way. We therefore return to the case where we begin with T2 on a rectangular mesh, setting ϕ_i = ±π.

Performing the artificially imposed analytical ("untwisting") deformation instead of relaxation, we see that derivatives calculated with the two new methods (fourth-order finite differences and the Stokes-based method) both give smaller maximum errors than with the original second-order scheme (see Table 2). For a less deformed mesh (ψ = π/2), the fourth-order finite differences perform better than the Stokes-based routine. However, for a more distorted mesh (ψ = π), the Stokes routine gives significantly lower errors than either of the finite-difference methods. It is particularly interesting to note that the errors for the Stokes-based method seem to scale relatively weakly with the mesh deformation, suggesting it to be a good choice for highly deformed meshes. Both new schemes, for reasonable resolution and levels of deformation (N ⩾ 41, ψ = π) give errors that are an order of magnitude lower than the original (second-order) method. This suggests these methods are worth pursuing, so we go on to perform relaxation simulations using them.

We now leave the artificially imposed deformation, and relax T2 using both the second- and fourth-order schemes. The relaxation is allowed to run until |J × B| as calculated by the relevant numerical scheme is reduced to 10⁻⁵. We compare this value with * (defined by Equation (19)), and also the maximum value of the Lorentz force obtained by calculating J via the Stokes-based routine, denoted J_s.

The results for two levels of initial twist $(\phi _i=\pm \frac{\pi }{2}, \phi _i=\pm \pi)$ are displayed in Tables 3 and 4. A number of points are immediately clear. First, in no simulation do we approach the apparent value of = 10⁻⁵. However, the use of fourth- rather than second-order finite differences improves the quality of the relaxed field by an order of magnitude in * when |ϕ| = π/2. Increasing the deformation in the final state (by increasing |ϕ| in the initial state T2) has a strong adverse effect on the relaxation process. This is found to be because spurious (unphysical) current concentrations arise where none should reasonably be expected. Examining the corresponding mesh, we find that these "false" current regions appear where the grid is most distorted—see Figure 3, and compare with Figure 1(b). The "current shards" shown in Figure 3(b) actually intensify as the relaxation proceeds. We find that this is possible since ∇ · J (approximated by interpolating J onto a uniform mesh) is not sufficiently close to zero. As a result there is no "return current" associated with these localized current regions, which might be expected to generate a Lorentz force that would act against the further intensification of the current shards (if they have no physical basis). In previous studies using such codes (e.g., Craig & Litvinenko 2005; Pontin & Craig 2005), intensification of |J| as |J × B| decreased in time was associated with current singularities, so at first sight it appears that these current shards could naively be interpreted as "current sheets", which would of course be unphysical. Note, however, that the most important signature of current singularity in previous studies was a (power-law) proportionality of the peak current with mesh resolution (for given ). We have found that in fact the current shards become less intense as N is increased (as required by the convergence of the scheme), so there is a clear distinction between the two phenomena.

Zoom In Zoom Out Reset image size

Table 3. Values of the Force-Free Quality Parameter * and the Lorentz Force Calculated via the Stokes-Based Routine, J_s × B/|B|², for Simulation Runs with Second- and Fourth-Order Finite Differences (Two n.n./Four n.n.) and Resolution N³, to Two Significant Figures

	Two n.n.		Four n.n.
N	*	Js × B/|B|2	*	Js × B/|B|2
21	0.11	0.053	0.048	0.063
41	0.054	0.046	0.0067	0.018
61	0.032	0.035	0.0042	0.0050
81	0.021	0.027	0.0015	0.0019

Note. Twist parameter ϕ_i = ±π/2.

Download table as:  ASCIITypeset image

Table 4. As Table 3, with Twist Parameter ϕ_i = ±π

	Two n.n.		Four n.n.
N	*	Js × B/|B|2	*	Js × B/|B|2
21	0.28	0.18	0.17	0.21
41	0.16	0.17	0.071	0.13
61	0.11	0.13	0.026	0.062
81	0.074	0.11	0.021	0.023

Download table as:  ASCIITypeset image

Finally, consider the values of J_s × B/|B|² we find for the relaxed fields (Tables 3 and 4). They are clearly of the same order as * (note that * based on J from the numerical scheme and * based on J_s are of the same order), and thus it seems that an implementation involving the Stokes-based routine has the capacity to yield a magnetic field that is much closer to being force-free (with lower *). This is illustrated in Figure 4. We consider the observed value of α* based on Equation (6) (dashed line), compared with a maximum allowable value for α*. Since α* is antisymmetric about z = 0, the maximum value allowed, based on Equation (9), is

$$\begin{equation} \alpha ^*_{\max } = \frac{L \epsilon }{d}, \end{equation} \tag{ 20 }$$

and we take |B| = 1 and $d=\sqrt{2}$ as before, and L = 20, the length of the domain. Then we obtain the maximum possible α* by taking to be the maximum value during the relaxation of J × B (solid line) or J_s × B (dot-dashed). We see that very early in the relaxation the actual value of α* becomes greater than the maximum allowed by J × B from the numerical scheme (second-order). Moreover, the discrepancy grows steadily. However, α* always remains less than the maximum allowed by the Stokes-based method, implying that this may be a more sound method to calculate the current and resulting Lorentz force.

Zoom In Zoom Out Reset image size