MAGNETO-FRICTIONAL MODELING OF CORONAL NONLINEAR FORCE-FREE FIELDS. I. TESTING WITH ANALYTIC SOLUTIONS

An ideal way to understand solar eruptive activities like flares and coronal mass ejections is to quantify the physical conditions such as density, temperature, velocity, and magnetic field in the most accurate way as a function of space and time. Knowing them all allows us to identify stability thresholds and the prevailing conditions during magnetic flux rope eruptions, magnetic reconnection, and particle acceleration. However, the low densities attained in the solar corona hamper our observational measurements of all physical parameters; in particular, it is challenging to get vector magnetic fields in a sufficiently large three-dimensional (3D) domain and time range. Since the magnetic field governs the structure and dynamics of the corona where the ratio of the gas pressure to the magnetic pressure (i.e., plasma β) is low (Gary 2001), knowing the vector magnetic field in the full 3D domain is prerequisite to reliably understanding the coronal atmosphere. The three components of the vector magnetic field are routinely observed in the photosphere (Lites 2000), although there are occasional direct measurements of the chromospheric and coronal magnetic field. For example, Lin et al. (2000) measured the line-of-sight coronal magnetic field directly by Zeeman splitting of the infrared Fe xiii 10747 Å forbidden spectral line. Lin et al. (2004) further obtained the orientation of the transverse magnetic field component using the same spectral line. Additionally, the Hanle effect is another promising technique to derive the chromospheric and coronal vector magnetic field (e.g., Stenflo 2013; Orozco Suárez et al. 2014; Schmieder et al. 2014; Supriya et al. 2014; Del Pino Alemán & Trujillo Bueno 2015; Li & Qu 2015).

Despite these successes, theoretical modeling to constrain the 3D magnetic field in large coronal volumes from observed boundary information will remain a necessary ingredient. If we consider static force-balanced equilibrium states where the Lorentz force dominates over the pressure gradient and gravity, the magnetic field must adopt the force-free state where the Lorentz force vanishes everywhere and magnetic pressure and tension balance each other. Including the divergence-free condition, any force-free magnetic field is described by the following equations:

$$\begin{eqnarray}&&({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0.\end{eqnarray} \tag{ 2 }$$

Equation (1) can be recast as ${\rm{\nabla }}\times {\boldsymbol{B}}=\alpha ({\boldsymbol{r}}){\boldsymbol{B}}$, where the torsional parameter $\alpha ({\boldsymbol{r}})$ is a scalar function of space ${\boldsymbol{r}}$ that is constant along any field line. If $\alpha ({\boldsymbol{r}})$ is zero or constant everywhere, which corresponds to a potential and linear force-free field respectively, the above equations are linear and can be solved by a Green function method or a Fourier transform method (Schmidt 1964; Nakagawa & Raadu 1972; Chiu & Hilton 1977; Seehafer 1978). When $\alpha ({\boldsymbol{r}})$ varies with position, there is no general analytic solution to the above equations because of nonlinearity.

Analytic or semi-analytic nonlinear force-free field (NLFFF) solutions can be derived under additional assumptions; for example, Low & Lou (1990) and Titov & Démoulin (1999) provide two classes of NLFFF in an axisymmetric setting. By considering only the resulting NLFFF field in a limited volume, both models can be transformed to seemingly 3D non-axisymmetric magnetic configurations, e.g., in Low & Lou (1990) by rotation and translation, or in Titov & Démoulin (1999) by considering only the part above a horizontal plane. The solutions of Low and Lou and the Titov–Démoulin models are commonly adopted as reference models to test numerical NLFFF extrapolations, which compute NLFFF solutions from given boundary and initial conditions.

Nowadays, various algorithms aim at reconstructing NLFFF models, which can be categorized as Grad–Rubin, upward integration, magneto-hydrodynamics (MHD) relaxation, optimization, or boundary integral equation methods. Grad & Rubin (1958) first proposed a method for computing force-free magnetic fields in fusion plasmas. Further applications and developments of the Grad–Rubin method have been made by Sakurai (1981), Amari et al. (1997, 2006, 2014), Wheatland (2004, 2006), Inhester & Wiegelmann (2006), and Gilchrist & Wheatland (2013, 2014). The vertical integration method computes the torsional parameter $\alpha ({\boldsymbol{r}})$ on the bottom boundary from observed vector magnetic fields. Then, it integrates over height to derive the vector magnetic field higher up (Nakagawa 1974). This suffers from the ill-posed nature of the governing equations, although various methods have been proposed to regularize the problem (Wu et al. 1985, 1990; Cuperman et al. 1989, 1990; Démoulin et al. 1992; Song et al. 2006). The MHD relaxation method in essence uses the momentum balance equation and the magnetic induction equation together, to reach a steady state (e.g., Chodura & Schlueter 1981). A variant of this method is adopted here, also known as the stress-and-relax (Roumeliotis 1996) or magneto-frictional method, as will be described in Section 2. The optimization method tries to minimize an objective functional, which sums over the force-free and divergence-free conditions. It was first proposed by Wheatland et al. (2000) and further implemented by Wiegelmann (2004, 2007) in Cartesian and spherical geometries, respectively. The spherical version has been tested and applied to observations by Tadesse et al. (2009, 2011). Another version of the optimization method has been implemented by J. McTiernan both in IDL⁴ and FORTRAN languages; its Cartesian version has been tested in Schrijver et al. (2006) and Metcalf et al. (2008), and its spherical version in Guo et al. (2012). The optimization method has also been implemented with a finite-element scheme by Inhester & Wiegelmann (2006), while most other codes adopt standard finite-difference schemes. Finite-element methodology is also commonly used in Grad–Shafranov solvers for (axisymmetric) tokamak equilibria, which balance Lorentz forces with pressure gradients within given outer boundaries. They allow for easy generalizations, e.g., including external gravity or equilibrium flows (e.g., Beliën et al. 2002; Petrie et al. 2007; Blokland & Keppens 2011). Finally, the boundary integral equation method was introduced by Yan & Sakurai (1997, 2000) to compute NLFFF models, and has been further developed and applied to observations (Yan et al. 2001; He & Wang 2006, 2008; Yan & Li 2006).

The advantages and disadvantages of the above algorithms have been discussed thoroughly in a series of tests using reference models. Schrijver et al. (2006) tested six different NLFFF algorithms in the Cartesian version with the solution of Low & Lou (1990). They found that the boundary conditions are crucial for the performance of the NLFFF algorithms. If all the six boundaries are provided by the reference model, the best-performing algorithm, which was the optimization method implemented by Wiegelmann (2004), could reconstruct the model with a mean vector error of about 2%. But if only the bottom boundary is provided by the reference model and the side and top boundaries use a potential field approximation or some other reasonable numerical extrapolation, only the magnetic field in an inner and lower volume can be reconstructed in agreement with the reference model. Metcalf et al. (2008) tested four NLFFF algorithms with a reference model derived by the flux rope insertion method (van Ballegooijen et al. 2000; van Ballegooijen 2004), which simulates forced photospheric and force-free chromospheric vector magnetic fields simultaneously. All tested codes could recover the magnetic flux rope and magnetic energy of the reference model when the force-free chromospheric boundary is provided. However, all NLFFF codes failed when the forced photospheric boundary data were provided, making preprocessing of the forced photospheric boundary data necessary, such as by the algorithm proposed in Wiegelmann et al. (2006b) and Fuhrmann et al. (2007).

There are ways to constrain the NLFFF modeling further by using additional observations other than the vector magnetic fields on the bottom boundary. The flux rope insertion method uses Hα, extreme-ultraviolet (EUV), or soft X-ray observations and the normal component of the magnetic field on the bottom boundary to construct an NLFFF model (van Ballegooijen 2004; Su et al. 2009). A magnetic flux rope partially constrained by the observations (for example, its footpoints, path, and sign of helicity) is inserted into the potential field computed from the normal magnetic field component. This configuration is relaxed to an NLFFF state with the magneto-frictional method. A series of such NLFFF models can be constructed by adjusting the parameters of the inserted magnetic flux rope. The final result is then determined by comparing selected field lines with EUV observations. Malanushenko et al. (2012, 2014) have developed a new NLFFF modeling method using a modified Grad–Rubin method (the so-called quasi-Grad–Rubin method). It only uses the normal magnetic field on the bottom boundary, but some $\alpha ({\boldsymbol{r}})$ distributions in the volume. In practice, these $\alpha ({\boldsymbol{r}})$ are derived along some trajectories of coronal loops fitted by the linear force-free field model.

Finally, we note that coronal magnetic field modeling in general also includes methods that do not adopt the force-free field assumption. For instance, non-force-free models have been developed with different assumptions (Zhao & Hoeksema 1993, 1994; Zhao et al. 2000; Hu & Dasgupta 2008; Hu et al. 2008). These models include the physical effects of pressure gradient and gravity, while the authors adopt additional assumptions. For example, Zhao & Hoeksema (1994) require that the direction of the electric current is perpendicular to gravity. Hu & Dasgupta (2008) require that the magnetic field can be expressed as the summation of three linear force-free fields with one of them being potential. Under those assumptions, the governing equations are simplified to linear problems. The solutions can be expressed as the summation involving spherical harmonics (Zhao & Hoeksema 1994) or the summation of linear force-free field solutions (Hu & Dasgupta 2008). Different from the above linear non-force-free models, full magneto-hydrostatic models try to balance all forces including the pressure gradient and gravity in a nonlinear way (Wiegelmann & Inhester 2003; Wiegelmann et al. 2007; Ruan et al. 2008).

In this paper, we document and thoroughly test the implementation of the magneto-frictional method in the open-source Message Passing Interface Adaptive Mesh Refinement (AMR) Versatile Advection Code (MPI-AMRVAC). MPI-AMRVAC has evolved into a versatile simulation code, which is equipped with a suite of spatial and temporal discretizations. It is designed especially for solving hyperbolic partial differential equations, which are common to multiple physics problems. The most relevant features of MPI-AMRVAC for the present effort are as follows. We take full advantage of (1) the block-based AMR technique; (2) the parallel computing ability (domain decomposition on a uniform grid or block-based decomposition on an AMR grid); (3) the options of different coordinate systems (Cartesian or spherical coordinates); (4) the potential and linear force-free models from magnetogram data in Cartesian coordinates; and (5) the potential field source surface (PFSS) model from full-disk synoptic magnetograms (Keppens et al. 2003, 2012; Porth et al. 2014). These are all useful for constructing realistic NLFFF models, which need to combine as large a field of view and as high a spatial resolution as possible. This requirement naturally introduces a big data problem, and can only be handled by the parallel and AMR abilities of MPI-AMRVAC. And the potential and linear force-free field models could serve as the initial condition for the NLFFF modeling.

The tests we perform revisit semi-analytic NLFFF solutions that have been used in previous NLFFF code comparisons, or have been tested with previous magneto-frictional implementations where their ability to obtain force-free volume-filling fields is well documented. Our implementation of the magneto-frictional method in a parallel block-adaptive framework is motivated by problems discussed in DeRosa et al. (2009), which indicate that the field of view for NLFFF modeling should be as large as possible to include all the possible magnetic field line connections. When the field of view is large and may include many active regions, the curvature of the Sun can no longer be neglected (Guo et al. 2012). Then, one must use spherical coordinates rather than Cartesian ones. Also, the spatial resolution should be preferentially high where electric current density enters the bottom boundary, most notably in active regions. In such a case, traditional uniform grids and single-processor executions will be challenging, which is where MPI-AMRVAC may offer a viable solution.

We first describe the implementation of the magneto-frictional module in MPI-AMRVAC in Section 2. Next, we test the performance of the code with the solution of Low & Lou (1990) in Cartesian coordinates in Section 3.1 and with the models of Titov & Démoulin (1999) in Section 3.2. Then, we test the NLFFF modeling in spherical coordinates in Section 3.3. The AMR ability of the magneto-frictional module in MPI-AMRVAC is demonstrated in Section 3.4. Finally, we summarize the test results and draw our conclusion in Section 4.

The magneto-frictional method can be seen as a special case of the more general MHD relaxation method, which tries to solve the momentum equation and magnetic induction equation to obtain an equilibrium state. Taking into account dissipative and resistive terms, we have in full

$$\begin{eqnarray}&&\rho \left(\displaystyle \frac{\partial }{\partial t}+{\boldsymbol{v}}\cdot {\rm{\nabla }}\right){\boldsymbol{v}}={\boldsymbol{J}}\times {\boldsymbol{B}}-{\rm{\nabla }}p+\rho {\boldsymbol{g}}+{\boldsymbol{D}}({\boldsymbol{v}}),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})+{\eta }_{{\rm{m}}}{{\rm{\nabla }}}^{2}{\boldsymbol{B}},\end{eqnarray} \tag{ 4 }$$

