Algorithms of the Potential Field Calculation in a Three Dimensional Box

TL;DR

This paper compares three methods for calculating the magnetic potential field inside a 3D box, assessing their accuracy and efficiency using artificial models with known solutions.

Contribution

It provides a detailed analysis and comparison of three potential field calculation methods in a 3D box, highlighting their performance and accuracy.

Findings

All three methods successfully compute the potential field.

Performance varies depending on the method used.

Accuracy is validated against known solutions.

Abstract

Calculation of the potential field inside a three-dimensional box with the normal magnetic field component given on all boundaries is needed for estimation of important quantities related to the magnetic field such as free energy and relative helicity. In this work we present an analysis of three methods for calculating potential field inside a three-dimensional box. The accuracy and performance of the methods are tested on artificial models with a priori known solutions.

Tables

Table 1. Table 1. : Model M charges subscript 𝑀 charges M_{\mathrm{charges}} : volume metrics of 𝐁 calc subscript 𝐁 calc {\mathbf{B}}_{\rm calc} and the scalar potential ϕ italic-ϕ \phi computed for the entire volume V 𝑉 V .

			${ϵ}_{\max }$
Case I	$9.734\times {10}^{-5}$	$7.850\times {10}^{-5}$	$8.244\times {10}^{-4}$	$9.734\times {10}^{-5}$	$9.745\times {10}^{-5}$
Case II	$1.111\times {10}^{-3}$	$1.091\times {10}^{-3}$	$6.401\times {10}^{-3}$	$1.038\times {10}^{-3}$	$2.925\times {10}^{-3}$
Case III	$8.197\times {10}^{-3}$	$6.935\times {10}^{-3}$	$5.524\times {10}^{-2}$	$5.343\times {10}^{-3}$	$1.084\times {10}^{-2}$

Table 2. Table 2. : Model M charges subscript 𝑀 charges M_{\rm charges} : surface metrics for 𝐁 calc subscript 𝐁 calc {\mathbf{B}}_{\rm calc} and the scalar potential ϕ italic-ϕ \phi computed for the whole boundary S 𝑆 S .

			${ϵ}_{\max }$
Case I	$1.043\times {10}^{-3}$	$4.679\times {10}^{-4}$	$2.594\times {10}^{-2}$	$1.278\times {10}^{-3}$	$1.419\times {10}^{-4}$
Case II	$1.999\times {10}^{-3}$	$1.487\times {10}^{-3}$	$1.918\times {10}^{-2}$	$2.331\times {10}^{-3}$	$4.184\times {10}^{-3}$
Case III	$1.293\times {10}^{-2}$	$9.985\times {10}^{-3}$	$1.847\times {10}^{-1}$	$1.311\times {10}^{-2}$	$1.445\times {10}^{-2}$

Table 3. Table 3. : Model M charges subscript 𝑀 charges M_{\rm charges} : volume metrics for ∇ × 𝐀 calc ∇ subscript 𝐀 calc {\bf\nabla}\times{\mathbf{A}}_{\rm calc} and the scalar potential ϕ italic-ϕ \phi computed for the entire volume V 𝑉 V .

			${ϵ}_{\max }$
Case I	$1.528\times {10}^{-3}$	$1.050\times {10}^{-3}$	$1.680\times {10}^{-2}$	$3.401\times {10}^{-3}$
Case II	$2.001\times {10}^{-3}$	$1.702\times {10}^{-3}$	$2.244\times {10}^{-2}$	$3.558\times {10}^{-3}$
Case III	$8.460\times {10}^{-3}$	$6.525\times {10}^{-3}$	$8.298\times {10}^{-2}$	$5.726\times {10}^{-3}$

Table 4. Table 4. : Model M charges subscript 𝑀 charges M_{\rm charges} : surface metrics for ∇ × 𝐀 calc ∇ subscript 𝐀 calc {\bf\nabla}\times{\mathbf{A}}_{\rm calc} and the scalar potential ϕ italic-ϕ \phi computed for the boundary surface S 𝑆 S .

			${ϵ}_{\max }$
Case I	$2.616\times {10}^{-3}$	$1.600\times {10}^{-3}$	$2.822\times {10}^{-2}$	$8.622\times {10}^{-3}$
Case II	$3.458\times {10}^{-3}$	$2.322\times {10}^{-3}$	$4.253\times {10}^{-2}$	$8.909\times {10}^{-3}$
Case III	$1.875\times {10}^{-2}$	$1.426\times {10}^{-2}$	$1.756\times {10}^{-1}$	$1.932\times {10}^{-2}$

