The PDFI_SS Electric Field Inversion Software

We start with the PTD for the magnetic field ${\boldsymbol{B}}$ in spherical coordinates (Chandrasekhar 1961; Backus 1986) in terms of the poloidal potential P and the toroidal potential T:

$$\begin{eqnarray}&&{\boldsymbol{B}}={\boldsymbol{\nabla }}\times {\boldsymbol{\nabla }}\times P\hat{{\boldsymbol{r}}}+{\boldsymbol{\nabla }}\times T\hat{{\boldsymbol{r}}},\end{eqnarray} \tag{ 4 }$$

where $\hat{{\boldsymbol{r}}}$ is the unit vector in the radial direction. Here, in a change from the notation used in KFW14, we use P for the poloidal potential, instead of ${ \mathcal B }$, and T for the toroidal potential, instead of ${ \mathcal J }$. This change was made for notational simplicity and also corresponds with the notation used by Lumme et al. (2017). We also note another useful form for Equation (4), namely,

$$\begin{eqnarray}&&{\boldsymbol{B}}=-\hat{{\boldsymbol{r}}}{{\rm{\nabla }}}_{h}^{2}P+{{\boldsymbol{\nabla }}}_{h}\left(\displaystyle \frac{\partial P}{\partial r}\right)+{\boldsymbol{\nabla }}\times T\hat{{\boldsymbol{r}}}.\end{eqnarray} \tag{ 5 }$$

The operator ${{\rm{\nabla }}}_{h}^{2}$ is the horizontal Laplacian, i.e., the full Laplacian but omitting the radial derivative contribution, and ${{\boldsymbol{\nabla }}}_{h}(\partial P/\partial r)$ represents the horizontal components of the gradient of the radial derivative of P. By uncurling Equation (4), it is clear that the vector potential ${\boldsymbol{A}}$ can be written in terms of P and T as

$$\begin{eqnarray}&&{\boldsymbol{A}}={\boldsymbol{\nabla }}\times P\hat{{\boldsymbol{r}}}+T\hat{{\boldsymbol{r}}},\end{eqnarray} \tag{ 6 }$$

where we have omitted an explicit gauge term.

The PTD, or "inductive" electric field ${{\boldsymbol{E}}}^{P}$, is related to the magnetic field ${\boldsymbol{B}}$ through Faraday's law:

$$\begin{eqnarray}&&\dot{{\boldsymbol{B}}}=-{\boldsymbol{\nabla }}\times c{{\boldsymbol{E}}}^{P},\end{eqnarray} \tag{ 7 }$$

where c is the speed of light, and where we use the overdot to denote a partial time derivative. Substituting Equation (4) into Faraday's law and uncurling, we find

$$\begin{eqnarray}&&c{{\boldsymbol{E}}}^{P}=-{\boldsymbol{\nabla }}\times \dot{P}\hat{{\boldsymbol{r}}}-\dot{T}\hat{{\boldsymbol{r}}},\end{eqnarray} \tag{ 8 }$$

where $\dot{P}$ and $\dot{T}$ are the partial time derivatives of P and T. The general description of the electric field will also include the gradient of a scalar potential in addition to the inductive solution in Equation (8), but we omit any explicit gradient contributions to $c{\boldsymbol{E}}$ here, and we discuss the gradient contributions in subsections further below.

By evaluating the radial component of Equation (7) when substituting Equation (8) for ${{\boldsymbol{E}}}^{P}$, we find that $\dot{P}$ obeys the 2D Poisson equation

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{h}^{2}\dot{P}=-{\dot{B}}_{r},\end{eqnarray} \tag{ 9 }$$

where B_r is the radial magnetic field component. Here, the right-hand side of Equation (9) is viewed as a source term that can be evaluated from the magnetogram data. By taking the radial component of the curl of Equation (7), we find that $\dot{T}$ obeys the Poisson equation

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{h}^{2}\dot{T}=-\hat{{\boldsymbol{r}}}\cdot \left({\boldsymbol{\nabla }}\times {\dot{{\boldsymbol{B}}}}_{h}\right),\end{eqnarray} \tag{ 10 }$$

where ${{\boldsymbol{B}}}_{h}$ are the horizontal components of ${\boldsymbol{B}}$. A third useful Poisson equation can be found by taking the divergence of the horizontal components of Equation (7):

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{h}^{2}\left(\displaystyle \frac{\partial \dot{P}}{\partial r}\right)={{\boldsymbol{\nabla }}}_{h}\cdot {\dot{{\boldsymbol{B}}}}_{h}.\end{eqnarray} \tag{ 11 }$$

The quantity $\partial \dot{P}/\partial r$ is important because it allows one to evaluate the radial derivative of the horizontal electric field components. To see this, one can evaluate the quantity

$$\begin{eqnarray}&&\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left({rc}{{\boldsymbol{E}}}_{h}^{P}\right),\end{eqnarray} \tag{ 12 }$$

where $c{{\boldsymbol{E}}}_{h}^{P}$ represents the horizontal components of $c{{\boldsymbol{E}}}^{P}$ from Equation (8):

$$\begin{eqnarray}&&c{{\boldsymbol{E}}}_{h}^{P}=-{\boldsymbol{\nabla }}\times \dot{P}\hat{{\boldsymbol{r}}}.\end{eqnarray} \tag{ 13 }$$

The θ and ϕ components of the curl in Equation (13) both contain leading factors of 1/r, as can be seen from the first term of Equation (1) and the second term of Equation (2). The radial derivative in Equation (12) therefore is applied directly to $\dot{P}$, resulting in

$$\begin{eqnarray}&&\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left({rc}{{\boldsymbol{E}}}_{h}^{P}\right)=-{\boldsymbol{\nabla }}\times \displaystyle \frac{\partial \dot{P}}{\partial r}\hat{{\boldsymbol{r}}}.\end{eqnarray} \tag{ 14 }$$

Expanding the radial derivative on the left-hand side of Equation (14), we then arrive at this expression for the radial derivative of ${{\boldsymbol{E}}}_{h}^{P}$:

$$\begin{eqnarray}&&c\displaystyle \frac{\partial {{\boldsymbol{E}}}_{h}^{P}}{\partial r}=-{\boldsymbol{\nabla }}\times \displaystyle \frac{\partial \dot{P}}{\partial r}\hat{{\boldsymbol{r}}}-c\displaystyle \frac{{{\boldsymbol{E}}}_{h}^{P}}{r}.\end{eqnarray} \tag{ 15 }$$

Here, the quantity r is the radius of the surface upon which the 2D Poisson equations are solved, which for nearly all of our purposes can be taken as the radius of the Sun R_⊙. Equation (15) was not given in KFW14 but, as shown here, is easy to derive. Note that the second term on the right-hand side of Equation (15) goes to zero as $r\to \infty$, meaning that in the Cartesian limit, this term vanishes. The radial derivative of the horizontal inductive electric field is useful because it allows one to compute the horizontal components of ${\boldsymbol{\nabla }}$ × ${\boldsymbol{E}}$.

The availability of time cadences of vector magnetic field measurements, such as from the HMI instrument on NASA's SDO Mission (Scherrer et al. 2012), enables the evaluation of the time derivatives as the source terms in the above Poisson equations, making such electric field solutions possible. With the data in hand, evaluation of ${{\boldsymbol{E}}}^{P}$ becomes a matter of solving the above Poisson equations on a region of the Sun's surface and then evaluating Equation (8).

The Doppler velocity, when combined with the magnetic field measurements, provides additional information about the electric field beyond the inductive contribution ${{\boldsymbol{E}}}^{P}$ (Ravindra et al. 2008). The relationship between the measured line-of-sight (LOS) velocity vector ${{\boldsymbol{V}}}_{\mathrm{LOS}}$ and the true plasma velocity ${\boldsymbol{V}}$ is given by

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{\mathrm{LOS}}=({\boldsymbol{V}}\cdot \hat{{\boldsymbol{\ell }}})\hat{{\boldsymbol{\ell }}},\end{eqnarray} \tag{ 16 }$$

where $\hat{{\boldsymbol{\ell }}}$ is the LOS unit vector pointing toward the observer from the surface of the Sun. Note that $\hat{{\boldsymbol{\ell }}}$ is a function of position on the solar surface, since the Sun's surface is curved. Near an LOS polarity inversion line (PIL), the ${{\boldsymbol{V}}}_{\mathrm{LOS}}$ flow carrying transverse components of the magnetic field perpendicular to $\hat{{\boldsymbol{\ell }}}$ results in an electric field contribution

$$\begin{eqnarray}&&c{{\boldsymbol{E}}}_{\mathrm{PIL}}=-{{\boldsymbol{V}}}_{\mathrm{LOS}}\times {{\boldsymbol{B}}}_{t},\end{eqnarray} \tag{ 17 }$$

where ${{\boldsymbol{B}}}_{t}={\boldsymbol{B}}-({\boldsymbol{B}}\cdot \hat{{\boldsymbol{\ell }}})\hat{{\boldsymbol{\ell }}}$ represents the components of ${\boldsymbol{B}}$ transverse to $\hat{{\boldsymbol{\ell }}}$.

When we are not near an LOS PIL, a nonzero Doppler velocity is less certain to be coming from a flow that transports ${{\boldsymbol{B}}}_{t}$, and instead could be a signature of flows parallel to ${\boldsymbol{B}}$, which have no electric field consequences. To account for this uncertainty away from LOS PILs, the electric field in Equation (17) is modulated by an empirical factor w_LOS given by the following expression:

$$\begin{eqnarray}&&{w}_{\mathrm{LOS}}=\exp \left(-\ \displaystyle \frac{1}{{\sigma }_{\mathrm{PIL}}^{2}}{\left|\displaystyle \frac{{B}_{\mathrm{LOS}}}{{{\boldsymbol{B}}}_{t}}\right|}^{2}\right),\end{eqnarray} \tag{ 18 }$$

where B_LOS is the LOS component of ${\boldsymbol{B}}$, and σ_PIL is an empirically adjustable parameter, commonly taken as unity.

To the extent that the Doppler contribution to the electric field contributes magnetic evolution, that contribution should already have been included in the inductive contribution described in Section 2.1. We therefore want to include any additional curl-free contribution to the electric field from the Doppler term. To do this, we will represent the Doppler contribution in Equation (17) by the gradient of a scalar potential, which we will call ψ^D:

$$\begin{eqnarray}&&c{{\boldsymbol{E}}}^{D}=-{\boldsymbol{\nabla }}{\psi }^{D},\end{eqnarray} \tag{ 19 }$$

where we note that this form automatically results in zero curl.

The details of how the equation defining ψ^D is derived are provided in Section 2.3.3 in KFW14. Here, we will simply write down the result, modified slightly to account for working in spherical, rather than Cartesian, coordinates:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\nabla }}}_{h}{\psi }^{D}\cdot {\hat{{\boldsymbol{q}}}}_{h}+{q}_{r}\displaystyle \frac{\partial {\psi }^{D}}{\partial r} & = & {w}_{\mathrm{LOS}}{\left({{\boldsymbol{V}}}_{\mathrm{LOS}}\times {{\boldsymbol{B}}}_{t}\right)}_{h}\cdot {\hat{{\boldsymbol{q}}}}_{h}\\ & & +\ {q}_{r}{w}_{\mathrm{LOS}}{\left({{\boldsymbol{V}}}_{\mathrm{LOS}}\times {{\boldsymbol{B}}}_{t}\right)}_{r},\end{array}\end{eqnarray} \tag{ 20 }$$

where $\hat{{\boldsymbol{q}}}$ is the unit vector pointing in the same direction as ${{\boldsymbol{V}}}_{\mathrm{LOS}}$ × ${{\boldsymbol{B}}}_{t}$, q_r is the radial component of $\hat{{\boldsymbol{q}}}$, and ${\hat{{\boldsymbol{q}}}}_{h}$ are the horizontal components of $\hat{{\boldsymbol{q}}}$. Equation (20) is solved using the "iterative" technique developed by coauthor Brian Welsch, initially described in Section 3.2 of Fisher et al. (2010), with subsequent changes discussed in Section 2.2 of KFW14. In this article, the current version of the iterative technique for PDFI_SS is described in Section 3.9.2. The iterative technique involves repeated solutions of a 2D Poisson equation that tries to best represent the Doppler electric field from the observed data by the gradient of ψ^D in the $\hat{{\boldsymbol{q}}}$-direction.

There are many techniques currently available to estimate velocities in the directions parallel to the solar surface by solving the "optical flow" problem on pairs of images closely adjacent in time to estimate these flows (see, e.g., reviews by Welsch et al. 2007; Schuck 2008; Tremblay et al. 2018). Here we estimate horizontal flow velocities ${{\boldsymbol{V}}}_{h}^{F}$ using the FLCT local correlation tracking code (Fisher & Welsch 2008) applied to images of B_r from the vector magnetogram sequence. The choice of FLCT is somewhat arbitrary; any other existing technique could be used as an alternative. We chose FLCT because we are very familiar with the algorithm and the code and have spent years making the code as computationally efficient as possible.

The use of FLCT in spherical geometry introduces some complications, since the FLCT algorithm is based strictly on assumptions of Cartesian geometry. We adopt the solution proposed in the appendix of Welsch et al. (2009), in which the B_r images are mapped to a Mercator projection, FLCT is run, and then the velocities are interpolated and rescaled back to spherical geometry. The FLCT code has been updated to perform this operation automatically if the input data are specified as being equally spaced in longitude and latitude, i.e., on a Plate Carrée grid. More details on recent versions of FLCT can be found in Section 4.4 and in the updated source code and documentation for FLCT, available at the software repository http://cgem.ssl.berkeley.edu/cgi-bin/cgem/FLCT/home.

Horizontal velocities ${{\boldsymbol{V}}}_{h}^{F}$ estimated from FLCT, and acting on the radial component of the magnetic field, provide an additional contribution to the electric field:

$$\begin{eqnarray}&&c{{\boldsymbol{E}}}_{\mathrm{FLCT}}=-{{\boldsymbol{V}}}_{h}^{F}\times {B}_{r}\hat{{\boldsymbol{r}}}.\end{eqnarray} \tag{ 21 }$$

We have neglected an additional term, contributing to E_r, coming from $-{{\boldsymbol{V}}}_{h}^{F}\times {{\boldsymbol{B}}}_{h}$ in Equation (21). In KFW14, we showed that in Cartesian coordinates this term is already accounted for by the inductive contribution to ${E}_{z}^{P}$. The same argument applies here, except the contribution is to ${E}_{r}^{P}$, the radial component.

In regions where the radial magnetic field component is small compared to the horizontal magnetic field components, we trust the FLCT velocities less because the radial magnetic field evolution is less likely to be due to advection by horizontal flows. Therefore, we introduce an empirical modulation function (1−w_r) that multiplies the right-hand side of Equation (21) and that reduces the amplitude of the electric field when the magnetic field is mostly horizontal. The quantity w_r is given by the expression

$$\begin{eqnarray}&&{w}_{r}=\exp \left(-\ \displaystyle \frac{1}{{\sigma }_{\mathrm{PIL}}^{2}}{\left|\displaystyle \frac{{B}_{r}}{{{\boldsymbol{B}}}_{h}}\right|}^{2}\right),\end{eqnarray} \tag{ 22 }$$

where B_r is the radial magnetic field component and ${{\boldsymbol{B}}}_{h}$ are the horizontal magnetic field components. The empirical factor σ_PIL is the same empirical factor (typically unity) used in the above discussion of the electric field from the Doppler velocity.

To avoid including inductive electric fields that are already accounted for by ${{\boldsymbol{E}}}^{P}$, we remove any inductive contributions by writing the FLCT-derived electric field in terms of the gradient of a scalar potential ψ^F:

$$\begin{eqnarray}&&c{{\boldsymbol{E}}}^{F}=-{\boldsymbol{\nabla }}{\psi }^{F}.\end{eqnarray} \tag{ 23 }$$

To derive an equation for ψ^F, we can take the horizontal divergence of $c{{\boldsymbol{E}}}^{F}$ and set it equal to the divergence of $(1-{w}_{r})c{{\boldsymbol{E}}}_{\mathrm{FLCT}}$, where $c{{\boldsymbol{E}}}_{\mathrm{FLCT}}$ is taken from Equation (21):

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}^{2}{\psi }^{F}={{\boldsymbol{\nabla }}}_{h}\cdot ((1-{w}_{r}){{\boldsymbol{V}}}_{h}^{F}\times {B}_{r}\hat{{\boldsymbol{r}}}).\end{eqnarray} \tag{ 24 }$$

Once this Poisson equation is solved for ψ^F, $c{{\boldsymbol{E}}}^{F}$ can be evaluated from Equation (23).

Most of the time, we expect that electric fields in the solar photosphere will be largely determined by the electric field in ideal MHD, namely, $c{\boldsymbol{E}}=-{\boldsymbol{V}}\times {\boldsymbol{B}}$, where ${\boldsymbol{V}}$ is the local plasma velocity. A consequence of this is that we expect ${\boldsymbol{E}}\cdot {\boldsymbol{B}}=0$. However, we found that if the PTD (inductive) contribution, Doppler contribution, and FLCT contribution are added together, the resulting electric field can have a significant component parallel to the direction of ${\boldsymbol{B}}$. We therefore want to find a way to add a scalar potential electric field

$$\begin{eqnarray}&&c{{\boldsymbol{E}}}^{I}=-{\boldsymbol{\nabla }}{\psi }^{I},\end{eqnarray} \tag{ 25 }$$

such that

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}{\psi }^{I}\cdot {\boldsymbol{B}}=c{{\boldsymbol{E}}}^{\mathrm{PDF}}\cdot {\boldsymbol{B}},\end{eqnarray} \tag{ 26 }$$

where ${{\boldsymbol{E}}}^{\mathrm{PDF}}$ is the sum of the PTD (inductive), Doppler, and FLCT electric field contributions. When ${{\boldsymbol{E}}}^{I}$ is added to the electric field, the result should have ${{\boldsymbol{E}}}^{\mathrm{PDFI}}\cdot {\boldsymbol{B}}\approx 0$, where ${{\boldsymbol{E}}}^{\mathrm{PDFI}}$ is now the complete PDFI electric field solution. In Section 2.2 of KFW14, the equation that ψ^I and its depth derivative obey is given in Equation (19) of that article. The form of the equation is changed only slightly in spherical coordinates and is

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}{\psi }^{I}\cdot {\hat{{\boldsymbol{b}}}}_{h}+{b}_{r}\displaystyle \frac{\partial {\psi }^{I}}{\partial r}=c{{\boldsymbol{E}}}_{h}^{\mathrm{PDF}}\cdot {\hat{{\boldsymbol{b}}}}_{h}+{{cE}}_{r}^{\mathrm{PDF}}{b}_{r}.\end{eqnarray} \tag{ 27 }$$

Here, $\hat{{\boldsymbol{b}}}$ is the unit vector pointing in the direction of ${\boldsymbol{B}}$, and b_r and ${\hat{{\boldsymbol{b}}}}_{h}$ are the radial and horizontal components of $\hat{{\boldsymbol{b}}}$, respectively. As described in Section 2.2 of KFW14, Equation (27) is solved using the "iterative" technique, also mentioned earlier in Section 2.2 of this article, with further details given in Section 3.9.2. Once ψ^I and ∂ψ^I/∂r have been found, the full PDFI solution is given by

$$\begin{eqnarray}&&c{{\boldsymbol{E}}}^{\mathrm{PDFI}}=c{{\boldsymbol{E}}}_{h}^{\mathrm{PDF}}-{{\boldsymbol{\nabla }}}_{h}{\psi }^{I}+\hat{{\boldsymbol{r}}}\left({{cE}}_{r}^{\mathrm{PDF}}-\displaystyle \frac{\partial {\psi }^{I}}{\partial r}\right).\end{eqnarray} \tag{ 28 }$$

It is important to note that this same procedure to "perpendicularize" ${\boldsymbol{E}}$ with respect to ${\boldsymbol{B}}$ can be performed with any combination of the other electric field contributions. For example, in cases where there are no Doppler or FLCT velocity flows available, one can substitute ${{\boldsymbol{E}}}^{P}$ for ${{\boldsymbol{E}}}^{\mathrm{PDF}}$ in Equations (27)–(28) to generate an electric field solution that should still minimize ${\boldsymbol{E}}\cdot {\boldsymbol{B}}$. This is described in detail in Section 2.3.4 of KFW14 (see also Table 1 of KFW14).

Once the PDFI electric fields have been computed, we can use them to estimate the Poynting flux of energy in the radial direction:

$$\begin{eqnarray}&&{S}_{r}=\hat{{\boldsymbol{r}}}\cdot \displaystyle \frac{1}{4\pi }c{{\boldsymbol{E}}}_{h}^{\mathrm{PDFI}}\times {{\boldsymbol{B}}}_{h}.\end{eqnarray} \tag{ 29 }$$

We can also compute the helicity injection rate contribution function, h_r:

$$\begin{eqnarray}&&{h}_{r}=\hat{{\boldsymbol{r}}}\cdot 2\,c{{\boldsymbol{E}}}_{h}^{\mathrm{PDFI}}\times {{\boldsymbol{A}}}_{P},\end{eqnarray} \tag{ 30 }$$

where ${{\boldsymbol{A}}}_{P}={\boldsymbol{\nabla }}\times P\hat{{\boldsymbol{r}}}$, and the poloidal potential P is found from Equation (9) but without the time derivative in the source term. The relative helicity injection rate was derived by Berger & Field (1984) in terms of a surface integral of Equation (30). Schuck & Antiochos (2019) argue that integrating over a finite area, as we do here, and neglecting the other surfaces of our spherical wedge volume domain may result in a loss of gauge invariance. For the time being, we ignore this possible complication and simply use Equation (30) as an integrand to estimate the helicity injection rate. Berger & Hornig (2018) have pointed out that using the PTD formalism for describing the magnetic field, which we employ for computing the inductive contribution to the electric field, has some useful properties. The magnetic helicity in a volume can be understood in terms of the linkage between the contribution to the magnetic field generated by the poloidal potential and that generated from the toroidal potential. That analysis also leads to the same surface integral of the quantity in Equation (30) for the relative helicity injection rate, assuming that the volume integral of ${\boldsymbol{E}}\cdot {\boldsymbol{B}}$ is zero if the coronal plasma is described by ideal MHD, and that the contributions from the other surfaces surrounding the volume can be neglected.

KFW14 studied the accuracy of the PDFI electric field solutions, as well as the Poynting flux and Helicity injection rates, for the case of an MHD simulation of magnetic flux emergence in a convecting medium, originally described by Welsch et al. (2007). They found that including the PTD, Doppler, FLCT, and ideal electric field contributions (in other words, the full PDFI electric field) resulted in the most accurate reconstruction of the electric field from the MHD simulation. The comparison is described in detail in Section 4 of KFW14. For details of tests of the PDFI solution from these simulation data using PDFI_SS, see Section 9.2.

The Helmholtz/Poisson equation subroutines for spherical coordinates in FISHPACK are named HSTSSP (the staggered grid case) and HWSSSP (the centered grid case). Important input arguments to these subroutines include the source term for the Poisson equation (a 2D array) and boundary conditions applied to the four edges of the problem domain (four 1D arrays). Note that the FISHPACK software assumes that the Poisson equation is multiplied by r² (where r is the radius), meaning that the source terms of all Poisson equations in PDFI_SS must also be multiplied by r² before calls to HSTSSP and HWSSSP. The boundary conditions most useful to PDFI_SS include the specification of derivatives of the solution in the directions normal to the domain boundary edges (Neumann boundary conditions). The problem domain is described further below. Both of these subroutines make certain assumptions about the geometry of the domain, the array dimensions, and how the arrays are ordered. To avoid confusion, we adopt exactly the same nomenclature that is used in FISHPACK throughout our software to describe the domain, its boundaries, and the grid spacing.

Spherical coordinates in FISHPACK assume spherical polar coordinates, with the first independent variable θ being colatitude and the second independent variable ϕ being azimuthal angle (longitude). See Section 3.2 and the figures in Section 3.4 for a comparison between spherical polar coordinates and the longitude–latitude coordinate system typically used to display magnetic field data.

Both colatitude and longitude are measured in radians. Colatitude θ ranges from 0 (North Pole) to π (South Pole). Longitude ϕ ranges from 0 to 2π, with negative values of longitude not allowed within these two subroutines. The finite-difference approximations to the Poisson equations are solved in a domain where the edges of the domain are defined by lines of constant colatitude and constant longitude. The northern edge of the domain is defined by θ = a and the southern edge of the domain by θ = b, where 0 < a < b < π. The leftmost edge of the domain is defined by ϕ = c, and the rightmost edge of the domain is defined by ϕ = d, where 0 < c < d < 2π.

In the colatitude direction, there are a finite number of cells, denoted m. In the longitude direction, there are a finite number of cells, denoted n. The angular extent of a cell in each direction is assumed constant, and these angular thicknesses are given by the following expressions:

$$\begin{eqnarray}&&{\rm{\Delta }}\theta =\displaystyle \frac{(b-a)}{m}\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&{\rm{\Delta }}\phi =\displaystyle \frac{(d-c)}{n}.\end{eqnarray} \tag{ 32 }$$

If we are using the staggered grid case, the number of variables in each direction is the number of cell centers; so in this case, solution arrays (without ghost cells) will have dimensions of (m, n). If we are assuming the centered grid case, the number of variables in each direction is the number of cell edges, which is one greater than the number of cell centers. In that case, the solution arrays (without ghost zones) will have dimensions of (m + 1, n + 1).

The quantities a, b, c, d, m, n, Δθ, and Δϕ will retain the meanings defined here throughout the rest of this article.

Given the implicit assumption in FISHPACK of spherical polar coordinates (colatitude and longitude) and the default assumption used nearly universally in solar physics of longitude–latitude array orientation, we are led immediately to the need for frequently transforming back and forth between longitude–latitude and colatitude–longitude array orientations. Thus, an important part of the PDFI_SS software consists of the ability to perform these transpose operations easily and routinely. Detailed discussion of the subroutines that perform these operations is described in Section 3.6 of this article; here we simply present a high-level view of where these transpose operations must be done.

First, if we are using vector magnetogram and Doppler data from HMI on SDO, these data are automatically provided in longitude–latitude orientation. Therefore, a first step is to transpose all the input data (vector magnetograms, velocity maps, LOS unit vectors) to colatitude–longitude orientation. Then, the PDFI solutions are obtained using FISHPACK software in colatitude–longitude order. Essentially all mathematical operations on the data and electric field solutions are done in colatitude–longitude orientation. Finally, because users expect the solutions to be in the same orientation as the HMI data, we must transpose computed results in the other direction to provide the electric field solutions and other related quantities in longitude–latitude order.

For further details, see Section 3.6.