where the electric current density ${\boldsymbol{J}}={\rm{\nabla }}\times {\boldsymbol{B}}/{\mu }_{0}$, the dissipative term ${\boldsymbol{D}}({\boldsymbol{v}})$ is a function of the (derivatives of) velocity ${\boldsymbol{v}}$, and ρ, p, ${\boldsymbol{g}}$, ${\eta }_{{\rm{m}}}$, and ${\mu }_{0}$ stand for the density, gas pressure, gravitational acceleration vector, magnetic diffusivity, and vacuum permeability, respectively. Note that ${\eta }_{{\rm{m}}}=1/(\sigma {\mu }_{0})$, and σ is the conductivity. The dissipative term can be written as a friction term proportional to the local speed such that ${\boldsymbol{D}}({\boldsymbol{v}})=-\nu {\boldsymbol{v}}$ (Yang et al. 1986; Roumeliotis 1996; Valori et al. 2005, 2007, 2010) or in a more physical viscosity form where ${\boldsymbol{D}}({\boldsymbol{v}})={\rm{\nabla }}\cdot (\nu \rho {\rm{\nabla }}{\boldsymbol{v}})$ (McClymont & Mikić 1994; Amari et al. 1996, 1997; McClymont et al. 1997), where ν is a free parameter to control the speed of the dissipative process. The MHD relaxation method can be used not only for NLFFF modeling but also for more general magneto-hydrostatic modeling (e.g., Chodura & Schlueter 1981).

Here, we consider the NLFFF model. Therefore, the inertia, pressure gradient, and gravity forces are omitted in Equation (3). We adopt the magneto-frictional form for the dissipative term, and the momentum equation is reduced to

$$\begin{eqnarray}&&{\boldsymbol{v}}=\displaystyle \frac{1}{\nu }{\boldsymbol{J}}\times {\boldsymbol{B}}.\end{eqnarray} \tag{ 5 }$$

In practice, we use the same expressions for the velocity ${\boldsymbol{v}}$ as in Valori et al. (2005, 2007) to optimize the relaxation process:

$$\begin{eqnarray}{\boldsymbol{v}} & = & \displaystyle \frac{{c}_{c}{c}_{y}}{\max (| \hat{{\boldsymbol{v}}}| )}\displaystyle \frac{{\rm{\Delta }}x}{{\rm{\Delta }}t}\hat{{\boldsymbol{v}}}{f}_{w}({\boldsymbol{x}}),\\ \hat{{\boldsymbol{v}}} & = & {\rm{\Delta }}x\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{| {\boldsymbol{B}}{| }^{2}},\end{eqnarray} \tag{ 6 }$$

where c_c and c_y are two free parameters to control the magnitude of the magneto-frictional velocity, the function max() evaluates the instantaneous maximum value of $\hat{{\boldsymbol{v}}}$ in the computational box at each iteration, and ${\rm{\Delta }}x$ and ${\rm{\Delta }}t$ are the spatial and temporal discretization steps. According to Equation (6), the magneto-frictional velocity is proportional to the Lorentz force, whose coefficient is a function of space and time. The coefficient ${c}_{c}/\max (| \hat{{\boldsymbol{v}}}| )$ is introduced as a heuristic Courant–Friedrichs–Lewy (CFL) criterion as explained in Valori et al. (2007). Provided the parameters ${c}_{c}\leqslant 1$ and ${c}_{y}\leqslant 1$, the magneto-frictional velocity ${\boldsymbol{v}}$ will always be less than or equal to the information propagation velocity, ${\rm{\Delta }}x/{\rm{\Delta }}t$, no matter what time and space discretization schemes are used. Considering the relaxation speed and the stability of the code simultaneously, we use c_c = 0.5 and c_y = 0.2 for most of the test and application cases unless specified otherwise. The coefficient $1/| {\boldsymbol{B}}{| }^{2}$ is introduced to speed up the relaxation process in weak-field regions (Yang et al. 1986). Additionally, the velocity field in the computational domain has been multiplied by a weighting function ${f}_{w}({\boldsymbol{x}})\in [0,1]$, which is equal to one throughout most of the box, but decreases from unity in a boundary layer: this expands over a distance L_j, which is 5% of the length of the computational box in direction j. Within that layer, ${f}_{w}({\boldsymbol{x}})$ drops to 0 in a parabolic profile toward the four side and top boundaries. The weighting function decreases the velocity, which minimizes changes in the four side and top boundaries by numerical extrapolations or symmetries, or limits large velocities to a lower value when the boundaries are provided by analytic data. Since the physical quantities are defined on the cell center while the weighting function equals zero on the cell edge of the outermost layer, the change in magnetic field is not exactly zero but is limited to a small value on the outermost layer in the physical domain.

The magnetic field in the computational domain is updated by iterating on the induction equation. The resistive term in this equation can introduce magnetic reconnection and changes in magnetic topology, which might be necessary for a successful relaxation, because the final state of an NLFFF usually has a different magnetic topology to its initial state, such as a potential field. On the other hand, the resistive term has to vanish when the relaxation reaches a steady state. Valori et al. (2005) found that the resistive term might pollute the desired magnetic configuration when there are complex mixtures of magnetic polarities on the bottom boundary. Therefore, we omit the resistive term. The desired change in magnetic topology is still guaranteed by two other factors, namely, the numerical dissipation term in the spatial discretization schemes we use (see below) and the fact that we face discontinuous changes between the field values on the bottom boundary from the magnetogram and the initial condition using a potential extrapolation.

The solenoidal condition of the magnetic field is controlled by adding a diffusive term in the induction equation (Keppens et al. 2003):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})+\delta {\rm{\nabla }}({\rm{\nabla }}\cdot {\boldsymbol{B}}),\end{eqnarray} \tag{ 7 }$$

where δ is a free parameter to control the diffusion speed. It is implemented as $\delta ={c}_{d}{\rm{\Delta }}{l}^{2}/{\rm{\Delta }}t$, where ${\rm{\Delta }}{l}^{2}$ is the harmonic mean (divided by the number of directions) of the square of the spatial steps in each direction, so ${\rm{\Delta }}{l}^{2}={\left(\tfrac{1}{{\rm{\Delta }}{x}^{2}}+\tfrac{1}{{\rm{\Delta }}{y}^{2}}+\tfrac{1}{{\rm{\Delta }}{z}^{2}}\right)}^{-1}$ in the Cartesian coordinate system and ${\rm{\Delta }}{l}^{2}={\left(\tfrac{1}{{\rm{\Delta }}{r}^{2}}+\tfrac{1}{{(r{\rm{\Delta }}\theta )}^{2}}+\tfrac{1}{{(r\sin \theta {\rm{\Delta }}\phi )}^{2}}\right)}^{-1}$ in the spherical coordinate system. The free parameter c_d has to satisfy $0\leqslant {c}_{d}\leqslant 2$ to conform with the CFL time-step limitation (van der Holst & Keppens 2007). We use c_d = 0.1 throughout the tests and applications unless specified otherwise.

For the magneto-frictional module in MPI-AMRVAC, we selected to compare two spatial discretization schemes in particular, one being the total variation diminishing Lax–Friedrichs (TVDLF) method that we often exploited in full dynamical MHD simulations, and the other a (Lax–Friedrichs stabilized) fourth-order central difference scheme. The TVDLF scheme in MPI-AMRVAC is described, e.g., in Keppens et al. (2012), and here in essence (1) uses a user-selected limited reconstruction strategy from cell center to cell edge to get a left and right cell-edge evaluation of the magnetic field ${{\boldsymbol{B}}}_{i+1/2}^{L,R}$ and velocity field ${{\boldsymbol{v}}}_{i+1/2}^{L,R};$ (2) uses them to evaluate the corresponding left and right physical flux ${\boldsymbol{vB}}-{\boldsymbol{Bv}}$ at the cell edges; and (3) takes as the numerical flux across the cell edge $i+1/2$ the expression

$$\begin{eqnarray}{F}_{i+1/2}^{\mathrm{LF}} & = & \displaystyle \frac{1}{2}[F({{\boldsymbol{B}}}_{i+1/2}^{L},{{\boldsymbol{v}}}_{i+1/2}^{L})+F({{\boldsymbol{B}}}_{i+1/2}^{R},{{\boldsymbol{v}}}_{i+1/2}^{R})\\ & & -{\epsilon }^{\mathrm{LF}}{c}_{i+1/2}^{\max }({U}_{i+1/2}^{R}-{U}_{i+1/2}^{L})].\end{eqnarray} \tag{ 8 }$$

In this expression, F denotes the (vector component of) the physical flux, U is the corresponding B_x, B_y, or B_z component, and we discuss the other parameters ${\epsilon }^{\mathrm{LF}}$ and ${c}^{\max }$ below. In the first reconstruction stage, many limited reconstructions can be used to get the cell interface from cell-center values. Here we will always adopt the "Koren" limiter (Koren 1993), which can achieve third-order accuracy for smooth regions with a stencil requiring two ghost cells only. We have verified that the more recently proposed third-order limiter "cada3" (Čada & Torrilhon 2009; Keppens & Porth 2014) yields essentially the same results for these magneto-frictional problems as the "Koren" limiter.

The second base spatial discretization scheme is a simple generalization of the (second-order) TVDLF scheme. This scheme combines a fourth-order central difference with the stabilizing local Lax–Friedrichs term (referred to as CD4LF hereafter), by taking