Table 5. Table 5. : Model M charges subscript 𝑀 charges M_{\rm charges} : magnetic energy calculation results.

	$E$	$E_{\varphi }$	$E_{virial}$
Model	$2.9780\times {10}^5$	$2.9710\times {10}^5$	$3.0072\times {10}^5$
Case I	$2.9776\times {10}^5$	$2.9705\times {10}^5$	$3.0069\times {10}^5$
Case II	$2.9727\times {10}^5$	$2.9654\times {10}^5$	$3.0020\times {10}^5$
Case III	$2.9684\times {10}^5$	$2.9673\times {10}^5$	$3.0334\times {10}^5$

Table 6. Table 6. : Model M ARs subscript 𝑀 ARs M_{\rm ARs} : volume metrics for 𝐁 calc subscript 𝐁 calc {\mathbf{B}}_{\rm calc} and scalar potential ϕ italic-ϕ \phi computed for the entire volume V 𝑉 V

			${ϵ}_{\max }$
Case I	$3.758\times {10}^{-4}$	$2.573\times {10}^{-4}$	$4.955\times {10}^{-1}$	$3.163\times {10}^{-4}$	$1.219\times {10}^{-3}$
Case II	$3.868\times {10}^{-4}$	$2.473\times {10}^{-4}$	$3.680\times {10}^{-1}$	$2.170\times {10}^{-4}$	$9.405\times {10}^{-4}$
Case III	$2.154\times {10}^{-2}$	$1.552\times {10}^{-2}$	$1.809\times {10}^0$	$6.932\times {10}^{-3}$	$6.209\times {10}^{-2}$

Table 7. Table 7. : Model M AR subscript 𝑀 AR M_{\rm AR} : surface metrics for 𝐁 calc subscript 𝐁 calc {\mathbf{B}}_{\rm calc} and the scalar potential ϕ italic-ϕ \phi computed for the boundary surface S 𝑆 S .

			${ϵ}_{\max }$
Case I	$2.475\times {10}^{-2}$	$1.404\times {10}^{-3}$	$8.130\times {10}^0$	$4.523\times {10}^{-2}$	$2.666\times {10}^{-3}$
Case II	$5.935\times {10}^{-2}$	$1.434\times {10}^{-3}$	$1.358\times {10}^0$	$4.002\times {10}^{-3}$	$2.115\times {10}^{-3}$
Case III	$7.682\times {10}^{-2}$	$4.792\times {10}^{-2}$	$1.308\times {10}^1$	$1.078\times {10}^{-1}$	$1.002\times {10}^{-1}$

Table 8. Table 8. : Model M AR subscript 𝑀 AR M_{\rm AR} : volume metrics for ∇ × 𝐀 calc ∇ subscript 𝐀 calc {\bf\nabla}\times{\mathbf{A}}_{\rm calc} and the scalar potential ϕ italic-ϕ \phi computed for the entire volume V 𝑉 V .

			${ϵ}_{\max }$
Case I	$5.776\times {10}^{-4}$	$2.816\times {10}^{-4}$	$1.452\times {10}^0$	$1.482\times {10}^{-3}$
Case II	$6.602\times {10}^{-4}$	$2.716\times {10}^{-4}$	$2.052\times {10}^0$	$1.643\times {10}^{-3}$
Case III	$2.152\times {10}^{-2}$	$1.550\times {10}^{-2}$	$1.494\times {10}^0$	$6.984\times {10}^{-3}$

Table 9. Table 9. : Model M AR subscript 𝑀 AR M_{\rm AR} : surface metrics for ∇ × 𝐀 calc ∇ subscript 𝐀 calc {\bf\nabla}\times{\mathbf{A}}_{\rm calc} and the scalar potential ϕ italic-ϕ \phi computed for the boundary surface S 𝑆 S .

			${ϵ}_{\max }$
Case I	$4.304\times {10}^{-2}$	$6.913\times {10}^{-4}$	$1.326\times {10}^1$	$9.649\times {10}^{-2}$
Case II	$3.555\times {10}^{-2}$	$8.725\times {10}^{-4}$	$1.017\times {10}^1$	$7.979\times {10}^{-2}$
Case III	$8.471\times {10}^{-2}$	$4.906\times {10}^{-2}$	$1.549\times {10}^1$	$1.276\times {10}^{-1}$

Table 10. Table 10. : Model M AR subscript 𝑀 AR M_{\rm AR} : magnetic energy calculation results.

	$E$	$E_{\varphi }$	$E_{\mathrm{virial}}$