In KFW14, we used a centered grid definition for finite-difference expressions for first derivatives (Equations (14)–(15) in KFW14) in the horizontal directions in our definitions for the curl and gradient. On the other hand, we also used a standard five-point expression for the Laplacian (Equation (16) in KFW14, also used by the Cartesian FISHPACK Poisson solver), which uses centered grid expressions for second derivatives. However, Equation (16) from KFW14 implicitly uses first-derivative finite-difference expressions that are centered half a grid point away from the central point. This means that there is an inconsistency between Equations (16) and (14)–(15) in KFW14. This inconsistency shows up when one uses the centered finite-difference expressions to evaluate $\hat{{\boldsymbol{z}}}\cdot ({\boldsymbol{\nabla }}\times {\boldsymbol{\nabla }}\times \dot{P}\hat{{\boldsymbol{z}}})$ and compares it to $-{{\boldsymbol{\nabla }}}_{h}^{2}\dot{P}$ using Equation (16) in KFW14. In the continuum limit, the two expressions should be identical, but the finite-difference approximations are not equal; the double curl expression using centered finite differences for first derivatives yields an expression like Equation (16) of KFW14, but with the grid points separated by 2Δx and 2Δy. If the grid resolves the solution well, the two different expressions will not differ greatly. This, in fact, is the case with the ANMHD simulation data analyzed in KFW14. But if the solution has structure on the same scale as the grid separation, the double curl expression and the horizontal Laplacian expression can differ significantly, rendering solutions to the PDFI equations quite inaccurate. This problem exists for both the Cartesian and spherical versions of the PDFI equations.

To illustrate the problem quantitatively, we have constructed a solution in spherical coordinates to the Poisson Equation (9) using a centered grid formulation of the finite differences, then evaluated $c{{\boldsymbol{E}}}_{h}$ from the horizontal components of Equation (8), and then evaluated $-\hat{{\boldsymbol{r}}}\cdot ({\boldsymbol{\nabla }}\times c{{\boldsymbol{E}}}_{h})$, which should be equal to the input source term, ${\dot{B}}_{r}$. In this case, the input source term is taken to be a field of random numbers ranging from −0.5 to 0.5, which has significant structure on the scale of the grid. Figure 1 shows an image of the original field of the assumed ${\dot{B}}_{r}$. Figure 2 shows a point-by-point scatterplot of the "recovered" versus original values of ${\dot{B}}_{r}$, showing about half the correct slope, and random errors of roughly 50%. The recovered values of ${\dot{B}}_{r}$ were computed as described above using the centered grid finite-difference expressions. Examination of Figure 2 shows that the centered grid finite-difference expressions do a poor job of describing the correct solution for this test case. The behavior of this test case is similar to what we might expect if the source term has significant levels of pixel-to-pixel noise, which is the case with real magnetogram data in weak-field regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By changing the definition of how finite-difference approximations to spatial derivatives are defined, and where different variables are located within the grid, we can improve this behavior dramatically. In the centered grid case, all variables are colocated at the same grid points. By defining the radial magnetic field and its time derivative to lie at the centers of cells, with the electric field components lying on the edges (or "rails") surrounding the cells, the finite-difference approximations to the derivatives can be made to obey Faraday's law to floating-point roundoff error. Figure 3 shows the analogous scatterplot shown in Figure 2, but using the staggered grid definition described above. The relationship is a straight line. Figure 4 shows the difference between the recovered and original values of ${\dot{B}}_{r}$. Note that the amplitude of the error is multiplied by 1 × 10¹², so that the error term is visible in the scatterplot. These two plots clearly show that solutions for the electric fields in a staggered grid formulation can do a far better job of representing the observed data than can the centered grid formulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These figures motivate the development of our more detailed staggered grid formalism, which is described in Section 3.4.

In three dimensions, Yee (1966) worked out a second-order-accurate finite-difference formulation for Maxwell's equations, pointing out that if one places different variables into different locations within the grid, the governing continuum equations (the curls in Maxwell's equations) become conservative when written down in a finite-difference form. In recent years, "mimetic methods" have been developed, which are higher-order analogs to the Yee grid, in that different variables are defined at different locations within a voxel (such as at interiors, faces, or edges), and with some internal structure in these voxel subdomains allowed. The locations of the variables depend on which integral conservation law is being applied (see, e.g., Candelaresi et al. 2014 and references therein). For our work, the second-order-accurate Yee grid is sufficient. The Yee grid is the basis for the numerical implementation of the MF code described by Cheung & DeRosa (2012).

In PDFI_SS, we have a slightly different situation, where most of the calculations are defined on a subdomain of a spherical surface, so that the domain is 2D, rather than 3D. Nevertheless, the exercise shown in Section 3.3 shows that we want to use the advantages of a staggered grid description of the finite-difference equations, which is inspired by the Yee grid. This is complicated by the fact that in addition to needing curl contributions to the electric field (see Section 2.1), we also need to represent gradient contributions to the electric field (Sections 2.2–2.4). This must be done in such a way that both contributions are colocated along the rails that surround a cell center. We found a way to satisfy these constraints with a specific staggered grid arrangement of the physical and mathematical variables within our 2D domain, summarized below.

First, we define six different grid locations for variables in PDFI_SS, illustrated schematically in Figures 5 and 6. For the moment, we use colatitude–longitude array index order in this discussion. We define the CE grid locations as being the centers of the 2D cells; the CE grid variables are dimensioned (m, n). Next, we define the interior corner grid locations, the "CO" grid, as residing at all the corners, or vertices of the cells, but specifically not including the vertices that lie along the domain edges. Variables lying on the CO grid will have dimension (m − 1, n − 1). Next, we define the "COE" grid, which is also located along corners of cells, but in this case the corners that lie along the domain edges are included. Variables lying on the COE grid will have dimension (m + 1, n + 1). Variables that lie along cell edges that have constant values of ϕ (or longitude) but are at midpoints in θ (or colatitude) lie on the "PE" (phi-edge) locations of the domain. Variables at PE locations have dimension (m, n + 1). Variables that lie along edges with constant θ (or colatitude) but are at midpoints in ϕ lie on the "TE" (theta-edge) grid locations. Variables at TE locations have dimension (m + 1, n). Finally, if we are describing these grid locations, but using longitude–latitude index order, the dimensions of the variables are just the reverse of the dimensions given above.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here are some examples of where different physical and mathematical variables are located using these grid definitions: B_r and $\dot{{B}_{r}}$ are located on the CE grid, B_θ and E_ϕ are located on the TE grid, B_ϕ and E_θ are located along the PE grid, E_r and $\dot{T}$ are located along the COE grid, and $\hat{{\boldsymbol{r}}}\cdot {{\boldsymbol{\nabla }}}_{h}\times {{\boldsymbol{B}}}_{h}$ is located along the CO grid. The scalar potentials defined in the PDFI equations (Sections 2.2–2.4) are located on the COE grid. The poloidal potential $\dot{P}$ is in principle located along the CE grid (with dimensions (m, n)), but we find it convenient to add ghost zones to $\dot{P}$ to implement Neumann (derivative-specified) boundary conditions. When that is done, we refer to this grid as a "CEG" grid (CE plus ghost zones). $\dot{P}$ is on the CEG grid and has dimensions of (m + 2, n + 2).

The placement of variables into the staggered grid locations described above is very similar to the placement used in the "constrained transport" MHD model of Stone & Norman (1992a, 1992b) and the filament construction model of van Ballegooijen (2004). Table 1 contains a list of variables in PDFI_SS and where they reside in terms of these grid locations.

Table 1.  PDFI_SS Variable Locations and Dimensions

Quantity	Grid	Dimension (tp)	Dimension (ll)
Input Data	COE	(m + 1, n + 1)	(n + 1, m + 1)
${B}_{r},{\dot{B}}_{r}$	CE	(m, n)	(n, m)
${B}_{\theta },{\dot{B}}_{\theta }$	TE	(m + 1, n)	(n, m + 1)
${B}_{\phi },{\dot{B}}_{\phi }$	PE	(m, n + 1)	(n + 1, m)
Vθ	TE	(m + 1, n)	(n, m + 1)
Vϕ	PE	(m, n + 1)	(n + 1, m)
VLOS	COE	(m + 1, n + 1)	(n + 1, m + 1)
${\hat{{\boldsymbol{\ell }}}}_{r,\theta ,\phi }$	COE	(m + 1, n + 1)	(n + 1, m + 1)
Er	COE	(m + 1, n + 1)	(n + 1, m + 1)
Eθ	PE	(m, n + 1)	(n + 1, m)
Eϕ	TE	(m + 1, n)	(n, m + 1)
$\dot{P}$	CEG	(m + 2, n + 2)	(n + 2, m + 2)
$\partial \dot{P}/\partial r$	CEG	(m + 2, n + 2)	(n + 2, m + 2)
$\dot{T}$	COE	(m + 1, n + 1)	(n + 1, m + 1)
ψ	COE	(m + 1, n + 1)	(n + 1, m + 1)
$\hat{{\boldsymbol{r}}}\cdot {\boldsymbol{\nabla }}\times {\boldsymbol{E}}$	CE	(m, n)	(n, m)
$\hat{{\boldsymbol{r}}}\cdot {\boldsymbol{\nabla }}\times \dot{{\boldsymbol{B}}}$	CO	(m − 1, n − 1)	(n − 1, m − 1)
${{\boldsymbol{\nabla }}}_{h}\cdot {\dot{{\boldsymbol{B}}}}_{h}$	CE	(m, n)	(n, m)
${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{E}}}_{h}$	CO	(m − 1, n − 1)	(n − 1, m − 1)
Sr	CE	(m, n)	(n, m)
Hm	CE	(m, n)	(n, m)
MCOE	COE	(m + 1, n + 1)	(n + 1, m + 1)
MCO	CO	(m − 1, n − 1)	(n − 1, m − 1)
MTE	TE	(m + 1, n)	(n, m + 1)
MPE	PE	(m, n + 1)	(n + 1, m)
MCE	CE	(m, n)	(n, m)

Download table as:  ASCIITypeset image

It is assumed by the PDFI_SS library that all magnetic field components on input to the library subroutines are in units of Gauss (G). Units of length are determined by the radius of the Sun, which is assumed to be given in kilometers (km). For solar calculations in spherical coordinates, we expect R_⊙ to be 6.96 × 10⁵ km, although in the software the radius of the Sun is an input parameter that can be set by the user. Units of time are assumed to be in seconds (s). Velocities are assumed to be expressed in units of km s⁻¹. For the "working" units of the electric field, the electric field is evaluated as $c{\boldsymbol{E}}$, i.e., the speed of light times the electric field vector, with each component having units of G km s⁻¹. The subroutines that compute the Poynting flux and the Helicity injection rate contribution function are exceptions to this rule and assume that electric field components on input are expressed in units of volts per cm (V cm⁻¹). To convert from G km s⁻¹ to V cm⁻¹, one can simply divide by 1000. To convert units in the opposite direction, one would multiply by 1000.

We mentioned in Section 3.2 that transpose operations from longitude–latitude array orientation to spherical polar coordinates (and the reverse) would need to be done frequently. Now that we have introduced our staggered grid definitions, we will describe in detail how these operations are done, as well as how the interpolation from the input data grid to the staggered grid locations is done. We will discuss how time derivatives are estimated, as well as the calculation of the strong magnetic field masks, designed to decrease the effects of noise from the magnetic field measurements in weak-field regions on the electric field solutions.

The source terms for the PTD contribution to the electric field (Section 2.1) consist of time derivatives of magnetic field components. To estimate these time derivatives from the data, we simply difference the magnetic field values at their staggered grid locations between two adjacent measurement times and divide by the cadence time period, Δt. Thus, if we have magnetic field measurements at times t₀ and t₁ = t₀ + Δt, then our electric field solution will be evaluated at time ${t}_{0}+\tfrac{1}{2}{\rm{\Delta }}t$ and will be assumed to apply over the entire time interval between t₀ and t₁. Furthermore, we assume that the magnetic field values needed to evaluate the other electric field contributions (Sections 2.2–2.4) will be the magnetic field values at ${t}_{0}+\tfrac{1}{2}{\rm{\Delta }}t$, which will be an average of the input values at the two times. Similarly, the other input variables that affect the calculation of ${\boldsymbol{E}}$ will also be an average of the variables at the two adjacent times. If our electric field solutions are conservative and accurately obey Faraday's law, then the computed electric field solutions should correctly evolve ${\boldsymbol{B}}$ from t₀ to t₀ + Δt = t₁ with minimal error.

Thus, for a single time step, the needed input data to evaluate the PDFI solutions are arrays of B_r, B_θ, B_ϕ, V_LOS, V_θ, V_ϕ, ℓ_r, ℓ_θ, and ℓ_ϕ at two adjacent measurement times, for a total of 18 input arrays. Because of the FISHPACK spherical coordinate solution constraints, the data will have to be evaluated using equally spaced colatitude and longitude grid separations, meaning constant spacing in Δθ and Δϕ, referred to as a "Plate Carrée" grid. In the case of HMI data from SDO, this is one of the standard mapping outputs for the magnetic field and Doppler measurements. For CGEM calculations of the electric field supported by the SDO JSOC, the values of Δθ and Δϕ are set to 003 in heliographic coordinates (converted to radians), coinciding closely with an HMI pixel size near disk center. The PDFI_SS software can accommodate values of Δθ and Δϕ that differ, but the FLCT software used upstream of PDFI_SS needs to have these values equal to one another. The JSOC software produces Plate Carrée data with Δθ = Δϕ.

We now briefly digress to describe the relationship between the mathematical coordinate system used by PDFI_SS, with angular domain limits a, b, c, and d, and the standard WCS keywords CRPIX1, CRPIX2, CRVAL1, CRVAL2, CDELT1, and CDELT2 that describe the position of the HMI data on the solar disk (Thompson 2006). We want the ability to concisely relate these two descriptions to each other. The quantities CRPIX1 and CRPIX2 denote longitude and latitude reference pixel locations (the center of the field of view measured from the lower left pixel at (1,1)), CRVAL1 and CRVAL2 denote the longitude and latitude (in degrees) of the reference pixel, and CDELT1 and CDELT2 denote the number of degrees in longitude and latitude between adjacent pixels. From the above description, we expect that CDELT1 and CDELT2 will be equal to 003 pixel^–1. We have written three subroutines, abcd2wcs_ss, wcs2mn_ss, and wcs2abcd_ss, the first of which converts a, b, c, d, m, and n to the WCS keywords CRPIX1, CRPIX2, CRVAL1, CRVAL2, CDELT1, and CDELT2; and in the reverse direction, wcs2mn_ss, which finds m and n from CRPIX1 and CRPIX2 for the COE grid, and wcs2abcd_ss, which converts the keywords CRVAL1, CRVAL2, CDELT1, and CDELT2 to a, b, c, and d. The subroutine abcd2wcs_ss computes the reference pixel locations CRPIX1 and CRPIX2 for all six grid cases, namely, the COE, CO, CE, CEG, TE, and PE grids. These results for the reference pixel locations are returned as six-element arrays, in the order given above.

Returning the discussion to how the input data arrays are processed, the data arrays, in longitude–latitude order, are assumed to be dimensioned (n + 1, m + 1), with all nine input arrays for each of the two times being colocated in space. The parameter a is the colatitude of the northernmost points in these arrays, and the parameter b is the colatitude of the southernmost points in the arrays. The parameters c and d are the leftmost and rightmost longitudes of the input arrays, respectively.

The first task is to transpose all 18 arrays from longitude–latitude to colatitude–longitude (spherical polar coordinates, or θ − ϕ order.) Basically, the transpose operation looks like

$$\begin{eqnarray}&&{A}_{\mathrm{tp}}(i,j)={A}_{{\ell }{\ell }}(j,m-i),\end{eqnarray} \tag{ 33 }$$

where j ∈ [0, n] and i ∈ [0, m], and where A_tp is the array in θ − ϕ index order and A_ℓℓ is the array in longitude–latitude order. Here "tp" in the subscript is meant as a shorthand for "theta-phi," and "ℓℓ" is meant as a shorthand for "longitude–latitude." An exception is for those arrays that represent the latitude components of a vector (like B_lat), in which case when transforming to B_θ the overall sign must also be changed since the unit vectors in latitude and colatitude directions point in opposite directions.

PDFI_SS has several subroutines to perform these transpose operations (and their reverse operations) on the COE grid, namely, brll2tp_ss bhll2tp_ss brtp2ll_ss bhtp2ll_ss.

Here the subroutines starting with "br" perform the transpose operation on scalar fields, while the subroutines starting with "bh" perform the transpose operations on pairs of arrays of the horizontal components of vectors. Subroutines containing the substring "ll2tp" perform the transpose operation going from longitude–latitude order to theta-phi (colatitude–longitude) order, while those with the substring "tp2ll" go in the reverse direction. When going from the input data to colatitude–longitude order, we use the subroutines containing ll2tp within their name. When examining the source code, the expressions will differ slightly from that in Equation (33) to conform with the default FORTRAN index range (where index numbering starts from 1).

Once the input data arrays on the COE grid have been transposed to colatitude–longitude order, we then interpolate the data to their staggered grid locations. B_r is interpolated to the CE grid, B_θ and V_θ to the TE grid, and B_ϕ and V_ϕ to the PE grid. In addition, to evaluate the FLCT electric field contribution, we also need to have B_r and $| {{\boldsymbol{B}}}_{h}|$ interpolated to both the TE and PE grids. Here, we use a simple linear interpolation, as given in these examples for the magnetic field components:

$$\begin{eqnarray}&&\begin{array}{l}{B}_{r}\left(i+\displaystyle \frac{1}{2},j+\displaystyle \frac{1}{2}\right)=\displaystyle \frac{1}{\sin {\theta }_{i}+\sin {\theta }_{i+1}}\\ \ \ \times \ \left(\displaystyle \frac{1}{2}\sin {\theta }_{i}({B}_{r}(i,j)+{B}_{r}(i,j+1))\right.\\ \ \ \left.+\ \displaystyle \frac{1}{2}\sin {\theta }_{i+1}({B}_{r}(i+1,j)+{B}_{r}(i+1,j+1))\right),\end{array}\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{B}_{\theta }\left(i,j+\displaystyle \frac{1}{2}\right)=\displaystyle \frac{1}{2}\left({B}_{\theta }(i,j)+{B}_{\theta }(i,j+1)\right),\end{eqnarray} \tag{ 35 }$$

and

$$\begin{eqnarray}&&{B}_{\phi }\left(i+\displaystyle \frac{1}{2},j\right)=\displaystyle \frac{1}{2}\left({B}_{\phi }(i,j)+{B}_{\phi }(i+1,j)\right).\end{eqnarray} \tag{ 36 }$$

The interpolations from the input data arrays on the COE grid to the staggered grid locations can be accomplished with the subroutines interp_data_ss interp_var_ss.

The linear interpolation is a conservative choice and results in a slight increase in signal-to-noise ratio if there is a high level of pixel-to-pixel noise variation. This interpolation slightly decouples the PDFI_SS electric field from the original input data on the COE grid: the near-perfect reproduction of ${\dot{B}}_{r}$ applies for the interpolations to cell center, but not necessarily for the original input B_r at COE locations.

The Doppler velocity and the LOS unit vector input data arrays are kept at the COE grid locations, so no interpolation of these data arrays is necessary.

In addition to interpolating the input data to the staggered grid locations, we must also construct masks, based on the input data, that reflect regions of the domain where we expect that noise in the magnetic field measurements will make the electric field calculation unreliable. In PDFI_SS the criterion for masks on the magnetic field variables is determined by a threshold on the absolute magnetic field strength, including radial and horizontal components. The mask value is set to unity if the absolute value of the magnetic field in the input data is greater than a chosen threshold for both of the time steps; otherwise, the mask value is set to zero. This calculation is done on the COE grid, after the transpose from longitude–latitude to theta-phi array order. The subroutine that does this is

find_mask_ss.

Subroutine find_mask_ss was originally written assuming that we were using data from three separate time steps, rather than the two time steps we now use. We now simply repeat the array inputs for one of the two time steps, which then results in the correct behavior. For HMI vector magnetogram data, we currently use a threshold value bmin of 250 G. The threshold value is a calling argument to the subroutine and thus can be controlled by the user.

We need to have mask arrays for all the staggered grid locations, not just the COE grid. To get mask arrays for the CE, TE, and PE locations and array sizes, we use a two-step process. First, we use the subroutine interp_var_ss to interpolate the COE mask array to the other staggered grid locations. Those interpolated points where input mask values transition between 0 and 1 will have mask values that are between 0 and 1. The subroutine fix_mask_ss can then be used to set intermediate mask values to either 0 or 1, depending on the value of a "flag" argument to the subroutine, which can be either 0 or 1. Setting flag to 0 is the more conservative choice, while setting flag to 1 is more trusting of the data near the mask edge values.

Once the strong magnetic field mask arrays have been computed, they can be used to multiply the corresponding magnetic field or magnetic field time derivative arrays on input to the subroutines that calculate electric field contributions. This can significantly reduce the impact of magnetogram noise on the electric field solutions in weak magnetic field regions of the domain.

We denote the mask arrays coinciding with different grid locations with the following notation: M^COE is the mask on the COE grid, M^CO denotes the mask on the CO grid, and M^TE, M^PE, and M^CE denote the masks for the TE, PE, and CE grid locations, respectively. The mask arrays are also shown in Table 1.

Once electric field solutions have been computed in spherical polar coordinates, we need the ability to transpose these arrays, as well as the staggered grid magnetic field arrays, back to longitude–latitude order. Because the array sizes are all slightly different depending on which grid is used for a given variable, we have written a series of subroutines designed to perform the transpose operations on our staggered grid, depending on variable type and grid location. There are subroutines to go from theta-phi (colatitude–longitude) order to longitude–latitude order, as well as those that go in the reverse direction. The subroutines that perform the transpose operations on staggered grid locations all have the substring "yee" in their name. As before, subroutine names that include a substring of "tp2ll" transpose the arrays from theta-phi to longitude–latitude array order, while those with "ll2tp" go in the reverse direction. The subroutines are bhyeell2tp_ss bryeell2tp_ss bhyeetp2ll_ss bryeetp2ll_ss ehyeell2tp_ss eryeell2tp_ss ehyeetp2ll_ss eryeetp2ll_ss.

Now that we have established how to generate input data on the staggered grid locations in spherical polar coordinates, we are ready to discuss how to perform vector calculus operations on those data. These operations are used inside the software that evaluates various electric field contributions, and they can also be used to perform other calculations using the electric field solutions.

The following expressions are the continuum vector calculus operations in spherical polar coordinates that are important in evaluating the PDFI equations of Section 2:

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}^{2}{\rm{\Psi }}=\displaystyle \frac{1}{{r}^{2}\sin \theta }\displaystyle \frac{\partial }{\partial \theta }\left(\sin \theta \displaystyle \frac{\partial {\rm{\Psi }}}{\partial \theta }\right)+\displaystyle \frac{1}{{r}^{2}{\sin }^{2}\theta }\displaystyle \frac{{\partial }^{2}{\rm{\Psi }}}{\partial {\phi }^{2}},\end{eqnarray} \tag{ 37 }$$

where Ψ is a scalar function defined in the θ, ϕ domain,

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\times {\rm{\Psi }}\hat{{\boldsymbol{r}}}=\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial {\rm{\Psi }}}{\partial \phi }\hat{{\boldsymbol{\theta }}}-\displaystyle \frac{1}{r}\displaystyle \frac{\partial {\rm{\Psi }}}{\partial \theta }\hat{{\boldsymbol{\phi }}},\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}{\rm{\Psi }}=\displaystyle \frac{1}{r}\displaystyle \frac{\partial {\rm{\Psi }}}{\partial \theta }\hat{{\boldsymbol{\theta }}}+\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial {\rm{\Psi }}}{\partial \phi }\hat{{\boldsymbol{\phi }}},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}\cdot {\boldsymbol{U}}=\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial }{\partial \theta }\left(\sin \theta \ {U}_{\theta }\right)+\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial {U}_{\phi }}{\partial \phi },\end{eqnarray} \tag{ 40 }$$

where ${\boldsymbol{U}}$ is an arbitrary vector field and ${{\boldsymbol{\nabla }}}_{h}\cdot$ represents the divergence in the horizontal directions, and finally

$$\begin{eqnarray}&&\hat{{\boldsymbol{r}}}\cdot {\boldsymbol{\nabla }}\times {\boldsymbol{U}}=\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial }{\partial \theta }\left(\sin \theta \ {U}_{\phi }\right)-\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial {U}_{\theta }}{\partial \phi }.\end{eqnarray} \tag{ 41 }$$

Here U_θ and U_ϕ are the θ and ϕ components of ${\boldsymbol{U}}$.

We now convert these differential expressions to finite-difference expressions, evaluated at various different grid locations in our staggered grid system. Many of these expressions must be evaluated separately depending on where the variables are located, or where we want the expression to be centered. For example, we will need to evaluate Equation (41) at both cell centers (the CE grid) and interior corners (the CO grid), and the exact expressions will differ depending on where the equations are centered.

The subroutines that evaluate the finite-difference expressions corresponding to the above equations are curl_psi_rhat_co_ss curl_psi_rhat_ce_ss gradh_co_ss gradh_ce_ss divh_co_ss divh_ce_ss curlh_co_ss curlh_ce_ss delh2_ce_ss delh2_co_ss.

These subroutines are discussed in more detail below.

In the following equations, we will use this notation to distinguish quantities lying along an edge, versus halfway between edges: an index denoted i or j denotes a location at a θ or ϕ edge, respectively, while an index denoted i + $\tfrac{1}{2}$ or j + $\tfrac{1}{2}$ denotes a location halfway between edges. For example, a quantity defined on the CO grid will have indices i, j, while a quantity defined on the CE grid will have indices $i+\tfrac{1}{2},j+\tfrac{1}{2}$. Similarly, a quantity defined on the TE grid will have the mixed index notation $i,j+\tfrac{1}{2}$, and one along the PE grid would have an index notation of $i+\tfrac{1}{2},j$.

We will start by evaluating Equation (38) assuming that the scalar field Ψ lies on the COE grid. An example of this case is evaluating the curl of the toroidal potential T times $\hat{{\boldsymbol{r}}}$ to find its contribution to the horizontal components of the magnetic field (see Equation (4)). Setting ${\boldsymbol{U}}={\boldsymbol{\nabla }}\times {\rm{\Psi }}\hat{{\boldsymbol{r}}}$, we find

$$\begin{eqnarray}&&{U}_{\theta }\left(i,j+\displaystyle \frac{1}{2}\right)=\displaystyle \frac{{{\rm{\Psi }}}^{\mathrm{COE}}(i,j+1)-{{\rm{\Psi }}}^{\mathrm{COE}}(i,j)}{r\ \sin {\theta }_{i}\ {\rm{\Delta }}\phi }\end{eqnarray} \tag{ 42 }$$

for i ∈ [0, m] and for $j+\tfrac{1}{2}\in \left[\tfrac{1}{2},n-\tfrac{1}{2}\right]$ and