$$\begin{eqnarray}{F}_{i+1/2}^{\mathrm{CD}} & = & \displaystyle \frac{1}{12}[-F({{\boldsymbol{B}}}_{i+2},{{\boldsymbol{v}}}_{i+2})+7F({{\boldsymbol{B}}}_{i+1},{{\boldsymbol{v}}}_{i+1})\\ & & +7F({{\boldsymbol{B}}}_{i},{{\boldsymbol{v}}}_{i})-F({{\boldsymbol{B}}}_{i-1},{{\boldsymbol{v}}}_{i-1})]\\ & & -\displaystyle \frac{1}{2}{\epsilon }^{\mathrm{LF}}{c}_{i+1/2}^{\max }({U}_{i+1/2}^{R}-{U}_{i+1/2}^{L}).\end{eqnarray} \tag{ 9 }$$

Without the dissipative local Lax–Friedrichs term, this would correspond to a cell i update using the familiar fourth-order centering, since then

$$\begin{eqnarray}{F}_{i+1/2}-{F}_{i-1/2} & = & \displaystyle \frac{1}{12}[-F({{\boldsymbol{B}}}_{i+2},{{\boldsymbol{v}}}_{i+2})+8F({{\boldsymbol{B}}}_{i+1},{{\boldsymbol{v}}}_{i+1})\\ & & +8F({{\boldsymbol{B}}}_{i-1},{{\boldsymbol{v}}}_{i-1})-F({{\boldsymbol{B}}}_{i-2},{{\boldsymbol{v}}}_{i-2})].\end{eqnarray} \tag{ 10 }$$

The local Lax–Friedrichs term not only acts to stabilize the scheme, but its corresponding numerical diffusion allows the magnetic topology to change in the process of time-advancing. As explained in Keppens et al. (2012), one might gradually reduce the dissipation term by multiplying by a factor, ${\epsilon }^{\mathrm{LF}}\leqslant 1$, to improve the local accuracy of steady-state solutions. Here, we kept the factor constant during the time-advancing in the model of Low and Lou and Titov–Démoulin, because we find that decreasing the dissipation term for matching these semi-analytic solutions does not improve the overall performance.

For the maximal physical propagation speed, we use ${c}_{i+1/2}^{\max }=\max ({c}_{\max }^{L},{c}_{\max }^{R})$, where we adopt in each direction j as a maximal propagation speed ${c}_{\max }=| {v}_{j}| +\tfrac{B}{\sqrt{{\mu }_{0}\rho }}$ where v_j is the magneto-frictional velocity in the jth direction, and where we use the local Alfvén speed obtained with a reference "density" value ρ. In this way, the dissipation term depends on the value of this density ρ. For the solar application, it makes sense to normalize all quantities with typical values in the solar atmosphere, such that the typical length, density, and magnetic field are chosen as 10⁹ cm, $2.34\times {10}^{-15}$ g cm⁻³, and 2 G, respectively. The normalized density is then uniform and equal to 1 in the full computational domain. As an auxiliary parameter that is used for determining the time step, the density does not change in the process of the magneto-frictional relaxation, since it does not appear in the governing equations. This may be changed later on, when we plan to use stratified densities when using NLFFF extrapolations as initial conditions for full MHD evolutions of specific solar events.

Equation (7) can be advanced by various explicit schemes, such as the single-step Euler, two-step predictor corrector, and the three-step Runge–Kutta methods, all available in MPI-AMRVAC. To determine the time step for advancing Equation (7), we adopt ${\rm{\Delta }}t={c}_{c}{\rm{\Delta }}x/{c}_{\max }$, with the same c_max (this is meant to consider all the directions j). The exact value of ${\rm{\Delta }}t$ plays only a minor role in the time-advancing. Combining Equations (6) and (7), the update of the magnetic field ${\boldsymbol{B}}$ does not depend on the time step if we do not consider the dissipation term in the spatial discretization schemes.

The initial condition for the magneto-frictional method is always provided by the potential field computed from the normal component on the bottom boundary. The Green function method proposed by Chiu & Hilton (1977) and the PFSS model (Altschuler & Newkirk 1969; Schatten et al. 1969; Schrijver & De Rosa 2003) in the Cartesian and spherical coordinate systems, respectively, are adopted to compute the potential field. Both methods have been implemented in MPI-AMRVAC as demonstrated in Porth et al. (2014). In principle, linear force-free field could also serve as initial condition for Cartesian cases. We do not consider it here, since the state of a linear force-free field with a constant torsional parameter α in the solar atmosphere is very rare, and observations show that α is highly nonlinear in active regions. Some tests also show that using a linear force-free field as the initial condition usually provides worse NLFFF models (e.g., Valori et al. 2010). At the same time, potential field has proven to be a good zero-order approximation for large-scale magnetic configurations.

Similar to previous studies, we will use and compare different ways to specify the boundary conditions, which in our code are always handled by introducing and prescribing (field and velocity) values in ghost cells. We typically use two ghost cell layers (although more cells can be used depending on the reconstruction strategy). In the test cases where all the information on the magnetic field is known, all the six boundaries can be provided for all edges at x_min, x_max, y_min, y_max, z_min, and z_max for Cartesian coordinates or r_min, r_max, θ_min, θ_max, ϕ_min, and ${\phi }_{\max }$ for spherical coordinates. This ideal setting excludes all the effects caused by the boundary conditions and is very useful for testing the best performance of different numerical schemes or for finding optimized free parameters. To simulate cases of solar-like application, where only the vector magnetic field on the bottom boundary is known, the boundary conditions at all other boundaries can be provided by the initial potential field or by adopting simple numerical extrapolations in both Cartesian and spherical coordinates, or by using periodicity in the spherical coordinate when $\theta \in [0,\pi ]$ or $\phi \in [0,2\pi ]$. Since there are many choices of the boundary conditions, we will explain more details in Section 3 and document performances for several combinations.

We note that Low (2013) questioned the magneto-frictional method when the velocity takes the form of Equation (5). The algebraic form of the velocity could bring unmeaningful results, such as a contradiction of the initial conditions with the boundary conditions, and can restrict the magnetic topology. However, not all arguments apply here, because the numerical implementation of the velocity as expressed in Equation (6) is not uniform in space and time, and decreases to zero toward the side and top boundaries in a prescribed (smoothed) fashion. The arguments in Low (2013) are valid when the MHD relaxation process is interpreted as a physical evolution process. However, we use the shock-capturing (hence discontinuity-capable) properties of the schemes to handle, e.g., the dramatic changes in boundary and initial conditions, which can be incompatible changeovers between NLFFF and potential field. The velocity computed at such boundaries is artificial, and the numerical diffusion can change the magnetic topology. Therefore, the relaxation process is meant as a numerical iteration process to converge from initial conditions to an NLFFF model that matches the boundary conditions. The intermediate iterations of the velocity computed by Equation (6) and the magnetic field computed by Equation (7) do not possess a physical meaning. The magnetic topology does change in the process of the magneto-frictional relaxation, as will be clear in the examples.

Finally, we remark that the MHD relaxation method (such as the magneto-frictional approach) has also been implemented by other authors more recently. Jiang et al. (2011) solved the full MHD equations with a new numerical method, namely, the conservation element and solution element method, to reconstruct an approximate NLFFF model. Jiang & Feng (2012) further implemented and tested the magneto-frictional method with the solution of Low & Lou (1990) and the reference model of van Ballegooijen (2004). They extended the magneto-frictional method to spherical coordinates (Jiang et al. 2012). In particular, they adopted a so-called Yin–Yang grid to avoid the singularity at the poles of the spherical coordinate system. Then, the full-sphere extrapolation is tested with the solution of Low and Lou. Finally, the MHD relaxation method has also been implemented and applied to observations by Inoue et al. (2011, 2012, 2014).

We adopt two sets of semi-analytic NLFFF solutions, namely the solutions of Low & Lou (1990) and the models of Titov & Démoulin (1999), to test the magneto-frictional module in MPI-AMRVAC. The solutions of Low and Lou can be obtained by solving a second-order ordinary differential equation with different eigenvalues and eigenfunctions. The solutions have a singular line at their symmetric axis, which can be avoided by restricting them to a region without this singular line. The solutions of Low and Lou are suitable for NLFFF code tests although electric currents are distributed relatively smoothly in the whole volume. Therefore, they do not necessarily represent realistic conditions in the solar atmosphere. We remark that there is a singular line along the symmetry axis of the solution of Low and Lou. The magnetic field is divergent and the divergence-free and force-free conditions are not satisfied along this line. It might therefore not be appropriate to use the solution of Low and Lou as a test for the full-sphere extrapolation if this singular line is not excluded from the computational volume. Besides, the magnetic topology of the solution of Low and Lou is relatively simple compared to a potential field with the same normal field on the boundary. On the contrary, the Titov–Démoulin models represent more realistic magnetic configurations in the solar atmosphere, where electric currents are more concentrated in active regions. The magnetic topology is very different from the corresponding potential field. Using these known NLFFF solutions, the performance of the code is tested with different coordinate systems (Cartesian or spherical), employed meshes (uniform or AMR grids), boundary conditions (fully prescribed or solar-like), spatial differentiation schemes (CD4LF or TVDLF), and time-advancing methods (one, two, or three steps). To distinguish the test cases with the aforementioned options, we adopt a set of labeling codes as listed in Table 1.

Table 1.  Explanation of the Labeling Codes

Option	Label	Explanation
Model	LL	Solution of Low and Lou
	TD-NO	Titov–Démoulin model without any hyperbolic flux tube (HFT)
	TD-LO	TD model with a low HFT
	TD-HI	TD model with a high HFT
Coordinate	CAR	Cartesian coordinates
	SPH	Spherical coordinates
Grid	UNI	Uniform grid
	AMR	Adaptive mesh refinement grid
Boundary	FUL	Fully prescribed boundary conditions. In the Cartesian coordinates, all the six boundaries for the magnetic field at ${x}_{\min ,\max },$ ${y}_{\min ,\max },$ ${z}_{\min ,\max }$ are fully provided by the semi-analytic data. In the spherical coordinates, the semi-analytic magnetic field is provided at ${\theta }_{\min ,\max },$ ${r}_{\min ,\max }$, and the boundary condition is periodic at ${\phi }_{\min ,\max }$ when $\phi \in [0,2\pi ]$. We note that each side, top, and bottom boundary includes two ghost layers.
	BOT	Solar-like boundary conditions. In the Cartesian coordinates, the boundaries for the magnetic field at ${x}_{\min ,\max },$ ${y}_{\min ,\max }$, zmax, and the outer zmin ghost layer are provided by the zero-gradient extrapolation, and the boundary at the inner zmin ghost layer is provided by the semi-analytic data. In the spherical coordinates, the boundaries for the magnetic field at ${\theta }_{\min ,\max },$ rmax, and the outer rmin ghost layer are provided by the zero-gradient extrapolation, those at ${\phi }_{\min ,\max },$ where $\phi \in [0,2\pi ]$, are periodic, and that at the inner rmin ghost layer is provided by the semi-analytic data.
	OTH	Other boundary conditions that will be explained for each case
Spatial differentiation	CD4LF	The fourth-order central difference scheme with a local Lax–Friedrichs dissipation term (default)
	TVDLF	The total variation diminishing Lax–Friedrichs scheme