Model	$7.7832\times {10}^{32}$	$7.7766\times {10}^{32}$	$7.7767\times {10}^{32}$
Case I	$7.7828\times {10}^{32}$	$7.7763\times {10}^{32}$	$7.7764\times {10}^{32}$
Case II	$7.7831\times {10}^{32}$	$7.7765\times {10}^{32}$	$7.7766\times {10}^{32}$
Case III	$7.7718\times {10}^{32}$	$7.7692\times {10}^{32}$	$7.7262\times {10}^{32}$

Equations

$$\Delta \varphi =0,\left(-\nabla \varphi \cdot \mathbf{n}=g\right){|}_S,$$

$$\left\{\begin{array}[]{ll}\ {\bf\nabla}\times{{\mathbf{B}}}={{\mathbf{0}}}\\ \ {\bf\nabla}\cdot{{\mathbf{B}}}=0\\ \ {{\mathbf{B}}}=-{\bf\nabla}\phi\\ \end{array}\right|_{\overline{V}},\ \ \ {\left.({\mathbf{B}}\cdot{\mathbf{n}}=g)\right|}_{S},$$

$$\left\{\begin{array}[]{ll}\ \Delta{\mathbf{A}}={\bf 0}\\ \ {\bf\nabla}\cdot{\mathbf{A}}=0\\ \ {\mathbf{B}}={\bf\nabla}\times{\mathbf{A}}\\ \ {\bf\nabla}\times{\mathbf{A}}=-{\bf\nabla}\phi\\ \end{array}\right|_{\overline{V}},\ \ \begin{array}[]{ll}\ \ \ \ {\left.({\bf\nabla}\times{\mathbf{A}}\cdot{\mathbf{n}}=g)\right|}_{S}\\ \mathrm{Option}\ \ {\left.({\mathbf{A}}\cdot{\mathbf{n}}=0)\right|}_{S}\\ \end{array}.$$

$$\underset{S}{\int }g\mathrm{d}s=0.$$

$$\int\limits_{S}\widetilde{g}_{i}\mathrm{d}s=0\ \ (i=1,2,...,6):\ \left\{\begin{array}[]{ll}\widetilde{g}_{i}({{\mathbf{r}}})\neq 0;{{\mathbf{r}}}\in S_{i}\\ \widetilde{g}_{i}({{\mathbf{r}}})=0;{{\mathbf{r}}}\notin S_{i}\end{array}\right..$$

$$g=g_a+\sum _{i=1}^6{\tilde{g}}_i,$$

$$\underset{S_i}{\int }{g}_a\mathrm{d}s=\underset{S_i}{\int }g\mathrm{d}s\equiv G_i.$$

$$\begin{array}[]{ll}\ \phi=\phi_{a}+\sum\limits_{i=1}^{6}\phi_{i},\\ \ {\mathbf{B}}={\mathbf{B}}_{a}+\sum\limits_{i=1}^{6}{\mathbf{B}}_{i},\\ \ {\mathbf{A}}={\mathbf{A}}_{a}+\sum\limits_{i=1}^{6}{\mathbf{A}}_{i},\end{array}$$

$$\begin{array}[]{ll}\ \phi_{a}=\sum\limits_{j=1}^{5}m_{j}\phi_{a}^{j},\\ \ {\mathbf{B}}_{a}=\sum\limits_{j=1}^{5}m_{j}{\mathbf{B}}_{a}^{j},\\ \ {\mathbf{A}}_{a}=\sum\limits_{j=1}^{5}m_{j}{\mathbf{A}}_{a}^{j}.\end{array}$$

$$\begin{array}[]{ll}\phi_{a}^{1}=-x,\\ \phi_{a}^{2}=-y,\\ \phi_{a}^{3}=-z,\\ \phi_{a}^{4}=-x^{2}+z^{2},\\ \phi_{a}^{5}=-y^{2}+z^{2};\end{array}$$

$$\begin{array}[]{ll}{\mathbf{B}}_{a}^{1}=(1,0,0),\\ {\mathbf{B}}_{a}^{2}=(0,1,0),\\ {\mathbf{B}}_{a}^{3}=(0,0,1),\\ {\mathbf{B}}_{a}^{4}=(2x,0,-2z),\\ {\mathbf{B}}_{a}^{5}=(0,2y,-2z);\end{array}$$

$$\begin{array}[]{ll}{\mathbf{A}}_{a}^{1}=(0,0,y),\\ {\mathbf{A}}_{a}^{2}=(z,0,0),\\ {\mathbf{A}}_{a}^{3}=(0,x,0),\\ {\mathbf{A}}_{a}^{4}=(0,-2xz,0),\\ {\mathbf{A}}_{a}^{5}=(2yz,0,0).\end{array}$$

$$\underset{S_i}{\int }{g}_a^j\mathrm{d}s=\underset{S_i}{\int }{\mathbf{B}}_a^j\cdot \mathbf{n}\mathrm{d}s=G_{ai}^j.$$

$$G_{a1}^3=\underset{0}{\overset{L_x}{\int }}\underset{0}{\overset{L_y}{\int }}{B}_{az}^3\times (-1)\mathrm{d}x\mathrm{d}y=-\underset{0}{\overset{L_x}{\int }}\underset{0}{\overset{L_y}{\int }}\mathrm{d}x\mathrm{d}y=L_x{L}_y.$$

$$G_{ai}^j{m}_j=G_i(i=1,...,5;j=1,...,5).$$

$${\tilde{g}}_i=g_i-g_{ai}(i=1,...,6).$$

$$\varphi =\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }{p}_{mn}{\varphi }_{mn}=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }{p}_{mn}{C}_x^m{C}_y^n{Z}^{mn},$$

$$\begin{array}[]{ll}B_{x}=\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}(B_{x})_{mn}=-\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}C_{x}^{{}^{\prime}m}C_{y}^{n}Z^{mn},\\ B_{y}=\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}(B_{y})_{mn}=-\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}C_{x}^{m}C_{y}^{{}^{\prime}n}Z^{mn},\\ B_{z}=\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}(B_{z})_{mn}=-\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}C_{x}^{m}C_{y}^{n}Z^{{}^{\prime}mn};\end{array}$$

