© ESO, 2013

1. Introduction

The calculation of three-dimensional non-linear magnetohydrostatic (MHS) equilibria is a non-trivial subject. In the solar physics context, progress has been made by considering linear subclasses of MHS equilibria. One such class, known as laminated equilibria, makes use of an Euler potential representation of the magnetic field (e.g. Low 1982). On each lamina, a 2D magnetic field is calculated and the 3D field comprises of the union of these laminas.

Another approach is to model the current density as a linear combination of field-aligned and cross-field currents (Low 1991, 1992). Here, the field-aligned current is that of a linear force-free field and the cross-field current depends on the variation of the magnetic field with height. Details of this will be presented in the next section. Neukirch & Rastätter (1999; NR99 hereafter) present a new formulation of Low’s model by writing the magnetic field in terms of poloidal and toroidal components. This has the advantage that the calculation of the magnetic field only involves one scalar function whereas, previously, one was forced to operate with all three components of the magnetic field independently. Petrie & Neukirch (2000) use the representation of NR99 and find closed-form solutions for MHS equilibria via Green’s functions. The price paid for finding closed-form solutions, however, is that they are forced to choose a simple term for the cross-field current. To date, the authors are unaware of any closed-form solutions using Green’s functions, other than those presented in Petrie & Neukirch (2000).

Although more detailed analytical solutions may prove difficult with the representation of NR99, it is, however, set up perfectly for an efficient numerical treatment. In this paper, we outline a simple and fast numerical procedure for calculating MHS equilibria based on the NR99 representation. The details of this are given in the next section. This is followed by some examples of force-free and linear MHS equilibria. For the linear MHS case, we investigate the effects of the cross-field currents on the equilibria. We then highlight how the scheme can be parallelized simply and effectively on a graphics processing unit (GPU). The paper concludes with a summary.

2. Model equations and solution method

Firstly, we shall outline where the equations to be solved come from. Fuller details can be found in NR99. This will be followed by an algorithm for the numerical solution.

2.1. The model equations

The MHS equations that we shall solve are

$\begin{array}{ccc}& & {\mathit{\mu }}_{\mathrm{0}}^{-1}\mathrm{(}\mathrm{\nabla }\mathrm{\times }{B}\mathrm{)}\mathrm{\times }{B}\mathrm{-}\mathrm{\nabla }\mathit{p}\mathrm{-}\mathit{\rho g}{ê}\mathit{z}\mathrm{=}{0}\mathit{,}\\ & & \mathrm{\nabla }\mathrm{\times }{B}\mathrm{=}{\mathit{\mu }}_{\mathrm{0}}{j}\mathit{,}\\ & & \mathrm{\nabla }\mathrm{·}{B}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$

where B is the magnetic induction (commonly referred to as the magnetic field), j is the current density, p is the plasma presure, ρ is the density, g is the (constant) gravitational acceleration and μ₀ is the permeability of free space. To complete the problem, boundary conditions must be specified. This choice is problem-dependent but will, at least, require B to be specified on the boundaries of the domain through Dirichlet or von Neuman conditions. In this paper we will consider a Cartesian domain. Following Low (1992), we assume the current density takes the following form ${\mathit{\mu }}_{\mathrm{0}}{j}\mathrm{=}\mathit{\alpha }{B}\mathrm{+}\mathrm{\nabla }\mathrm{\left(}\mathit{gF}\mathrm{\right)}\mathrm{\times }{ê}\mathit{z}\mathit{,}$(4)where α is a constant and F is an arbitrary function. The first term on the right-hand-side is the aligned current associated with a linear force-free field. The second term is the non-force-free cross current due to the gravitational stratification. The magnetic induction can be written as ${B}\mathrm{=}\mathrm{\nabla }\mathrm{\times }\mathrm{[}\mathrm{\nabla }\mathrm{\times }\mathrm{(}\mathit{P}{ê}\mathit{z}\mathrm{)}\mathrm{+}\mathit{T}{ê}\mathit{z}\mathrm{]}\mathit{,}$(5)where P and T are scalar functions corresponding to the poloidal and toroidal components respectively. This form satisfies Eq. (3). Using this in Eq. (2) with Eq. (4) gives