Time-advancing	ONE	One-step Euler type (default)
	TWO	Two-step predictor corrector method
	THR	Three-step Runge–Kutta method

Download table as:  ASCIITypeset image

To quantitatively evaluate the basic properties of an NLFFF, which is supposed to be divergence-free and force-free, we adopt the following two metrics proposed by Wheatland et al. (2000). The divergence-free condition is described by the volume-weighted average of the absolute value of the fractional magnetic flux change:

$$\begin{eqnarray}&&\langle | {f}_{i}| \rangle =\displaystyle \frac{{\displaystyle \sum }_{i}| {f}_{i}| {\rm{\Delta }}{V}_{i}}{{\displaystyle \sum }_{i}{\rm{\Delta }}{V}_{i}},\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{f}_{i}=\displaystyle \frac{{({\rm{\nabla }}\cdot {\boldsymbol{B}})}_{i}{\rm{\Delta }}{V}_{i}}{{B}_{i}{A}_{i}},\end{eqnarray} \tag{ 12 }$$

and ${\rm{\Delta }}{V}_{i}$, A_i, and B_i are the cell volume, cell surface area, and the magnitude of the magnetic field ${\boldsymbol{B}}$ at the grid point i. The force-free condition is described by the current- and volume-weighted average of the sine of the angle between the magnetic field ${\boldsymbol{B}}$ and the current density ${\boldsymbol{J}}$:

$$\begin{eqnarray}&&{\sigma }_{J}=\displaystyle \frac{\sin {\theta }_{i}{\displaystyle \sum }_{i}{J}_{i}{\rm{\Delta }}{V}_{i}}{{\displaystyle \sum }_{i}{J}_{i}{\rm{\Delta }}{V}_{i}},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&\sin {\theta }_{i}=\displaystyle \frac{| {{\boldsymbol{J}}}_{i}\times {{\boldsymbol{B}}}_{i}| }{{J}_{i}{B}_{i}},\end{eqnarray} \tag{ 14 }$$

and J_i is the magnitude of ${\boldsymbol{J}}$ at grid point i. The ghost layers are always excluded in computing the metrics; therefore, only the physical computational domain is considered. We note that $| {f}_{i}|$ is weighted by ${\rm{\Delta }}{V}_{i}$ in Equation (11) and $\sin {\theta }_{i}$ is weighted not only by J_i but also by ${\rm{\Delta }}{V}_{i}$ in Equation (13), since the cell volume changes with position in spherical coordinate or AMR grids. In the case of Cartesian coordinate and uniform girds, Equations (11) and (13) are equivalent to the expression in Wheatland et al. (2000). The magnetic field divergence, ${\rm{\nabla }}\cdot {\boldsymbol{B}}$, and the electric current density, ${\boldsymbol{J}}$, are evaluated numerically with the second-order central difference scheme. We note that ${\rm{\nabla }}\cdot {\boldsymbol{B}}$ in spherical geometry is computed by an integral numerical scheme, $\tfrac{{\int }_{V}{\rm{\nabla }}\cdot {\boldsymbol{B}}{dV}}{V}$, in each cell. This scheme differs little from the central difference scheme except in the polar regions. We also adopt the following quantity to describe the average angle between ${\boldsymbol{B}}$ and ${\boldsymbol{J}}$:

$$\begin{eqnarray}&&{\theta }_{J}={\sin }^{-1}{\sigma }_{J}.\end{eqnarray} \tag{ 15 }$$

If the magnetic field is perfectly divergence-free and force-free, all the metrics $\langle | {f}_{i}| \rangle$, ${\sigma }_{J}$, and ${\theta }_{J}$ will be equal to zero.

In the test cases, it is necessary to quantify how the reconstructed magnetic field ${\boldsymbol{b}}$ is in agreement with the reference (known) magnetic field ${\boldsymbol{B}}$. The following metrics have been used in Schrijver et al. (2006) and Metcalf et al. (2008) for such a purpose, namely the vector correlation metric:

$$\begin{eqnarray}&&{C}_{\mathrm{vec}}=\displaystyle \frac{{\displaystyle \sum }_{i}{{\boldsymbol{B}}}_{i}\cdot {{\boldsymbol{b}}}_{i}}{{\left({\displaystyle \sum }_{i}{B}_{i}^{2}{\displaystyle \sum }_{i}{b}_{i}^{2}\right)}^{1/2}},\end{eqnarray} \tag{ 16 }$$

the Cauchy–Schwartz metric:

$$\begin{eqnarray}&&{C}_{\mathrm{CS}}=\displaystyle \frac{1}{M}{\displaystyle \sum }_{i}\displaystyle \frac{{{\boldsymbol{B}}}_{i}\cdot {{\boldsymbol{b}}}_{i}}{{B}_{i}{b}_{i}},\end{eqnarray} \tag{ 17 }$$

the normalized vector error:

$$\begin{eqnarray}&&{E}_{{\rm{n}}}=\displaystyle \frac{{\displaystyle \sum }_{i}| {{\boldsymbol{b}}}_{i}-{{\boldsymbol{B}}}_{i}| }{{\displaystyle \sum }_{i}{B}_{i}},\end{eqnarray} \tag{ 18 }$$

and the mean vector error:

$$\begin{eqnarray}&&{E}_{{\rm{m}}}=\displaystyle \frac{1}{M}{\displaystyle \sum }_{i}\displaystyle \frac{| {{\boldsymbol{b}}}_{i}-{{\boldsymbol{B}}}_{i}| }{{B}_{i}}.\end{eqnarray} \tag{ 19 }$$

When the reconstructed magnetic field ${\boldsymbol{b}}$ is identical to the reference magnetic field ${\boldsymbol{B}}$, ${C}_{\mathrm{vec}}=1$, ${C}_{\mathrm{CS}}=1$, ${E}_{{\rm{n}}}=0$, and ${E}_{{\rm{m}}}=0$. When ${\boldsymbol{b}}$ deviates from ${\boldsymbol{B}}$, the metric C_vec decreases to 0, C_CS decreases to −1, while E_n and E_m could increase to infinity. To be compared with each other easily, $1-{E}_{{\rm{n}}}$ and $1-{E}_{{\rm{m}}}$ will be used. Then, the closer to 1 all these four metrics are, the better the two magnetic fields match each other. Finally, the magnetic energy fraction of magnetic field ${\boldsymbol{b}}$ compared to ${\boldsymbol{B}}$ is defined as

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{\displaystyle \sum }_{i}{b}_{i}^{2}{\rm{\Delta }}{V}_{i}}{{\displaystyle \sum }_{i}{B}_{i}^{2}{\rm{\Delta }}{V}_{i}}.\end{eqnarray} \tag{ 20 }$$

If the two fields are identical, the magnetic energy fraction would be 1.

3.1. Solution of Low and Lou in Cartesian Coordinates

The solution of Low and Lou adopted in the following test cases is obtained by using the solution to Equation (5) in Low & Lou (1990) corresponding to their parameters n = 1 and eigenvalue order m = 1. This solution is essentially two-dimensional and can be extended to 3D by the assumption of rotational symmetry. Therefore, the 3D solution of Low and Lou in the local coordinate system is axisymmetric. The symmetry property is reduced after the 3D local solution gets transformed to a physical coordinate system by rotation and translation. In that process, we employ a rotation about the y-axis counterclockwise by an angle ${\rm{\Phi }}=\pi /10$. This is followed by a shift of the origin of the local coordinate system to $(-{l}_{x},0,0)$ in the physical coordinate system with ${l}_{x}=-0.25$. We then select a subdomain of the solution of Low and Lou in a rectangular box with $x\in [1.0,3.0]$, $y\in [-1.0,1.0]$, and $z\in [-0.8,1.2]$ in the physical coordinate. The box is resolved into 100³ grid points (not including ghost cells). The units of length and magnetic field strength in the various steps are arbitrary and we normalize and scale the solution of Low and Lou in the physical domain by a constant factor to ensure that the maximum value of $| {B}_{z}|$ is 500 G while our grid cell size corresponds to ${\rm{\Delta }}x={\rm{\Delta }}y={\rm{\Delta }}z={10}^{8}\,{\rm{cm}}$. The parameters chosen in this study are not the same as in other works with a similar goal. This makes the detailed comparison of the metrics between different articles slightly less significant. However, we believe that the metrics presented here are representative of the overall performance of our code when applied to the family of equilibria of Low and Lou.

First, we verify the "convergence" of the magneto-frictional method with the following settings, hereafter referred to as model LL-CAR-UNI-FUL. Here, convergence is verified by showing the evolution of the divergence-free and force-free metrics during the iterations. The initial condition is provided by the potential magnetic field computed from the bottom boundary of the solution of Low and Lou. The three components of the magnetic field are prescribed on all the six boundaries by the solution of Low and Lou. The CD4LF scheme and one-step time-advancing method are adopted. Figure 1(a) shows the evolution of $\langle | {f}_{i}| \rangle$ and ${\sigma }_{J}$ during the iteration. The divergence-free metric $\langle | {f}_{i}| \rangle$ increases in the first $2\times {10}^{3}$ steps and then decreases until we reach the enforced maximum iteration step number of $6\times {10}^{4}$. The force-free metric ${\sigma }_{J}$ decreases the whole way and reaches ${ \mathcal O }(0.01)$ at the end. Note that the reference solutions themselves are not perfect either, in the sense that using the same metrics for quantifying their force-freeness yields similar values. However, these metrics do not always decrease monotonically, as will be shown clearly in the following cases, although the aim is to decrease during the iteration process. This is one reason to fix the maximum step number as a termination condition, instead of an automatic criterion based on relative changes of $\langle | {f}_{i}| \rangle$ or ${\sigma }_{J}$.

Zoom In Zoom Out Reset image size

Figures 2(a) and (b) show the comparison of the initial condition (the potential field) with the solution of Low and Lou. They have obvious differences: the potential field lines are nearly perpendicular to the polarity inversion line, while the field lines in the solution of Low and Lou are highly sheared. After the magneto-frictional iteration, the magnetic field lines of the reconstructed field are almost identical to the original semi-analytic solution (see Figures 2(c) and (d)). The divergence-free and force-free metrics of the reference solution of Low and Lou and for the reconstructed field are listed in Table 2 as models LL-CAR-UNI and LL-CAR-UNI-FUL, respectively. The metrics are computed in two different volumes of the computational box, i.e., the entire 100³ one and an inner 50³ one excluding the outermost quarter in all four sides and the top half in the vertical direction. From Table 2, we find that the solution of Low and Lou is well reconstructed from the potential field. The five metrics C_vec, C_CS, $1-{E}_{{\rm{n}}}$, $1-{E}_{{\rm{m}}}$, and for model LL-CAR-UNI-FUL in both the entire and inner volume are comparable to the best results as listed in Tables 1 and 2 in Schrijver et al. (2006).

Zoom In Zoom Out Reset image size

Table 2.  Metrics for the Solution of Low and Lou in the Cartesian Coordinates and Uniform Grid and the Nonlinear Force-free Fields Derived with Different Boundary Conditions