$$\begin{array}[]{ll}A_{x}=\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}(A_{x})_{mn}=-\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}C_{x}^{m}C_{y}^{{}^{\prime}n}Z^{{}^{\prime}mn}\frac{\pi n}{q^{2}L_{y}},\\ A_{y}=\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}(A_{y})_{mn}=\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}p_{mn}C_{x}^{{}^{\prime}m}C_{y}^{n}Z^{{}^{\prime}mn}\frac{\pi m}{q^{2}L_{x}},\\ A_{z}=(A_{z})_{mn}=0,\end{array}$$

$$q_{mn}=\sqrt{{\left(\frac{\pi m}{L_x}\right)}^2+{\left(\frac{\pi n}{L_y}\right)}^2},$$

$$\begin{array}[]{ll}C_{x}^{m}=\cos\left(\frac{\pi m}{L_{x}}x\right),\ \ C_{x}^{{}^{\prime}m}=-\frac{\pi m}{L_{x}}\sin\left(\frac{\pi m}{L_{x}}x\right),\\ C_{y}^{n}=\cos\left(\frac{\pi n}{L_{y}}y\right),\ \ C_{y}^{{}^{\prime}n}=-\frac{\pi n}{L_{y}}\sin\left(\frac{\pi n}{L_{y}}y\right),\\ Z^{mn}=\left(1-\frac{1}{2}\delta_{0m}-\frac{1}{2}\delta_{0n}\right)\frac{4}{L_{x}L_{y}q_{mn}}e^{-q_{mn}z}\frac{1+e^{-2q_{mn}(L_{z}-z)}}{1-e^{-2q_{mn}L_{z}}},\\ Z^{{}^{\prime}mn}=-\left(1-\frac{1}{2}\delta_{0m}-\frac{1}{2}\delta_{0n}\right)\frac{4}{L_{x}L_{y}}e^{-q_{mn}z}\frac{1-e^{-2q_{mn}(L_{z}-z)}}{1-e^{-2q_{mn}L_{z}}}.\end{array}$$

$$p_{mn}=-\underset{0}{\overset{L_x}{\int }}\underset{0}{\overset{L_y}{\int }}(B_z{)}_{mn}{\tilde{g}}_i\mathrm{d}x\mathrm{d}y,(z=0).$$

$$(\varphi ,\mathbf{B},\mathbf{A})=({\varphi }_a,{\mathbf{B}}_a,{\mathbf{A}}_a)+\sum _{i=1}^6(\tilde{\varphi },\tilde{\mathbf{B}},\tilde{\mathbf{A}}).$$

$$ϵ=\frac{\left|{\mathbf{B}}_{\mathrm{calc}}-{\mathbf{B}}_{\mathrm{model}}\right|}{\left|{\mathbf{B}}_{\mathrm{model}}\right|}.$$

$${⟨ϵ⟩}_{\mathrm{mean}}=\mathrm{mean}(ϵ).$$

$${⟨ϵ⟩}_{\mathrm{median}}=\mathrm{median}(ϵ).$$

$${⟨ϵ⟩}_{\mathrm{w}}=\frac{\mathrm{mean}(\left|{\mathbf{B}}_{\mathrm{model}}\right|ϵ)}{\mathrm{mean}(\left|{\mathbf{B}}_{\mathrm{model}}\right|)}.$$