$$\begin{eqnarray}&&{U}_{\phi }\left(i+\displaystyle \frac{1}{2},j\right)=-\ \displaystyle \frac{{{\rm{\Psi }}}^{\mathrm{COE}}(i+1,j)-{{\rm{\Psi }}}^{\mathrm{COE}}(i,j)}{r{\rm{\Delta }}\theta },\end{eqnarray} \tag{ 43 }$$

for $i+\tfrac{1}{2}\in \left[\tfrac{1}{2},m-\tfrac{1}{2}\right]$ and for j ∈ [0, n]. Here, U_θ is dimensioned (m + 1, n) and is defined on the TE grid, while U_ϕ is dimensioned (m, n + 1) and lies on the PE grid. To evaluate Equations (42) and (43), one can use subroutine curl_psi_rhat_co_ss from the PDFI_SS library.

If Ψ is defined on the CEG grid, this has an array size of (m + 2, n + 2), and we then have

$$\begin{eqnarray}&&\begin{array}{l}{U}_{\theta }\left(i+\displaystyle \frac{1}{2},j\right)\\ \ =\ \displaystyle \frac{{{\rm{\Psi }}}^{\mathrm{CEG}}\left(i+\tfrac{1}{2},j+\tfrac{1}{2}\right)-{{\rm{\Psi }}}^{\mathrm{CEG}}\left(i+\tfrac{1}{2},j-\tfrac{1}{2}\right)}{r\ \sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}\phi }\end{array}\end{eqnarray} \tag{ 44 }$$

for $i+\displaystyle \frac{1}{2}\in \left[\displaystyle \frac{1}{2},m-\displaystyle \frac{1}{2}\right]$ and for j ∈ [0, n] and

$$\begin{eqnarray}&&\begin{array}{l}{U}_{\phi }\left(i,j+\tfrac{1}{2}\right)\\ \ =\ -\ \displaystyle \frac{{{\rm{\Psi }}}^{\mathrm{CEG}}\left(i+\tfrac{1}{2},j+\tfrac{1}{2}\right)-{{\rm{\Psi }}}^{\mathrm{CEG}}\left(i-\tfrac{1}{2},j+\tfrac{1}{2}\right)}{r{\rm{\Delta }}\theta }\end{array}\end{eqnarray} \tag{ 45 }$$

for i ∈ [0, m] and for $j+\tfrac{1}{2}\in \left[\tfrac{1}{2},n-\tfrac{1}{2}\right]$.

Here, U_θ is dimensioned (m, n + 1) and lies on the PE grid, while U_ϕ is dimensioned (m + 1, n) and lies on the TE grid. In these expressions, array values of Ψ^CEG at $j+\tfrac{1}{2}=-\tfrac{1}{2}$ and $j+\tfrac{1}{2}=n+\tfrac{1}{2}$ refer to the ghost zone values (with similar expressions for i + $\tfrac{1}{2}$). To evaluate Equations (44) and (45), one can use subroutine curl_psi_rhat_ce_ss.

Note that these expressions require the evaluation of $\sin \theta$ at both edge locations and at cell centers in θ. The need for these geometric factors is ubiquitous in PDFI_SS, so we have written a subroutine sinthta_ss

to precompute these array values before calling many of the vector calculus subroutines.

Turning now to the discretization of Equation (39), namely, evaluating ${\boldsymbol{U}}={{\boldsymbol{\nabla }}}_{h}{\rm{\Psi }}$, the finite-difference expressions for Ψ lying on the COE grid are

$$\begin{eqnarray}&&{U}_{\theta }\left(i+\tfrac{1}{2},j\right)=\displaystyle \frac{{{\rm{\Psi }}}^{\mathrm{COE}}(i+1,j)-{{\rm{\Psi }}}^{\mathrm{COE}}(i,j)}{r{\rm{\Delta }}\theta },\end{eqnarray} \tag{ 46 }$$

where $i+\tfrac{1}{2}\in \left[\tfrac{1}{2},m-\tfrac{1}{2}\right]$ and j ∈ [0, n], and

$$\begin{eqnarray}&&{U}_{\phi }\left(i,j+\tfrac{1}{2}\right)=\displaystyle \frac{{{\rm{\Psi }}}^{\mathrm{COE}}(i,j+1)-{{\rm{\Psi }}}^{\mathrm{COE}}(i,j)}{r\ \sin {\theta }_{i}\ {\rm{\Delta }}\phi },\end{eqnarray} \tag{ 47 }$$

where i ∈ [0, m] and $j+\tfrac{1}{2}\in \left[\tfrac{1}{2},n-\tfrac{1}{2}\right]$. In this case U_θ is dimensioned (m, n + 1) and lies on the PE grid, while U_ϕ is dimensioned (m + 1, n) and lies on the TE grid. These two equations are relevant for computing electric field contributions from the gradients of scalar potentials; subroutine gradh_co_ss can be used to compute these arrays.

When Ψ lies on the CEG grid, we have

$$\begin{eqnarray}&&\begin{array}{l}{U}_{\theta }\left(i,j+\tfrac{1}{2}\right)\\ \ =\ \displaystyle \frac{{{\rm{\Psi }}}^{\mathrm{CEG}}\left(i+\tfrac{1}{2},j+\tfrac{1}{2}\right)-{{\rm{\Psi }}}^{\mathrm{CEG}}\left(i-\tfrac{1}{2},j+\tfrac{1}{2}\right)}{r{\rm{\Delta }}\theta }\end{array}\end{eqnarray} \tag{ 48 }$$

for i ∈ [0, m] and $j+\tfrac{1}{2}\in \left[\tfrac{1}{2},n-\tfrac{1}{2}\right]$ and

$$\begin{eqnarray}&&\begin{array}{l}{U}_{\phi }\left(i+\tfrac{1}{2},j\right)\\ \ =\ \displaystyle \frac{{{\rm{\Psi }}}^{\mathrm{CEG}}\left(i+\tfrac{1}{2},j+\tfrac{1}{2}\right)-{{\rm{\Psi }}}^{\mathrm{CEG}}\left(i+\tfrac{1}{2},j-\tfrac{1}{2}\right)}{r\ \sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}\phi }\end{array}\end{eqnarray} \tag{ 49 }$$

for $i+\tfrac{1}{2}\in \left[\tfrac{1}{2},m-\tfrac{1}{2}\right]$ and j ∈ [0, n]. Here U_θ is dimensioned (m + 1, n) and lies on the TE grid, while U_ϕ is dimensioned (m, n + 1) and lies on the PE grid. To compute U_θ and U_ϕ from Ψ lying on the CEG grid, one can use subroutine gradh_ce_ss.

Moving on now to the discretization of Equation (40), the horizontal divergence of a vector field, this can be evaluated on either the CO grid or the CE grid. Setting ${\rm{\Phi }}={{\boldsymbol{\nabla }}}_{h}\cdot {\boldsymbol{U}}$, we find for the CO grid locations

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}(i,j)=\tfrac{1}{r\sin {\theta }_{i}{\rm{\Delta }}\theta }\\ \ \ \times \ \left({U}_{\theta }\left(i+\tfrac{1}{2},j\right)\sin {\theta }_{i+\tfrac{1}{2}}-{U}_{\theta }\left(i-\tfrac{1}{2},j\right)\sin {\theta }_{i-\tfrac{1}{2}}\right)\\ \ \ +\ \tfrac{1}{r\sin {\theta }_{i}{\rm{\Delta }}\phi }\left({U}_{\phi }\left(i,j+\tfrac{1}{2}\right)-{U}_{\phi }\left(i,j-\tfrac{1}{2}\right)\right),\end{array}\end{eqnarray} \tag{ 50 }$$

where i ∈ [0, m − 2] and j ∈ [0, n − 2]. Here, Φ lies on the CO grid and is dimensioned (m − 1, n − 1), the input array U_θ lies on the PE grid and is dimensioned (m, n + 1), and U_ϕ lies on the TE grid and is dimensioned (m + 1, n). Subroutine divh_co_ss can be used to evaluate the horizontal divergence on the CO grid.

To evaluate the horizontal divergence on the CE grid, we have

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}\left(i+\tfrac{1}{2},j+\tfrac{1}{2}\right)=\tfrac{1}{r\sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}\theta }\\ \ \times \ \left({U}_{\theta }\left(i+1,j+\tfrac{1}{2}\right)\sin {\theta }_{i+1}-{U}_{\theta }\left(i,j+\tfrac{1}{2}\right)\sin {\theta }_{i}\right)\\ \ +\ \tfrac{1}{r\sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}\phi }\times \left({U}_{\phi }\left(i+\tfrac{1}{2},j+1\right)-{U}_{\phi }\left(i+\tfrac{1}{2},j\right)\right),\end{array}\end{eqnarray} \tag{ 51 }$$

where $i+\tfrac{1}{2}\in \left[\tfrac{1}{2},m-\tfrac{1}{2}\right]$ and $j+\tfrac{1}{2}\in \left[\tfrac{1}{2},n-\tfrac{1}{2}\right]$. Here, Φ lies on the CE grid and has dimensions of (m, n). The arrays U_θ and U_ϕ lie on the TE and PE grids, respectively, with dimensions (m + 1, n) and (m, n + 1). Subroutine divh_ce_ss can be used to evaluate the horizontal divergence on the CE grid.

Finally, we address the discretization of Equation (41). Setting ${\rm{\Phi }}=\hat{{\boldsymbol{r}}}\cdot {\boldsymbol{\nabla }}\times {\boldsymbol{U}}$, we can evaluate Φ on the CO or the CE grid. For the CO grid, we have

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}(i,j)=\tfrac{1}{r\sin {\theta }_{i}{\rm{\Delta }}\theta }\\ \ \times \ \left({U}_{\phi }\left(i+\tfrac{1}{2},j\right)\sin {\theta }_{i+\tfrac{1}{2}}-{U}_{\phi }\left(i-\tfrac{1}{2},j\right)\sin {\theta }_{i-\tfrac{1}{2}}\right)\\ \ -\ \tfrac{1}{r\sin {\theta }_{i}{\rm{\Delta }}\phi }\left({U}_{\theta }\left(i,j+\tfrac{1}{2}\right)-{U}_{\theta }\left(i,j-\tfrac{1}{2}\right)\right)\end{array}\end{eqnarray} \tag{ 52 }$$

for i ∈ [0, m − 2] and j ∈ [0, n − 2], with the output dimensioned (m − 1, n − 1). The input array U_θ is dimensioned (m + 1, n) and lies on the TE grid, and U_ϕ is dimensioned (m, n + 1) and lies on the PE grid. Subroutine curlh_co_ss can be used to evaluate Equation (52).

Evaluating Φ on the CE grid, we have

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}\left(i+\tfrac{1}{2},j+\tfrac{1}{2}\right)=\tfrac{1}{r\sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}\theta }\\ \ \times \ \left({U}_{\phi }\left(i+1,j+\tfrac{1}{2}\right)\sin {\theta }_{i+1}-{U}_{\phi }\left(i,j+\tfrac{1}{2}\right)\sin {\theta }_{i}\right)\\ \ -\ \tfrac{1}{r\sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}\phi }\times \left({U}_{\theta }\left(i+\tfrac{1}{2},j+1\right)-{U}_{\theta }\left(i+\tfrac{1}{2},j\right)\right)\end{array}\end{eqnarray} \tag{ 53 }$$

for $i+\tfrac{1}{2}\in \left[\tfrac{1}{2},m-\tfrac{1}{2}\right]$ and $j+\tfrac{1}{2}\in \left[\tfrac{1}{2},n-\tfrac{1}{2}\right]$, with the output dimensioned (m, n). The input array U_θ is dimensioned (m, n + 1) and lies on the PE grid, and U_ϕ is dimensioned (m + 1, n) and lies on the TE grid. Subroutine curlh_ce_ss can be used to evaluate Equation (53).

Note that we have not written down a finite-difference expression for the horizontal Laplacian, Equation (37). The finite-difference expression can be found in FISHPACK documentation, in Swarztrauber & Sweet (1975), and in Section 6. The horizontal Laplacian of Ψ, when Ψ is on the COE grid, can be computed with subroutine delh2_co_ss.

This subroutine just uses a call to gradh_co_ss, followed by a call to divh_co_ss. The result is the Laplacian of Ψ evaluated on the CO grid. Similarly, the horizontal Laplacian of P, which lies on the CEG grid, can be computed with subroutine delh2_ce_ss.

This just uses a call to gradh_ce_ss, followed by a call to divh_ce_ss. The result is the Laplacian of P evaluated on the CE grid. Solutions of Poisson's equation found from FISHPACK can be tested with these subroutines in PDFI_SS and then the results compared with the source terms used on input to the Poisson equation. In all cases tested, we find agreement that is close to floating-point roundoff error.

When closely examining the source code for the above subroutines, one may find that the index range differs from the ranges mentioned above; this is done to adhere to default array index ranges in FORTRAN. However, the array dimensions and grid locations will be consistent with those described above.

When normal derivative (Neumann) boundary conditions are input into FISHPACK subroutines, the solution is computed only on the "active" part of the grid. To ensure that the finite-difference expressions given in Section 3.7 are consistent with the boundary conditions, we add extra "ghost zones" to the solutions that ensure that the boundary conditions are obeyed within these expressions. The clearest example of how this is done are the two extra ghost zones added in θ and in ϕ to get values of $\dot{P}$ on the CEG grid from the solution returned by FISHPACK on the CE grid. In this case, we are using subroutine HSTSSP, for which the first active point in ϕ is at $c+\tfrac{1}{2}{\rm{\Delta }}\phi$, while the boundary is at ϕ = c. Thus, for $\dot{P}$, we add a row of m cells centered at $\phi =c-\tfrac{1}{2}{\rm{\Delta }}\phi$, for which

$$\begin{eqnarray}&&\begin{array}{l}\dot{P}\left({\theta }_{i+\tfrac{1}{2}},c-\tfrac{1}{2}{\rm{\Delta }}\phi \right)=\dot{P}({\theta }_{i+\tfrac{1}{2}},c+\tfrac{1}{2}{\rm{\Delta }}\phi )\\ \ \ -\ {\rm{\Delta }}\phi {\left.\left(\displaystyle \frac{\partial \dot{P}}{\partial \phi }\right)\right|}_{\phi =c}({\theta }_{i+\tfrac{1}{2}}),\end{array}\end{eqnarray} \tag{ 54 }$$

where ${\left.\left(\displaystyle \frac{\partial \dot{P}}{\partial \phi }\right)\right|}_{\phi =c}({\theta }_{i+\tfrac{1}{2}})$ is the derivative of $\dot{P}$ specified at ϕ = c for the m points along the left boundary in the call to HSTSSP. There is a similar operation to determine ghost cell values for the other three domain boundaries. See Section 3.9.1 for a discussion of the Neumann boundary conditions for $\dot{P}$ and $\partial \dot{P}/\partial r$.

Because the CEG grid is dimensioned (m + 2, n + 2), there are four unused "corner" values for the $\dot{P}$ array at CEG locations. These four values could be set to any value, but for display purposes we set the corner values to be the average of the two closest neighbor points, so that the $\dot{P}$ array can be viewed as a continuous function when visualized.

In this section, we describe in detail how the four different electric field contributions to the PDFI electric field can be computed with various subroutines within PDFI_SS. We first describe the calculation of the PTD (inductive) electric field. Next, we discuss the software for performing the "iterative" method, necessary to evaluate the Doppler and ideal contributions to the electric field. Following this, we discuss the calculation of the Doppler electric field, followed by the contribution from FLCT (or other "optical flow" derived horizontal velocities) to the electric field. Finally, we discuss the calculation of the "ideal" contribution to ${\boldsymbol{E}}$.

3.9.1. Numerical Solution for E^P in PDFI_SS

The calculation of ${{\boldsymbol{E}}}^{P}$ (the PTD, or inductive contribution to ${\boldsymbol{E}}$) depends exclusively on time derivatives of ${\boldsymbol{B}}$. In Section 3.6 we described the estimation of time derivatives in terms of simple differences in the magnetic field components that take place between two successive times in an assumed time cadence. Once the data have been interpolated to the staggered grid locations, we will simply define ${\dot{B}}_{r}$ from the data as

$$\begin{eqnarray}&&{\dot{B}}_{r}\left({t}_{0}+\tfrac{1}{2}{\rm{\Delta }}t\right)=\displaystyle \frac{{B}_{r}({t}_{0}+{\rm{\Delta }}t)-{B}_{r}({t}_{0})}{{\rm{\Delta }}t},\end{eqnarray} \tag{ 55 }$$

where Δt is the assumed time separation of the cadence. Here, the subtraction is understood to apply in a whole array sense, i.e., to all the points in the CE grid locations for both arrays of B_r, with similar expressions for ${\dot{B}}_{\theta }$ and ${\dot{B}}_{\phi }$, defined for their own array sizes and locations (see Table 1). If no masking for weak magnetic field regions is desired, this definition for the time derivatives can then be used to derive the PTD solution. However, we have found that in weak-field regions the PTD solution ${{\boldsymbol{E}}}^{P}$ can be strongly affected by noise. One can suppress much of this noise by multiplying ${\dot{B}}_{r}$, ${\dot{B}}_{\theta }$, and ${\dot{B}}_{\phi }$ by their respective strong magnetic field mask arrays M^CE, M^TE, and M^PE as defined in Section 3.6. Whether the input is masked or not, the resulting electric field contributions for ${{\boldsymbol{E}}}^{P}$ are computed by two subroutines, called in succession: ptdsolve_ss e_ptd_ss.

Subroutine ptdsolve_ss solves the three Poisson Equations (9)–(11) for $\dot{P}$, $\dot{T}$, and $(\partial \dot{P}/\partial r)$. Once these have been computed, subroutine e_ptd_ss uses $\dot{P}$ and $\dot{T}$ to compute the three components of ${{\boldsymbol{E}}}^{P}$. If desired, one can then use a third subroutine dehdr_ss to use $(\partial \dot{P}/\partial r)$ to compute the radial derivatives of ${{\boldsymbol{E}}}_{h}^{P}$.

There are a lot of assumptions made about boundary conditions and equation centering in subroutine ptdsolve_ss that need to be mentioned. First, the Poisson equations for $\dot{P}$ and $\partial \dot{P}/\partial r$ are centered on the CE grid, since their source terms ${\dot{B}}_{r}$ and ${{\boldsymbol{\nabla }}}_{h}\cdot {\dot{{\boldsymbol{B}}}}_{h}$ are both centered on the CE grid, meaning that FISHPACK subroutine HSTSSP will be used for the solution. Second, the Poisson equation for $\dot{T}$ is centered on the CO grid, since its source term $-\hat{{\boldsymbol{r}}}\cdot {{\boldsymbol{\nabla }}}_{h}\times {\dot{{\boldsymbol{ \mathcal B }}}}_{h}$ is located on the CO grid. This means that FISHPACK subroutine HWSSSP will be used for its solution.

A physically meaningful boundary condition for ${{\boldsymbol{E}}}_{h}^{P}$ is to specify the electric field component tangential to the domain boundary. The tangential electric field is directly related to the derivatives of $\dot{P}$ in the directions normal to the boundary. During the development phase of the PDFI_SS software, we initially assumed zero tangential electric field along the domain boundary. However, this assumption was in conflict with the actual HMI data: for many regions, the net radial magnetic flux is not balanced and can increase or decrease in time, which was inconsistent with our assumptions. Therefore, we have modified our assumed boundary conditions for $\dot{P}$, and we now specify a boundary condition on E_t (the electric field tangential to the boundary) that is consistent with the observed increase or decrease in the net radial magnetic flux.

We first evaluate the time rate of change of the net radial magnetic flux in the domain:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Phi }}}{\partial t}={\rm{\Delta }}\phi {\rm{\Delta }}\theta \ {r}^{2}\sum _{i+\tfrac{1}{2},j+\tfrac{1}{2}}{\dot{B}}_{r}\left(i+\tfrac{1}{2},j+\tfrac{1}{2}\right)\ \sin {\theta }_{i+\tfrac{1}{2}},\end{eqnarray} \tag{ 56 }$$

where the sum is over the cell centers of the domain (i.e., over the CE grid). Next, we evaluate the perimeter length of the domain along the north and south edges:

$$\begin{eqnarray}&&{L}_{\mathrm{perim}}=r\left[(d-c)(\sin \,a+\sin \,b)\right].\end{eqnarray} \tag{ 57 }$$

Assuming that the tangential electric field is zero along the left and right edges of the domain, we can then use Stokes's theorem to integrate Faraday's law over the domain to find a constant amplitude of the electric field on the north and south edges, ${{cE}}_{\mathrm{perim}}$:

$$\begin{eqnarray}&&{{cE}}_{\mathrm{perim}}=-\ \displaystyle \frac{\partial {\rm{\Phi }}}{\partial t}{/L}_{\mathrm{perim}},\end{eqnarray} \tag{ 58 }$$

where it is understood that ${{cE}}_{\mathrm{perim}}$ along the domain edges points in the counterclockwise direction if it is positive. The minus sign in Equation (58) comes from the minus sign in Faraday's law. We assign the normal derivatives of $\dot{P}$ to either zero (left and right boundaries) or ${{cE}}_{\mathrm{perim}}$ (north and south boundaries) in FISHPACK subroutine HSTSSP after accounting for the spherical geometry factors (see Equation (38)) and the sign of ${{cE}}_{\mathrm{perim}}$ relative to the ϕ unit vector along the domain boundaries:

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial \dot{P}}{\partial \theta }\right|}_{\theta =a}\left(j+\tfrac{1}{2}\right)=-\,{{cE}}_{\mathrm{perim}}\,r,\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial \dot{P}}{\partial \theta }\right|}_{\theta =b}\left(j+\tfrac{1}{2}\right)={{cE}}_{\mathrm{perim}}\,r,\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial \dot{P}}{\partial \phi }\right|}_{\phi =c}\left(i+\tfrac{1}{2}\right)=0,\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial \dot{P}}{\partial \phi }\right|}_{\phi =d}\left(i+\tfrac{1}{2}\right)=0.\end{eqnarray} \tag{ 62 }$$

We still have, in the PDFI_SS library, the subroutine that assumes zero tangential electric field along boundary edges if the user wants to use this assumption. That subroutine is named ptdsolve_eb0_ss.

Assuming that the boundary is far away from the most rapid evolution of the data, and that there is no significant change in the net radial current density in the domain, we set the electric field ${E}_{r}^{P}$ to zero at the boundary. Since this is proportional to $\dot{T}$, we use homogeneous Dirichlet boundary conditions for $\dot{T}$ ($\dot{T}$ is set to zero at θ = a, θ = b, ϕ = c, and ϕ = d).

The physical boundary condition for $\partial \dot{P}/\partial r$ is to specify the component of ${\dot{{\boldsymbol{B}}}}_{h}$ normal to the domain boundary, since the gradient of $\partial \dot{P}/\partial r$ is equal to ${\dot{{\boldsymbol{B}}}}_{h}$ when $\dot{T}$ is zero on the boundary, as can be seen from Equation (5). This leads to these Neumann boundary conditions (see Equation (39)) for this Poisson equation:

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial }{\partial \theta }\right|}_{\theta =a}\displaystyle \frac{\partial \dot{P}}{\partial r}=r{\dot{B}}_{\theta }\left(\theta =a,{\phi }_{j+\tfrac{1}{2}}\right),\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial }{\partial \theta }\right|}_{\theta =b}\displaystyle \frac{\partial \dot{P}}{\partial r}=r{\dot{B}}_{\theta }\left(\theta =b,{\phi }_{j+\tfrac{1}{2}}\right),\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial }{\partial \phi }\right|}_{\phi =c}\displaystyle \frac{\partial \dot{P}}{\partial r}=r\sin {\theta }_{i+\tfrac{1}{2}}\,{\dot{B}}_{\phi }\left(\phi =c,{\theta }_{i+\tfrac{1}{2}}\right),\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial }{\partial \phi }\right|}_{\phi =d}\displaystyle \frac{\partial \dot{P}}{\partial r}=r\sin {\theta }_{i+\tfrac{1}{2}}\,{\dot{B}}_{\phi }\left(\phi =d,{\theta }_{i+\tfrac{1}{2}}\right).\end{eqnarray} \tag{ 66 }$$

Note that the use of a staggered grid decouples the boundary conditions for $\dot{T}$ and $\partial \dot{P}/\partial r$ that existed for the centered grid, as described in Fisher et al. (2010) and KFW14.

The solutions for $\dot{P}$ and $\partial \dot{P}/\partial r$ are returned by ptdsolve_ss on the CEG grid, which means that there is an extra row of ghost zones returned along each of the four sides of the boundary (see discussion above in Section 3.8).

Subroutine ptdsolve_ss can also be used to find the poloidal potential P, the toroidal potential T, and the radial derivative of the poloidal potential ∂P/∂r, if, instead of inputting the time derivatives of the magnetic field components, one inputs the magnetic field components themselves. The vector potential ${\boldsymbol{A}}$ can be computed from subroutine e_ptd_ss, but a minus sign must be applied to the output.

3.9.2. Implementation of the "Iterative" Method in PDFI_SS

The "iterative" method for finding a scalar potential whose gradient is designed to closely match a given vector field was developed by coauthor Brian Welsch and initially described in Section 3.2 of Fisher et al. (2010). In KFW14, the method was used to derive the scalar potential representing the Doppler electric field, as well as the ideal electric field contribution, which has the goal of setting ${\boldsymbol{E}}\cdot {\boldsymbol{B}}=0$. Applying the technique was relatively simple because, using the centered grid formalism, ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$ were colocated. In PDFI_SS, however, the various components of ${\boldsymbol{B}}$ and ${\boldsymbol{E}}$ all lie in different locations, making the algorithm less straightforward to implement directly. We now describe our current solution to the problem.

We describe the iterative technique in terms of generating the ideal component of the electric field (rather than the Doppler contribution), because we think that the logic is easier to follow in that case, but the solution method applies equally well to both the Doppler and ideal electric field contributions.

In PDFI_SS the subroutine that computes solutions using the iterative method is called relax_psi_3d_ss.