Modela	$\langle | {f}_{i}| \rangle$	${\sigma }_{J}$	${\theta }_{J}$	Cvec	CCS	$1-{E}_{{\rm{n}}}$	$1-{E}_{{\rm{m}}}$
	$({10}^{-4})$	$({10}^{-2})$	(deg)
The metrics are computed in the entire volume:
LL-CAR-UNI	0.04	0.19	0.11	⋯	⋯	⋯	⋯	1
LL-CAR-UNI-FUL	0.63	1.36	0.78	1.00	1.00	0.99	0.98	1.00
LL-CAR-UNI-OTH1	1.23	4.83	2.77	1.00	1.00	0.97	0.97	1.05
LL-CAR-UNI-OTH2	0.96	2.21	1.27	1.00	1.00	0.97	0.96	1.00
LL-CAR-UNI-OTH3	10.38	15.92	9.16	0.98	0.88	0.75	0.57	0.88
LL-CAR-UNI-OTH4	9.69	17.72	10.21	0.98	0.87	0.73	0.55	0.93
LL-CAR-UNI-OTH5	8.67	16.39	9.43	0.97	0.86	0.72	0.54	0.87
LL-CAR-UNI-BOT	0.98	6.16	3.53	0.99	0.92	0.81	0.64	0.99
LL-CAR-UNI-OTH6	0.90	5.38	3.08	0.99	0.91	0.80	0.63	0.98
The metrics are computed in the inner volume:
LL-CAR-UNI	0.10	0.22	0.12	⋯	⋯	⋯	⋯	1
LL-CAR-UNI-FUL	0.57	1.61	0.92	1.00	1.00	0.99	0.98	1.00
LL-CAR-UNI-OTH1	2.89	5.39	3.09	1.00	1.00	0.95	0.94	1.07
LL-CAR-UNI-OTH2	0.92	2.43	1.39	1.00	1.00	0.97	0.95	1.01
LL-CAR-UNI-OTH3	3.55	3.63	2.08	1.00	0.97	0.90	0.82	0.97
LL-CAR-UNI-OTH4	4.26	6.70	3.84	0.99	0.96	0.88	0.79	1.04
LL-CAR-UNI-OTH5	3.37	4.57	2.62	0.99	0.95	0.86	0.76	0.97
LL-CAR-UNI-BOT	3.00	5.69	3.26	1.00	0.99	0.92	0.88	1.04
LL-CAR-UNI-OTH6	2.83	5.67	3.25	1.00	0.99	0.92	0.87	1.03

Note.

^aRefer to Table 1 for the explanation of the models. The boundary conditions other than the fully prescribed and the solar-like are explained as follows: LL-CAR-UNI-OTH1: fixed solution of Low and Lou at ${x}_{\min ,\max },$ ${y}_{\min ,\max }$, z_max, inner z_min + zero-gradient at outer z_min. LL-CAR-UNI-OTH2: fixed solution of Low and Lou at ${x}_{\min ,\max },$ ${y}_{\min ,\max }$, z_max, inner z_min + linear at outer z_min. LL-CAR-UNI-OTH3: fixed potential field at ${x}_{\min ,\max },$ ${y}_{\min ,\max },$ z_max + fixed solution of Low and Lou at z_min. LL-CAR-UNI-OTH4: fixed potential field at ${x}_{\min ,\max },$ ${y}_{\min ,\max }$, z_max + fixed solution of Low and Lou at inner z_min + zero-gradient at outer z_min. LL-CAR-UNI-OTH5: fixed potential field at ${x}_{\min ,\max },$ ${y}_{\min ,\max },$ z_max + fixed solution of Low and Lou at inner z_min + linear at outer z_min. LL-CAR-UNI-OTH6: linear at ${x}_{\min ,\max },$ ${y}_{\min ,\max },$ z_max + zero-gradient at outer z_min + fixed solution of Low and Lou at inner z_min.

Download table as:  ASCIITypeset image

Second, we test the TVDLF scheme with one-step time-advancing, again using the fully prescribed boundary condition and the potential field as the initial condition (named as model LL-CAR-UNI-FUL-TVDLF). The evolution of the metrics $\langle | {f}_{i}| \rangle$ and ${\sigma }_{J}$ as shown in Figure 1(b) indicates that neither the divergence-free nor the force-free condition decreases to values as small as the CD4LF scheme, as shown in Figure 1(a). Besides, the metrics show an oscillatory behavior in the iteration process. The main difference in the CD4LF scheme compared to TVDLF is the raised order of central difference in the non-dissipative physical flux evaluation, where the CD4LF scheme uses cell-center values directly to avoid the nonlinear limited reconstruction process in the TVDLF scheme. Valori et al. (2007) also found that the fourth-order spatial discretization is essential to derive better divergence-free and force-free conditions than flux-limited second-order schemes. Any numerical treatment for NLFFF modeling should avoid the decoupling between even and odd grid layers generated by imposing boundary conditions only on the bottom boundary, while relying on outward extrapolation for all non-photospheric ones. In a sense, the diffusive part of the fluxes improves the coupling, but it also limits the level of force-freeness that can be attained. In this sense, the simple fourth-order central difference scheme is able to improve the propagation of the photospheric information to all layers, without adding terms that are incompatible with the force-free equations. Therefore, a high-order scheme is clearly a desired feature for the convergence of the magneto-frictional method. We note that the ${\sigma }_{J}$ curve shows a plateau in the later iteration steps with the TVDLF scheme (Figure 1(b)), but clearly the magnetic field does not converge to the NLFFF state as well as the CD4LF scheme.

Third, we tried different time-advancing schemes with the fully prescribed boundary condition, the potential field as the initial condition, and the CD4LF scheme. Two time-advancing schemes are adopted, namely the two-step (LL-CAR-UNI-FUL-TWO) and three-step (LL-CAR-UNI-FUL-THR) methods. Since the magnetic field is updated twice and three times in the two- and three-step methods, we correspondingly reduce the maximum iteration step to $3\times {10}^{4}$ steps (Figure 1(c)) and $2\times {10}^{4}$ steps (Figure 1(d)), respectively. The multi-step methods are applied only to the advection term in Equation (7). Therefore, the magnetic field divergence cleaner proceeds at the same pace as the one-step scheme, and the divergence-free metric $\langle | {f}_{i}| \rangle$ does not decrease to the same value as that in the one-step case, which is iterated with more steps. The force-free metric ${\sigma }_{J}$ in the three-step case decreases to a value similar to that in the one-step case (Figure 1(d)). We find that ${\sigma }_{J}$ in the two-step case has large oscillations in the late iteration steps (Figure 1(c)). As expected, the run times for the three time-advancing schemes are very similar to each other due to the selected maximum iteration steps. Unlike in a physical evolution process, multi-step time-advancing does not provide higher accuracy, and we therefore use the simple one-step method in the following tests and applications.

Fourthly, we continue the relaxation up to $2\times {10}^{5}$ steps using the CD4LF scheme, one-step time-advancing method, the fully prescribed boundary condition, and the potential field as the initial condition. This is also labeled as model LL-CAR-UNI-FUL as shown in Figure 1(e), but the relaxation step is larger than that in Figure 1(a) and Table 2, which is $6\times {10}^{4}$. Figure 1(e) shows that the force-free metric ${\sigma }_{J}$ decreases much more slowly in the later steps, to reach a "plateau." The magneto-frictional velocity in Equation (5) becomes smaller when a force-free state is reached. This is why the relaxation process slows down and the plateau is reached in the later steps. The metrics for model LL-CAR-UNI-FUL in Table 2 show that the magnetic field has already been relaxed to values that are one order of magnitude away from the metrics of the reference model LL-CAR-UNI. This may be hard to further improve upon, so we do not compute all test cases to $2\times {10}^{5}$ steps but find the chosen $6\times {10}^{4}$ steps to be a good compromise between execution time and accuracy reached. Quantities in Figure 1 are displayed on logarithmic scales, but we note that the force-free metric ${\sigma }_{J}$ has reached values at $6\times {10}^{4}$ steps that are not far away from the one finally attained within $2\times {10}^{5}$ steps. The divergence-free metric $\langle | {f}_{i}| \rangle$ even eventually increases as shown in Figure 1(e). This illustrates the difficulty of defining a clear stopping criterion; e.g., in some cases as in Figures 1(b) and (c), the solution is no longer converging. The computation time for the case with $6\times {10}^{4}$ steps is about 9.8 hr on a local desktop, with four Intel® Xeon® E5-2630 processors with 2.40 GHz central processing unit (CPU).

Finally, we simulate solar-like cases where only one layer of vector magnetic field on the bottom boundary is available. This still leaves quite some freedom, and we could loosen these boundary conditions step by step. These models are labeled as LL-CAR-UNI-OTH and LL-CAR-UNI-BOT depending on the chosen boundary conditions. Since we always use two ghost layers for each boundary, we always prescribe the known vector magnetic field to the inner ghost layer at the bottom boundary. As an intermediate step, the four side and one top boundaries are prescribed by their local solution of Low and Lou, but the outer ghost layer at the bottom boundary is handled by the following different numerical extrapolations. There, we compare two different numerical schemes, which are loosely referred to as zero-gradient and linear extrapolations. The zero-gradient extrapolation means using ${u}_{i}^{\prime }=0$ in its one-sided finite-difference expression, making ${u}_{i}=\tfrac{1}{25}(-3{u}_{i+4}+16{u}_{i+3}-36{u}_{i+2}+48{u}_{i+1})$ with the fourth-order accuracy and ${u}_{i}=\tfrac{1}{3}(-{u}_{i+2}+4{u}_{i+1})$ with the second-order accuracy. The linear extrapolation uses ${u}_{i}^{\prime\prime }=0$ from a one-sided second-order difference expression, making ${u}_{i}=-\tfrac{1}{2}(-{u}_{i+3}+4{u}_{i+2}-5{u}_{i+1})$. The tangential magnetic vector components on the outer bottom ghost layer are computed either by the fourth-order accuracy zero-gradient extrapolation or the second-order accuracy linear extrapolation, while the normal component is computed from ${\rm{\nabla }}\cdot {\boldsymbol{B}}=0$ in a second-order central difference manner, which is consistent with its use in handling divergence or in the evaluation of the metrics. We impose a boundary condition for the velocity, which is simply using an asymmetric prescription in both ghost cells for all the six boundaries. A linear extrapolation for the velocity boundary condition has also been tested, which yields very similar results to the asymmetric boundary condition. The models LL-CAR-UNI-OTH1 and LL-CAR-UNI-OTH2 in Table 2 show that the linear extrapolations provide slightly better divergence-free and force-free conditions than the zero-gradient extrapolation, but the mean vector error of the latter is better than that of the former in the entire computational volume.

Next, we further loosen the boundary conditions at the four side and one top boundaries, which are now fixed to the initial potential field values instead of the local values of Low and Lou. When both ghost layers of the bottom boundary are given by the local solution of Low and Lou, we have model LL-CAR-UNI-OTH3 in Table 2, which shows that the solution of Low and Lou cannot be recovered very well in the entire volume, but the solution in the inner volume is better. We also vary the outer layer of the bottom boundary as prescribed by the zero-gradient (LL-CAR-UNI-OTH4) and linear (LL-CAR-UNI-OTH5) extrapolations. The results are very similar to that of model LL-CAR-UNI-OTH3. The zero-gradient and linear extrapolation boundaries provide similar results for NLFFF extrapolation.