$\begin{array}{ccc}& & \mathit{T}\mathrm{=}\mathit{\alpha P,}\\ & & {\mathrm{\nabla }}^{\mathrm{2}}\mathit{P}\mathrm{+}{\mathit{\alpha }}^{\mathrm{2}}\mathit{P}\mathrm{+}\mathit{gF}\mathrm{=}\mathrm{0.}\end{array}$

Following NR99 and Low (1992), we set

$\mathit{F}\mathrm{=}{\mathit{g}}^{-1}\mathit{\xi }\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{B}}_{\mathit{z}}\mathit{.}$With this identity, Eq. (7) becomes $\mathrm{[}\mathrm{1}\mathrm{-}\mathit{\xi }\mathrm{(}\mathit{z}\mathrm{\left)}\mathrm{\right]}\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right)\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{+}{\mathit{\alpha }}^{\mathrm{2}}\mathit{P}\mathrm{=}\mathrm{0.}$(8)For ξ(z) < 1, Eq. (8) is elliptic. Now the magnetic field can be found by solving this for the poloidal scalar potential P and then using Eqs. (5) and (6).

Once the magnetic field is found, the plasma pressure and density are given by $\mathit{p}\mathrm{=}{\mathit{p}}_{\mathrm{b}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{-}\mathit{\xi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\frac{{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}$(9)and $\mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{\mathit{g}}\left(\frac{\mathrm{d}\mathit{\xi }}{\mathrm{d}\mathit{z}}\frac{{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\mu }}_{\mathrm{0}}}\mathit{\xi }{B}\mathrm{·}\mathrm{\nabla }{\mathit{B}}_{\mathit{z}}\right)\mathit{.}$(10)p_b and ρ_b are the background equilibrium plasma pressure and density respectively. A derivation of these formulae is given in NR99.

2.2. The numerical procedure

An efficient numerical solution of this class of the MHS equations rests on reducing the equations to a set of linear, scalar elliptic problems. I.e. those of the form

$\mathit{a}\mathrm{\left(}{x}\mathrm{\right)}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{u}\mathrm{+}\mathit{b}\mathrm{\left(}{x}\mathrm{\right)}\mathit{u}\mathrm{=}\mathit{c}\mathrm{\left(}{x}\mathrm{\right)}\mathit{,} {x}\mathrm{\in }\mathrm{𝒟}\mathit{,}$with u defined on . Here the Laplacian operator can be two-dimensional (for boundaries) or three-dimensional (for the interior) and a, b and c are known scalar functions.

There exists a wide variety of techniques to solve such problems efficiently, from Krylov subspace methods to multigrid methods. Although solution methods are problem-dependent, the multigrid method is often an optimal solver for discrete Poisson problems (Elman et al. 2005), i.e. its convergence rate is independent of the problem mesh size. It can also be parallelized efficiently, and for these reasons we adopt it in this work. The multigrid method, in a different representation to the one of this paper, was used successfully for MHS models of flux tubes bounded by current sheets (Henning & Cally 2001).

As mentioned above, the MHS equations are reduced to a set of linear elliptic problems defined on the boundaries and in the interior. Equation (8) is the three-dimensional elliptic problem for the poloidal scalar potential P defined in the interior of the domain. The other elliptic equations required to be solved depend on the boundary conditions. On the lower boundary, the vertical component of the magnetic field, B_z, is prescribed. From Eq. (5), it is clear that, on the lower boundary, $\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right)\mathrm{=}\mathrm{-}{\mathit{B}}_{\mathit{z}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathit{.}$(11)In principle other elliptic equtions can be defined on the other boundaries of the Cartesian domain. For simplicity, however, we shall assume the B = 0 on these boundaries. For the scalar elliptic problem, this translates to P = 0 on top and side boundaries. This boundary condition also applies to the sides of the 2D domain (the lower boundary) where Eq. (11) is solved. In short, the numerical procedure to determine the poloidal scalar potential P involves:

• 1.

  Set P = 0 everywhere, initially, on the grid.

• 2.

  Use the multigrid method to solve Eq. (11) on the lower boundary with P = 0 on the sides.

• 3.

  Use the multigrid method to solve Eq. (8) within the domain, using P = 0 on the top and side boundaries and the distribution of P on the lower boundary determined from the previous step.

After these steps are completed, the magnetic field B can be calculated. From Eq. (5), the other components of the magnetic field are given by

${\mathit{B}}_{\mathit{x}}\mathrm{=}\mathit{\alpha }\frac{\mathit{\partial P}}{\mathit{\partial y}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial x\partial z}}$and

${\mathit{B}}_{\mathit{y}}\mathrm{=}\mathrm{-}\mathit{\alpha }\frac{\mathit{\partial P}}{\mathit{\partial x}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial y\partial z}}\mathrm{·}$These, together with Eq. (11), can be approximated using finite differences, giving B_x, B_y and B_z on the grid to a required accuracy. In this paper, we use second-order accurate central finite differences.