Our approach is to perform the iterative procedure on a temporary, centered grid, coinciding with the CO grid, but computed using centered grid finite-difference expressions, centered on cell vertices. On input to relax_psi_3d_ss, we need values of ${{\boldsymbol{E}}}^{\mathrm{PDF}}$ and ${\boldsymbol{B}}$ lying on the CO grid. To construct values of ${{\boldsymbol{E}}}^{\mathrm{PDF}}$ on that grid from their native staggered grid locations, we can use linear interpolation, computed from subroutine interp_eh_ss, to interpolate the horizontal electric field contributions to the CO grid. The solutions for E_r are already given on the COE grid, from which the CO grid is just a subset. For magnetic field values on the CO grid, we note that the original input magnetogram data were initially provided on the COE grid, so no interpolation to CO is necessary. Subroutine relax_psi_3d_ss assumes a domain size that is smaller in extent than our overall domain boundary, with its boundaries given by a' = a + Δθ, b' = b − Δθ, c' = c + Δϕ, and d' = d − Δϕ. Internal to relax_psi_3d_ss, boundary conditions assumed at a', b', c', and d' are that the normal derivatives of the scalar potential ψ are zero (homogeneous Neumann boundary conditions). This is implemented by assigning ghost zone values for ψ that enforce this boundary condition under the centered grid assumption, resulting in the array for ψ, including the ghost zones, lying on the COE grid. Once the iterative procedure has been completed, but before exiting the subroutine, a final set of boundary conditions are applied using ghost zones for ψ implemented on the edges of the COE grid, in which the homogeneous Neumann boundary conditions for E_n (the normal components of ${{\boldsymbol{\nabla }}}_{h}\psi$) are applied using the staggered grid formalism, at a boundary halfway between the edges of the CO and COE grids, at $a+\tfrac{1}{2}{\rm{\Delta }}\theta$, $b-\tfrac{1}{2}{\rm{\Delta }}\theta$, $c+\tfrac{1}{2}{\rm{\Delta }}\phi$, and $d-\tfrac{1}{2}{\rm{\Delta }}\phi$. This results in slight changes to the ghost zone values of ψ located at the edges of the COE grid, compared to their values computed as ghost zones using the centered grid formalism. We then also set boundary conditions for ∂ψ/∂r = 0 at the regular domain boundaries, a, b, c, and d. After exiting relax_psi_3d_ss, gradients of ψ are then computed using the staggered grid formalism of Section 3.7, rather than using the centered grid description that is used internally within that subroutine.

We now summarize the details of how ψ and ∂ψ/∂r are computed in the subroutine, combining material originally in Section 3.2 of Fisher et al. (2010) and Section 2.2 of KFW14 and using spherical coordinates. The procedures described here are slight updates of the original iterative method, as the code has evolved from its original formulation:

Step 1:

Decompose ${\boldsymbol{\nabla }}\psi$ as

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\psi ={s}_{1}\hat{{\boldsymbol{b}}}+{s}_{2}\ \hat{{\boldsymbol{r}}}\times \hat{{\boldsymbol{b}}}+{s}_{3}\ \hat{{\boldsymbol{b}}}\times \left(\hat{{\boldsymbol{r}}}\times \hat{{\boldsymbol{b}}}\right),\end{eqnarray} \tag{ 67 }$$

where s₁, s₂, and s₃ are understood to be functions of θ and ϕ. Here $\hat{{\boldsymbol{b}}}$ is the unit vector pointing in the same direction as ${\boldsymbol{B}}$. Quantities ${\hat{{\boldsymbol{b}}}}_{h}$ and b_r, used below, denote the horizontal and radial components of $\hat{{\boldsymbol{b}}}$, respectively, and ${b}_{h}^{2}$ is the square of the amplitude of ${\hat{{\boldsymbol{b}}}}_{h}$.

Step 2:

Set

$$\begin{eqnarray}&&{s}_{1}={{\boldsymbol{E}}}^{\mathrm{PDF}}\cdot \hat{{\boldsymbol{b}}},\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&{s}_{2}\,=\,0,\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}&&{s}_{3}\,=\,0,\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\partial \psi }{\partial r}\right)}_{0}={s}_{1}\ {b}_{r},\end{eqnarray} \tag{ 71 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}{\psi }_{0}={s}_{1}{\hat{{\boldsymbol{b}}}}_{h}.\end{eqnarray} \tag{ 72 }$$

The quantities ψ and ∂ψ/∂r are both regarded as functions of θ and ϕ, and subscript 0 denotes the "zeroth" iterative approximation to ψ. Equation (72) should result in cancellation of the component of ${{\boldsymbol{E}}}^{\mathrm{PDF}}$ parallel to ${\boldsymbol{B}}$ if $-{\boldsymbol{\nabla }}\psi$ is added to it. To obtain the guess for ψ₀, we can take the divergence of Equation (72):

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}^{2}{\psi }_{0}={{\boldsymbol{\nabla }}}_{h}\cdot ({s}_{1}\ {\hat{{\boldsymbol{b}}}}_{h}).\end{eqnarray} \tag{ 73 }$$

The horizontal divergence operation on the right-hand side of Equation (73) is computed using a centered grid formalism using subroutine divh_sc.

We can solve this Poisson equation for ψ₀ using FISHPACK subroutine HWSSSP, subject to homogeneous Neumann boundary conditions at the "primed" values for a, b, c, and d noted above. The quantity s₁ will be held fixed throughout the iterative sequence. Once we have a solution for ψ₀, we can evaluate ${{\boldsymbol{\nabla }}}_{h}{\psi }_{0}$ using the centered grid subroutine gradh_sc.

Step 3 (the beginning of the iterative sequence):

Given the current guess for ψ and ${\boldsymbol{\nabla }}\psi$, evaluate

$$\begin{eqnarray}&&{s}_{2}=\displaystyle \frac{\hat{{\boldsymbol{r}}}\cdot \left({\hat{{\boldsymbol{b}}}}_{h}\times {{\boldsymbol{\nabla }}}_{h}\psi \right)}{{b}_{h}^{2}}\end{eqnarray} \tag{ 74 }$$

and

$$\begin{eqnarray}&&{s}_{3}=\displaystyle \frac{\partial \psi }{\partial r}-\displaystyle \frac{{b}_{r}({{\boldsymbol{\nabla }}}_{h}\psi \cdot {\hat{{\boldsymbol{b}}}}_{h})}{{b}_{h}^{2}}.\end{eqnarray} \tag{ 75 }$$

Equations (74) and (75) can be derived by dotting both sides of Equation (67) with the vectors $\hat{{\boldsymbol{r}}}\times \hat{{\boldsymbol{b}}}$ and $\hat{{\boldsymbol{b}}}\times (\hat{{\boldsymbol{r}}}\times \hat{{\boldsymbol{b}}})$, respectively.

Step 4:

Given s₂ and s₃ from Step 3, evaluate the horizontal divergence of Equation (67)

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}^{2}\psi ={{\boldsymbol{\nabla }}}_{h}\cdot \left({s}_{1}{\hat{{\boldsymbol{b}}}}_{h}+{s}_{2}\ \hat{{\boldsymbol{r}}}\times \hat{{\boldsymbol{b}}}-{s}_{3}\ {b}_{r}{\hat{{\boldsymbol{b}}}}_{h}\right)\end{eqnarray} \tag{ 76 }$$

and then solve this Poisson equation for the next iterative solution for ψ. The update for ∂ψ/∂r is given by the following expression:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \psi }{\partial r}={s}_{1}{b}_{r}+{s}_{3}{b}_{h}^{2}.\end{eqnarray} \tag{ 77 }$$

Step 5:

If the number of iterations is less than the maximum number (current default value is 25), go back to Step 3. If the maximum number is exceeded, then exit the iteration procedure. The resulting arrays of ψ and ∂ψ/∂r on the CO grid are the final arrays for the iterative solution. We note again that the ghost zone values, located on the edges of the COE grid, are adjusted from the values computed from the centered grid finite differences to make them consistent with the use of the staggered grid finite differences.

We now remark on several properties of the iterative technique described above. First, as noted in Section 3.2 of Fisher et al. (2010) and Section 2.2 of KFW14, the mathematical problem that the iterative method is designed to solve has no unique solution; nevertheless, this method appears to find a unique solution, meaning that most likely the method imposes other hidden constraints, which makes it behave like a unique solution. See Section 2.2 of KFW14 for further discussion.

Second, the iterative improvement in the solution is rapid for the first few iterations, and then improvement slows dramatically. We have found that implementing an error convergence criterion, as originally suggested in Fisher et al. (2010), has proven unreliable and difficult. We follow the suggestion of KFW14 that setting a fixed number of iterations is a better implementation. We adopt the suggestion from KFW14 of 25 iterations. Experiments we have done have shown that using much larger numbers of iterations (100 or 1000, for example) actually makes the solution worse.

Third, in contrast to the error in, e.g., Faraday's law, which is near the floating-point roundoff error, the ability of the iterative technique to exactly cancel the component of ${\boldsymbol{E}}$ in the direction of ${\boldsymbol{B}}$ is much less precise. We find in PDFI_SS that the typical angle between ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$ once $-{\boldsymbol{\nabla }}\psi$ has been added to ${{\boldsymbol{E}}}^{\mathrm{PDF}}$ for a number of different test cases is within 2° of 90°. Histograms of the cosine of the angle between the two vectors are sharply peaked at zero (see Figure 7), but with significant tails in the distribution. Examination of where the outlier points are located shows a concentration near the boundaries of the strong magnetic field mask, suggesting that the iterative procedure has its worst performance near the mask boundaries.

Zoom In Zoom Out Reset image size

Fourth, the iterative technique is sensitive to noisy input data in ${{\boldsymbol{E}}}^{\mathrm{PDF}}$. For example, if we compute a solution to ${{\boldsymbol{E}}}^{P}$ using unmasked magnetic field time derivatives and then try to find an electric field contribution that attempts to make ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$ perpendicular, the iterative technique can diverge, rather than converge. For this reason, we strongly recommend using masking for the magnetic field time derivative arrays on input to ptdsolve_ss, if the iterative method will then be used to compute the "ideal" contribution.

Finally, we note again that once the solutions for ψ and ∂ψ/∂r are obtained after exiting relax_psi_3d_ss, the resulting contribution from $-{\boldsymbol{\nabla }}\psi$ to ${\boldsymbol{E}}$ is evaluated at the staggered grid locations. This has the advantage of producing a curl-free contribution to ${\boldsymbol{E}}$, but at the cost of losing the direct connection to the angle between ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$, since the vectors are no longer colocated. To evaluate the angle between ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$, the electric fields must once again be interpolated to the CO grid before this comparison can be done. Subroutine angle_be_ss can be used to interpolate electric field components to the CO grid and then evaluate the cosine of the angle between the ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$ vectors on the CO grid.

3.9.3. Implementing the Doppler Electric Field Contribution

3.9.4. Computing the FLCT Contribution

3.9.5. Computing the Ideal Contribution to the PDFI Electric Field

Once the PDFI electric field has been computed, there are a number of other useful quantities that can be computed with it, including the radial component of the Poynting flux, as well as the contribution function to the relative helicity injection rate.

These quantities are computed by the subroutines sr_ss and hm_ss.

The subroutine sr_ss takes as input the horizontal components of both the electric field and magnetic field in their staggered grid locations on the TE and PE grids and computes the radial component of the Poynting flux at cell centers (the CE grid). While most of the PDFI_SS software assumes that electric fields are computed as $c{\boldsymbol{E}}$, in units of G km s⁻¹, subroutine sr_ss assumes that the input electric fields do not include the factor of c and are given in units of V cm⁻¹. To convert from $c{\boldsymbol{E}}$ in units of G km s⁻¹ to ${\boldsymbol{E}}$ in units of V cm⁻¹, one can simply divide by a factor of 1000.

We find that in the weak-field regions, the Poynting flux can be quite unreliable, so we recommend that after output from subroutine sr_ss the resulting Poynting flux array should be multiplied by the strong magnetic field mask for the CE grid, M^CE. If the strong-field masks have been used to compute the electric field contributions, then, for consistency, the masks should also be applied on either the input horizontal magnetic fields or the Poynting flux output (which is what we do). The assumed units on output from sr_ss for the Poynting flux are erg cm⁻² s⁻¹.

To compute the total magnetic energy input rate from the radial component of the Poynting flux, one can use subroutine srtot_ss

to integrate the radial Poynting flux contribution over area. The output is a single value, computed in units of erg s⁻¹.

The subroutine hm_ss is used to compute the contribution function for the relative helicity injection rate. On input, it uses the poloidal potential P computed from the radial component of the magnetic field using subroutine ptdsolve_ss and the horizontal components of ${\boldsymbol{E}}$ from the PDFI solution. We typically compute P using arrays of B_r that are multiplied by the strong-field mask M^CE before ptdsolve_ss is called, so that the vector potential for the potential magnetic field ${{\boldsymbol{A}}}_{P}$ does not contain contributions from the weak-field regions. The vector potential ${{\boldsymbol{A}}}_{P}$ is computed from P using subroutine curl_psi_rhat_ce_ss within hm_ss. Values of ${{\boldsymbol{A}}}_{P}$ and ${\boldsymbol{E}}$ components are then interpolated to the CE grid, and then the quantity $\hat{{\boldsymbol{r}}}\cdot 2c{{\boldsymbol{E}}}_{h}\times {{\boldsymbol{A}}}_{P}$ is computed on the CE grid. The input units for the horizontal electric field are assumed to be in units of V cm⁻¹, and the units for P are assumed to be G km². The output array is computed in units of Mx² cm⁻² s⁻¹.

One can integrate the contribution function to get a total relative helicity injection rate by calling subroutine hmtot_ss using the output from hm_ss as input. The output is a single value, given in units of Mx² s⁻¹.

The preceding parts of this section have described in detail how the HMI input data are transposed, interpolated, and then used to compute the various contributions to the PDFI electric field, and how that can then be used to create maps of the Poynting flux and the contribution function for the relative helicity injection rate. We have written a subroutine in the PDFI_SS library, pdfi_wrapper4jsoc_ss, that combines all of these pieces together. The SDO JSOC calls this subroutine to compute the electric field and related variables to create the CGEM data series that is distributed by the JSOC. The subroutine is also useful, in that it can serve as a template for a customized calculation of the electric field, allowing a user to eliminate unwanted terms, experiment with various masking strategies for input data, or experiment with new electric field contributions.

The list of major tasks performed by pdfi_wrapper4jsoc_ss, along with the subroutines used for these tasks, is given in order below:

• 1.  
  Transpose the 18 input arrays from longitude–latitude order to colatitude–longitude (theta-phi) order (brll2tp_ss, bhll2tp_ss).
• 2.  
  Convert Doppler velocities from m s⁻¹ to km s⁻¹ and change sign convention to positive for upflows.
• 3.  
  Compute strong-field mask arrays for staggered grid locations from arrays of input magnetic field arrays on the COE grid (find_mask_ss, fix_mask_ss).
• 4.  
  Interpolate input data to staggered grid locations (interp_data_ss, interp_var_ss).
• 5.  
  Compute $\sin \theta$ arrays at colatitude edges and cell centers (sinthta_ss).
• 6.  
  Compute $\dot{P}$, $\dot{T}$, $\partial \dot{P}/\partial r$ and P, T, ∂P/∂r (ptdsolve_ss).
• 7.  
  Compute PTD electric field contribution ${{\boldsymbol{E}}}^{P}$ (e_ptd_ss).
• 8.  
  Compute Doppler electric field contribution ${{\boldsymbol{E}}}^{D}$ (e_doppler_ss, relax_psi_3d_ss).
• 9.  
  Compute FLCT electric field contribution ${{\boldsymbol{E}}}^{F}$ (e_flct_ss).
• 10.  
  Compute ideal electric field contribution ${{\boldsymbol{E}}}^{I}$ (e_ideal_ss, relax_psi_3d_ss).
• 11.  
  Add all four contributions for ${{\boldsymbol{E}}}^{\mathrm{PDFI}}$, convert units to V cm⁻¹.
• 12.  
  Compute radial derivatives of horizontal components of electric field (dehdr_ss).
• 13.  
  Compute Poynting flux and its area integral (sr_ss, srtot_ss).
• 14.  
  Compute contribution function for helicity injection and its area integral (hm_ss, hmtot_ss).
• 15.  
  Transpose all output arrays to longitude–latitude array order (bhyeetp2ll_ss,bryeetp2ll_ss, ehyeetp2ll_ss, eryeetp2ll_ss).
• 16.  
  Return as calling arguments the staggered grid arrays of all three magnetic field components, all three electric field components, the radial derivative of the horizontal electric field components, the radial component of the Poynting flux, the relative helicity injection contribution function, the energy input rate into the upper atmosphere, and the relative helicity injection rate. Note that for the radial electric field component, we output both the total radial electric field and the purely inductive component. The inductive component is used when computing the horizontal components of the curl of ${\boldsymbol{E}}$, whereas the total radial electric field would be used for the evaluation of, e.g., ${\boldsymbol{E}}$ × ${\boldsymbol{B}}$, or for evaluating the angle between ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$ (subroutine angle_be_ss). The strong-field mask arrays for the COE, CO, CE, TE, and PE grids are also returned. All returned arrays are oriented in longitude–latitude index order.

The input data sets to and the output data sets from pdfi_wrapper4jsoc_ss are archived and publicly available through the SDO data center with the series names cgem.pdfi_input and cgem.pdfi_output, respectively. They can be directly accessed through the SDO JSOC website http://jsoc.stanford.edu, as are all SDO/HMI and AIA data, or through a variety of other means, including the Solarsoft IDL packages or the SunPy Python package. Users are referred to the SDO data analysis guides for data query and retrieval methods, such as http://jsoc.stanford.edu/How_toget_data.html and https://www.lmsal.com/sdouserguide.html. Each record in these two data series can be uniquely identified via two keywords, CGEMNUM and T_REC, which indicates the CGEM identification number and the nominal observation time, respectively. The CGEMNUM is currently defined to be identical to the NOAA AR number when the CGEM region coincides with a single named active region, and 100,000 plus the SHARP number (Bobra et al. 2014), when the CGEM region consists of multiple active regions or no named active regions. For cgem.pdfi_output, the nominal T_REC is designated at 06, 18, 30, 42, and 54 minutes after the hour. For example, users can find a pair of input records for AR 11158 at the beginning of 2011 February 15 with cgem.pdfi_input[11158][2011.02.15_00:00-2011.02.15_00:12], which includes vector magnetic field, the FLCT velocity field, the Doppler velocity, and the local unit normal vectors. The corresponding PDFI output can be found with cgem.pdfi_output[11158]. The processing necessary to prepare the input data (cgem.pdfi_input) is described in Section 4.

There is currently no formal way for deriving errors in the electric fields within the PDFI_SS software. The fact that the electric field solutions are derived from solutions of elliptic equations means that any magnetic field or Doppler velocity errors result in nonlocalized errors in the resulting electric fields, making analytic error propagation studies difficult. The effects of random errors in the magnetic field measurements and how these propagate into the PDFI electric field inversions in HMI data have been studied by Kazachenko et al. (2015) and Lumme et al. (2019). In Kazachenko et al. (2015), given estimated errors in the radial (30 G) and the two horizontal components (100 G) of the magnetic field determined from the width of distribution functions in the weak-field regions of NOAA AR 11158, this resulted in estimated relative errors of 15%–20% in the three electric field components at a given pixel location for a given pair of AR magnetograms. These results were derived by applying Monte Carlo techniques. Lumme et al. (2019) performed a more detailed error analysis on the PDFI electric fields that was focused primarily on global quantities such as the spatially and/or temporally integrated Poynting flux and helicity injection rate contribution functions. They showed that spatial averaging and temporal integration resulted in significantly lower relative errors than one obtained for individual pixel values for a pair of magnetograms.

Neither of these studies addresses another source of error, the systematic effects to the velocity and magnetic field signals that are due to incomplete corrections for the daily orbital motion of the SDO spacecraft around Earth. These effects appear to generate a false temporal signal at the first few harmonics of the orbital frequency in the magnetic and velocity signals. A false temporal signal in the magnetic field will generate a false electric field through Faraday's law. These systematic errors in the observed quantities from orbital artifacts are characterized by Hoeksema et al. (2014) and Schuck et al. (2016). Schuck et al. (2016) provide a suggested correction for the Doppler velocity that appears to remove much of the artificial temporal signal, but as of yet, no similar correction for the magnetic field components is available. While these systematic errors can affect the electric field solutions over several-hour timescales, short-term variations are small, and the work of Lumme et al. (2019) indicates that they do not greatly affect the time evolution on longer timescales. Nevertheless, the results of the PDFI_SS electric field solutions would be improved if these artifacts could be removed.

It is possible that the user may wish to obtain electric field solutions at a different resolution than the 003 resolution provided by the JSOC upstream processing (described in Section 4). One might be tempted to simply interpolate the output electric field results to a different resolution, but doing so will generally destroy the adherence of the solutions to Faraday's law (an exception to this rule is the flux-preserving "downsampling" subroutines, described in Section 5.1, but these only work for certain specified cases to decrease the resolution). We have found that if one wants electric field inversions with an arbitrary change of the resolution, the best solution is to interpolate the input data to the desired resolution and then compute the solutions from scratch from, e.g., subroutine pdfi_wrapper4jsoc_ss.

The interpolation technique we have used for this process is the ninth-order B-spline, a subset of interpolation solutions described by Thevenaz et al. (2000). The low-level source code for this interpolation procedure was written by coauthor Dave Bercik, inspired by Thevenaz et al. (2000) and the accompanying C source code at http://bigwww.epfl.ch/thevenaz/interpolation/. It is implemented in subroutine bspline_ss.

To interpolate a single one of the 18 input data arrays to a different resolution (either coarser or finer), one can use subroutine interp_hmidata_ll,

where the original and desired array dimensions can be specified. In this subroutine, the degree of the B-spline can be specified, but we recommend setting degree to 9. To interpolate the entire 18-level stack of arrays, input as a 3D array, interpolated to a new 18-level 3D array, one can use subroutine interp_hmidata_3d_ll.

This subroutine assumes degree=9. The latter two subroutines assume that the domain boundaries a, b, c, and d remain the same in the output interpolated data arrays as those values for the input arrays.

The Doppler velocity calibration software that computes the convective blueshift (Welsch et al. 2013) was initially written in FORTRAN by coauthor Brian Welsch and then modified by coauthor Xudong Sun to be called from an HMI module written in C. The module uses full-disk vector magnetograms and Doppler data as input and estimates a "bias" that we later subtract from the Doppler shift measurement. Additional output includes both LOS and radial PIL masks for the LOS magnetic field B_ℓ and the radial field B_r. The source code for this module can be seen by clicking on the "view" link at http://jsoc.stanford.edu/cvs/JSOC/proj/cgem/prep/apps/doppcal_estimate.f90. We have chosen to work with the Doppler data derived from the full spectral inversion rather than the traditional Doppler data derived from the LOS field pipeline, following the recommendation of the HMI Team. We have performed a comparison between the "vector Doppler" and "LOS Doppler" data. The comparison was done in a cutout that tracked NOAA AR 11158 in full-disk Doppler velocity maps, and it revealed that the two types of raw uncalibrated Doppler maps have systematic differences, with median difference oscillating in phase with the radial velocity of the SDO spacecraft. The removal of the convective blueshift using the method of Welsch et al. (2013) reduces the median difference between the two velocities significantly, particularly in strong-field pixels ($| {\boldsymbol{B}}| \gt 300\,{\rm{G}}$). Subsequent tests of the impact of the differences between Doppler velocities from the two different data sets on the calculation of the integrated Poynting flux and helicity injection rate showed only a modest difference. We conclude that while there are differences in the results using the two different data sets, our processing reduces these differences, and there is not a substantial difference in the final results.

Once the convective blueshift has been computed, it is used to correct the Doppler velocity measurements during the step described in Section 4.3 below.

This module extracts a series of AR vector field patches in native coordinates from full-disk data, with constant center latitude (rounded to the nearest pixel), and tracks them at a constant rotation rate. These patches are given a unique "CGEM number" as an identifier and are used as input for the subsequent modules.

This module takes a series of AR patches from Section 4.2; flips ambiguity choices that create large, transient changes in azimuth (see Welsch et al. 2013 for a detailed description); corrects the Doppler velocity with the bias computed in Section 4.1, computes the LOS unit vector; and maps these quantities onto a Plate Carrée (uniformly spaced in longitude and latitude) coordinate system with a pixel spacing of 0.03 heliographic degrees (coinciding closely with the HMI pixel size near disk center). We remove differential rotation based on the fit of Snodgrass (1984), remove the spacecraft velocity, and then correct for the Doppler bias computed from Section 4.1.

We currently use the local correlation tracking code FLCT (Fisher & Welsch 2008) to estimate horizontal flow velocities, which are then used to compute the noninductive contribution to the horizontal electric field described in Sections 2.3 and 3.9.4. The original idea for local correlation tracking was first described by November & Simon (1988).

The basic idea of the FLCT code is to link small changes in two images taken at two closely spaced times to a 2D flow velocity that moves features in the first image toward the corresponding features in the second image. To compute the "optical flow" velocity at a given pixel location, both images are multiplied by a windowing function, assumed to be a Gaussian of width σ_FLCT, centered at that given pixel location, to de-emphasize parts of the two images that are far away from the given location. The cross-correlation function of the resulting subimages is computed using Fourier transform techniques, and the location of the peak of the cross-correlation function is found to subpixel accuracy. The difference between the location of the peak and the original pixel location is assigned to be the distance of the pixel shift (in both horizontal directions), and this shift, divided by the time difference between images, is identified with the horizontal flow velocity at that pixel. This procedure is then repeated for all pixel locations in the two images. To compensate for noisy data in the images, the algorithm allows one to select a threshold parameter thr. If the average image value has an absolute value less than thr, no velocity is computed, and a mask value for that pixel is set to zero, to indicate that no value was computed. The velocity itself is then set to zero as well at that pixel. The code also allows the user to filter the images with a low-pass filter before computing the cross-correlation function, if there is a large-degree small-scale noise.

The FLCT algorithm as originally conceived was described in Welsch et al. (2004), with major improvements to the algorithm described in Fisher & Welsch (2008). Since the publication of that article, coauthors Fisher and Welsch have made a number of improvements to the algorithm and the code to increase the accuracy and speed of FLCT. The developer site for the FLCT code is a fossil repository, located at http://cgem.ssl.berkeley.edu/cgi-bin/cgem/FLCT/index. The latest version can always be downloaded there. We have also published a recent snapshot of the FLCT software from the above repository as an archive on Zenodo (Fisher & Welsch 2020).

First, the original C code as described in Fisher & Welsch (2008) was written as a stand-alone executable, intended to be used while running in an IDL session. To read in the image data, and to write out the resulting velocity fields, the information was communicated with IDL using disk I/O. While this works fine for an IDL session, it is inefficient and does not allow the FLCT method to be easily incorporated into other software. Therefore, the current version of FLCT has been rewritten as a library of functions, easily callable from C, FORTRAN, or Python programs. There is still also a stand-alone FLCT executable that has the same user interface as the original version, but this stand-alone code now consists mainly of I/O tasks and calls functions from the FLCT library to perform the main computation. The construction of the library was done in consultation with coauthors Erkka Lumme and Xudong Sun to make sure it could be used from the ELECTRICIT (Lumme et al. 2017, 2019) Python software, from the JSOC's HMI software, and from other FORTRAN test programs.

Second, the FLCT algorithm was rewritten so that the means of the subimages described above are subtracted from the subimages before the cross-correlation function was computed. We found that this resulted in more accurate results.