The fixed potential field at the four side and one top boundaries can also be loosened. These boundaries of the test cases LL-CAR-UNI-BOT and LL-CAR-UNI-OTH6 as listed in Table 2 are then also prescribed by the second-order accuracy zero-gradient and linear extrapolations, respectively. Model LL-CAR-UNI-BOT provides the best mean vector error among all the solar-like cases (the last four models in Table 2). Since this metric is the most sensitive one to describe how the original magnetic field is recovered, we select this setting for the boundaries in all following tests and applications. Figure 1(f) displays the evolution of the divergence-free and force-free metrics in the iteration process. Figures 2(e) and (f) compare the extrapolated magnetic field and the original solution of Low and Lou. It is found that higher extrapolated magnetic field lines, which are connected to the side or top boundaries, are not recovered by the extrapolation method. The lower-lying extrapolated magnetic field lines coincide with the original solution to a good degree. This result is understandable, since the information on the four side and one top boundaries is unknown and the potential field is far from the solution of Low and Lou in this case. The metrics listed in Table 2 for model LL-CAR-UNI-BOT in both the entire and inner volumes are comparable to the best results as listed in Tables 1 and 2 in Schrijver et al. (2006). Those metrics for model LL-CAR-UNI-BOT in the inner volume are also very similar to that derived by Valori et al. (2007) in their Table 1.

3.2. Titov–Démoulin Model in Cartesian Coordinates

The models of Titov & Démoulin (1999) are approximate semi-analytic solutions with different parameters for force-free magnetic field configurations with a twisted magnetic flux rope. The models are a superposition of three parts: a current channel in a torus with a minor radius a and major radius R, two magnetic charges with strength q on the symmetry axis of the torus, and a line current with magnitude of I₀ along the symmetry axis. The two magnetic charges are located at a distance L on each side of the torus plane. To avoid the singularity of the magnetic charges, the symmetry axis of the torus is placed at depth d in a selected photosphere plane that is perpendicular to the torus plane. The magnetic field above the photosphere plane is obtained by tuning parameters to reach force equilibrium of the torus current channel both inside and outside. The external force balance is between the torus hoop force and the magnetic tension force of the external (to the torus) potential field (Kliem & Török 2006). The electric current inside the torus is determined by this external force balance. The internal force balance of the torus is between magnetic tension and magnetic pressure gradients, and sets the internal magnetic field along the direction of the torus (the so-called toroidal component). Toroidal magnetic field is provided by the line electric current. The poloidal magnetic field is build up from the magnetic charges and the torus current.

Here, we adopt the same reference Titov–Démoulin models as in Valori et al. (2010). Three major modifications have been made to these models compared to the original models. The line current is replaced with two pairs of dipoles to simulate a more realistic coronal field. The current profile in the torus is changed from uniform to parabolic. The semi-analytic solutions are further relaxed to a numerical equilibrium with a separate MHD code (Török & Kliem 2003). This is needed since the solutions are derived by omitting contributions from terms of order a/R. The model parameters are listed in Table 1 of Valori et al. (2010). We test only the first three cases with or without a hyperbolic flux tube (HFT). As listed in Table 1 of Valori et al. (2010), a = 0.67 and R = 1.83 for all the three models that are adopted here. Therefore, a/R has finite values to induce non-force-free terms, which need to be relaxed with a numerical method. The reference models have a slightly larger field of view than what we use. In practice, we use them also to provide values for ghost cells, so our computational box is selected as $x\in [-2.70,2.70]$, $y\in [-4.50,4.50]$, and $z\in [-0.06,1.74]$. These ranges are resolved into $90\times 150\times 60$ grid points. The data are again scaled such that the maximum value of $| {B}_{z}|$ is 500 G, and we take ${\rm{\Delta }}x={\rm{\Delta }}y={\rm{\Delta }}z={10}^{8}\,{\rm{cm}}$.

By adjusting the free parameters, especially the distance L at which the two magnetic charges are located, three different Titov–Démoulin models—without an HFT, with a low HFT, or with a high HFT—are constructed. An HFT is a three-dimensional magnetic structure where two quasi-separatrix layers (QSLs) intersect (Titov et al. 2002). Due to the sharp changes in connectivity occurring there, QSLs and HFTs are preferred locations for the formation of current layers and for the occurrence of magnetic reconnection (e.g., Aulanier et al. 2005). Hence, it is important to verify that the extrapolation correctly reproduces such detailed magnetic field structures. The magnetic topology of these models is distinct from that of the potential field models derived from the normal component of the vector magnetic field on the bottom boundary. Therefore, magnetic topology must change in a magneto-frictional relaxation of the potential field to the original NLFFF. In such a sense, these tests are more challenging than the solution of Low and Lou.

We prepared two different boundary conditions for each of the three Titov–Démoulin models. One is the fully prescribed boundary condition where all the six boundaries (x_min, x_max, y_min, y_max, z_min, and z_max), each with two ghost layers, are fixed by the Titov–Démoulin models. The other is the solar-like boundary condition whereby only the inner ghost layer of the bottom boundary is fixed by the Titov–Démoulin models, and all the others are computed by the zero-gradient extrapolation. Figure 3 displays the evolution of the divergence-free and force-free metrics during the magneto-frictional process for all the six models (TD-NO-CAR-UNI-FUL, TD-NO-CAR-UNI-BOT, TD-LO-CAR-UNI-FUL, TD-LO-CAR-UNI-BOT, TD-HI-CAR-UNI-FUL, TD-HI-CAR-UNI-BOT). The three Titov–Démoulin models with the same boundary conditions show very similar evolutions. We could focus on the first model without an HFT. With the fully prescribed boundary as shown in Figure 3(a), the divergence-free metric first increases and then decreases to reach a value much smaller than its initial value. The general trend of the force-free metric always decreases until it reaches a plateau after about $4\times {10}^{4}$ steps. For the solar-like boundary condition as shown in Figure 3(b), the divergence-free metric increases first, then decreases, and finally increases to reach a plateau. The divergence-free metric does not decrease below the value of the initial state because we use an average metric, $\langle | {f}_{i}| \rangle$, but not the maximum value. The largest magnetic divergences in the initial state are on the bottom boundaries, where the analytic vector magnetic field is different from the initial potential field. The diffusive term to clean the divergence diffuses them into the computational box. Therefore, $\langle | {f}_{i}| \rangle$ for the final state might be larger than that for the initial one, but the maximum value of $| {f}_{i}|$ decreases as is shown in Figure 8 of Valori et al. (2007). The level of $\langle | {f}_{i}| \rangle$ for model TD-NO-CAR-UNI-BOT, which is $1.44\times {10}^{-4}$ computed in the entire volume, is acceptable. From Table 1 of Valori et al. (2013), this level of finite magnetic divergence would cause an error in magnetic energy less than a few per cent. For comparison, $\langle | {f}_{i}| \rangle$ is $1.03\times {10}^{-4}$ and $3.52\times {10}^{-5}$ for Valori et al. (2010) and Jiang & Feng (2016), respectively. The force-free metric always decreases, but to a value that is larger (and therefore worse) than the fully prescribed boundary condition. The computation time for all the cases is about 7.9 hr with four Intel® Xeon® E5-2630 processors with 2.40 GHz CPU.

Zoom In Zoom Out Reset image size

Figure 4 compares different extrapolated models with the Titov–Démoulin model without HFT. The potential field as shown by the white lines in Figures 4(a) and (b) obviously deviates from the original NLFFF. If all the boundaries are provided by the Titov–Démoulin model, the reconstructed NLFFF coincides with the original one both for higher field lines (Figure 4(c)) and for lower field lines (Figure 4(d)). When only the inner layer on the bottom boundary is provided, the NLFFF can still be reconstructed very close to the original model. Unlike the test of Low and Lou as shown in Figure 2(e), the difference between the reconstructed field and the original Titov–Démoulin model is very subtle even for higher magnetic field lines as shown in Figure 4(e). This is an encouraging result, indicating that if the magnetic field on the boundaries is close to the potential field derived from the bottom boundary, the magneto-frictional method could recover the NLFFF reliably. The Titov–Démoulin models are intended to simulate the magnetic field configurations in the solar atmosphere, where the magnetic field far from active regions might be close to potential.

Zoom In Zoom Out Reset image size

We list all the metrics for the three Titov–Démoulin models with two ways of providing the boundary conditions, and both in the entire computational domain and in an inner volume, where 22, 37, and 30 grid points are subtracted from ${x}_{\min ,\max }$, ${y}_{\min ,\max }$, and z_max, respectively. If a fully prescribed boundary condition is provided, all the three models can be reconstructed to a very reliable state, where the complement of the mean vector error ($1-{E}_{{\rm{m}}}$) reaches 0.99 for the best cases (TD-NO-CAR-UNI-FUL and TD-LO-CAR-UNI-FUL in the entire volume) and 0.96 for the worst case (TD-HI-CAR-UNI-FUL in the inner volume). With the solar-like boundary condition, the best $1-{E}_{{\rm{m}}}$ metric reaches 0.95 for the TD-NO-CAR-UNI-BOT and TD-LO-CAR-UNI-BOT cases, and the worst is 0.94 for TD-HI-CAR-UNI-BOT case, all in the inner volume. These metrics can be compared with previous studies. For example, Wiegelmann et al. (2006a), Valori et al. (2010), and Jiang & Feng (2016) have also tested different NLFFF extrapolation algorithms with the Titov–Démoulin models. The complement of the mean vector error ($1-{E}_{{\rm{m}}}$) in our tests is better than that derived by Wiegelmann et al. (2006a), comparable to that of Valori et al. (2010), and slightly worse than that of Jiang & Feng (2016). The force-free and divergence-free metrics are larger (and therefore worse) than those derived by Valori et al. (2010) and Jiang & Feng (2016), but the good match of field lines shows that their influence on the connection of the Titov–Démoulin equilibrium is limited.

3.3. Solution of Low and Lou in Spherical Coordinates

There are several reasons to demonstrate the capability of computing NLFFF models in the spherical coordinate system. First, the solar surface is approximately spherical in nature. The magnetic field observed on the photosphere, which is usually adopted as the boundary condition, is located on such a spherical surface, where it is natural to use spherical coordinates. Next, NLFFF modeling requires as large a field of view as possible to include more magnetic field connections, and to decrease the contributions of the electric current in active regions on the side and top boundaries. Guo et al. (2012) have shown that when the field of view is large, projection of the photospheric data to a plane causes large errors. The spherical nature of the solar surface cannot be neglected in such cases. Additionally, full-disk vector magnetic field data are available from the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012; Schou et al. 2012) on board the Solar Dynamics Observatory (SDO), which allows us to consider NLFFF modeling in large fields of view. Finally, there are many large-scale phenomena in the solar atmosphere, such as long filaments, coronal holes, trans-equatorial coronal loops, and so on, that require large fields of view, or even a full spherical modeling. Magnetic field modeling in spherical geometry is also required to connect NLFFF extrapolations with full dynamical models, where global-scale solar activities such as coronal mass ejections and solar winds need to be considered.