$${⟨ϵ⟩}_{\mathrm{w}}=\mathrm{mean}\left(2\frac{\left|{\varphi }_{\mathrm{calc}}-{\varphi }_{\mathrm{model}}\right|}{\left|{\varphi }_{\mathrm{calc}}\right|+\left|{\varphi }_{\mathrm{model}}\right|}\right).$$

$${⟨{ϵ}_z⟩}_{\mathrm{mean}}=\mathrm{mean}(ϵ):(z=\mathrm{const}).$$

$${⟨{ϵ}_z⟩}_{\mathrm{median}}=\mathrm{median}(ϵ):(z=\mathrm{const}).$$

Full text

{article}

{opening}

1 Introduction

Currently, many studies of the activity in solar active regions rely on the photospheric magnetic field vector measurements and the results of their subsequent processing. The latter includes resolving the 180-degree ambiguity, non-linear force free field (NLFFF) extrapolation, and estimating macroscopic magnetic characteristics, such as free energy and relative helicity. Research activities in these areas are aimed at finding fast and reliable methods for solving these problems on the basis of various physical and mathematical concepts. During the last decades, a significant progress has been made in these areas \inlineciteciteM1;\inlineciteciteM2; \inlineciteciteM3; \inlineciteciteM4; \inlineciteciteM5; \inlineciteciteM6; \inlineciteciteM7; \inlineciteciteM8; \inlineciteciteM9; \inlineciteciteM10; \inlineciteciteM11; \inlineciteciteM12; \inlineciteciteMVal; \inlineciteciteM13, providing us with a reliable (taking into account limitations of the available measurements) description of the three-dimensional magnetic structure and dynamics of real magnetic regions in some cases \inlineciteciteA1; \inlineciteciteA2; \inlineciteciteA3; \inlineciteciteA4; \inlineciteciteA5; \inlineciteciteSun; \inlineciteciteKost.

Despite the significant progress, the accuracy of existing magnetic data processing methods requires improvement. For instance, free energy estimations should have relative error of less than $10^{-2}$ to make tracking of flare energy release possible \inlineciteciteSun;\inlineciteciteKost. Similar accuracy, presumably, is required for relative magnetic helicity calculations. Current methods provide an accuracy of $10^{-1}$ – $10^{-2}$ in relative error [Sun et al. (2012)], which is not enough for a reliable description of energy dynamics in real active regions.

In this paper we present the analysis of three methods for potential field calculation in a 3D rectangular box. This problem is one of key elements of the magnetic free energy and relative helicity estimation. Seeking a precise computation of the potential field, we, at least partially, solve the problem of the accuracy of free energy estimation for solar active regions. However, the potential field calculation is not the only source of uncertainty in this problem. A part is caused by the systematic errors in the measurements of the photospheric magnetic field and NLFFF extrapolation issues. These errors are hard to estimate and, in principle, they can significantly exceed the required threshold.

We focus on two aspects of the considered numerical solutions: accuracy and performance. We demonstrate that a balance between the accuracy and the execution speed is achieved for one the methods, which can calculate both the potential magnetic field and its vector potential.

2 Derivation of the Potential Field, Scalar Potential, and Normalised Vector Potential for the Neuman BVP in a Rectangular Box

Let us consider the boundary value problem (BVP) in a rectangular box $\overline{V}$ with Neuman boundary conditions on its sides. In the following, we describe the algorithm for calculating three representation of this problem: scalar potential $\phi$, magnetic field ${\mathbf{B}}$ and vector potential ${\mathbf{A}}$. The problem consist of solving the following equations for $\phi$, ${\mathbf{B}}$ and ${\mathbf{A}}$:

[TABLE]

[TABLE]

[TABLE]

Here ${\mathbf{n}}$ is the outward normal to the boundary of $\overline{V}$, and $g$ is a normal field component given on the boundary surface $S$: $V\cup S=\overline{V}(x,y,z)$, $\left(x\in[0,L_{x}];\ y\in[0,L_{y}];\ z\in[0,L_{z}]\right)$. We assume that the divergence free condition for $g$ is satisfied:

[TABLE]

In the following, we consider the above BVP (Equations (1) – (3)) without imposing the optional boundary condition in Equation (3). However this condition can be satisfied by adding to the calculated vector potential a harmonic function, found using the algorithm for solving the BVP formulated in Equation (2).

Our algorithm is based on the reduction of the original problem to a series of BVPs with additional to Equation (4) the requirement of a flux balance:

[TABLE]

Thus, the magnetic flux must be balanced for each side of the rectangular volume $\overline{V}$ individually. In the general case, the boundary conditions do not satisfy this requirement (Equation (5)), but we always can represent an arbitrary boundary condition in the following form:

[TABLE]

where $g_{a}$ is the boundary condition for a potential field ${\mathbf{B}}_{a}$ that compensates the unbalanced flux $G_{i}$ on each side $i$:

[TABLE]

Construction of such a compensating field $(\phi_{a},{\mathbf{B}}_{a},{\mathbf{A}}_{a})$ in analytical form is the first essential step of our algorithm. When the compensating field is found, one can represent the general solution in the following form:

[TABLE]

where $(\phi_{i},{\mathbf{B}}_{i},{\mathbf{A}}_{i})$ are the solutions of the BVPi (Equations (1) – (4)) for six sides of the volume $\overline{V}$.

2.1 Construction of a Compensating Potential Field

Let us find the compensating potential field $(\phi_{a},{\mathbf{B}}_{a},{\mathbf{A}}_{a})$ as a linear combination of five harmonic solutions.

[TABLE]

The essential requirement to the harmonic fields ${\mathbf{B}}_{a}^{j}$ is linear independence. Hence, one can select an arbitrary set of linearly independent harmonic functions. For the case of simplicity, we have selected them to be simple quadratic linearly independent polynomials:

[TABLE]

[TABLE]

[TABLE]

Since the selection of Equation (10) is arbitrary, one can freely select another set. For example, $-x^{2}+z^{2}$ and $-y^{2}+z^{2}$ can be replaced by $xz$ and $yz$.

The fluxes through the box sides can be calculated analytically for each basis solution:

[TABLE]

For instance, let us calculate the flux of the basis field ${\mathbf{B}}_{a}^{3}$ through the bottom side of the box $(x,y,0)$. Applying Equation (13), we get

[TABLE]

After calculation of the $G_{ai}^{j}$ values, we find the unknown coefficients $m_{j}$ by solving a linear system of algebraic equations, representing the flux balance conditions on five sides (Equation (7)):

[TABLE]

The condition in Equation (14) in combination with the divergence-free condition (Equation (4)) automatically provides flux balance (Equation (7)) on all six sides.

The found solution $(\phi_{a},{\mathbf{B}}_{a},{\mathbf{A}}_{a})$ allows us to obtain modified boundary conditions $\widetilde{g}_{i}$ balanced on each side of the volume:

[TABLE]

2.2 Solution of a BVP with Non-zero Boundary Conditions on One Side of a Rectangular Box

We present two different approaches for solving this problem: an analytic and a numerical one. In the first case, we derive exact analytic equations for the magnetic field components, and its scalar and vector potentials. In practical implementations the integrals which appear are computed directly with the use of 2D local quadratic splines allowing for analytical integration of their convolution with trigonometric functions. The second approach solves the problem in the discrete Fourier space and uses the fast Fourier transform (FFT) to speed up the computation.

2.2.1 Analytical Solution

Let us consider in detail the solution of one of the sub-problems $($BVP$)_{i}$ $(g=\widetilde{g}_{i})$ with non-zero boundary conditions $\widetilde{g}_{i}$ given on the bottom side of the computational box $(x,y,0)$. For other sides, the solution can be found in the same way. We search for the solution in the form of a decomposition into a set of harmonic basis functions of the following form:

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Each partial harmonic solution $(\phi_{mn},{\mathbf{B}}_{mn},{\bf A}_{mn})$ satisfies Equations (1) – (4) with $g\neq 0$ on the side $(x,y,0)$ and $g=0$ on the other sides. The coefficients $p_{mn}$ are determined by the following formula:

[TABLE]

The solutions $(\widetilde{\phi}_{i},\widetilde{{\mathbf{B}}}_{i},\widetilde{{\mathbf{A}}}_{i})$ found in this way for each side $i$ allow us to obtain the final solution of the problem (Equations (1) – (4)) as the following sum:

[TABLE]

2.3 Variants of the BVP Numerical Implementation

In a practical implementation of this scheme, the number of terms in the expansion Equations (16) – (18) can be limited by the size of the grid $(x,y)$. In the following, we benchmark two different implementations of our method. In the first implementation (CASE I), the coefficients $p_{mn}$ are calculated by the direct integration of Equation (21). To improve the accuracy, we present the boundary condition $\widetilde{g}_{i}$ in the form of 2D local quadratic splines allowing for analytical integration of their convolution with trigonometric functions. In our second implementation, we calculate $p_{mn}$ coefficients with the use of the FFT (CASE II). Then the components of the field and vector potential are computed in terms of FFT coefficients in a form similar to Equations (17) and (18).