Third, while the FLCT algorithm as written strictly only applies in Cartesian coordinates, Welsch et al. (2009) described in an appendix of that article how data on a spherical surface can be mapped into a conformal Mercator projection. FLCT can then be run in this projection, and once the velocities are derived, they can be scaled and mapped back onto the spherical surface. We have now modified the FLCT code so that if the input images are given on a Plate Carrée grid, the code itself handles the mapping to the Mercator projection, runs the FLCT algorithm to find the velocities on the Mercator map, and then rescales and remaps the data back to the Plate Carrée grid.

Fourth, we have we performed a study of biases in the calculation of velocities using the FLCT code. A number of published studies have shown that FLCT tends to underestimate flow velocities in cases where the flow velocities are known. Two especially insightful articles on this topic are Freed et al. (2016) and Löptien et al. (2016). The appendix of Freed et al. (2016) quantifies this behavior as a function of FLCT input parameters. Our own study identifies a likely reason for the systematic velocity underestimates, in that the Gaussian windowing function at the heart of the algorithm is centered at the same pixel location in both images, even though the second image has been slightly shifted. We have developed an experimental technique to correct for this bias, which is an input option to the FLCT library functions.

Further details regarding these changes can be viewed in the README file in the latest FLCT distribution, along with documentation files in the doc folder within the distribution. A more complete discussion of the updated FLCT code will be described in a future article.

To compute the FLCT flow velocities in the Plate Carrée data for input to PDFI_SS, for each time, we use pairs of images of B_r that are one time step behind and one time step ahead of the current time. So, for the nominal HMI cadence of 12 minutes, the images are 24 minutes apart. The parameter σ_FLCT is chosen to be 5 pixels. The value of the threshold thr is set to 200 G. We also have chosen not to apply any low-pass filtering of the images in the FLCT code, as we find we get better results overall. For now, we have not implemented the experimental bias correction, but we may apply it in the future.

We have found that the properties of the electric field solutions are improved by adding a region of "padding" around the input data, in which a ribbon of data with a width of approximately 50–60 pixels is added to each of the four boundaries, with the values of the padded data for all 18 input arrays set to zero. The exact width for each padded region varies slightly, such that the resulting values of m and n are each divisible by 12. This property of the resulting data arrays facilitates the use of the electric field data by the CGEM magnetofrictional model of Cheung & DeRosa (2012), because this property of m and n makes it easier to set up computational runs that use many processors.

Adding the padding is done as the last step before defining the input data for the electric field inversions, and it is performed as part of the HMI magnetic pipeline. To mimic the padding operation within the PDFI_SS library, we have written several FORTRAN subroutines that do the same thing as the JSOC padding process. The subroutine pad_int_gen_ss

takes as input the unpadded values of m and n and "first-guess" values of the amounts of latitude and longitude padding, mpad0 and npad0, and computes output values of m and n and also outputs the exact amounts of padding that will be applied along each of the four boundaries, such that the output values of m and n are divisible by 12.

The adjusted values of a, b, c, and d are computed from the original values of a, b, c, and d, plus the four padding amounts returned by pad_int_gen_ss, by subroutine pad_abcd_as_ss. The padded arrays themselves can then have their interiors filled with the original, unpadded input data, by calling subroutine add_padding_as_ss.

In the test_wrapper.f test program (see Section 9.1), which mimics the call of pdfi_wrapper4jsoc_ss from the JSOC software, these three padding subroutines are called to mimic the same padding procedure performed by the JSOC software. The trial padding values mpad0 and npad0 are set to 50 pixels.

One of the important components of the CGEM project is an electric-field-based Surface Flux Transport Model (SFTM), developed by coauthors DeRosa and Cheung, the details of which will be described in a future publication. A summary can be found in the CGEM final report at http://cgem.ssl.berkeley.edu. The SFTM computes the global horizontal electric field in spherical coordinates based on differential rotation and meridional flows acting on the radial component of the magnetic field, along with a term that describes the dispersal of magnetic flux by supergranular motions. The electric field in the two horizontal directions is then used to evolve the radial magnetic field at the photosphere. The SFTM is used in regions of the Sun for which no PDFI electric fields have been computed, mainly outside of ARs. Where PDFI solutions are computed with PDFI_SS, the model inserts the PDFI solutions into the global domain and evolves B_r by using the PDFI solutions, rather than the SFTM solutions. There are two complications to doing this: first, the SFTM generally uses a coarser grid than is used by the PDFI solutions, and second, there will generally be a solution mismatch at the boundary between the PDFI_SS domain and the global SFTM model. Such a mismatch, if not corrected, results in a large, spurious curl of ${\boldsymbol{E}}$ at the boundaries, which will then result in a spurious evolution of B_r at the boundaries or "seams" where the PDFI electric fields are inserted. We now describe how we cope with these two complications.

In general, the SFTM is run with considerably coarser resolution than the 003 resolution computed by default with, e.g., pdfi_wrapper4jsoc_ss in PDFI_SS. Before the PDFI_SS solutions can be inserted into the SFTM model, both solutions must have the same grid resolution. Our approach is to perform a flux-preserving "downsampling" of the higher-resolution electric field results to the same grid resolution that is used by the SFTM. This must be done in such a way that the magnetic flux evolution in the coarser grid is physically consistent with that in the finer grid.

Our solution is to define "macropixels" for the coarser grid in terms of the fine grid, such that there is a whole integer number of fine grid edges fitting within the macropixel edges, in both horizontal directions, and that the line integral of the electric field around the edges of a macropixel is equal to the line integral of the electric field along those fine grid pixels that touch the macropixel boundary. This condition is illustrated schematically in Figure 8.

Zoom In Zoom Out Reset image size

The downsampling, considering only horizontal components of the electric field, can be accomplished with subroutine downsample_ss

when using colatitude–longitude array orientation, or subroutine downsample_ll

when using longitude–latitude array orientation, which is the relevant case for the SFTM model. Once downsample_ll has been called, then the coarse-resolution PDFI_SS horizontal electric fields can be inserted into the SFTM model results.

For completeness, we have also written two additional subroutines, downsample3d_ss and downsample3d_ll, which downsample not only the horizontal components of the electric field but also the radial component E_r and the radial derivatives of the horizontal components of the electric field. This additional information is needed to create a downsampled three-component electric field that can be used to compute all three components of Faraday's law in the coarser grid in a way that is consistent with the solutions on the original finer grid.

Now we discuss the problem of the mismatch between electric fields in the SFTM and the PDFI solutions, once the latter have been downsampled to the same grid resolution in SFTM. The idea is to add a solution to the PDFI results that has zero curl, but which then matches the SFTM results at the PDFI domain boundaries.

For a curl-free electric field with specified values of the tangential electric field on its boundaries, $\dot{P}$ obeys the Laplace equation

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{h}^{2}\dot{P}\,=\,0,\end{eqnarray} \tag{ 78 }$$

where the tangential component of the horizontal electric field on the boundaries is related to $\dot{P}$ by Equation (13). It thus follows that the Neumann boundary conditions needed by the FISHPACK subroutine HSTSSP are given by

$$\begin{eqnarray}&&{\left.{\left.\displaystyle \frac{\partial \dot{P}}{\partial \theta }\right|}_{\theta =a,b}={{rcE}}_{\phi }\right|}_{\theta =a,b}\end{eqnarray} \tag{ 79 }$$

for the n points along the north and south boundaries at θ = a and θ = b, respectively, and

$$\begin{eqnarray}&&{\left.{\left.\displaystyle \frac{\partial \dot{P}}{\partial \phi }\right|}_{\phi =c,d}=-r\sin {\theta }_{i+\tfrac{1}{2}}{{cE}}_{\theta }\right|}_{\phi =c,d}\end{eqnarray} \tag{ 80 }$$

for the m points along the left and right boundaries at ϕ = c and ϕ = d, respectively. Once the Laplace equation for $\dot{P}$ is solved with these boundary conditions, the horizontal components of the electric field within the domain are evaluated by taking minus the curl of $\dot{P}\hat{{\boldsymbol{r}}}$. These operations are performed by subroutine e_laplace_ss

for arrays in colatitude–longitude orientation and by e_laplace_ll,

where the input electric field components at the boundaries and the output electric fields within the domain are computed in longitude–latitude orientation. The latter case is the one relevant to SFTM, which uses longitude–latitude orientation exclusively.

In the SFTM model, the electric field components at the boundaries on input to these subroutines in Equations (79) and (80) are defined by the difference between the initial SFTM electric field values and the downsampled PDFI electric field values at the boundary locations.

Another useful application of our curl-free electric field solutions is to match boundary conditions assumed in computational models for the solar atmosphere. The PDFI electric field solutions computed by subroutine pdfi_wrapper4jsoc_ss can have nonzero electric field components parallel to the domain boundaries, originating from the contributions from gradients in the scalar potentials. If a computational model requires that the electric field parallel to the boundary is zero, and if the radial magnetic field data are flux balanced, then subroutines e_laplace_ss or e_laplace_ll can be used to compute a curl-free electric field solution that exactly matches the PDFI solution for the tangential component of ${\boldsymbol{E}}$ on the boundary. That solution can then be subtracted from the PDFI solution, yielding solutions for ${\boldsymbol{E}}$ that have zero tangential electric field on the boundaries but still obey all three components of Faraday's law.

Imagine that we have a computational model for the temporal evolution of ${\boldsymbol{B}}$ in a volume, with the lower boundary surface of the volume coinciding with the photosphere, where we have evaluated B_r at the centers of cells in a Plate Carrée grid, and for which we have computed electric field solutions on the edges or rails that surround the cells, using PDFI_SS solutions. The HMI data and the electric field solutions together define a time sequence of magnetic field and electric field solutions that are consistent with one another, at least in terms of Faraday's law. However, the computational model will in general be based on an additional set of physical or mathematical assumptions that can contain far more constraints on how ${\boldsymbol{B}}$ behaves in the model. Given some initial condition for ${\boldsymbol{B}}$ at t = 0 that matches B_r at the photosphere, is there any guarantee that the model's evolution for ${\boldsymbol{B}}$ will be consistent with how the HMI magnetic field behaves at the photosphere? In general, the answer to this is no. Given that, sooner or later, the computational model will "go off the rails" as compared to how the observed magnetic field changes over time, what can we do to "nudge" the model to get back on track?

In PDFI_SS, we have developed a series of subroutines that are designed to compute a nudging electric field, in effect giving the computational model a "kick" to make its evolution behave more consistently with the observed magnetic field data. The idea is to use the mismatch between the computational model and the data to compute an electric field that is designed to return the model's magnetic field evolution to match the photospheric magnetic field evolution.

To illustrate this in the simplest way, we consider the computational model to be the spherical version of the magnetofrictional coronal model developed by Cheung & DeRosa (2012). In this model, the observed values of the horizontal components of the magnetic field are not used, and the model is constrained to match the observed evolution of B_r at the centers of the photospheric cells, i.e., on the CE grid at the photosphere:

$$\begin{eqnarray}&&\displaystyle \frac{\delta {B}_{r}(t+{\rm{\Delta }}t)-\delta {B}_{r}(t)}{{\rm{\Delta }}t}=-{\boldsymbol{\nabla }}\times \delta c{{\boldsymbol{E}}}_{h},\end{eqnarray} \tag{ 81 }$$

where $\delta {B}_{r}(t)={B}_{r}^{\mathrm{target}}(t)-{B}_{r}^{\mathrm{model}}(t)$. We can use the PTD approximation to compute $\delta c{{\boldsymbol{E}}}_{h}$, where

$$\begin{eqnarray}&&\delta c{{\boldsymbol{E}}}_{h}=-{\boldsymbol{\nabla }}\times \dot{P}\hat{{\boldsymbol{r}}},\end{eqnarray} \tag{ 82 }$$

and $\dot{P}$ obeys the 2D Poisson Equation (9), with the source term equal to the left-hand side of Equation (81). The boundary conditions assumed are the same as those employed in determining $\dot{P}$ in subroutine ptdsolve_ss (Section 3.9.1).

A useful way to think about this is to imagine what happens over a single time step Δt taken by the model, assuming that both the target and model magnetic field values are equal at time t:

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{r}^{\mathrm{target}}(t+{\rm{\Delta }}t)-{B}_{r}^{\mathrm{model}}(t+{\rm{\Delta }}t)}{{\rm{\Delta }}t}=-{\boldsymbol{\nabla }}\times \delta c{{\boldsymbol{E}}}_{h}.\end{eqnarray} \tag{ 83 }$$

The vector quantity $\delta c{{\boldsymbol{E}}}_{h}$ is the electric field that must be added to the model's electric field to return B_r to its observed value at time t + Δt.

Depending on the details of the computational model, such a nudging step could be taken within the model's own time-advancing algorithm, or alternatively, if the error is small, it can just be added to the model's electric field on output and applied to the calculation of B_r for the next time step evolution. The latter case is how the nudging electric field is used within CGEM's spherical magnetofrictional model.

We must add an additional comment on the use of nudging when using electric fields determined from PDFI solutions, as described in Section 3, to drive numerical models. The procedure described in that section recommends using strong-field masks for computing solutions for the electric field. From the perspective of deriving electric fields from the data that are physically meaningful in the presence of magnetic field noise, this is the correct thing to do. However, using these electric fields to drive a numerical model without also using nudging can result in inconsistencies when particular regions of the domain move from being within the strong-field mask region to being in the weak-field region, as time evolves: a given pixel initially within the strong-field region that has a nonzero curl of ${\boldsymbol{E}}$ will suddenly have zero curl, meaning that the magnetic field at that point will no longer evolve forward in time if only the PDFI solutions of Section 3 are used to drive the model. The use of an additional nudging electric field step, with no strong-field masks applied, will then allow regions of the domain that move between strong-field and weak-field regions to evolve in a way that is consistent with the magnetic field data in both regions. This is how the CGEM magnetofrictional model uses nudging electric fields to address this particular problem.

The nudging electric field in Equation (81) is formally just the horizontal components of the PTD (inductive) electric field from $\dot{P}$. It can be computed with subroutine enudge_ss

for input B_r error terms in Equation (83), and with outputs E_θ and E_ϕ. If the B_r error term is oriented in longitude–latitude array order, one can use subroutine enudge_ll

to compute the corresponding components E_lon and E_lat in longitude–latitude array order.

Another useful application of enudge_ll is within the CGEM SFTM for handling the magnetic field evolution of ARs rotating onto the disk from the east limb of the Sun. The magnetic field observations, if incorporated directly into the SFTM, have unwanted impacts on global magnetic flux balance and other distortions due to the extreme viewing angle. Once the magnetic structure of the rotating AR becomes clearer a couple of days after the AR has rotated onto the disk, the nudging software can be used to reconstruct an artificial but physically reasonable evolution, by using as input to enudge_ll the term on the left-hand side of Equation (81) equal to the difference of the magnetic field 2 days after rotating onto the disk and an initial δB_r of 0. The value of Δt in Equation (81) is then set to 2 days. This allows the AR to grow in a natural way within the SFTM without having the unwanted global impacts on the SFTM solution. At that point, the PDFI_SS solutions can begin to be inserted directly into the SFTM as described in Section 5.1.

The concept of nudging can be generalized to include the calculation of all three components of ${\boldsymbol{E}}$, in response to evolution in a computational model that computes all three components of ${\boldsymbol{B}}(t)$, instead of just the evolution of B_r(t). The same general concept is used: differences between a model's temporal evolution of ${\boldsymbol{B}}$ versus a "target" observed evolution can be used to derive corrective values for all three components of ${\boldsymbol{E}}$ using the PTD solutions. The primary difference is that in the latter case the electric field components are computed on all the edges (rails) in a 3D layer of voxels bisected by the photosphere, in contrast with the 2D case, in which the horizontal electric fields are computed along the edges surrounding the photospheric face with B_r computed on the CE grid. The subroutine enudge3d_ss

can be used to compute the electric fields on all the edges of the voxels, given source term time derivatives for B_r, B_θ, and B_ϕ and the radial thickness of the voxels. This subroutine assumes that all arrays are in colatitude–longitude order.

The emphasis of most of the software in PDFI_SS is for spherical wedge domains that subtend only a subset of the full spherical domain. However, for completeness, we have written some electric field software for the global domain (4π sr), including PTD solutions that could also be used for computing nudging solutions in a global domain. This software takes advantage of the special case of "global" boundary conditions available in some of the FISHPACK Helmholtz/Poisson equation subroutines, plus some additional constraints to be applied at the north and south poles in PDFI_SS subroutines. A good discussion of the "global" boundary conditions at the poles can be found in the description of the FISHPACK subroutine PWSSSP in Swarztrauber & Sweet (1975).

The global versions of enudge_ss and enudge_ll, which compute horizontal electric field components from a global distribution of ${\dot{B}}_{r}$, are computed by subroutines enudge_gl_ss, and enudge_gl_ll,

for arrays oriented in colatitude–longitude and longitude–latitude order, respectively. Note that with ${\dot{B}}_{r}$ defined on the CE grid, there are no values of B_r defined at the north or south poles. On the other hand, the output arrays of the azimuthal or longitudinal component of the electric field are defined at the poles. Note, however, that physical considerations demand that this component of ${\boldsymbol{E}}$ must be zero at the poles, or else the behavior of B_r would become singular. The colatitudinal (or latitudinal) component of ${\boldsymbol{E}}$ is not defined at the poles. The global solutions for the poloidal potential $\dot{P}$ within these two subroutines do not include ghost zones, in contrast to the spherical wedge solutions.

While localized spherical wedge solutions can have a flux imbalance, the global solutions must be flux balanced, to avoid a monopole term in B_r or ${\dot{B}}_{r}$. Any existing monopole term in the input data is removed before the electric fields are computed. The subroutines fluxbal_ss and fluxbal_ll are used to compute a corrected input field that is flux balanced. The flux balance is corrected in such a way that the locations of preexisting PILs are not moved. The algorithm can be summarized as follows: The positive and negative magnetic fluxes within the domain are summed separately. Whichever polarity is the minority polarity then has the flux in each of its constituent pixels enhanced by a constant relative factor such that the total net radial flux is zero. This is similar to a technique proposed by Yeates (2017), except that in the latter case both polarity regions are adjusted.

In an analogy with the subroutine enudge3d_ss, we can also compute global PTD solutions for all three components of ${\boldsymbol{E}}$, given the time derivatives ${\dot{B}}_{r}$, ${\dot{B}}_{\theta }$, and ${\dot{B}}_{\phi }$ given on a global, staggered grid. The electric field components are computed on all the edges (rails) of a global set of spherical voxels, similar to the geometry assumed in enudge3d_ss. The solutions are computed by subroutines enudge3d_gl_ss and enudge3d_gl_ll,

for arrays oriented in colatitude–longitude and longitude–latitude order, respectively. The poloidal and toroidal potentials are solved using the "global" FISHPACK boundary conditions mentioned earlier.

There are a number of specific considerations for the north and south poles and the left and right boundaries that must be mentioned. First, the toroidal potential $\dot{T}$ is defined at the north and south poles. Physically, $\dot{T}$ is related to ${\dot{J}}_{r}$ at the poles through the Poisson Equation (10), so we need to evaluate the radial current density at the poles. We estimate this quantity by using Ampere's law for ${\dot{B}}_{\phi }$ along the highest latitudes and then dividing by the area subtended by this small disk to estimate ${\dot{J}}_{r}$ at the poles. Similarly, we can use periodic boundary conditions for ${\dot{B}}_{\phi }$ to evaluate ${\dot{J}}_{r}$ at the left and right boundaries in ϕ. Second, the quantity ${\dot{B}}_{\theta }$ can be defined at the north and south poles, but its value has no effect on the calculation of electric fields from Faraday's law, because the amount of magnetic flux across the θ faces at the poles is zero. Therefore, we assume that ${\dot{B}}_{\theta }$ is zero at the north and south poles, for simplicity. Internal to these subroutines, there are no ghost zones used in the solutions for the poloidal or toroidal potentials.

In order to test the accuracy with which electric field solutions obey Faraday's law, we need to be able to calculate the curl of ${\boldsymbol{E}}$. Here we describe a number of subroutines we have written to do this.

One simple and common example is taking the radial component of the curl of the horizontal components of ${\boldsymbol{E}}$. Given arrays of E_θ and E_ϕ on the PE and TE grids, respectively, subroutine curlehr_ss

will compute $\hat{{\boldsymbol{r}}}\cdot {\boldsymbol{\nabla }}\times c{{\boldsymbol{E}}}_{h}$ evaluated on the CE grid. This can be compared directly to ${\dot{B}}_{r}$, the radial time derivative of B_r (they should be equal and opposite). This subroutine assumes that the arrays are all in colatitude–longitude order.

There are several approaches to computing the curl of ${\boldsymbol{E}}$ for all three components. If the components of ${\boldsymbol{E}}$ are computed on all the rails of a layer of voxels, subroutines curle3d_ss and curle3d_ll

can compute all three components of the curl of ${\boldsymbol{E}}$. The quantities returned by these subroutines are actually minus the curl of ${\boldsymbol{E}}$, so they can be compared directly with the time derivative of the three components of the magnetic field and should be equal. It is important that the radial component of ${\boldsymbol{E}}$ contain only the inductive contribution to E_r. The subroutine curle3d_ss assumes that arrays are in colatitude–longitude order, while curle3d_ll assumes that arrays are in longitude–latitude order. Both of these subroutines can handle either spherical wedge electric field solutions or global electric field solutions.

If one is dealing strictly with electric fields defined at the photosphere, and not in a layer of spherical voxels bisected by the photosphere, there is a different approach that can be used. First, if radial derivatives of the horizontal electric field components have been computed at the photosphere, as is the case when using, e.g., subroutine pdfi_wrapper4jsoc_ss, or when downloading the electric field solutions from the JSOC, one can use subroutine curle3dphot_ss to compute all three components of the curl of ${\boldsymbol{E}}$, evaluated at the photosphere. If using quantities downloaded from the JSOC, you will need to (1) transpose the arrays from longitude–latitude to colatitude–longitude array order, (2) convert the units of the radial and horizontal E-fields from V cm⁻¹ to G km s⁻¹ by multiplying by 1000, and (3) convert the units of ∂E_θ/∂r and ∂E_ϕ/∂r from V cm⁻² to G s⁻¹ by multiplying by 10⁸. It is essential that E_r include only the inductive contribution to E_r when evaluating the curl. The three components returned by the subroutine are actually minus the curl of ${\boldsymbol{E}}$, so they can be compared directly with the time derivatives of the three magnetic field components.

If the radial derivatives of E_θ and E_ϕ are not available, they can be computed with subroutine dehdr_ss,

which uses as input the solutions for the poloidal potential and its radial derivative, as returned from, e.g., ptdsolve_ss.

The broad outline of the procedure for finding the poloidal potential P is as follows: (1) convert the continuum version of Bercik's Equation (88) to a second-order-accurate finite-difference equation for P; (2) convert the finite-difference version of the azimuthal second derivative term in the horizontal Laplacian to an eigenvalue problem, by Fourier transforming the second-order finite-difference contribution in the azimuth (longitude) direction; (3) derive a 2D finite-difference expression as a function of colatitude and radius for the amplitude of each Fourier mode, with each mode obeying specified boundary conditions in r and θ; (4) solve each one of these resulting 2D elliptic problems using the BLKTRI subroutine in FISHPACK, and (5) inverse transform the resulting solution back to a function in 3D space.

These tasks are all performed within subroutine scrbpot_ss,

which returns the solution P as a 3D array. We now describe these steps in detail.

6.1.1. The Continuum Equation for P and Defining the Finite-difference Grid

The model is based on a solution for P in a 3D domain, where P obeys Bercik's Equation (88). We first multiply Equation (88) by r² and can then write the resulting equation as

$$\begin{eqnarray}&&{L}_{\phi }(P)+{L}_{\theta }(P)+{L}_{r}(P)=0,\end{eqnarray} \tag{ 93 }$$

where we have decomposed the left-hand side of the equation into three operators acting on the poloidal potential P. Using Equation (37) for ${{\boldsymbol{\nabla }}}_{h}^{2}P$, we can write these operators as

$$\begin{eqnarray}&&{L}_{\phi }(P)=\displaystyle \frac{1}{{\sin }^{2}\theta }\displaystyle \frac{{\partial }^{2}P}{\partial {\phi }^{2}},\end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray}&&{L}_{\theta }(P)=\displaystyle \frac{1}{\sin \theta }\displaystyle \frac{\partial }{\partial \theta }\left(\sin \theta \displaystyle \frac{\partial P}{\partial \theta }\right),\end{eqnarray} \tag{ 95 }$$

and

$$\begin{eqnarray}&&{L}_{r}(P)={r}^{2}\left(\displaystyle \frac{{\partial }^{2}P}{\partial {r}^{2}}\right).\end{eqnarray} \tag{ 96 }$$

The locations where P is defined must be consistent with the 2D staggered grid locations defined in Section 3.4. In three dimensions, we have 3D voxels instead of 2D cells. P must be defined on the radial faces of the voxels and in the center of these faces, to be consistent with the location of P on the CE grid at the surface of the photosphere. We will therefore denote the indices for P with i + $\tfrac{1}{2}$ as the θ index, j + $\tfrac{1}{2}$ as the ϕ index, and q for the index of radial faces, in keeping with the index notation described in Section 3.7. We will place the source surface as the last radial face of the active part of the domain. If there are p + 1 radial faces in the r-direction, it means that the radial spacing Δr is given by

$$\begin{eqnarray}&&{\rm{\Delta }}r=({R}_{\mathrm{SS}}-{R}_{\odot })/p,\end{eqnarray} \tag{ 97 }$$

where p is the number of voxels between R_⊙ and R_SS. The dimension of P is therefore (m, n, p + 1) in the θ-, ϕ-, and r-directions, respectively. Figure 9 shows a diagram of a voxel in this 3D grid.

Zoom In Zoom Out Reset image size

In this section of the article, where describing finite-difference expressions, we assume index ranges that start at 0 and go to p in the radial direction for P, from $\tfrac{1}{2}$ to m − $\tfrac{1}{2}$ in the θ-direction, and from $\tfrac{1}{2}$ to n − $\tfrac{1}{2}$ in the ϕ-direction. But keep in mind that when examining these expressions in the source code we have used default index ranges in FORTRAN, where the first index starts from 1.

6.1.2. Fourier Transform P in Azimuth and Derive Finite-difference Equations for Each Fourier Mode

We now make the assumption that the solution for P can be separated into a product of eigenfunctions in the azimuthal direction multiplied by coefficients that are a function of colatitude θ and radius r, where each eigenfunction can be enumerated by a wavenumber index, j'.

Let

$$\begin{eqnarray}&&{P}_{i+\tfrac{1}{2},j+\tfrac{1}{2},q}=\sum _{j^{\prime} =0}^{n-1}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}{{\rm{\Phi }}}_{j^{\prime} }({\phi }_{j+\tfrac{1}{2}}^{{\prime} }),\end{eqnarray} \tag{ 98 }$$

where ${{\rm{\Phi }}}_{j^{\prime} }(\phi ^{\prime} )$ are a series of orthogonal basis functions in ϕ' and ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}$ are the amplitudes for each one of these n basis functions. Here, the azimuthal variable ϕ' has its range normalized to be from 0 to 2π, instead of from c to d:

$$\begin{eqnarray}&&\phi ^{\prime} =\displaystyle \frac{(\phi -c)}{(d-c)}\times 2\pi .\end{eqnarray} \tag{ 99 }$$

The expression for L_ϕ operating on P when using second-order-accurate finite differences in ϕ becomes, for each Fourier mode j',

$$\begin{eqnarray}&&\begin{array}{l}{L}_{\phi }({Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}{{\rm{\Phi }}}_{j^{\prime} }({\phi }_{j+\tfrac{1}{2}}^{{\prime} }))={\left(\displaystyle \frac{2\pi }{d-c}\right)}^{2}\displaystyle \frac{{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}}{{\sin }^{2}\,{\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}\phi {{\prime} }^{2}}\\ \ \ \times \ \left({{\rm{\Phi }}}_{j^{\prime} }({\phi }_{j+\tfrac{1}{2}}^{{\prime} }+{\rm{\Delta }}\phi ^{\prime} )+{{\rm{\Phi }}}_{j^{\prime} }({\phi }_{j+\tfrac{1}{2}}^{{\prime} }-{\rm{\Delta }}\phi ^{\prime} )\right.\\ \ \ \left.-\ 2{{\rm{\Phi }}}_{j^{\prime} }({\phi }_{j+\tfrac{1}{2}}^{{\prime} }\right).\end{array}\end{eqnarray} \tag{ 100 }$$

Here, the factor of (2π/(d−c))² accounts for the scaling of the second derivative between ϕ and ϕ', and the quantity Δϕ^'2 is the square of the corresponding spacing between grid points in ϕ' (Δϕ' = 2π/n).

If we let the basis functions ${{\rm{\Phi }}}_{j^{\prime} }(\phi ^{\prime} )$ be complex exponentials

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{j^{\prime} }(\phi ^{\prime} )=\exp ({ik}(j^{\prime} )\phi ^{\prime} )\end{eqnarray} \tag{ 101 }$$

(or sines and cosines over the same range of ϕ'), then it is straightforward to show that the above finite-difference expression becomes

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\phi }\left({Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}{{\rm{\Phi }}}_{j^{\prime} }({\phi }_{j+\tfrac{1}{2}}^{{\prime} }\right) & = & \displaystyle \frac{-2\left(1-\cos (k(j^{\prime} ){\rm{\Delta }}\phi ^{\prime} \right)}{{\sin }^{2}{\theta }_{i+\tfrac{1}{2}}\,{\rm{\Delta }}{\phi }^{2}}\\ & & \times \ {Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}{{\rm{\Phi }}}_{j^{\prime} }({\phi }_{j+\tfrac{1}{2}}^{{\prime} }).\end{array}\end{eqnarray} \tag{ 102 }$$

Note that the factors of (2π/(d−c))² and the expression for Δϕ^'2 = (2π/n)² occuring in Equation (100) result simply in division by Δϕ² in Equation (102). Equation (102) depends explicitly on wavenumbers k(j'), whose values depend on the details of the Fourier transform implementation. In the limit of low wavenumber, the cosine expression in Equation (102) results in the second derivative term being proportional to −k(j')², as one would expect, but as the wavenumber increases, the behavior deviates from this, also as one might expect since the finite-difference expression begins to deviate from the spectral derivative result.

The most important point is that the result of applying the L_ϕ operator to ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}{{\rm{\Phi }}}_{j^{\prime} }(\phi ^{\prime} )$ is simply multiplication by a factor (the eigenvalue of the operator) times that same function. Since the other two operators L_θ and L_r that define Bercik's equation do not depend on ϕ at all, the result will be a common factor of ${{\rm{\Phi }}}_{j^{\prime} }$ for all three operators, which can then be factored out. Furthermore, since the ${{\rm{\Phi }}}_{j^{\prime} }$ are all orthogonal to each other, the sum of the three operators for the Bercik equation acting on the solution must be zero not only for the entire solution but also for each individual term in the expansion (98). Therefore, for each value of the Fourier mode j', we need to determine only the coefficients ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}$. When evaluating the finite-difference versions of L_θ and L_r, we will therefore consider their action only on ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}$.

Evaluating Equation (95) using second-order-accurate finite differences, applied to ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}$, we find

$$\begin{eqnarray}&&\begin{array}{l}{L}_{\theta }({Q}_{j^{\prime} }^{i+\tfrac{1}{2},q})=\displaystyle \frac{\sin {\theta }_{i}}{\sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}{\theta }^{2}}{Q}_{j^{\prime} }^{i-\tfrac{1}{2},q}-\,\displaystyle \frac{\sin {\theta }_{i+1}+\sin {\theta }_{i}}{\sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}{\theta }^{2}}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}\\ \ \ \ +\ \displaystyle \frac{\sin {\theta }_{i+1}}{\sin {\theta }_{i+\tfrac{1}{2}}{\rm{\Delta }}{\theta }^{2}}{Q}_{j^{\prime} }^{i+\displaystyle \frac{3}{2},q}\end{array}\end{eqnarray} \tag{ 103 }$$

for i + $\tfrac{1}{2}$ that is not adjacent to the θ = a or θ = b boundaries. For $i+\tfrac{1}{2}=\tfrac{1}{2}$, adjacent to the θ = a boundary, the homogeneous Neumann boundary condition on P means that the ghost zone value ${Q}_{j^{\prime} }^{-\tfrac{1}{2},q}$ must be equal to ${Q}_{j^{\prime} }^{\tfrac{1}{2},q}$. Since the expression for the operator that will be input into BLKTRI cannot involve ghost zones, we can make that substitution into the operator equation to eliminate ${Q}_{j^{\prime} }^{-\tfrac{1}{2},q}$ and then find

$$\begin{eqnarray}&&{L}_{\theta }({Q}_{j^{\prime} }^{\tfrac{1}{2},q})=-\,\displaystyle \frac{\sin {\theta }_{1}}{\sin {\theta }_{\tfrac{1}{2}}{\rm{\Delta }}{\theta }^{2}}{Q}_{j^{\prime} }^{\tfrac{1}{2},q}+\displaystyle \frac{\sin {\theta }_{1}}{\sin {\theta }_{\tfrac{1}{2}}{\rm{\Delta }}{\theta }^{2}}{Q}_{j^{\prime} }^{\displaystyle \frac{3}{2},q}.\end{eqnarray} \tag{ 104 }$$

Doing a similar exercise for $i+\tfrac{1}{2}=m-\tfrac{1}{2}$, adjacent to the θ = b boundary, after applying the homogeneous Neumann boundary condition we have

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\theta }({Q}_{j^{\prime} }^{m-\tfrac{1}{2},q}) & = & \displaystyle \frac{\sin {\theta }_{m-2}}{\sin {\theta }_{m-\tfrac{1}{2}}{\rm{\Delta }}{\theta }^{2}}{Q}_{j^{\prime} }^{m-\displaystyle \frac{3}{2},q}\\ & & -\ \displaystyle \frac{\sin {\theta }_{m-2}}{\sin {\theta }_{m-\tfrac{1}{2}}{\rm{\Delta }}{\theta }^{2}}{Q}_{j^{\prime} }^{m-\tfrac{1}{2},q}.\end{array}\end{eqnarray} \tag{ 105 }$$

For the finite-difference version of the L_r operator Equation (96) acting on ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}$ we have

$$\begin{eqnarray}&&{L}_{r}({Q}_{j^{\prime} }^{i+\tfrac{1}{2},q})=\displaystyle \frac{{r}_{q}^{2}}{{\rm{\Delta }}{r}^{2}}\left({Q}_{j^{\prime} }^{i+\tfrac{1}{2},q-1}-2{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}+{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q+1}\right).\end{eqnarray} \tag{ 106 }$$

The boundary condition at the last radial point, r_p = R_SS, is determined by the outer boundary condition that Ψ is a constant we can set to 0, meaning that the radial derivative of ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},p}$ is zero. The ghost zone value, ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},p+1}$, must therefore be equal to ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},p}$, which then results in the equation for the operator acting on the last radial point

$$\begin{eqnarray}&&{L}_{r}({Q}_{j^{\prime} }^{i+\tfrac{1}{2},p})=\displaystyle \frac{{r}_{p}^{2}}{{\rm{\Delta }}{r}^{2}}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},p-1}-\displaystyle \frac{{r}_{p}^{2}}{{\rm{\Delta }}{r}^{2}}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},p}.\end{eqnarray} \tag{ 107 }$$

6.1.3. Getting the Finite-difference Equations into Block Tri-diagonal Form

6.1.4. Fourier Transform Details: Applying Azimuthal Boundary Conditions and Determining Wavenumbers

6.1.5. Matching the Solution to Observed B_r at Photosphere, and Photospheric Boundary Conditions for Fourier Coefficients

At the photospheric layer (q = 0, or r = R_⊙) we specify the arrays an, bn, cn so that at the first array elements all three array values are set to 0. This means that at this layer we ignore the radial variation of ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}$ and instead will set the horizontal Laplacian (determined by the am, bm, cm arrays) to match the observed values of B_r. To do this, we must determine from the observed data what the value of each Fourier coefficient ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q=0}$ is at the photosphere.

The procedure is straightforward. Given the observed photospheric array of B_r on the CE grid $\left(\mathrm{at}\ i+\tfrac{1}{2},j+\tfrac{1}{2}\right)$, we first solve the horizontal Poisson equation

$$\begin{eqnarray}&&{R}_{\odot }^{2}{{\boldsymbol{\nabla }}}_{h}^{2}P(q=0)=-{R}_{\odot }^{2}{B}_{r}\end{eqnarray} \tag{ 108 }$$

using FISHPACK subroutine HSTSSP. Once we have the solution, we then Fourier transform the solution in the ϕ-direction. The Fourier transform then results in the values of ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q=0}$. Applying the horizontal Laplacian operator (from arrays am, bm, cm) will result in the Fourier transform of −r²B_r at the photosphere, $(-{r}^{2}{B}_{r}{)}_{j^{\prime} }^{i+\tfrac{1}{2},q=0}$. This can be used to specify a 2D source term array expected by BLKTRI, y, which is dimensioned (m, p + 1). For each Fourier mode j', we can set y(0:m-1,0)=$(-{r}^{2}{B}_{r}{)}_{j^{\prime} }^{i+\tfrac{1}{2},q=0}$. All the other values of y for q > 0 are set to 0, consistent with the right-hand side of Bercik's equation being set to 0 for all radial layers above the photosphere. In principle, one should be able to Fourier transform the observed array B_r directly and do the same thing, but in practice this produces significant artifacts mainly due to the effects of flux imbalance in the input data. The procedure as described above, on the other hand, appears to be robust and accurate.

6.1.6. Assembling the 3D Solution from BLKTRI

6.1.7. Testing the Accuracy of scrbpot_ss

We first make some general comments about the properties of the solution for P determined from subroutine scrbpot_ss: (1) for well-resolved solutions, the resulting 3D array P can be huge; (2) if the input radial magnetic field at the photosphere has a net flux imbalance, the P solution will not reflect the flux imbalance (i.e., it is consistent with a net radial magnetic flux of zero); and (3) the solution for P includes no ghost zones. Part of the reason for not constructing ghost zones is because the array is already so large. In addition, we find that where ghost zones are needed to evaluate curls or gradients, we can add them on an as-needed, layer-by-layer basis. All of the subroutines we discuss here perform this operation internally when necessary.

While the solution for P computed by scrbpot_ss has no net radial flux, we feel that it is important for the potential field solutions to include a net flux when the user desires it. We therefore include a subroutine, mflux_ss, which can be used to compute the net radial flux from the input radial magnetic field data. The resulting net radial flux is then used to augment the solution for P to result in a potential magnetic field solution that is consistent with the data. The output from mflux_ss is a single value of the net radial flux Φ_M over the photospheric domain defined by the values of R_⊙, a, b, c, and d.

The vector potential ${{\boldsymbol{A}}}^{P}$ within the 3D volume can be computed in a straightforward way from P and from Φ_M. The radial component of ${{\boldsymbol{A}}}^{P}$ is zero, so only the horizontal components of ${{\boldsymbol{A}}}^{P}$ are computed.

The vector potential ${{\boldsymbol{A}}}^{P}$ is computed by subroutine ahpot_ss, using P and Φ_M on input, and on output computing the two components of ${{\boldsymbol{A}}}^{P}$, A_θ and A_ϕ, each of which are 3D arrays of dimension (m, n + 1, p + 1) and (m + 1, n, p + 1), respectively. The quantity A_θ is computed along radial faces of the voxels, on the PE (phi-edge) grid locations in the horizontal directions, and A_ϕ is also on radial faces but on the TE (theta-edge) grid locations in the horizontal directions. To compute ${{\boldsymbol{A}}}^{P}$, there is an outer loop over the radial index q. Then, for each radial layer, ghost zones are added to P that correspond to the boundary conditions assumed at the θ and ϕ edges of the domain. Then, ${\boldsymbol{\nabla }}\times \hat{{\boldsymbol{r}}}P$ is computed for that given radius using subroutine curl_psi_rhat_ce_ss, populating the A_θ and A_ϕ arrays at that radius. After that, an additional term is added to A_ϕ:

$$\begin{eqnarray}&&{A}_{\phi }^{{{\rm{\Phi }}}_{M}}=-{B}_{0}\displaystyle \frac{{R}_{\odot }^{2}}{{r}_{q}}\cot ({\theta }_{i}),\end{eqnarray} \tag{ 109 }$$

where B₀ = Φ_m/A_phot and θ_i are the colatitude values of the cell edges in the θ-direction. Here, the photospheric area of the domain is given by

$$\begin{eqnarray}&&{A}_{\mathrm{phot}}={R}_{\odot }^{2}(\cos (a)-\cos (b))\times (d-c).\end{eqnarray} \tag{ 110 }$$

After adding ${A}_{\phi }^{{{\rm{\Phi }}}_{M}}$ to A_ϕ, the vector potential preserves any radial net flux that is included with the observed radial magnetic field data. If zero net radial flux is desired, one can simply set Φ_M = 0 on input to ahpot_ss.

Once the vector potential has been computed with ahpot_ss, it can be used to compute all three components of the potential magnetic field by using subroutine curlahpot_ss.

The outputs from this subroutine are B_θ, B_ϕ, and B_r, with each component of ${\boldsymbol{B}}$ computed at the corresponding face centers of each voxel: B_θ is computed at the θ face centers, B_ϕ is computed at the ϕ face centers, and B_r is computed at radial face centers. B_θ is computed from radial derivatives of A_ϕ, B_ϕ from radial derivatives of A_θ, and B_r using subroutine curlh_ce_ss acting on A_θ and A_ϕ. The dimensions of these arrays are (m + 1, n, p) for B_θ, (m, n + 1, p) for B_ϕ, and (m, n, p + 1) for B_r.

We also have the ability to compute the magnetic field components directly from P and Φ_M, if desired, with subroutines brpot_ss and bhpot_ss.

Subroutine brpot_ss does an outer loop over radial index q. For each radial layer, ghost zones are added to P to make the solution consistent with the applied boundary conditions on the θ and ϕ edges of the domain. Then, the pair of subroutines curl_psi_rhat_ce and curlh_ce_ss are called in succession to compute the horizontal Laplacian of P, which then results in the values of B_r within that radial layer. An additional term is then added to the solution, ${B}_{0}\times {R}_{\odot }^{2}/{r}_{q}^{2}$, where, as before, B₀ = Φ_M/A_phot, to account for any net radial magnetic flux. The radial magnetic field component lies in the center of radial voxel faces.

Subroutine bhpot_ss does an outer loop over the radial index q, and first differences P between two adjacent levels in r to evaluate ∂P/∂r, after ghost zones have been added to each of the two layers. This derivative is evaluated at voxel centers in radius and also in θ and ϕ. Then, both B_θ and B_ϕ are evaluated by using subroutine gradh_ce_ss to take the horizontal gradient of ∂P/∂r. Here, a net radial magnetic flux plays no role in the result and so is ignored. The output arrays B_θ and B_ϕ are computed at θ and ϕ face centers, respectively.

The array dimensions of B_θ, B_ϕ, and B_r when computed by brpot_ss and bhpot_ss are identical to those computed by curlahpot_ss. The values of all three magnetic field components computed using the two different methods agree with each other to a high degree of accuracy, with an error level just slightly worse than roundoff error.

One defect of the staggered grid formalism used here is that the horizontal magnetic field components computed with the model at the lowest radial layer do not lie on the photosphere, where we have the magnetic field measurements, but instead lie half a voxel above the photosphere. We would like to compare and contrast horizontal magnetic field components from the data with their potential field counterparts lying on the photosphere. Fortunately, we can use the finite-difference form of Bercik's equation to infer the photospheric values of the horizontal magnetic field components. Equation (91) shows that

$$\begin{eqnarray}&&{B}_{r}^{\mathrm{phot}}=\displaystyle \frac{{\left(\partial P/\partial r\right)}_{\tfrac{1}{2}{\rm{\Delta }}r}-{\left(\partial P/\partial r\right)}_{-\tfrac{1}{2}{\rm{\Delta }}r}}{{\rm{\Delta }}r},\end{eqnarray} \tag{ 111 }$$

where ${(\partial P/\partial r)}_{-\tfrac{1}{2}{\rm{\Delta }}r}$ would be the ghost zone value for ∂P/∂r just below the photosphere. We know B_r at the photosphere from the data, and from the solution for P, we can difference P to find ${(\partial P/\partial r)}_{\tfrac{1}{2}{\rm{\Delta }}r}$, so we can solve for the ghost zone value below the photosphere:

$$\begin{eqnarray}&&{\left(\partial P/\partial r\right)}_{-\tfrac{1}{2}{\rm{\Delta }}r}={\left(\partial P/\partial r\right)}_{\tfrac{1}{2}{\rm{\Delta }}r}-{\rm{\Delta }}{{rB}}_{r}^{\mathrm{phot}}.\end{eqnarray} \tag{ 112 }$$

Once this has been done, we can then interpolate (average) ∂P/∂r between values at $r=-\tfrac{1}{2}{\rm{\Delta }}r$ and $r=+\tfrac{1}{2}{\rm{\Delta }}r$ to get the photospheric value of ∂P/∂r,

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\partial P}{\partial r}\right)}_{\mathrm{phot}}={\left(\displaystyle \frac{\partial P}{\partial r}\right)}_{\tfrac{1}{2}{\rm{\Delta }}r}-\left(\displaystyle \frac{{\rm{\Delta }}r}{2}\right){B}_{r}^{\mathrm{phot}}.\end{eqnarray} \tag{ 113 }$$

Then, the horizontal gradient of this quantity yields the potential field values of B_θ and B_ϕ at the photosphere. These operations are carried out by subroutine bhpot_phot_ss.

The subroutine uses the solution P computed by scrbpot_ss and observed photospheric values of B_r on input and computes B_θ and B_ϕ at the photosphere on output. B_θ lies along the TE grid, and B_ϕ lies along the PE grid. These arrays can be compared directly to the staggered grid values of the observed data, for comparisons of the differences and similarities between the observations and what the potential field model predicts. Figures 10–13 show test results of the potential field software, using data from a vector magnetogram of AR 11158 taken on 2011 February 15, 22:47 UT. The values of m, n, and p for this calculation are 632, 654, and 1000, respectively. The value of the assumed source surface was R_SS = 2 R_⊙. Here, the homogeneous Neumann boundary condition in ϕ was used to compute the solution. On an Apple MacBook Pro (early 2015), with 16 GB of RAM and an SSD disk drive, these solutions can be derived and written to disk on timescales of roughly 10 minutes, using a single processor, with the compute time noticeably shorter for periodic boundary conditions in ϕ as compared to homogeneous Neumann boundary conditions. The compute time is dominated by subroutine scrbpot_ss. Once the solution for P has been obtained, evaluating the magnetic field components takes a small fraction of the total compute time. In the horizontal directions, the angular resolution is close to that of HMI; in the radial direction, Δr is roughly 700 km, about twice the horizontal spacing as that at the photosphere.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While our potential field software was designed for deriving solutions on AR-sized portions of the Sun, it can also be used for deriving solutions that lie above a very large fraction of the solar disk. First, the range of ϕ can be extended to the entire circumference of the Sun by simply choosing c = 0 and d = 2π and then choosing the periodic boundary condition option in ϕ (bcn=0) when calling scrbpot_ss. Second, we have tested the software by choosing very small values of a and values of b that approach π, with no major ill effects or artifacts near the poles. In particular, we have chosen a and b such that their values differ from 0 and π by only 001 without difficulty. We then compared the morphology of the solutions at various radii computed with the spherical-harmonic-based PFSS model of D. J. Bercik & J. G. Luhmann (2019, in preparation), using moderate numbers for maximum spherical harmonic degree, with solutions from the software in PDFI_SS described here, using compatible resolution for the finite-difference equations presented here. The solutions seemed compatible overall at several different radii between the photosphere and the source surface.

If, instead of specifying the radial magnetic field at the photosphere, one specifies the partial time derivative of the radial magnetic field at the photosphere, subroutine scrbpot_ss will find $\dot{P}$ instead of P. In that case, if one then calls subroutine ahpot_ss with $\dot{P}$ as input, the output will be the electric field components cE_θ and cE_ϕ (both with a minus sign) throughout the coronal volume. These are the electric fields that correspond to the time derivative of the corresponding potential magnetic fields in the volume. Calling subroutine curlahpot_ss will then compute the time derivative of all the magnetic field components in the volume defined by the potential field software. If homogeneous Neumann boundary conditions in ϕ are chosen (bcn=3), these solutions will be compatible with the PTD solution for $c{{\boldsymbol{E}}}_{h}$ at the photosphere computed from ptdsolve_ss and e_ptd_ss, apart from the minus sign. It then becomes possible to perform detailed investigations of how the electric field corresponding to the changing potential magnetic field distributions behaves in the coronal volume.

Electric fields computed in this way appear to be generally consistent with one of the contributions to the electric field ${{\boldsymbol{E}}}_{P}$ identified with the changing potential magnetic field in the analysis of Schuck & Antiochos (2019). They identify that electric field as coming from a solenoidal contribution, denoted ${{\boldsymbol{\Sigma }}}_{P}$, plus that from the gradient of a scalar potential, denoted ${\boldsymbol{\nabla }}{{\rm{\Lambda }}}_{P}$. The calculation described above will compute a solenoidal (inductive) version of ${{\boldsymbol{\Sigma }}}_{P}$. However, looking at trial test cases indicates that along the $\hat{{\boldsymbol{\phi }}}$ and $\hat{{\boldsymbol{\theta }}}$ normal faces of the spherical wedge volume, the component of the electric field normal to these surfaces is not zero, in contrast to the boundary condition (32c) assumed for ${{\boldsymbol{\Sigma }}}_{P}$ in Schuck & Antiochos (2019). This is true for either bcn=3 or bcn=0, neither of which constrains the behavior of the component of ${\boldsymbol{E}}$ normal to the faces. Thus, the PDFI_SS solutions do not appear to be consistent with the assumed boundary conditions for ${{\boldsymbol{\Sigma }}}_{P}$ in Schuck & Antiochos (2019).

Apart from this, we have not yet pursued any further detailed studies using this possible application of the PDFI_SS potential field software at this time, but simply point out this possibility for future work.

Once the 3D distributions of the magnetic fields have been computed, it is straightforward to estimate the energy in the potential magnetic field, by either performing an integral of B²/(8π) over the computational volume or estimating this quantity through a photospheric surface integral, using Gauss's theorem, and ignoring side and top boundaries. We have written three subroutines to provide such estimates. These subroutines are emagpot_ss, emagpot_srf_ss, and emagpot_psi_ss.

Subroutine emagpot_ss takes as input the three 3D arrays B_θ, B_ϕ, and B_r; interpolates the magnetic field components from voxel faces to the voxel centers; evaluates B²/(8π) at the center of each voxel; and then sums up the magnetic energy density from each individual voxel. The advantage of this subroutine is that no assumptions about side-wall boundary conditions are made; a possible disadvantage is that magnetic field energy outside the volume is not computed and therefore underestimated. Energies are computed in units of ergs.

Subroutine emagpot_srf_ss uses the fact that the volume integral of B²/(8π) can also be written, using Gauss's theorem, as an area integral of $(1/(8\pi )){\boldsymbol{A}}\times {\boldsymbol{B}}\cdot \hat{{\boldsymbol{n}}}$ over all the surfaces surrounding the volume, where $\hat{{\boldsymbol{n}}}$ is the outward surface normal vector for each surface. However, the subroutine ignores all the surfaces except the photosphere and simply integrates $(-1/(8\pi ))\hat{{\boldsymbol{r}}}\cdot ({{\boldsymbol{A}}}_{h}\times {{\boldsymbol{B}}}_{h})$ over the photospheric domain. On input, the subroutine takes the photospheric values of A_θ, A_ϕ, and the potential field values (not the observed values!) of B_θ and B_ϕ as, e.g., computed from bhpot_phot_ss. We find that emagpot_srf_ss tends to overestimate magnetic energies to a modest degree when compared to the results from the volumetric integral computed by emagpot_ss. Most likely this is due to the effects of ignoring all the nonphotospheric surfaces in the area integral.

A third subroutine for computing the magnetic energy is emagpot_psi_ss. Here, the photospheric values of the potential magnetic field B_r and the scalar potential Ψ are used on input to compute the magnetic energy by integrating ΨB_r/(8π) over the photospheric surface. This equation also results from a use of Gauss's theorem, when ${\boldsymbol{B}}$ is expressed as minus the gradient of the scalar potential Ψ. With our solutions derived in terms of P, this is less convenient to use than emagpot_srf_ss, since Ψ = −∂P/∂r is normally not evaluated at the photosphere, but half a voxel above it. However, Equation (113) shows a strategy to evaluate Ψ = −∂P/∂r at the photosphere, which is implemented in PDFI_SS by calling subroutine psipot_phot_ss.

We also find that emag_psi_ss tends to overestimate magnetic energies compared to the volumetric integral computed by emagpot_ss, probably for similar reasons to emagpot_srf_ss.

The potential magnetic field solutions are computed in colatitude–longitude order for computational purposes within the PDFI_SS software, but for display purposes, as well as other applications that use longitude–latitude array orientation, we want an efficient capability to transform the 3D solutions from colatitude–longitude orientation to longitude–latitude orientation.

We have written several subroutines to perform these transpose operations, ahpottp2ll_ss, bhpottp2ll_ss, and brpottp2ll_ss.

The subroutine ahpottp2ll_ss converts the two components of the vector potential from colatitude–longitude order to longitude–latitude order and also changes the sign from A_θ to A_lat. The subroutine bhpottp2ll_ss does the same operation on the horizonal components of the potential magnetic field, also changing the sign of B_θ when converting to B_lat. The subroutine brpottp2ll_ss transposes the radial magnetic field array, B_r. It can also be used to transpose the poloidal potential itself, P.