We show here that the spherical version of the magneto-frictional method in MPI-AMRVAC can be used in a wedge-shaped space with a range of $r\in [{r}_{\min },{r}_{\max }]$, $\theta \in [{\theta }_{\min },{\theta }_{\max }]$, and $\phi \in [{\phi }_{\min },{\phi }_{\max }]$. Using the solution of Low and Lou as a test, two factors prevent us from using the code with $\theta \in [0,\pi ]$. One is that the solution of Low and Lou has a singular line along its symmetry axis. The magnetic field is not divergence-free and force-free along this singular line. The other is that we did encounter numerical problems at the physical poles where $\theta =0$ or $\theta =\pi$. Although we can rely on the standard π-periodicity boundary condition at $\theta =0$ and $\theta =\pi$, we noticed that the Lorentz force becomes too large at the poles and errors spread into the computational volume. For the actual test, we prepare a different solution of Low and Lou compared to the Cartesian case: the same solution from the local coordinate space is now transformed by the rotational parameter ${\rm{\Phi }}=-\pi /20$ and translation parameter l_x = 0.02. This solution allows us to use a large range for θ and we scale the magnetic field such that the largest $| {B}_{r}|$ at r = 1 is 500 G. Since the PFSS model requires B_r on the full global surface as a boundary condition, this scaling is performed on the full data cube with $r\in [1,2.5]$, $\theta \in [0,\pi ]$, and $\phi \in [0,2\pi ]$. The computational box for the NLFFF, however, is selected as $r\in [1,2.5]$, $\theta \in [0.15\pi ,0.85\pi ]$, and $\phi \in [0,2\pi ]$, which are resolved into $100\times 140\times 200$ grid points. The length is scaled such that $r\in [1{R}_{\odot },2.5{R}_{\odot }]$, where ${R}_{\odot }=6.955\times {10}^{5}$ km.

Similar to previous cases, two different boundary conditions are prepared. In the fully prescribed boundary case (model LL-SPH-UNI-FUL), all the boundary conditions at r_min, r_max, ${\theta }_{\min }$, and ${\theta }_{\max }$ are provided by the solution of Low and Lou and the boundary conditions at ${\phi }_{\min }$ and ${\phi }_{\max }$ are periodic. In the solar-like case (model LL-SPH-UNI-BOT), only the inner ghost layer of the bottom boundary is provided by the solution of Low and Lou, and the others use the second-order zero-gradient prescription. Figures 5(a) and (b) display the evolution of divergence-free and force-free metrics in these two ways of providing boundary conditions. All the metrics decrease to a very small value, although the fully prescribed boundary case yields smaller, and therefore better, metrics. The computation time is about 3 hr with 40 Intel® Xeon® E5-2680v2 processors with 2.80 GHz CPU.

Zoom In Zoom Out Reset image size

Figure 6 shows the modeling results of the magneto-frictional method with the solution of Low and Lou in the spherical coordinate system. The initial condition is computed with the PFSS model as shown by the white field lines in Figures 6(a) and (b). They are obviously different from the original solution of Low and Lou, which is a NLFFF. With the fully prescribed boundary conditions, the initial potential field can be relaxed to the NLFFF with the magneto-frictional method in the whole computational volume (Figures 6(c) and (d)). But with the solar-like boundary conditions, the magnetic field at higher altitudes cannot be relaxed to the original one (Figure 6(e)), while lower magnetic field lines that do not connect to the side and top boundaries can be relaxed (Figure 6(f)).

Zoom In Zoom Out Reset image size

Table 3.  Metrics for the Titov–Démoulin Models and the Nonlinear Force-free Fields Derived with Different Boundary Conditions

Modela	$\langle | {f}_{i}| \rangle$	${\sigma }_{J}$	${\theta }_{J}$	Cvec	CCS	$1-{E}_{{\rm{n}}}$	$1-{E}_{{\rm{m}}}$
	$({10}^{-4})$	$({10}^{-2})$	(deg)
The metrics are computed in the entire volume:
TD-NO-CAR-UNI	0.31	0.75	0.43	⋯	⋯	⋯	⋯	1
TD-NO-CAR-UNI-FUL	0.50	4.17	2.39	1.00	1.00	0.98	0.99	0.99
TD-NO-CAR-UNI-BOT	1.44	12.5	7.18	1.00	1.00	0.93	0.90	1.01
TD-LO-CAR-UNI	0.32	0.85	0.48	⋯	⋯	⋯	⋯	1
TD-LO-CAR-UNI-FUL	0.50	3.31	1.90	1.00	1.00	0.98	0.99	0.99
TD-LO-CAR-UNI-BOT	1.43	12.3	7.05	1.00	1.00	0.93	0.90	1.01
TD-HI-CAR-UNI	0.32	1.01	0.58	⋯	⋯	⋯	⋯	1
TD-HI-CAR-UNI-FUL	0.60	4.00	2.29	1.00	1.00	0.98	0.98	1.00
TD-HI-CAR-UNI-BOT	1.47	11.6	6.65	1.00	1.00	0.93	0.90	1.02
The metrics are computed in the inner volume:
TD-NO-CAR-UNI	0.62	0.44	0.25	⋯	⋯	⋯	⋯	1
TD-NO-CAR-UNI-FUL	0.65	3.61	2.07	1.00	1.00	0.98	0.98	0.98
TD-NO-CAR-UNI-BOT	3.06	9.32	5.35	1.00	1.00	0.95	0.95	1.00
TD-LO-CAR-UNI	0.63	0.52	0.30	⋯	⋯	⋯	⋯	1
TD-LO-CAR-UNI-FUL	0.64	2.94	1.68	1.00	1.00	0.98	0.98	0.99
TD-LO-CAR-UNI-BOT	3.12	9.35	5.37	1.00	1.00	0.95	0.95	1.01
TD-HI-CAR-UNI	0.78	0.65	0.37	⋯	⋯	⋯	⋯	1
TD-HI-CAR-UNI-FUL	0.73	3.36	1.92	1.00	1.00	0.96	0.96	0.99
TD-HI-CAR-UNI-BOT	3.22	9.00	5.17	1.00	1.00	0.94	0.94	1.01

Note.

^aRefer to Table 1 for the explanation of the models.

Download table as:  ASCIITypeset image

The metrics listed in Table 4 for the test cases of LL-SPH-UNI-FUL and LL-SPH-UNI-BOT confirm the previous findings. If we focus on the entire volume, only the fully prescribed boundary condition could yield good force-free and vector comparison metrics. Compared to the LL-SPH-UNI-BOT case in the entire volume, the force-free and vector comparison metrics of the solar-like case are much better in the inner volume, which is defined as the inner region excluding 50 grid points in the r_max boundary and 35 grid points in the ${\theta }_{\min ,\max }$ boundaries. These results can also be compared with previous tests with the solution of Low and Lou in the spherical geometry. Only the solar-like case in an inner volume is considered. We find that the complements of the normal vector error ($1-{E}_{{\rm{n}}}$) and mean vector error ($1-{E}_{{\rm{m}}}$) in this paper are better than those derived by Wiegelmann (2007) and Guo et al. (2012). And $1-{E}_{{\rm{n}}}$ in this paper is better than that in Tadesse et al. (2009), $1-{E}_{{\rm{m}}}$ is comparable in the two studies. We note that the models of Low and Lou and the NLFFF extrapolation algorithms from these studies differ from our current setting; still, these comparisons provide an idea of the performance.

Table 4.  Metrics for the Solution of Low and Lou and the Nonlinear Force-free Fields Derived with Different Boundary Conditions in the Spherical Geometry with or without AMR Grid, and in the Cartesian Geometry with AMR Grid

Modela	$\langle | {f}_{i}| \rangle$	${\sigma }_{J}$	${\theta }_{J}$	Cvec	CCS	$1-{E}_{{\rm{n}}}$	$1-{E}_{{\rm{m}}}$
	$({10}^{-4})$	$({10}^{-2})$	(°)
The metrics are computed in the entire volume:
LL-SPH-UNI	0.01	0.06	0.04	⋯	⋯	⋯	⋯	1
LL-SPH-UNI-FUL	0.70	0.82	0.47	1.00	1.00	0.95	0.94	0.95
LL-SPH-UNI-BOT	2.06	12.70	7.29	0.98	0.91	0.83	0.70	0.95
LL-CAR-AMR	0.02	0.22	0.12	⋯	⋯	⋯	⋯	1
LL-CAR-AMR-FUL	0.77	3.35	1.92	1.00	1.00	0.98	0.96	0.98
LL-CAR-AMR-BOT	1.01	5.93	3.40	1.00	0.93	0.93	0.78	0.85
LL-SPH-AMR	0.01	0.04	0.02	⋯	⋯	⋯	⋯	1
LL-SPH-AMR-FUL	0.58	2.44	1.40	1.00	1.00	0.98	0.97	0.97
LL-SPH-AMR-BOT	1.07	2.58	1.48	0.99	0.99	0.93	0.91	0.97
The metrics are computed in the inner volume:
LL-SPH-UNI	0.01	0.08	0.04	⋯	⋯	⋯	⋯	1
LL-SPH-UNI-FUL	0.28	0.33	0.19	1.00	1.00	0.95	0.95	0.95
LL-SPH-UNI-BOT	2.94	5.80	3.33	1.00	0.99	0.94	0.92	1.00
LL-CAR-AMR	0.04	0.21	0.12	⋯	⋯	⋯	⋯	1
LL-CAR-AMR-FUL	0.37	1.65	0.94	1.00	1.00	0.99	0.98	0.98
LL-CAR-AMR-BOT	1.17	4.04	2.32	1.00	1.00	0.97	0.95	0.95
LL-SPH-AMR	0.01	0.05	0.03	⋯	⋯	⋯	⋯	1
LL-SPH-AMR-FUL	0.28	0.98	0.56	1.00	1.00	0.99	0.99	0.98
LL-SPH-AMR-BOT	1.11	2.12	1.21	1.00	1.00	0.98	0.97	0.99

Note.

^aRefer to Table 1 for the explanation of the models.

Download table as:  ASCIITypeset image

3.4. Test of AMR

The AMR technique allows us to model a physical space with as large a field of view and as high a spatial resolution as possible at the same time. One could exploit high AMR levels only in the regions of interest, such as the lower solar atmosphere with high magnetic field strength. Other places are resolved with coarser grid points. Hence, the goal of using high spatial resolution in a big field of view is achieved with a reasonable data volume. The implementation of AMR in MPI-AMRVAC has been described in detail in Keppens et al. (2012). Here, we test whether the magneto-frictional method works with the AMR technique. The solution of Low and Lou is adopted as the reference model. The initial condition is the potential field derived from the normal component of the magnetic field on the bottom boundary. Four test cases are detailed as below with fully prescribed or solar-like boundary conditions and in Cartesian or spherical coordinates (models LL-CAR-AMR-FUL, LL-CAR-AMR-BOT, LL-SPH-AMR-FUL, and LL-SPH-AMR-BOT). Figures 5(c)–(f) display the evolution of the divergence-free and force-free metrics for these four test cases. The curves indicate that the initial conditions converge to NLFFF models with low magnetic field divergence and Lorentz force.