As a reference, we also include in our test the results of a proprietary Poisson solver built into the Intel MKL library (CASE III). Since the above package allows us to calculate only the scalar potential $\phi$, the magnetic field ${\mathbf{B}}$ is computed as the gradient of ${\phi}$ (Equation (2), the third equation of the system) using the finite difference scheme. The vector potential ${\mathbf{A}}$ is calculated from the obtained field ${\mathbf{B}}$ with the use of the algorithm developed by [Valori et al. (2012)]. Note that the vector potential ${\mathbf{A}}$ computed by this algorithm generally does not satisfy the conditions Equation (3) (the second equation of the system) and Equation (3) (the first equation of the system).

3 Numerical Tests

We test our method on two models: the field of two magnetic charges situated outside the box and the magnetic field of a typical solar active region obtained by potential field extrapolation from the normal component measured at the photospheric level under the condition of a finite field at infinite distance.

3.1 Notations and Metrics

Let us introduce the following notations and metrics:

• •

Numerical solution: $(\phi_{\rm calc},{\mathbf{B}}_{\rm calc},{\mathbf{A}}_{\rm calc})$

• •

Model solution: $(\phi_{\rm model},{\mathbf{B}}_{\rm model},{\mathbf{A}}_{\rm model})$

• •

Local relative error:

[TABLE]

• •

Average relative error:

[TABLE]

• •

Median relative error:

[TABLE]

• •

Weighted mean relative error:

[TABLE]

• •

The average error of the calculated scalar potential:

[TABLE]

• •

Average relative error:

[TABLE]

• •

Median relative error:

[TABLE]

• •

Weighted mean relative error:

[TABLE]

• •

Relative amount of nodes with an error is less than the $\epsilon$:

[TABLE]

• •

Magnetic field energy inside the computational box $V$ (without boundaries):

[TABLE]

• •

Energy, calculated by the virial theorem. Here $S^{\prime}$ is the boundary of the inner region $V^{\prime}$ and ${\mathbf{n}}$ is the outward normal to the boundary $S^{\prime}$:

[TABLE]

• •

Energy, calculated from the normal component of the magnetic field and scalar potential given on the Surface $S^{\prime}$:

[TABLE]

3.2 Model of Two Charges

The model of two charges ($M_{\rm charges}$) is defined analytically as follows:

[TABLE]

[TABLE]

whith $q_{1}=1000$, $q_{2}=-1500$, ${{\mathbf{r}}}_{1}=(-0.2,-0.18,-0.22)$, ${{\mathbf{r}}}_{2}=(0.2,0.18,-0.22)$, $\overline{V}=[-1,0.233]\times[-1,0.953]\times[0,1.610]$, and grid size $\left[n_{x},n_{y},n_{z}\right]=[37,58,48].$

The model parameters are selected to give a magnetic field configuration with large values on one of the side faces of the box, as demonstrated in Figure 1.

For the first assessment of the method’s accuracy, we present the relative errors in the $B_{x}$ component extracted from two arbitrarily chosen one-dimensional cross-sections of the computational box (See Figures 2 – 5). All methods show good agreement with the model, both inside the box, and on its boundaries. The best accuracy is demonstrated by Case I, followed by the Case II. As expected, the accuracy of the field calculated from the vector potential $({\bf\nabla}\times{\mathbf{A}})_{{\rm calc},x}$ using the finite difference method is similar for all 3 cases.

The global metrics Equations (24) – (27) characterizing the accuracy of the three BVP solutions in the volume $V$ are presented in Tables 1 and 2. The metrics for the boundary $S$ are given in Tables 3 and 4.

Tables 1 and 2 confirm that Case I provides the most accurate solution. Both volume and surface metrics are an order of magnitude better than in Case II. The Case II method is more accurate than Case III by an order of magnitude for most of the metrics. The metrics for the magnetic field reconstructed from the vector potential ${\bf\nabla}\times{\mathbf{A}}_{\rm calc}$ (Tables 3 and 4) are less precise due to the usage of the finite difference approximation. However, the relative errors of all methods are small both in the volume and on the boundary, demonstrating sufficiently high quality for practical application.

The additional metrics Equations (28) – (31) shown in Figures 6 – 8 demonstrate that the solution accuracy increases with height ($z$). Hence, all algorithms tend to produce more accurate solutions in areas of smooth and low magnitude fields. Figure 8 shows that the fraction of grid points with large relative errors ($>1\%$) is extremely low in the entire volume $V$.

In Table 5, we compare the values of the total magnetic energy calculated in different ways: by integration over the volume $V$ (Equation (32)) and over the boundary $S$ (Equations (33) and (34)). Also, Table 5 gives an idea of how accurate are the different algorithms for the calculation of the magnetic field energy that is an important macroscopic parameter of a solar active region.