Since this is a very memory-intensive operation, in the potential field documentation file Potential-Fields-Spherical.txt located in the fossil repository http://cgem.ssl.berkeley.edu/cgi-bin/cgem/PDFI_SS in the doc folder, there is a discussion of how one can use a combination of C and FORTRAN pointers to perform the transpose operations "in place" if desired, which uses significantly less memory. No change to the existing source code for the transpose subroutines in the PDFI_SS library is necessary to do this; this memory-sharing operation is done entirely in the calling program. This same principle is also outlined in Section 9.4, which describes a test program for doing these transpose operations.

Welsch & Fisher (2016) showed an alternative method for deriving potential magnetic fields from vector magnetogram data, where, instead of matching B_r at the photosphere, one could match ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$ as measured from the data. They found that these solutions could result in quite different values of B_r at the photosphere as compared to the observations, just as potential field models based on B_r can have horizontal magnetic fields that differ considerably from the observed values (see, e.g., Figures 11 and 12). Welsch & Fisher (2016) found that potential field solutions matching ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$ can have substantially smaller magnetic energies than those matching B_r. We have implemented a technique for finding potential field solutions that match ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$ using much of the same potential field framework described above. We now outline this technique, which is included in the PDFI_SS software.

We first note that relating the poloidal potential P to ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$ in a way that can use BLKTRI is not as straightforward as it was for relating P to B_r. However, if we solve for the scalar potential Ψ instead of P, then much of the mathematical and numerical framework we use in scrbpot_ss can be adapted to solve for Ψ.

Writing ${\boldsymbol{B}}=-{\boldsymbol{\nabla }}{\rm{\Psi }}$ and then taking the horizontal divergence of ${\boldsymbol{B}}$ at the photosphere, we have

$$\begin{eqnarray}&&{R}_{\odot }^{2}{{\rm{\nabla }}}_{h}^{2}{\rm{\Psi }}=-{R}_{\odot }^{2}{{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h},\end{eqnarray} \tag{ 114 }$$

where in the volume above the photosphere Ψ obeys the Laplace equation

$$\begin{eqnarray}&&{r}^{2}{{\rm{\nabla }}}_{h}^{2}{\rm{\Psi }}+\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\displaystyle \frac{\partial {\rm{\Psi }}}{\partial r}\right)=0.\end{eqnarray} \tag{ 115 }$$

Equation (114) is of exactly the same form as Equation (108), but involving Ψ and ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$ instead of P and B_r. Equation (115) has a somewhat different radial term than does Bercik's Equation (88), but it is compatible with the use of BLKTRI. The source surface boundary condition at r = R_SS will differ between P and Ψ. We now describe the details of how the equation for Ψ is solved, particularly where the details differ from those in Section 6.1.

6.7.1. Finite-difference Expressions for the Laplace Equation for Ψ and the Solution Procedure

Our strategy for solving the Laplace Equation (115) is identical to that for solving Bercik's equation. We will convert the azimuthal, colatitude, and radial derivatives to finite differences and then Fourier transform the equations in the azimuthal direction and derive finite-difference equations for the amplitude of each Fourier mode as a function of θ and r. Using the same notation of Section 6.1.2, the operators L_ϕ and L_θ, when converted to finite-difference form, will be identical to the operators in Section 6.1.2, where we also assume homogeneous Neumann boundary conditions at θ = a and θ = b.

As in Section 6.1.2, we write the solution Ψ as

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{i+\tfrac{1}{2},j+\tfrac{1}{2},q}=\sum _{j^{\prime} =0}^{n-1}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}{{\rm{\Phi }}}_{j^{\prime} }({\phi }_{j+\tfrac{1}{2}}^{{\prime} }),\end{eqnarray} \tag{ 116 }$$

where ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}$ is the amplitude of the coefficient of ${{\rm{\Phi }}}_{j^{\prime} }$ as a function of colatitude and radius indices. The actions of the L_ϕ and L_θ operators on ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}$ are identical to those in Section 6.1.2.

Since the L_r operator differs from that in Section 6.1.2, we write down the result:

$$\begin{eqnarray}&&\begin{array}{l}{L}_{r}({Q}_{j^{\prime} }^{i+\tfrac{1}{2},q})=\displaystyle \frac{{\left({r}_{q}-\tfrac{1}{2}{\rm{\Delta }}r\right)}^{2}}{{\rm{\Delta }}{r}^{2}}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q-1}\\ \ \ \ -\ \displaystyle \frac{{\left({r}_{q}-\tfrac{1}{2}{\rm{\Delta }}r\right)}^{2}+{\left({r}_{q}+\tfrac{1}{2}{\rm{\Delta }}r\right)}^{2}}{{\rm{\Delta }}{r}^{2}}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q}\\ \ \ \ +\ \displaystyle \frac{{\left({r}_{q}+\tfrac{1}{2}{\rm{\Delta }}r\right)}^{2}}{{\rm{\Delta }}{r}^{2}}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},q\,+\,1},\end{array}\end{eqnarray} \tag{ 117 }$$

for values of q that are in the interior of the problem. For the outermost radial position q = p, we want to impose the boundary condition that ${\rm{\Psi }}\to 0$ as $r\to {R}_{\mathrm{SS}}$, so that ${{\boldsymbol{B}}}_{h}\to 0$. This means we have a ghost zone value of ${Q}_{j^{\prime} }^{i+\tfrac{1}{2},p+1}=-{Q}_{j^{\prime} }^{i+\tfrac{1}{2},p}$. Since BLKTRI does not use ghost zones, we make this substitution into Equation (117), resulting in

$$\begin{eqnarray}&&\begin{array}{l}{L}_{r}({Q}_{j^{\prime} }^{i+\tfrac{1}{2},p})=\displaystyle \frac{{\left({r}_{p}-\tfrac{1}{2}{\rm{\Delta }}r\right)}^{2}}{{\rm{\Delta }}{r}^{2}}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},p-1}\\ \ \ -\ \displaystyle \frac{{\left({r}_{p}-\tfrac{1}{2}{\rm{\Delta }}r\right)}^{2}\,+\,2{\left({r}_{p}+\tfrac{1}{2}{\rm{\Delta }}r\right)}^{2}}{{\rm{\Delta }}{r}^{2}}{Q}_{j^{\prime} }^{i+\tfrac{1}{2},p}.\end{array}\end{eqnarray} \tag{ 118 }$$

From Equations (117) and (118) we can easily determine the values of the arrays an, bn, and cn for radial positions above the photosphere in subroutine BLKTRI.

As in Section 6.1.5, the first (photospheric) values of the an, bn, and cn arrays are set to zero.

The solution procedure for the finite-difference equations for Ψ is otherwise identical to that described in Section 6.1 for P. Getting the solution for Ψ is carried out by subroutine psipot_ss.

To test the accuracy with which Laplace's equation is obeyed by Ψ, we have written the subroutine laplacetest_ss, which outputs separately the horizontal and radial contributions to the Laplacian. We found that this was useful in debugging psipot_ss.

Finally, there is an issue that the solution to the Laplace equation can contain a spurious artifact in Ψ that is proportional to r⁻¹. The origin of this artifact appears to be an interaction between the source surface boundary condition for Ψ, which essentially sets Ψ = 0 at r = R_SS, and the homogeneous Neumann boundary conditions in θ for the k = 0 mode at the photosphere, when BLKTRI is called for this particular mode. This, in effect, allows Ψ at the photosphere to "float," i.e., to have an arbitrary constant added to it. At radii in between the photosphere and source surface, the solution is connected by an r⁻¹ dependence that satisfies the Laplace equation. This solution, while mathematically legitimate, has no physical basis and results in a sometimes large, horizontally uniform radial component B_r that must be removed from the solution. We have written a subroutine psi_fix_ss, which evaluates and removes this artifact, by evaluating the B_r term at the photosphere and then removing it from the entire solution volume. This subroutine is called just before exiting subroutine psipot_ss. Subroutine psi_fix_ss can also be used to impose an observed nonzero net radial flux to the solution for Ψ, if desired.

6.7.2. Getting the Potential Magnetic Field Components from Ψ: Subroutine bpot_psi_ss

In principle, once Ψ is known, the magnetic field components can be found by simply taking minus the gradient of Ψ. In practice, there is a major challenge to deriving the magnetic field components from Ψ: gradients of Ψ put the magnetic field values in a different location in the grid from where we need it. We now describe how we evaluate the magnetic field components at the grid locations where we need them.

The scalar potential Ψ lies at the centers of the radial faces of our voxels. In the horizontal directions, Ψ lies on the CE grid (see Figure 9). The grid location of Ψ in the horizontal directions is different from the grid locations of the scalar potentials ψ for computing the electric field contributions at the photosphere, where these various scalar potentials all lie on the COE grid. This is why we have used the notation of the uppercase Ψ to distinguish the scalar potential for the potential magnetic field derived from ${{\boldsymbol{B}}}_{h}$ from the notation of the lowercase ψ contributions to the photospheric electric field.

When we take horizontal gradients of Ψ, the results lie on the θ and ϕ edges of the radial faces, but we need them at midpoints in radius (midway in r between radial faces) at the centers of the horizontal faces of the voxels in θ and ϕ). Similarly, when we take the radial component of the gradient of Ψ to get B_r, it is evaluated at midpoints in r within the voxels, but we need B_r on radial face centers.

We now describe how we interpolate the horizontal magnetic field components from the θ and ϕ edges of radial faces to the θ and ϕ face centers of our voxels.

If we imagine that we have another set of voxels ("offset voxels") that are offset by $\tfrac{1}{2}{\rm{\Delta }}r$ from our grid voxels, then we also imagine that each of our voxels contains the upper half of an offset voxel in the bottom half of our given voxel, and the bottom half of the next highest offset voxel in the top half of our voxel. We want the θ and ϕ magnetic fluxes from our voxel to match the flux from the top half of the lower offset voxel, plus the flux from the bottom half our the upper offset voxel. These considerations result in the following expression for the interpolated horizontal magnetic field components:

$$\begin{eqnarray}\begin{array}{c}\begin{array}{c}\begin{array}{rcl}{B}_{\theta }^{i,j+\tfrac{1}{2},q+\tfrac{1}{2}} & = & \tfrac{1}{{r}_{q+1}^{2}-{r}_{q}^{2}}\times \left(({r}_{q+\tfrac{1}{2}}^{2}-{r}_{q}^{2}\right){B}_{\theta }^{i,j+\tfrac{1}{2},q}\\ & & +\,\left.\left({r}_{q+1}^{2}-{r}_{q+\tfrac{1}{2}}^{2}\right){B}_{\theta }^{i,j+\tfrac{1}{2},q+1}\right)\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 119 }$$

and

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{B}_{\phi }^{i+{\textstyle \tfrac{1}{2}},j,q+{\textstyle \tfrac{1}{2}}} & = & {\textstyle \tfrac{1}{{r}_{q+1}^{2}-{r}_{q}^{2}}}\times \left((,{r}_{q+{\textstyle \tfrac{1}{2}}}^{2}-{r}_{q}^{2}\right){B}_{\phi }^{i+{\textstyle \tfrac{1}{2}},j,q}\\ & & +\,\left.\left({r}_{q+1}^{2}-{r}_{q+{\textstyle \tfrac{1}{2}}}^{2}\right){B}_{\phi }^{i+{\textstyle \tfrac{1}{2}},j,q+1}\right),\end{array}\end{array}\end{eqnarray} \tag{ 120 }$$

where we use the fact that the area of the side faces of a voxel are proportional to ${r}_{q+1}^{2}-{r}_{q}^{2}$ if the bottom and top radial faces of the voxel are located at r_q and r_q+1.