All that remains now is the calculation of the pressure p and the density ρ. Once B_z is determined on the grid, it is clear from Eq. (9) that p can be evaluated directly. In Eq. (10) the derivatives must be dealt with using finite differences. After this, ρ is determined on the grid.

3. Examples

Here we present examples of equilibria to demonstrate the effectiveness of our numerical scheme. Details of computational aspects and parallelization are discussed in Sect. 4.

For the cases under consideration we non-dimensionalize Eqs. (1) to (3) with respect to photospheric values. The variables are plasma pressure, p₀ = 1.4 × 10⁴ Pa; density, ρ₀ = 3 × 10^-4 kg/m³; scale height H₀ = 340 km; magnetic induction B₀ = (2μ₀p₀)^1/2 = 1.3 × 10³ G and temperature, T₀ = p₀/(Rρ₀) = 5.6 × 10³ K. Here, R is the gas constant.

We assume a background temperature profile of the form

${\mathit{T}}_{\mathrm{b}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathrm{1}\mathit{,}& \mathrm{0}\mathrm{\le }\mathit{z}\mathrm{\le }\mathrm{5}\mathit{,}\\ {\mathit{T}}_{\mathrm{cor}}^{\mathrm{(}\mathit{z}\mathrm{-}\mathrm{5}\mathrm{)}\mathit{/}\mathrm{5}}\mathit{,}& \mathrm{5}\mathit{<}\mathit{z}\mathrm{\le }\mathrm{10}\mathit{,}\\ {\mathit{T}}_{\mathrm{cor}}\mathit{,}& \mathit{z}\mathit{>}\mathrm{10.}\end{array}$T_cor = 150 is the non-dimensional coronal temperature. The model photosphere/chromosphere ranges from 0 ≤ z ≤ 5, the transition region ranges from 5 < z ≤ 10 and the corona is in z > 10. With an ideal gas equation

${\mathit{p}}_{\mathrm{b}}\mathrm{=}{\mathit{\rho }}_{\mathrm{b}}{\mathit{T}}_{\mathrm{b}}\mathit{,}$hydrostatic pressure balance can be written as

$\frac{\mathrm{d}{\mathit{p}}_{\mathrm{b}}}{\mathrm{d}\mathit{z}}\mathrm{=}\mathrm{-}\frac{{\mathit{p}}_{\mathrm{b}}\mathit{g}}{{\mathit{T}}_{\mathrm{b}}}\mathrm{·}$The equation is solved for p_b and then ρ_b follows from the ideal gas law. The non-dimensionalization and atmospheric model presented here are similar to those in flux emergence studies (e.g. MacTaggart 2011; Hood et al. 2012; McLaughlin et al. 2012). For the following cases, we set B = 0 on the side and top boundaries. The size of the domain is [−5.5,5.5] ² ×  [0,12] and we use a resolution of 128³.

3.1. Force-free case

For a force-free field, ξ = 0. This, of course, means that p = p_b and ρ = ρ_b. In this example we take α = 0.4. We model a magnetic configuration with three sources, similar to that in Régnier et al. (2005). For each source, we define a Gaussian profile for B_z, on the lower boundary, of the form

${\mathit{B}}_{\mathit{z}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{0}}\exp \mathrm{(}\mathrm{-}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{l}}^{\mathrm{2}}\mathrm{)}\mathit{.}$B₀ is the field strength at the centre of the source, l is source width and r² = (x − x₀)² + (y − y₀)² with source centre (x₀,y₀). Here we take l = 0.3 for all the sources. There is one positive source with B_z = 1 and two negative sources with B_z = −0.5. The position of the positive source is (1.5, 1.5). The two negative sources have positions (−1.5, −1.5) and (1.5, −1.5). Figure 1 shows some field lines of the calculated region, with a magnetogram of B_z displayed at the base.

Fig. 1 Magnetic field lines of the linear force-free region are displayed in orange. A magnetogram of Bz is given at the base to indicate the positions of the sources.