For the Cartesian test case, we use a similar model of Low and Lou to that in the uniform case. Here, three AMR levels are adopted. The first level is resolved into $50\times 50\times 50$ grid points. Each higher level has a spatial resolution twice that of the lower one. The AMR level is block-based, and the size of each block in this case is $10\times 10\times 10$ grid points. We let the local level be determined by a Löhner-type estimator (Keppens et al. 2012), which evaluates a weighted second-order derivative from instantaneous values. The mesh is then hierarchically organized automatically, but we only use this functionality at the beginning and keep the hierarchical mesh unchanged during the magneto-frictional iteration. As an example, Figures 7(a) and (c) display the meshes for the fully prescribed (LL-CAR-AMR-FUL) and solar-like (LL-CAR-AMR-BOT) boundary cases, respectively. The meshes are refined to the highest level at most places of the six boundaries for the fully prescribed case; while for solar-like cases, we consistently find the mesh blocks of the highest level near the bottom boundary. The reason for this difference is that magnetic field in the ghost layers and in the computational domain in the fully prescribed case jumps from the solution of Low and Lou to the potential field at boundaries, which is detected by the error estimator that determines the local refinement level. On the other hand, in the solar-like cases, the magnetic field in the lateral and top ghost layers and in the computational domain is potential, and the magnetic field changes smoothly across these boundaries. So, the meshes at these boundaries are not refined. The three AMR levels for the fully prescribed boundary case have 27, 298, and 3888 blocks from the lowest level to the highest one. Since each block has 10³ cells, the data volume is equivalent to a uniform mesh with about 162³ cells (and has an effective resolution of 200³). The computation time is 4.7 hr with 40 Intel® Xeon® E5-2680v2 processors with 2.80 GHz CPU. Similarly, there are 56, 387, and 1320 blocks on each AMR level for the solar-like boundary case, which is equivalent to a uniform mesh with about 121³ cells. The computation time is about 1.9 hr with the same resource.

Zoom In Zoom Out Reset image size

The relaxed magnetic field is very close to the semi-analytic solution of Low and Lou for both boundary conditions as shown in Figure 7. The relaxed magnetic field in the horizontal direction with fully prescribed boundaries is closer to the original solution than the case with solar-like conditions as shown in Figures 7(b) and (d). We also provide a close-up side view of the magnetic field lines for both cases in Figure 8. Unlike the top view, the magnetic field with solar-like boundaries is closer to the original solution than the case with fully prescribed ones in this vertical direction. The divergence-free, force-free, and vector comparison metrics for these two cases are listed as models LL-CAR-AMR-FUL and LL-CAR-AMR-BOT in Table 4. The metrics for LL-CAR-AMR-FUL in the entire region are slightly worse than for model LL-CAR-UNI-FUL as listed in Table 2, but they are very similar in the inner domain. This suggests that AMR with a base resolution of 50³ cells and three levels could perform as well as the uniform-mesh case with a resolution of 100³ cells, at least in the inner region. Although the total data volume for the AMR (162³ cells) is still larger than in the uniform-mesh case, the advantage of AMR is that it allocates finer meshes only where a higher spatial resolution is needed, while the total physical space can be enlarged. Most of the metrics with solar-like boundaries and AMR (model LL-CAR-AMR-BOT as listed in Table 4) are better than or similar to those of the uniform case (model LL-CAR-UNI-BOT as listed in Table 2), except for the magnetic energy in the entire volume.

Zoom In Zoom Out Reset image size

Finally, we test the magneto-frictional method with AMR and in spherical geometry. The solution of Low and Lou similar to the uniform mesh is also adopted for this purpose. We select a different field of view with $r\in [1,2.5]$, $\theta \in [0.125\pi ,0.875\pi ]$, and $\phi \in [0,2\pi ]$, which is resolved into $32\times 48\times 64$ cells at the base AMR level. Three AMR levels are used here, and each block contains 8³ cells. The initial condition is provided by the PFSS model. Two different boundary conditions are considered. In the fully prescribed case (LL-SPH-AMR-FUL), all the ghost layers for ${r}_{\min ,\max }$ and ${\theta }_{\min ,\max }$ are provided by the solution of Low and Lou. In the solar-like case (LL-SPH-AMR-BOT), the innermost ghost layer for the r_min boundary is provided by the solution of Low and Lou, and those for r_max and ${\theta }_{\min ,\max }$ are provided by the second-order zero-gradient extrapolation. Here, we must take care with the different sizes of the meshes. The vector magnetic field used as the bottom boundary condition is always fixed at the highest AMR level in the innermost ghost layer. The boundary conditions for ${\phi }_{\min ,\max }$ boundaries for both cases are periodic. A slice of the mesh for each case is displayed in Figure 9. The three AMR levels consist of 48, 556, and 4768 blocks for the fully prescribed boundary case, and 128, 300, 1696 blocks for the solar-like boundary case. Since each block consists of 8³ cells, the two cases have data volumes equivalent to 140³ and 103³ cells, respectively. The computation times are about 3.9 and 1.6 hr for the two cases, respectively, with 40 Intel® Xeon® E5-2680v2 processors with 2.80 GHz CPU.

Zoom In Zoom Out Reset image size

The magneto-frictional method relaxes the initial PFSS model to an NLFFF model for both boundary conditions. If all the boundaries are known, the relaxed magnetic field coincides with the original one, not only for the lower altitudes (Figure 9(b)) but also for the higher ones (Figure 9(a)). If only the vector magnetic field on the bottom is known, only the lower magnetic field can relax to the original solution (Figure 9(d)). Higher magnetic field lines still deviate from the original ones, since the information on the side and top boundaries is provided by a field inconsistent with the solution of Low and Lou. The divergence-free, force-free, and vector comparison metrics listed for models LL-SPH-AMR-FUL and LL-SPH-AMR-BOT in Table 4 confirm these findings. We could also compare the results with the uniform-mesh case (models LL-SPH-UNI-FUL and LL-SPH-UNI-BOT in Table 4). In both the entire and inner volumes and for the fully prescribed and solar-like boundary cases, AMR grids provide similar or slightly better metrics compared to the uniform-mesh case, except for the force-free metric in the fully prescribed boundary case.

We have implemented the magneto-frictional module in MPI-AMRVAC to reconstruct NLFFF models with vector magnetic field provided on the boundaries. The distinctive feature of the present code is its versatile abilities combined with AMR for NLFFF modeling. It works in both Cartesian and spherical coordinates, and can exploit the parallelism of MPI-AMRVAC in both uniform (domain-decomposed) and AMR (block-based) grids. These features allow us to model the magnetic field in the corona in more realistic conditions. On the one hand, from both the practical needs and the NLFFF modeling itself, one needs to consider as large a field of view of the vector magnetic field boundary as possible. Then, we have to consider the curvature of the photosphere by employing spherical coordinates. On the other hand, much useful information about the magnetic field is in small-scale structures, such as elementary magnetic flux tubes and electric current channels. We need high spatial resolution to resolve such localized features or strong magnetic structures such as sunspots. To achieve these two goals—a large field of view and high spatial resolution—at the same time it becomes natural to consider AMR grids to reduce the data volume and the computation time.

As the code is parallelized with MPI, we can reduce the computation time by employing more resources. We find that the magneto-frictional module shows good scaling properties when the number of processors is increased. The tests in Section 3 can be sped up by a factor of 8.8–9.8 when the number of processors is increased from 4 to 40. The case reported using AMR in spherical coordinates with solar-like boundary conditions is 6.7 times faster when there are 10 times as many processors.

We test the magneto-frictional module in MPI-AMRVAC with the solution of Low & Lou (1990) and the models of Titov & Démoulin (1999). It is found that the CD4LF scheme and one-step time-advancing perform better than the TVDLF scheme and multi-step time-advancing. Thanks to the modular structure of MPI-AMRVAC, alternative spatial differentiation and time-advancing schemes are available besides the tested schemes. So, it is easy to implement extensions that allow for future developments of the code. When all the boundary conditions are known, the magneto-frictional method could reconstruct the NLFFF from a potential field in the whole computational box. If only the bottom boundary is provided, only the magnetic field in an inner region could be relaxed to NLFFF. We have listed all the divergence-free, force-free, and vector comparison metrics for the test cases with different boundaries (fully prescribed or solar-like), models (Low and Lou or Titov–Démoulin), coordinates (Cartesian or spherical), and grids (uniform or AMR). All of the tests can be compared with previous studies, and our results compare favorably. The magneto-frictional module in MPI-AMRVAC performs as well as previous studies in the cases of uniform grids in both Cartesian and spherical coordinates. Tests with AMR have been demonstrated here, which pave the way for future applications to cases with boundaries constrained by observational data (documented in an accompanying paper). We find that even with AMR, the initial potential field can be made to converge to the solution of Low and Lou in Cartesian and spherical coordinates, with fully prescribed and solar-like boundary conditions. And we note that even solar-like boundary conditions could deliver good results in an inner computational volume. The convergence in the outer volume may require further improvements and information, but the global topology is already reproduced.

There are four important free parameters that would affect the performance of the magneto-frictional method, namely, c_c and c_y to control the magneto-frictional velocity, c_d to control the diffusion of the magnetic divergence, and ${\epsilon }^{\mathrm{LF}}$ to control the dissipation term. They influence the divergence-free, force-free, vector comparison metrics, and the convergence speed of the code. By some trial-and-error experiments, we adopt c_c = 0.5, c_y = 0.2, c_d = 0.1, and ${\epsilon }^{\mathrm{LF}}=1$ for all the test cases. We have also made a series of tests for the three Titov–Démoulin models by decreasing the parameter ${\epsilon }^{\mathrm{LF}}$ during the iteration process. Specifically, we decrease ${\epsilon }^{\mathrm{LF}}$ from 1 by multiplying by a factor of 0.9998 at each iteration step. It turns out that the divergence-free and force-free metrics can be improved but the vector correlation metrics become worse than in the present cases as listed in Table 3.

The magneto-frictional module developed in MPI-AMRVAC enables us to apply it to observations simultaneously with high spatial resolution and large field of view, such as those data provided by SDO/HMI. We demonstrate such applications and discuss issues specifically related to observations in the second paper of this series. The magneto-frictional method also builds initial conditions of the magnetic field for MHD simulations. It is a big advantage to prepare the NLFFF in the same grid as that used for the dynamic evolution in time. Such a procedure avoids the artificial magnetic divergence and Lorentz force caused by the interpolation in grid conversions.

Y.G. is supported by NSFC (11203014 and 11533005), NKBRSF (2011CB811402 and 2014CB744203), and the mobility grant from the Belgian Federal Science Policy Office (BELSPO). This research was supported by the Interuniversity Attraction Poles Programme (initiated by the Belgian Science Policy Office, IAP P7/08 CHARM) and by the KU Leuven GOA/2015-014. Computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center) funded by the Hercules foundation and the Flemish government, department EWI. C.X. acknowledges support from FWO. G.V. acknowledges the support of the Leverhulme Trust Research Project Grant 2014-051. Y.G. and G.V. thank the ISSI international team on magnetic helicity estimations in models and observations of the solar magnetic field.