Now that we have values of B_θ and B_ϕ interpolated to horizontal face centers, we can evaluate ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$ at voxel centers, using subroutine divh_ce_ss, where the result projected onto the horizontal directions lies on the CE grid. We can then use the constraint ${\boldsymbol{\nabla }}\cdot {\boldsymbol{B}}=0$ to derive the radial derivative of r²B_r:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial r}({r}^{2}{B}_{r})=-{r}^{2}{{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}.\end{eqnarray} \tag{ 121 }$$

Since B_r evaluated from $-{\boldsymbol{\nabla }}{\rm{\Psi }}$ is also colocated with ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$, we can use Equation (121) to extrapolate B_r to the upper and lower radial faces, where we want the values. The evaluation of ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$ and the extrapolation to radial faces are done within subroutine br_voxels3d_ss.

All of these tasks are accomplished within subroutine bpot_psi_ss.

If a nonzero net radial flux is desired for the potential field, one can specify its value in the input variable mflux in the call to bpot_psi_ss.

If one wants to compute a solution for Ψ itself that is consistent with an imposed net radial flux, one can call subroutine psi_fix_ss with a nonzero value of the net radial flux mflux. Doing this is advisable if using Ψ to compute magnetic energies with subroutine emagpot_psi_ss.

6.7.3. Testing Potential Field Models that Use B_h at Photosphere

How can we characterize the accuracy of solutions to the potential field models that use photospheric values of ${{\boldsymbol{B}}}_{h}$? The interpolation and extrapolation steps described in Section 6.7.2 will introduce some amount of error into the magnetic field solutions, as compared to the solutions based on B_r, where these steps are not needed. Our objective here is to provide a method for estimating errors in the solution obtained from ${{\boldsymbol{B}}}_{h}$.

The test described here is based on the following procedure: (1) First, obtain the potential field solution that matches observed photospheric values of B_r using subroutines scrbpot_ss, ahpot_ss and curlahpot_ss, along with subroutine bhpot_phot_ss to compute the horizontal potential field components at the photosphere. (2) Using the photospheric potential field components B_θ and B_ϕ computed from the above, compute the scalar potential Ψ to match ${{\boldsymbol{\nabla }}}_{h}\cdot {{\boldsymbol{B}}}_{h}$ at the photosphere by calling subroutine psipot_ss. (3) Compute the magnetic field components by calling subroutine bpot_psi_ss. (4) Compare the original solutions for B_θ, B_ϕ, and B_r with those computed from Step 3, and evaluate the discrepancies. Figure 14 shows the recovered values of B_r versus the original values of B_r from the data, showing an rms difference of ∼18 G. For comparison, the scatter plots of the recovered values of B_θ and B_ϕ (not shown) look like straight lines and have much smaller errors. Similarly, the original potential field model based on B_r shows errors of ∼10⁻⁶ G, in recovered versus observed values of B_r, close to roundoff error. Thus, the largest source of error seems to be the interpolation and extrapolation procedures in subroutine bpot_psi_ss that were needed to compute B_r at radial face centers. While these errors are visible here, they are smaller than the quoted HMI errors in B_r, so we feel that the solutions are accurate enough for many scientific studies.

Zoom In Zoom Out Reset image size

6.7.4. Applications of Potential Field Solutions from B_h

Welsch & Fisher (2016) proposed the idea that potential field models derived from vector magnetograms can include observed data from ${{\boldsymbol{B}}}_{h}$ and B_r and proposed composite models where both solutions can be used, with weights for each based on measurement errors for the different components of ${\boldsymbol{B}}$ considered separately. Because our potential field software in PDFI_SS includes the ability to compute both solutions, this can be done in a straightforward way.

We also note that Welsch & Fisher (2016) proposed that differences in B_r between the observations and the ${{\boldsymbol{B}}}_{h}$ potential field solutions may provide a diagnostic for the existence of horizontal currents. For AR 11158, we show such a difference image of B_r in Figure 15. It is interesting that in the sunspots there is only a slight difference in B_r, but there are also large-scale patterns elsewhere in the AR showing a significant difference. This potential field software makes more detailed studies a practical possibility.

Zoom In Zoom Out Reset image size

The first published study in which time-dependent vector magnetic fields were used to derive electric fields was that of Fisher et al. (2010; earlier, Mikić et al. 1999 described electric field solutions determined from time derivatives of the radial component of ${\boldsymbol{B}}$). In this case, ANMHD magnetoconvection simulation data (Welsch et al. 2007), which had known electric field solutions, were used to create a vector magnetogram data sequence, which could then be analyzed by computing PTD solutions for ${\boldsymbol{E}}$. While there was a broad resemblance between the inverted PTD solutions and the actual electric fields, there were also a number of artifacts. The authors developed an "iterative" technique to compute an additional scalar potential, whose gradient could be added to the PTD solutions to make ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$ perpendicular to each other, resulting in a moderately better agreement between the inverted and actual electric fields. A further investigation in that article used a variational approach, to impose a "smallness" constraint to the electric field solutions, which resulted in a poor match with the actual ANMHD electric fields. We subsequently gave up on using the variational approach.

The PTD Poisson equations were solved in Fisher et al. (2010) using a FORTRAN version of the Newton–Krylov technique, originally developed for the first version of RADMHD (Abbett 2007), since the required boundary conditions were inconsistent with the use of FFTs. Solutions obtained with the iterative method, which used repeated solutions of a Poisson equation, were performed in IDL using FFTs.

A great improvement in the accuracy of the electric field inversions of the ANMHD simulations was made in Fisher et al. (2012), in which it was realized that adding information about Doppler shifts, which can be measured, resulted in dramatically better solutions for the electric field. They derived Poisson equations for contributions from both Doppler shifts and horizontal flows derived from Local Correlation Tracking, which were then solved in IDL using FFTs, with the solutions added to the PTD solutions obtained with the Newton–Krylov software.

In 2011–2012, coauthors Maria Kazachenko, Brian Welsch, and George Fisher realized they needed more efficient software for solving the Poisson equations, in which many more types of boundary conditions could be applied, and which would be faster than their existing Newton–Krylov code. They tested several numerical techniques and concluded that the elliptic equation package FISHPACK was ideally suited to these tasks. They proceeded to write a very general executable program in FORTRAN, which could be spawned from IDL, which would read input data, compute the solutions using FISHPACK, and then write the solutions to a file, which could then be read back into the IDL session. This software model for PDFI existed from roughly 2012 through 2015, during which the centered, Cartesian version of the PDFI software was developed. We now refer to this version as PDFI_CC, where "CC" refers to "Cartesian centered." This is the version of the software that was used to perform the research described in Kazachenko et al. (2014, 2015).

Starting in 2013, the above coauthors received funding for the CGEM project (Fisher et al. 2015), in which they proposed to take the existing PDFI software and (1) convert it to spherical coordinates and (2) rewrite it in an efficient computer language that could be run automatically from the SDO JSOC. This process happened in several stages. First, the Cartesian IDL source code had to be converted from Cartesian to spherical coordinates. This process took roughly 6 months and maintained the use of a centered grid. (This version was called PDFI_SC, where "SC" denotes "spherical centered.") In the meantime, by studying MuRAM MHD simulation results obtained from Matthias Rempel at HAO/NCAR, which had turbulent structures at the scale of the grid, the authors realized that the centered grid finite-difference formulation was simply unable to obey Faraday's law accurately when the solutions were so highly structured. They realized they needed to convert their finite-difference equations into a conservative, staggered grid coordinate system. After investigating several different formulations of staggered grid systems, they finally arrived at the system described in Section 3.4.

Once the PDFI_SC version was written and working in IDL, the next step was to convert the IDL code from the spherical centered grid to the spherical-staggered grid scheme. This process occurred during the first half of 2015. By July of 2015, an IDL version of PDFI_SS was operational and had successfully undergone a number of tests.

To deliver the software in a form that could be run automatically at the SDO JSOC, the coauthors knew that the IDL code would have to be converted to FORTRAN, since it relies so heavily on the FISHPACK FORTRAN library. The conversion of the code from IDL to FORTRAN was done during the last half of 2015 and early 2016. Since that time, nearly all development effort has been on the FORTRAN version of the software. The FORTRAN version of PDFI_SS now contains a much broader spectrum of capabilities than the original IDL version did. While we continue to keep the IDL legacy version within the PDFI_SS developer site, we no longer actively maintain the IDL branch of the software. We do find that the existing IDL code, particularly the procedures for performing vector calculus operations, can still be useful when analyzing output from PDFI_SS.

There were a number of choices made in how the PDFI_SS FORTRAN code was written. Here we briefly comment on these and discuss the motivations for these choices.

First, the calling arguments for all the subroutines include input and output arrays, as well as other important information provided as single real scalar values or as integers. All quantities defined as arrays have their dimensions defined in terms of integer values passed into the subroutine by the user. Modern FORTRAN allows one to determine array sizes and shapes by querying the attributes of these arrays, potentially reducing the number of necessary calling arguments. However, we found that these advanced features did not work when the PDFI_SS subroutines were invoked from other C and Python software. Thus, we define all array dimensions from other calling arguments in the subroutines.

Second, we have avoided any use of "common-block" variables, or other global parameters or variables, which can obscure the dependencies of output variables on input arguments. All input data are passed explicitly as calling arguments into the FORTRAN subroutines. This constraint eases the ability to use the library from languages other than FORTRAN.

Third, all floating-point operations are performed using 64 bit reals. All reals, either scalars or arrays, are declared as real*8 variables in the source code, a choice that seems to work correctly with all FORTRAN compilers attempted thus far. All integer arguments to PDFI_SS subroutines are assumed to be default integers in FORTRAN, which are 32 bit integers.

Fourth, the source code assumes that all input and output arrays are dimensioned or allocated (and deallocated) by the user in the calling programs. This is essential for the software to be used from languages other than FORTRAN. Thus, very few of the arrays in PDFI_SS are dynamically allocated within the source code. The one exception to this rule are the work arrays needed by FISHPACK subroutines. In this case, for each PDFI_SS subroutine that calls a FISHPACK subroutine, the work array is both allocated and then deallocated within that same subroutine.

Fifth, to facilitate ease of interoperability with C code, character string arguments have been completely avoided in the PDFI_SS software. Character strings have a different representation in memory between FORTRAN and C.

Finally, the FORTRAN syntax is implemented using the older .f suffix for the source code file names, rather than the more modern .f90 suffix. While the latter choice results in more flexible syntax for, e.g., line continuation, the former choice helps enforce 80-character line limits, which makes viewing the source code much easier from the default 80-character width of a terminal window.

The first step in compiling PDFI_SS is to download, compile, and install the FISHPACK FORTRAN library. Links for the FISHPACK version 4.1 source code are given in the introduction to Section 3.

After unpacking the tarball, we recommend that you replace the contents of file make.inc in the top folder of the FISHPACK distribution, with the contents of the file fishpackmake.inc located in the doc folder in the PDFI_SS distribution, then replace the file Makefile in the top folder of the FISHPACK distribution with the contents of the file fishpackmake in the doc folder in PDFI_SS, and finally replace the file Makefile in the src folder of the FISHPACK distribution with the contents of the file fishpackmakesrc in the doc folder in PDFI_SS. You may need to edit the file make.inc to (1) make sure that the name of the FORTRAN compiler coincides with the name of the FORTRAN compiler you have. We have specifically included lines for the gfortran and Intel compilers in the Linux part of the make.inc file, and the gfortran compiler for the Mac (Darwin) portion of make.inc. If you are using another compiler, you will need to edit the compiler definition F90 so that it reflects your compiler. (2) Make sure that the options included in defining the FORTRAN compiler also ensure that all reals are set to 64 bit reals. (3) Check that compiler options for compiling position-independent code, needed if FISHPACK will be used for languages other than FORTRAN, are invoked. (4) Edit definitions for make and ar, if they are different from what is defined in this file.

We cannot overemphasize how important it is to invoke the compiler option that all reals are treated as 64 bit reals. If this is not done, the attempted use of FISHPACK with PDFI_SS is doomed to fail, in ways that are not always easy to diagnose.

To compile the FISHPACK library, type "make." Once the FISHPACK library is compiled, you can install it into a location of your choosing by typing "make install" (or "sudo make install" if this requires root priviledge). Alternatively, you will need to remember the exact path to the location of the library file libfishpack.a.

Once FISHPACK has been compiled and installed, we are ready to compile PDFI_SS. As noted earlier, the PDFI_SS software developer site is http://cgem.ssl.berkeley.edu/cgi-bin/cgem/PDFI_SS/index. By clicking on "Login," one can log in as anonymous, and then by clicking on "Files," one should find a blue hexidecimal link, the ID for the latest software release. By clicking on that link, one should then be led to links for tarball or zip archives for the software.

Once the tarball has been downloaded and unpacked, you should see three subfolders: IDL, doc, and fortran. Descend into the fortran folder with a terminal window.

The next step is to open the file "Makefile" with an ASCII text editor, such as vi or emacs. You will most likely need to edit this file before you can compile the library. Currently, the Makefile is set up assuming you will be running on either a Mac or a Linux machine of some kind, and that you have access to a FORTRAN compiler. The file assumes you will have access to either the gnu/gfortran compiler or the Intel compiler, ifort. If you plan to use a different FORTRAN compiler, you will need to edit Makefile to add the name of that compiler and to add compiler options for it that coincide with the meanings of the compiler options for gfortran or ifort.

To compile the PDFI_SS library file, libpdfi_ss.a, type "make." To install the library into a specified location, edit the definition of INSTALLDIR, and then type "make install" (or possibly "sudo make install," if you need root privilege for the specified location). The default value of INSTALLDIR is /usr/local/lib.

Once the PDFI_SS library has been installed, linking to other FORTRAN programs is straightforward. For the gfortran and ifort compilers, linking to the library is invoked with the -lpdfi_ss -lfishpack (in that order) linking commands. If the libraries are not stored in "standard" locations, you may need to specify the location of each library with the -L<dir> directive. Specific examples can be found in the test programs, described in further detail in Section 9.

If there are no character string arguments, calling a FORTRAN subroutine from a C function is very straightforward, if one just remembers some basic rules: (1) From C, a FORTRAN subroutine is a function of type void (i.e., the function returns nothing). All input and output are handled through the calling arguments. (2) The name of a FORTRAN subroutine is changed by the FORTRAN compiler ("FORTRAN name mangling"); typically this is done by adding a trailing underscore. This practice is observed by both the gfortran and ifort compilers. In other words, in C, if one wants to call ahpot_ss, the corresponding function name in C is ahpot_ss_. (3) In FORTRAN, all arguments are called by reference, not by value. This means that when calling a FORTRAN subroutine from C, all arguments must be passed by reference, i.e., as pointers. For example, if in a C function calling a FORTRAN subroutine the variables m and n are declared as integers, their pointers &m and &n would be used in the call to the subroutine. (4) FORTRAN is a column-major language. For multidimensional arrays in FORTRAN, the first index always varies in memory the fastest. For example, a 2D array brll, dimensioned (n + 1, m + 1) in FORTRAN, assumed to be in longitude–latitude orientation, is ordered such that we start with the smallest latitude value, increase the longitude index from the smallest to the maximum value, then repeat the process with the next-lowest latitude index value, etc. C is considered to be a row-major language, so that given a 2D array in C, the second index varies the fastest in memory.

From our experience, the easiest way to deal with this possible source of confusion is, first, to stick with using 1D arrays in C of length (n + 1) ∗ (m + 1) using the above example and, second, to make sure that all input 1D arrays are arranged in column-major order before calling the FORTRAN subroutine. On output, we also recommend defining 1D arrays in C, keeping in mind that the output data will be ordered by the FORTRAN subroutine into column-major order. If you need the data arranged in a different order, you will need to do that rearrangement after the subroutine call. Fortunately, 1D arrays in C map neatly onto multidimensional arrays in FORTRAN, provided that one keeps in mind the assumed column-major order. For FORTRAN arrays of three or more dimensions, the same principle works: define an array in C of length equal to the product of the FORTRAN dimensions, and make sure that the first index varies the fastest, followed by the second index, followed by the next index.

The size of default integers in FORTRAN is 32 bits, so the C calling program should be sure to not use 16 bit or 64 bit integers when calling the subroutines. For nearly all systems, a declaration of int in C should be compatible with FORTRAN integer arguments to PDFI_SS subroutines. Similarly, all real variables in PDFI_SS are 64 bit reals, compatible with the double precision (double) declaration in C. The PDFI_SS subroutines assume that the calling program has already allocated memory for both input and output arrays. All memory management for the calling arguments to PDFI_SS subroutines is assumed to be handled by the calling program and is not done within PDFI_SS itself.

We have written as one of our test programs (see Section 9.8) a simple C program that calls the brll2tp_ss subroutine. The program shows explicitly how the input array is constructed and ordered into column-major order before calling the subroutine, and when the output array is printed, one sees that the output array is also arranged in column-major order, using the transposed dimensions.

We strongly recommend defining function prototype statements for any PDFI_SS subroutines you call from a C program (in C99, these statements are required). This reduces the chance of making errors in calling the subroutine from C and can be helpful in debugging the code by warning the user when calling arguments disagree with those of the function prototype. We have written an include file (pdfi_ss.h) that contains the function prototypes in C for all of the user-callable subroutines in PDFI_SS. This file can be included in any C code that calls PDFI_SS FORTRAN suboutines. In our test C program (Section 9.8), our test program C source code includes this file.

There is little difference in calling PDFI_SS subroutines from a C++ program compared to a C program. The same rules about passing arguments (all arguments are passed by reference, i.e., as pointers) apply. The main difference is that (1) in the C++ program you will need to set the lang=C option, and (2) the compiler options for position-independent code must be invoked when compiling PDFI_SS. In our Makefile, we have endeavored to make sure this option is chosen for the gfortran and ifort compilers.

Linking the PDFI_SS library into the Python programming language allows effective use of the software with solar-physics-related open-source Python packages, such as SunPy (Mumford et al. 2015) and astropy (Astropy Collaboration et al. 2013), as well as easy manipulation, analysis, and plotting of the input and output data using the basic Python modules NumPy, SciPy (Jones et al. 2001), and matplotlib (Hunter 2007). The PDFI_SS-Python linking has also been used to implement the PDFI_SS electric field inversion into ELECTRIC field Inversion Toolkit, ELECTRICIT (Lumme et al. 2017, 2019). ELECTRICIT is an easy-to-use Python software toolkit for downloading and processing of SDO/HMI data and inverting the photospheric electric field from the data using a range of state-of-the-art methods.

We have successfully created a working Python interface for several PDFI_SS functions using the F2PY FORTRAN to Python interface generator https://docs.scipy.org/doc/numpy/f2py/, which is a part of the NumPy package. F2PY is compatible with FORTRAN 77/90/95 languages and allows partly automated creation and compilation of Python interfaces for FORTRAN routines and functions. The generator includes several methods of creating the interface, from which we have chosen to use the method based on signature files. The process has the following steps: (1) The F2PY package is used to automatically create a signature file (e.g., pdfi_ss.pyf) from the FORTRAN source code. The signature file specifies the Python wrapping of the PDFI_SS routines of interest (e.g., pdfi_wrapper4jsoc_ss). (2) The automatically created signature file is then modified to ensure working wrapping of the FORTRAN routines (usually only modest changes are required). (3) Finally, F2PY is used to compile an extension module (e.g., pdfi_ss.so) from the modified signature file and FORTRAN source code and/or compiled libraries. The extension module and its functions are then importable and callable in Python (import pdfi_ss, output=pdfi_ss.pdfi_wrapper4jsoc_ss(arg1,arg2,...)).

We have retained the original IDL procedures that we used in the early phases of the development of the PDFI_SS library, although this software is no longer maintained. We find that the software is sometimes useful in the analysis of magnetic and electric field data generated by the library.

In this version of the software, there is still the need for a FORTRAN executable to solve the PDFI Poisson equations, but this executable is spawned from the IDL code when needed, and nearly all of the computational results apart from the solutions themselves are performed in IDL. The source code for the FORTRAN executable xpoisson is contained within the file poisson_arguments_stag.f and is compiled and installed with the Makefile that is in the IDL folder. The xpoisson executable is a very general wrapper for the FISHPACK subroutines HWSCRT, HSTCRT, HWSSSP, and HSTSSP and allows one to select either Cartesian or spherical coordinates and either centered or staggered grid solutions. The xpoisson executable does extensive error checking on all the input parameters for the FISHPACK subroutines before solving the Helmholtz or Poisson equation. To communicate the input data to xpoisson, and to read the output solutions from xpoisson, the "Simple Data Format" or sdf binary data format is used, and this library must be compiled and installed before xpoisson can be compiled and run.

The sdf library is written in C but is designed to be used from either C or FORTRAN. The objective of the sdf library is to read and write binary files containing both simple variables and large arrays, by calling simple subroutines or functions from FORTRAN or C. It was developed to aid in the debugging of numerical codes, making it easy to output and examine the contents of large arrays. The sdf library also has a set of IDL procedures to read and write sdf files, making this a convenient way of communicating between an IDL session and the FORTRAN executable. Coauthor Fisher developed and maintains the sdf library. The source code for sdf can be downloaded from http://solarmuri.ssl.berkeley.edu/~fisher/public/software/SDF/. Use the latest version. An archive of the latest version of this software at the time this article was published can also be downloaded from Zenodo (Fisher et al. 2020a).

The PDFI_SS IDL source code contains the core abilities to compute the PTD and FLCT terms, the relaxation procedure (needed for the Doppler and ideal contributions) to the PDFI electric field, but lacks much of the additional capabilities of the FORTRAN library. Computing solutions with the IDL code is also much more time-consuming than using the FORTRAN library software. Nevertheless, we sometimes find that the vector calculus procedures, which have nearly the same names as the corresponding FORTRAN subroutines, can be useful in analyzing the results from PDFI_SS solutions.

The source code in test_wrapper.f is designed to mimic the JSOC's call to subroutine pdfi_wrapper4jsoc_ss. In a nutshell, it reads in test magnetogram and Doppler data from HMI, along with stored FLCT estimates of horizontal flows, then adds padding to the data in a manner consistent with how this process is done upstream of the PDFI_SS call by the JSOC, and then calls the subroutine pdfi_wrapper4jsoc_ss. The output from the subroutine is written to an output file. The test program also independently computes each of the four electric field contributions and also writes these to the output file, so that a detailed and independent comparison can be made. The program computes and writes to the output file many other diagnostic quantities.

We must caution that the input data used here are taken from a preliminary test data series for NOAA AR 11158 generated several years ago and do not reflect a number of improvements to the data analysis that have occurred since that time. The FLCT flow velocities, in particular, were generated with nonoptimal parameter choices. Nevertheless, since these data are fixed here for the purpose of testing the PDFI_SS code, not the data analysis procedures, we have retained their use in this test program.

The documentation at the front of test_wrapper.f includes a list of the variables from the output file, if read into an IDL session with the procedure sdf_read_varlist. Here, we will discuss only a summary of the overall PDFI_SS electric field results for this particular test case. Particularly important diagnostics one can examine with the output data include comparison of the curl of ${\boldsymbol{E}}$ with the temporal difference in magnetic field components between the two adjacent vector magnetic field measurements.

One can use quantities in the output file to examine a detailed breakdown of the PDFI solutions into their four contributions, which are computed independently within test_wrapper.f. For further details, see the documentation at the head of the test_wrapper.f source code file.

We end our discussion of test_wrapper.f by displaying in a series of grayscale figures (Figures 16–22) the three magnetic field components for the test data, along with computed PDFI electric field inversions for the three components of ${\boldsymbol{E}}$. Both the inductive and total contributions to E_r are shown. We also show in Figure 23 the radial component of the Poynting flux computed for the pair of magnetograms. These figures provide an overall picture for the electric field and Poynting flux morphology that can be compared to the magnetic field components for context.

The purpose of the test_anmhd.f program is to use the PDFI_SS software to compute the electric fields for the ANMHD test case using vector magnetic field and Doppler data from a horizontal slice of an ANMHD simulation of magnetoconvection (Welsch et al. 2007; Kazachenko et al. 2014) and then compare that solution with the $-{\boldsymbol{V}}\times {\boldsymbol{B}}$ electric field computed from the simulation itself. This provides a good independent test of the PDFI_SS solution technique, and it can also be compared with the results obtained by Kazachenko et al. (2014) in Section 4 of that article using PDFI solutions that assume a centered grid formalism in Cartesian coordinates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because the ANMHD simulation was performed in Cartesian coordinates, the first task is to use the formalism described in Section 7 to map the Cartesian domain onto a small surface patch bisected by the equator on a very large sphere. The resulting radius of the sphere in this case is 9.998 × 10⁸ km, well over 1000 times larger than R_⊙. The colatitude range b–a = 10⁻⁴ rad, and it is also equal to the longitude range d–c, since the domain is a square in Cartesian coordinates.

There is a subroutine in the PDFI_SS library, pdfi_wrapper4anmhd_ss, which closely mimics the functionality of subroutine pdfi_wrapper4jsoc_ss, but with some differences needed to accommodate this special case (for example, no "zero padding" is done by the latter subroutine, since padding was also not done in Kazachenko et al. 2014). test_anmhd.f calls this subroutine and then writes the results to an output file. When the resulting electric field solutions are compared with those from the ANMHD simulation itself, our use of the staggered grid means that the comparison is a little more complicated than it was in Kazachenko et al. (2014). First, so that we can compare quantities directly, we must interpolate the simulation magnetic and electric fields from a centered grid (which the simulation used) to our staggered grid locations; this step was not necessary for the comparison in Section 4 of Kazachenko et al. (2014). Second, because the magnetic field time derivatives were not masked in Kazachenko et al. (2014), we also do not mask them here (in contrast to what is done with HMI data in pdfi_wrapper4jsoc_ss). Similarly, we also did not mask the FLCT electric field or the ideal electric field. However, we found that if we did not mask the Doppler contribution to the electric field, noise in the definition of the $\hat{{\boldsymbol{q}}}$ unit vector in the weak-field regions would wreak havoc on the solutions; therefore, we did use the strong-field mask on the Doppler electric field solutions. The threshold for the strong magnetic field mask is 370 G, chosen to be consistent with the threshold used in Kazachenko et al. (2014).

The output file anmhd_output_file_pdfi_ss.sdf from test_anmhd.f can be read into an IDL session with sdf_read_varlist. There is an extensive amount of diagnostic data that can be analyzed, as detailed in the documentation near the front of the file test_anmhd.f. Here, we display just a few aspects of this output data, where the results can be compared with those of Kazachenko et al. (2014). Figures 24–26 show side-by-side comparisons of the longitudinal, latitudinal, and radial electric field images, with the ANMHD simulation results on the left side of the figures, while the PDFI_SS inversion results are shown on the right side of the figures. Both the ANMHD simulation results and the PDFI_SS results have been multiplied by the strong magnetic field masks, as was done for similar figures in Kazachenko et al. (2014). Figure 27 shows the side-by-side comparison of the radial Poynting flux. Finally, Figure 28 shows a scatterplot of the radial component of the Poynting flux from the inversion versus that from the ANMHD simulation and provides a good indication of the resulting error levels from the inversion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall, while we find that the ANMHD electric field results are recovered well, we find that the quality of the inversion is not as good as it was for the centered grid case used in Kazachenko et al. (2014). There are a number of possible reasons for this, including the fact that the simulation data must be interpolated to the staggered grid locations and the fact that in the PDFI_SS inversions Faraday's law is obeyed to roundoff error, whereas in the simulation data it is not, as one can see in the bottom right panel of Figure 1 of Welsch et al. (2007). The latter is a consequence of the ANMHD simulations being run with spectral techniques used to compute spatial derivatives, whereas the curl of the simulation data was computed using finite differences. In spite of these differences, these figures show clearly that the PDFI_SS technique is able to reproduce the main morphological features and amplitudes of the ANMHD electric fields.

The potential field software described in Section 6 is tested by test_bpot.f, in which the radial magnetic field on the CE grid at the photosphere is used to compute a potential magnetic field distribution in the volume above the photosphere. In the test program, the angular resolution of the solution is the same as that for the photospheric data. In radius, the test program assumes 1000 voxels in the radial direction, with a source surface height of 2 R_⊙ (i.e., one solar radius above the photosphere). The user can choose whether to assume periodic boundary conditions in ϕ by setting the variable bcn to 0 or homogeneous Neumann boundary conditions in ϕ on the poloidal potential P by setting bcn to 3. The test program takes at least several minutes to run (about 10 minutes on a MacBook Pro with 16 GB of memory and a solid-state disk). With limited amounts of memory, it will likely take considerably longer.

The executable xbpot produces an output file, test_bpot_output.sdf, with a series of 2D and 3D arrays. The file can be read in using IDL procedure sdf_read_varlist. Reading in the file can take a long time, as the output file is very large. Once the file is read in, one can type the "help" command in IDL to see the list of output arrays. The solution for the poloidal potential P is stored in the 3D array scrb3d. The three magnetic field components for the solution are computed two different ways: First, subroutines brpot_ss and bhpot_ss are called to compute B_r and ${{\boldsymbol{B}}}_{h}$ using the poloidal potential P as input. The resulting 3D magnetic field arrays are brpot3d, btpot3d, and bppot3d. Second, subroutine ahpot_ss is used to compute the vector potential ${{\boldsymbol{A}}}^{P}$ from P. The theta and phi components of the vector potential are returned into the 3D arrays atpot and appot. Then, the three magnetic field components can be computed from the vector potential with subroutine curlahpot_ss. The resulting three 3D magnetic field arrays are brpotvp, btpotvp, and bppotvp. We find that the magnetic field components computed in the two different ways differ very little.

It is important to note that the potential field solution is computed on the faces of the voxels, as described in Section 6, and that the horizontal magnetic field components at the bottom layer in the 3D arrays are located half a voxel above the photosphere. The photospheric values of the potential field components are computed with subroutine bhpot_phot_ss and output into the 2D arrays btpotphot and bppotphot.

Note that all of these output arrays are stored in colatitude–longitude-radius index order. The following test program will take the output file from this test program and rotate some of the 3D arrays into longitude–latitude-radius index order.

The test program test_bptrans.f uses the 3D transpose subroutines ahpottp2ll_ss, bhpottp2ll_ss, and brpottp2ll_ss to transpose the 3D arrays computed by test_bpot.f from colatitude–longitude index order to longitude–latitude index order. The program uses the output file from test_bpot.f as an input file and writes transposed output arrays into the file test-inplace-transpose.sdf. The contents of the file can be read into an IDL session using sdf_read_varlist. The arrays in the file include the 3D arrays alon and alat, the longitudinal and latitudinal components of the vector potential, respectively, and blonpot, blatpot, and brllpot, the 3D arrays of the longitudinal, latitudinal, and radial components of ${{\boldsymbol{B}}}^{P}$, respectively. The transposed 2D photospheric magnetic field components are in the arrays blonphot and blatphot.

When working with large 3D arrays such as these, we used a memory-saving technique in FORTRAN that is worth describing. Instead of having to create two separate arrays of the same size, but with different shapes, it is possible, using a combination of C and FORTRAN pointers, to create the two separate arrays such that they occupy the same locations in memory but still have different shapes. The two different array names can then be successfully passed as arguments into our PDFI_SS 3D transpose subroutines. We will now illustrate how we do this using the radial magnetic field component of the potential field solution as an example.

The ability to combine the usage of C and FORTRAN pointers can be invoked by the declaration use, intrinsic :: iso_c_binding

within a FORTRAN program. The original array brpot in colatitude–longitude index order is declared as a third-rank allocatable array: real*8, allocatable,target :: brpot(:,:,:). The transposed array, brllpot, is declared as a third-rank FORTRAN pointer: real*8, pointer :: brllpot(:,:,:).

We also define the two shape arrays integer :: shapebrpot(3), shll(3).

Then, one defines and reads in the integers m,n,p. The array brpot is then allocated as allocate brpot(m,n,p+1).

Once the array is allocated, its shape can be determined with the FORTRAN shape function: shapebrpot=shape(brpot).

The shape of the transpose array can then be determined by the following three statements: shll(1)=shapebrpot(2) shll(2)=shapebrpot(1) shll(3)=shapebrpot(3).

The next step is to read in the brpot array from the input file. Once that is done, we can assign the FORTRAN pointer for the transposed array: call c_f_pointer(c_loc(brpot),brllpot,shll).

What this statement does is define the brllpot array to occupy the same location in memory as the brpot array, but with the shape of the transpose array. One can then call the PDFI_SS subroutine brpottp2ll_ss, with the brpot array as input and brllpot as the output array: call brpottp2ll_ss(m,n,p,brpot,brllpot).

The source code for brpottp2ll_ss is written in such a way that the transpose is done sequentially for each horizontal layer of the 3D arrays, with a copy of the 2D array made before any transpose operation has been done. Thus, the 3D transpose can be done "in place," without destroying any data. It is very important to note that once the 3D transpose subroutine has been called, the original untransposed array is essentially destroyed by scrambling it into the new shape and is hence useless. Thus, the 3D transpose operation should be the last operation done using the original array.

This concept is used for all the 3D transpose operations in test_bptrans.f. The contents of the output file were used to generate the potential field figures shown in Section 6.2.

The purpose of the test_psipot.f test program is to test the ability to compute the potential magnetic field by using the observed horizontal components of the magnetic field at the photosphere as input. These solutions are computed by subroutines psipot_ss and bpot_psi_ss, as described in Section 6.7. On input, the test program reads in arrays of the observed horizontal magnetic fields at the photosphere at the staggered grid locations from file test-fortran-input.sdf, and on output it computes the scalar potential Ψ and the 3D magnetic field components, all in colatitude–longitude-radius index order. The results of the potential field calculation are written to the output file test-psipot-output.sdf.

The solution for the scalar potential Ψ is returned as two separate 3D arrays, psi3d, which assumes zero net radial flux, and psi3df, in which a net radial flux is imposed with subroutine psi_fix_ss that coincides with the net radial flux from the observed radial components of the magnetic field. The potential magnetic field components are computed with bpot_psi_ss and are given in the arrays btpot3d, bppot3d, and brpot3d. Photospheric values of the horizontal potential magnetic field components are given in the 2D arrays btpot and bppot. The horizontal divergences of the observed photospheric horizontal magnetic fields are given in the 2D array divbh, and the divergences of the potential photospheric horizontal magnetic fields are given in the 2D array divbhpot.

As an alternative to the observed horizontal magnetic field components in test-fortran-input.sdf, one can instead change the input file name to test-fortran-inputpotphot.sdf, where in this case the horizontal components of ${\boldsymbol{B}}$ are computed from the potential magnetic field solution chosen to match B_r. In this case, the radial magnetic field component computed by test_psipot.f should be close to the observed radial magnetic field component. The error between the two was shown earlier in Figure 14.

One can choose to use periodic boundary conditions in ϕ by setting the variable bcn to 0 in test_psipot.f, or homogeneous Neumann boundary conditions by setting bcn to 3. For this test case, homogeneous Neumann boundary conditions require more compute time.

The test program test_global.f is designed to test the ability of subroutine enudge3d_gl_ll to compute a global (4π sr) PTD solution covering the entire surface of the Sun. On input, arrays of ${\dot{B}}_{\mathrm{lon}}$, ${\dot{B}}_{\mathrm{lat}}$, and ${\dot{B}}_{r}$ defined on the PE, TE, and CE grids, respectively, spanning the global Sun, are computed from two sets of vector magnetogram data, using a temporal finite difference to define the time derivative. In this case, we have used AR 11158 vector magnetogram data and mapped them onto the global Sun geometry, setting a = 0, b = π, c = 0, and d = 2π, i.e., we have allowed AR 11158 to take over the entire surface of the Sun. Because in a global geometry we cannot have a nonzero net radial magnetic flux, or its time derivative, we use subroutine fluxbal_ll to zero out the average before calling enudge3d_gl_ll, which computes ${\boldsymbol{E}}$ on all the rails of the spherical voxels, which are bisected in radius by the photosphere.

To test whether this global solution for ${\boldsymbol{E}}$ obeys all three components of Faraday's law, we compute the curl of ${\boldsymbol{E}}$ with subroutine curle3d_ll, which can accommodate both spherical wedge and global geometries. The resulting components of the curl of ${\boldsymbol{E}}$ can then be compared against the input time derivatives ${\dot{B}}_{\mathrm{lon}}$, ${\dot{B}}_{\mathrm{lat}}$, and ${\dot{B}}_{r}$. The documentation near the top of test_global.f provides the names of the appropriate arrays for the time derivatives and the three components of the curl of ${\boldsymbol{E}}$, as well as the arrays for ${\boldsymbol{E}}$ itself. Scatter plots of the curl of ${\boldsymbol{E}}$ versus the magnetic field time derivatives should show straight lines for each component. The output file for test_global.f is global_test_out_ar11158.sdf, and the contents of the file can be read in with IDL procedure sdf_read_varlist.

In Section 3.13, we described several subroutines designed to use a B-spline to interpolate the input data for computing PDFI electric field inversions to a different resolution. Here, the test program test_interp.f performs a test of this B-spline interpolation, by first interpolating an observed HMI image of B_r, which includes regions of zero padding, to a new resolution that is about 20% higher than the original resolution, and then it reinterpolates that image back to the original resolution. The original B_r data array can then be compared directly to the twice-interpolated array of the same size, and the quality of the twice-interpolated image can be evaluated. The output file, brint.sdf, contains the original image as the array brdat, the interpolation to the higher resolution as the array brint, and the twice-interpolated array as brback. The test program uses the subroutine interp_hmidata_ll to perform the interpolations. The default value of the degree of the spline, deg, is set to 9, but the user can experiment with other values. Allowed values are 3, 5, 7, and 9. Thus far we have found that setting deg=9 seems to produce the most accurate results. The output file can be read in with the IDL procedure sdf_read_varlist to study the results.

In Section 8.5, we described the principles of calling PDFI_SS subroutines from C/C++ programs. In the test program test_pdfi_c.c, we illustrate many of the important points made in that section with this very simple test program, which calls a single PDFI_SS subroutine, brll2tp_ss. This test program (1) illustrates how to define a 1D array in C that maps onto 2D arrays in FORTRAN, (2) shows how to define the input array in column-major order, consistent with its use in FORTRAN, (3) shows that calling arguments are all called by reference, and (4) uses the output array to demonstrate the column-major nature of the output generated by the FORTRAN subroutine.

In the C test program, pointers for the arrays brll (the input array) and br (the output array) are defined and initialized to NULL. Next, integer variables np1 and mp1 (representing n + 1 and m + 1, respectively) are defined and set to 12 and 10, respectively. The integers m and n are then defined by decrementing np1 and mp1 by one. Next, the input and output arrays brll and br are each allocated to have the size (n + 1) ∗ (m + 1), i.e., the product of the FORTRAN dimensions. Next, the array brll is defined so that regarded as a FORTRAN array, the value of brll(j,i)=(i-1)+(j-1). This is done by using an outer loop over the latitude index i that goes from 0 to m and an inner loop over the longitude index j that goes from 0 to n. The calculation of the array value is brll[i∗np1+j]=(double) i + j;. This is the essence of constructing the array in C to have column-major order, before the FORTRAN subroutine is called. Note that the j index is the most rapidly varying. (The values of i and j differ between C and FORTRAN because in FORTRAN the default index values start from 1 whereas in C they start from 0.)

Next, the FORTRAN subroutine brll2tp_ss is called: brll2tp_ss_ (&m, &n, brll, br);. Note the trailing underscore in the subroutine name and the fact that the pointers to m and n are used in the call (the arrays brll and br are already defined as pointers).

Following the subroutine call, both the brll and br arrays are printed to the screen. In the loop that prints the values of brll, the outer loop is over the latitude index and the inner loop is over the longitude index, consistent with column-major order in FORTRAN. When printing out the output array br to the screen, column-major ordering results in the colatitude index i varying most rapidly, and the longitude index j varies more slowly: the outer loop is over j, and the inner loop is over i. The FORTRAN array br(i,j) is indexed in C as br[j∗mp1+i].

These same principles in setting up 1D arrays in C that map onto multidimensional FORTRAN arrays should work for any of the subroutines in PDFI_SS, as long as one pays careful attention to the dimensions of the arrays defined in the FORTRAN subroutines and remembers that FORTRAN assumes that the arrays are defined in column-major index order.

When running the executable xctest, the array brll is printed to the screen, followed by its colatitude–longitude transpose, br.

We have prepared an example of the PDFI_SS-Python linking described in Section 8.6 in the subfolder fortran/test-programs/python-linking. The package contains a working signature file pdfi_ss.pyf created for routines add_padding_as_ss, pad_abcd_as_ss, pad_int_gen_ss, and pdfi_wrapper4jsoc_ss in the PDFI_SS library. The F2PY extension module pdfi_ss.so can be compiled by typing "make" in the terminal while within the subfolder python-linking. Typing "make clean" removes the extension module. The fishpack library location must be correctly specified in the Makefile, and the FISHPACK library must be compiled with the -fPIC flag (see discussion in Section 8.3).

The Python script pdfi_wrapper4jsoc_script.py imports the Python interfaces of the PDFI_SS subroutines from the compiled extension module and calls them to process the input data and to estimate the electric field by calling the pdfi_wrapper4jsoc_ss subroutine (Section 3.11). Thus, the script reproduces the first part of the test program test_wrapper.f (Section 9.1). The input data required by the pdfi_wrapper4jsoc_ss can be downloaded in Python-compatible .sav format from the URL described in the introduction to Section 9. In pdfi_wrapper4jsoc_script.py, the output arrays of pdfi_wrapper4jsoc_ss are returned as NumPy arrays, and the script saves them to a NumPy .npz file. This output can be then compared to the output of the pure FORTRAN version of pdfi_wrapper4jsoc_ss executed within test_wrapper.f. Example output files created by executing test_wrapper.f can be downloaded from the above URL in Python-compatible .sav format, and their content can be compared to the output of the pdfi_wrapper4jsoc_script.py using the compare_wrapper_outputs.py script. The differences printed by the script should be small and close to floating-point precision.

If the user wishes to do the comparison from scratch on his/her local system, the package includes also IDL helper scripts for transforming the .sdf input and output files of test_wrapper.f program into Python-compatible .sav format. The README.txt file in this subfolder contains also step-by-step instructions for creating the F2PY interface from scratch.