Between the two negative polarities, there is a rapid divergence in the field lines as they connect down to the positive polarity. This indicates the presence of a null point at the base of the domain. Further evidence for this can found by looking at the distribution of | B | near the base of the domain. Figure 2 displays this at z = 0.1.

Fig. 2 A map of | B |, within the range given on the scale, at z = 0.1. The null point is clearly highlighted between the two negative polarities.

3.2. Non-force-free case

Having demonstrated that our multigrid scheme can successfully compute linear force-free equilibria with rapidly changing field lines and null points, we now turn our attention to the linear MHS case. For this we examine two profiles for ξ: $\begin{array}{ccc}{\mathit{\xi }}_{\mathrm{1}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& \mathrm{0.7}\exp \mathrm{(}\mathrm{-}\mathrm{0.2}\mathit{z}\mathrm{)}\mathit{,}\\ {\mathit{\xi }}_{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& \mathrm{[}\mathrm{0.7}\mathrm{+}\mathrm{0.3}\sin \mathrm{(}\mathit{\pi z}\mathrm{\left)}\mathrm{\right]}\exp \mathrm{(}\mathrm{-}\mathrm{0.1}\mathit{z}\mathrm{)}\mathit{.}\end{array}$These profiles are displayed in Fig. 3. ξ₁(z) is a simple exponential decay which models the fact the magnetic field becomes more force-free as one moves from the photosphere to the corona. ξ₂(z) is also exponentially decaying but has the additional complexity of a sine wave superimposed on it. Note that both of these profiles satisfy the condition ξ(z) < 1 to ensure that Eq. (8) is elliptic.

Fig. 3 The profiles of ξ(z) from the cross-field current.

3.2.1. Magnetic field and current

The sigmoidal magnetic field for profile ξ₁ is displayed in Fig. 4. That for ξ₂ is in Fig. 5. In both figures, a magnetogram of B_z is placed at z = 0 and the field lines are shaded with | ∇ × B |.

Fig. 4 Magnetic field lines for the MHS case with ξ1. The field lines are coloured with | ∇ × B |.

Fig. 5 Magnetic field lines for the MHS case with ξ2. The field lines are coloured with | ∇ × B |.

The magnetic field of ξ₂ is the more twisted of the two. The current is stronger at the footpoints and this is due to the initial increase in the ξ₂ profile. Although the magnetic fields for both profiles look similar, superficially at least, the structure of the current density is different. Figure 6 displays an isosurface of | j | for the ξ₁ profile field. Here, the current density is concentrated at the footpoints. Figure 7 displays the corresponding isosurface for the ξ₂-profile field. The oscillating ξ₂-profile allows for additional structure within the current density, as evidenced by an additional bridging arch in Fig. 7. As one moves higher into the corona, the current density becomes weaker and the magnetic field becomes close to potential.

Fig. 6 An isosurface of | j| = 0.03 for the MHS solution with ξ1.

Fig. 7 An isosurface of | j| = 0.03 for the MHS solution with ξ2.

3.2.2. Pressure and density

As is evident from Eqs. (9) and (10), the plasma pressure and density are highly dependent on the magnetic field and the ξ-profile. Although the ξ-profiles are decreasing with height (on average for ξ₂), they are always positive. From Eq. (9) it is clear that the effect of a positive ξ-profile is to introduce a pressure deficit where there is strong vertical magnetic field. Since the temperature is constant in the region near the photosphere, the geometrical depression, as measured from the pressure p, corresponds approximately to the Wilson depression (Spruit 1976). For the linear MHS equilibria of this paper, a source of size ≈300 km produces a Wilson depression of ≈100 km. These values lie within the same order of magnitude as those given in Spruit (1976).

Fig. 8 The non-dimensional density along a diagonal cut through the domain at z = 0.46875. The two troughs correspond to the locations of the magnetic footpoints. This profile is for the linear MHS equilibrium with ξ2.

Equation (10) demonstrates that the density enhancement or depletion is due to the competition of a magnetic pressure term combined with ξ′(z) and a magnetic tension term. Considering the ξ₁-profile first, ${\mathit{\xi }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{<}\mathrm{0}$ for all z. Hence, the “pressure” term is negative. The “tension” term is also negative due to the fast drop in field strength of B_z with height. For the ξ₂-profile the situation is slightly more complex as ${\mathit{\xi }}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{z}\mathrm{\right)}$ is zero, positive and negative at different heights. This means that the “pressure” term changes sign with height. For the numerical values chosen in this paper, however, the “tension” term dominates in magnitude. This results in two density deficits over the locations of the footpoints, as shown in Fig. 8. If one chooses a “pressure” term that is greater in magnitude than the “tension” term, then one could produce a linear MHS model with density enhancements at different heights.

4. Computational aspects of the implementation

The scheme outlined in this paper reduces the problem of solving the full linear MHS equations to a finite set of scalar elliptic problems. We have implemented a multigrid method as an efficient solver of such problems. To reduce the time taken for calculations, we have made use of the massive parallelization capabilities of graphical processing units (GPUs). The development of GPUs has been driven by the computer games industry. In the last few years, however, general purpose GPU computing (GPGPU) has grown due to the development of APIs such as CUDA (for NVIDIA GPUs) and OpenCL. The use of these languages is now an important part of scientific computing. Within solar physics, GPUs have not been exploited greatly yet. This is likely to change, however, as they give scientists the capability of parallelizing their codes with hundreds or thousands of threads without having to completely redesign them.

For the problem in hand, the multigrid method involves various stages that can be treated in parallel. The relaxation steps (Gauss-Jacobi), the calculation of the residual and the interpolations between grids (reduction and prolongation) can each be treated in a parallel fashion. To perform this, one has to map the computational grid onto the GPU’s grid of threads that can be run in parallel. We demonstrate this below with an example of how to perform Gauss-Jacobi iteration on a GPU. For an introduction to GPUs and programming in CUDA, we recommend Sanders & Kandrot (2011).

Subroutines that are run on the GPU are known as kernels. Below is pseudocode showing the key elements of a Gauss-Jacobi kernel used in the multigrid evaluation of Eq. (11).

ix = threadIdx.x + blockIdx.x*blockDim.x;

iy = threadIdx.y + blockIdx.y*blockDim.y;

off = ix + iy * blockDim.x * gridDim.x;

io = 1;

jo = DIM;

invdx2 = 1.0/(dx*dx);

invdy2 = 1.0/(dy*dy);

if ( on boundary )

{

set d_P2[off] here;

}

else {

d_P2[off] = ((d_P1[off+io] + d_P1[off-io])*invdx2

+(d_P1[off+jo] + d_P1[off-jo])*invdy2

+Bz[off]) / (2.0*(invdx2+invdy2));

}

The first three lines of pseudocode link the node positions of the computational grid to that of the GPU, which consists of blocks that are divided into threads. Here we have one thread per node. In the x-direction, threadIdx.x is the thread index, blockIdx.x is the block index and blockDim.x is the number of threads per block. This combination gives the ix position. The iy position is found similarly. The third line of pseudocode calculates an offset that allows one to write all the positions on a 2D grid within a 1D array. Here, gridDim.x is the number of blocks in the x-direction. In this example, it is the same number as in the y-direction.

Now that a thread is linked to every node, we must be able to perform calculations involving nodes at different positions. To move to a neighbouring node in the x-direction, io=1 is added to or taken from off. Similarly, to move to a neighbouring node in the y-direction, jo=DIM is added to or taken from off. Here, DIM is the dimension of the grid in the x-direction.

With the grid set up, and the ability to move through it, the finite difference implementation of one iteration of Gauss-Jacobi relaxation is simple to implement. In the pseudocode, there is one command for points on the boundary and another for points in the interior. Values from the previous iteration are read from d_P1[]. The new values of the Gauss-Jacobi iteration are written to d_P2[]. Both variables begin with d to highlight that they are on the device (GPU) as opposed to the host (CPU).

This example demonstrates how one can parallelize a finite difference scheme simply. Edits to codes can be made function by function, rather than having to redesign from scratch. The purpose of this endeavour is, of course, to achieve a sizeable speed-up in the running of a code. Figure 9 compares the run times of the serial and parallel versions of the multigrid code for different resolutions. For these calculations, we use the same parameters as for the MHS case described in the paper but set ξ = 0. i.e. we solve the force-free case. As the resolution increases, the benefits of the parallel code run on the GPU become obvious. For resolutions of 128³ and 256³, impressive speed-ups of ×32.418 and ×29.632 are achieved respectively.

Fig. 9 The run times for the serial (blue) and parallel (green) versions of the multigrid code plotted against resolution. Each run time corresponds to the combined time for a 2D multigrid method calculation on the base of the domain and a 3D multigrid method calculation within the domain.

5. Summary

In this paper we calculate linear MHS equilibria with an efficient numerical scheme based on the representation of NR99. The outline algorithm of this scheme is as follows:

• Calculate the background equilibrium plasma plasma pressurep_b and density ρ_b.

• Prescribe B_z on the lower boundary and find the poloidal scalar potential P on this boundary by solving, using a multigrid method, the Poisson equation $\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\mathrm{=}\mathrm{-}{\mathit{B}}_{\mathit{z}}\mathit{,}$with P defined on the side boundaries.

• Define the other boundary conditions for P, choose ξ < 1 and α, and solve, using a multigrid method, $\mathrm{[}\mathrm{1}\mathrm{-}\mathit{\xi }\mathrm{(}\mathit{z}\mathrm{\left)}\mathrm{\right]}\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right)\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{P}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{+}{\mathit{\alpha }}^{\mathrm{2}}\mathit{P}\mathrm{=}\mathrm{0.}$

• Calculate the magnetic field by applying finite differences to ${B}\mathrm{=}\mathrm{\nabla }\mathrm{\times }\mathrm{[}\mathrm{\nabla }\mathrm{\times }\mathrm{(}\mathit{P}{ê}\mathit{z}\mathrm{)}\mathrm{+}\mathit{\alpha P}{ê}\mathit{z}\mathrm{]}\mathit{.}$

• Calculate the pressure directly using $\mathit{p}\mathrm{=}{\mathit{p}}_{\mathrm{b}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{-}\mathit{\xi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\frac{{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}\mathrm{·}$

• Use finite differences to calculate the density with $\mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{\mathit{g}}\left(\frac{\mathrm{d}\mathit{\xi }}{\mathrm{d}\mathit{z}}\frac{{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\mu }}_{\mathrm{0}}}\mathit{\xi }{B}\mathrm{·}\mathrm{\nabla }{\mathit{B}}_{\mathit{z}}\right)\mathit{.}$

The elegance of the NR99 formulation allows for an efficient numerical solution. We demonstrate the above scheme by calculating linear force-free and linear MHS equilibria. For the force-free case, we calculate an equilibrium with three footpoints and a null point. For the linear MHS case, we consider a region with two footpoints and investigate the effects of changing ξ. We consider profiles that cannot be treated, in the NR99 formalism, analytically and demonstrate how these can change the structure of the current distribution within the domain. We also show that the linear MHS equilibria are consistent with the Wilson effect for the size of the regions considered.

As well as demonstrating the scheme to be accurate, we show that the calculation time can be significantly reduced through parallelization on a GPU. This is achieved function by function and speed-ups of ×30 are realized. This result is significant as GPUs are inexpensive and readily available, unlike compute clusters and supercomputers. This fast and accurate scheme could be used for calculating equilibria in their own right or as input to MHD codes.

Acknowledgments

D.M. would like to thank the Royal Astronomical Society for financial support. J.M. acknowledges IDL support provided by STFC.

References

All Figures

	Fig. 1 Magnetic field lines of the linear force-free region are displayed in orange. A magnetogram of Bz is given at the base to indicate the positions of the sources.
In the text

	Fig. 2 A map of | B |, within the range given on the scale, at z = 0.1. The null point is clearly highlighted between the two negative polarities.
In the text

	Fig. 3 The profiles of ξ(z) from the cross-field current.
In the text

	Fig. 4 Magnetic field lines for the MHS case with ξ1. The field lines are coloured with | ∇ × B |.
In the text

	Fig. 5 Magnetic field lines for the MHS case with ξ2. The field lines are coloured with | ∇ × B |.
In the text

	Fig. 6 An isosurface of | j| = 0.03 for the MHS solution with ξ1.
In the text

	Fig. 7 An isosurface of | j| = 0.03 for the MHS solution with ξ2.
In the text

	Fig. 8 The non-dimensional density along a diagonal cut through the domain at z = 0.46875. The two troughs correspond to the locations of the magnetic footpoints. This profile is for the linear MHS equilibrium with ξ2.
In the text

	Fig. 9 The run times for the serial (blue) and parallel (green) versions of the multigrid code plotted against resolution. Each run time corresponds to the combined time for a 2D multigrid method calculation on the base of the domain and a 3D multigrid method calculation within the domain.
In the text