3.3 Model of a Real Active Region

The realistic $M_{\rm AR}$ is built by the potential extrapolation of the magnetic field of AR 11158 from the normal magnetic field component observed at the photospheric level by SDO/HMI at 2011-02-15 01:36:00 UT. The potential field calculation is performed using the FFT based algorithm proposed by [Alissandrakis (1981)]. The size of the computational box is 300$\times$300$\times$256 voxels. Results of the analysis are presented in the form of the same tables and figures as in Section 3.2.

Case I and Case II demonstrate similar accuracy for most of the metrics calculated for the $M_{\rm AR}$ model while the metrics for Case III are approximately two orders of magnitude worse. The magnetic energy calculations are also much more accurate for Case I and Case II in comparison with the Case III results.

3.4 Computing times for the three algorithms

All calculations are performed on a computer with an Intel(R) Core(TM) i5-4460 3.20GHz processor. The computation times for different algorithms are the following:

• •

For $M_{\rm charges}$

• –

Case I -16.48 sec (without ${\bf\nabla}\times{\mathbf{A}}$), 21.18 sec (with ${\bf\nabla}\times{\mathbf{A}}$);

• –

Case II - 0.23 sec (without ${\bf\nabla}\times{\mathbf{A}}$), 1.24 sec (with ${\bf\nabla}\times{\mathbf{A}}$);

• –

Case III- 0.065 sec (without ${\bf\nabla}\times{\mathbf{A}}$), 1.02 sec (with ${\bf\nabla}\times{\mathbf{A}}$).

• •

For $M_{\rm AR}$

• –

Case I -13.6 h (without ${\bf\nabla}\times{\mathbf{A}}$), 27.4 h (with ${\bf\nabla}\times{\mathbf{A}}$);

• –

Case II - 202.91 sec (without ${\bf\nabla}\times{\mathbf{A}}$), 424.25 sec (with ${\bf\nabla}\times{\mathbf{A}}$);

• –

Case III - 0.93 sec (without ${\bf\nabla}\times{\mathbf{A}}$), 253.05 sec (with ${\bf\nabla}\times{\mathbf{A}}$).

4 Conclusions

We consider three different algorithms for the numerical solution of the Neumann BVP for the Laplace equation in a 3D rectangular box. The first two options are based on dividing the whole problem into six sub-problems corresponding to six faces of the box. The solution of each sub-problem uses the ordinary Fourier decomposition of the solution into a set of harmonic functions satisfying the Laplace equation. In one variant (Case I) the solution is obtained through direct integration while in the other case (Case II) the FFT is applied to speed up the calculation. The third algorithm (Case III) uses the proprietary Poisson solver from the Intel MKL library. The solutions are presented as a set of three quantities, namely the magnetic field itself, and its scalar and vector potentials. In Cases I and II all of the quantities are expressed and calculated as a superposition of the harmonic functions. The Intel MKL library used in Case III allows for computation of the scalar potential only. The magnetic field vector is then computed using the finite difference scheme. The vector potential is calculated from the field using the method proposed by [Valori et al. (2012)].

All three methods provide high accuracy of the potential field calculation. However, the Case I and Case II methods proposed here provide much more precise solutions with metrics which are around two orders of magnitude better than the Poisson solver supplied with the Intel MKL library. The Case II method based on FFT provides the possibility of fast and accurate calculation of the potential field energy that is required for the estimation of the free magnetic in energy in solar active regions. Although the direct spline based integration (Case II gives the best results among all considered methods, it is a very computationally expensive approach especially for high resolution grids. Hover, this algorithm can be used for small and moderate-size computational domains. It also allows for potential field computation at arbitrary location, providing a continuous and precise solution. Although the Poison solver from the Intel MKL library (Case III) is not as precise as other approaches, it is very fast and gives the result in less than a second even for the high resolution grids.

In conclusion, we should note the quality of the algorithm for the vector potential calculation [Valori et al. (2012)] that is used in the Case III method to compute the vector potential of a given potential field. Our tests show that the reverse reconstruction of the potential field from the vector potential does not lead to significant worsening of the solution. It confirms the efficiency of the method, that also can be applied in the more general case of an arbitrary divergence-free field.

Acknowledgements

This study was supported by the Russian Foundation of Basic Research under grants 15-02-01077, 16-32-00315_mol_a, 15-02-03835_a, and 15-02-01089_a; by the ISSI International Team on Magnetic Helicity estimations in models and observations of the solar magnetic field. This study was supported by the Program of basic research of the RAS Presidium No. 7.

Disclosure of Potential Conflicts of Interest

The authors declare that they have no conflicts of interest.