Hyperbolic Cell-centered Finite Volume Method for Obtaining Potential Magnetic Field Solutions

Coronal magnetic fields play a significant role in the existence and variability of a wide variety of solar activities ranging from slowly evolving features such as sunspots (Solanki 2003), solar prominences (Mackay et al. 2010), and coronal loops (Reale 2010), to highly dynamic phenomena such as solar flares (Benz 2017) and coronal mass ejections (CMEs) (Forbes et al. 2006; Chen 2011). These solar activities are closely connected with some adverse phenomena in the Earth's environment that can affect human activities, space, and ground based technologies. For instance, the collapse of the Hydro-Quebec power grid in 1989 was caused by a geomagnetic storm triggered by a CME (Allen et al. 1989). Therefore, the study on the coronal magnetic field strength and topology is one of major goals of solar physics. Currently, coronal magnetic fields cannot be directly measured. Their quantitative information has mainly been provided by some models extrapolated from accessible measurements. The major models currently in use for deducing the global magnetic structures of solar corona include the potential field source surface (PFSS) model, nonlinear force-free field model, and multidimensional MHD model. Among these models, the PFSS model is basic in that it can not only reproduce the large-scale geometry of the solar coronal magnetic fields but is also simple to develop and implement (Riley et al. 2006; Arden et al. 2014; Wiegelmann et al. 2014). Other models are sophisticated and numerically expensive. For the application of the PFSS model, it can provide an initial magnetic field distribution for coronal and solar wind MHD simulations. Besides, it is also an important component of the Wang–Sheeley–Arge model that is regularly used to predict the background solar wind speed and the interplanetary magnetic field (IMF) polarity from photospheric observation (Arge et al. 2003, 2004; Sun & Hoeksema 2007; Yang et al. 2018). These two parameters, velocity and IMF polarity, can be applied to forecast the geomagnetic Ap index (Detman & Vassiliadis 1997). In a word, the PFSS model is a very useful tool in the space physics research field.

The PFSS model was first developed in the late 1960s by Schatten et al. (1969) and Altschuler & Newkirk (1969). It begins with the synoptic magnetograms of the Sun and extrapolates the magnetic field to a sphere called the "source surface" where the coronal magnetic field is purely radial. Usually, the source surface height R_ss is taken to be 2.5 R_s, where R_s is the solar radius. In this model, a key assumption is that there are no currents in the region $1{R}_{s}\leqslant r\leqslant {R}_{{ss}}$. Mathematically the problem is described as follows: given the magnetogram data ${B}_{r}(\theta ,\phi )$ at $r=1{R}_{s}$, find the scalar potential Φ so that

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {\rm{\Phi }}}{\partial r}\right|}_{r=1{R}_{s}}={B}_{r}(\theta ,\phi ),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{\Phi }}{| }_{r={R}_{{ss}}}=0.\end{eqnarray} \tag{ 3 }$$

Here, θ ∈ [0, π] and ϕ ∈ [0, 2π] are the co-latitude and longitude coordinates, respectively. Once the solution of Equation (1) is found, the potential magnetic field solution is obtained as

$$\begin{eqnarray*}&&{\boldsymbol{B}}={\rm{\nabla }}{\rm{\Phi }},\end{eqnarray*}$$

which trivially satisfies both the divergence-free and the current-free properties

$$\begin{eqnarray*}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}={\rm{\nabla }}\cdot ({\rm{\nabla }}{\rm{\Phi }})=0,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm{\nabla }}\times {\boldsymbol{B}}={\rm{\nabla }}\times ({\rm{\nabla }}{\rm{\Phi }})=0.\end{eqnarray*}$$

In order to obtain the potential magnetic field solutions to the Poisson Equation (1) with boundary conditions (2) and (3), there are currently two different approaches: the spherical harmonic expansion (SHE) approach (e.g., Altschuler et al. 1977; Mackay & Yeates 2012; Nikolić & Trichtchenko 2012) and the iterative finite difference method (e.g., Tóth et al. 2011). The SHE approach is an analytical solver. It works reasonably well when the order of spherical harmonics is limited to be small relative to the resolution of the magnetogram, although some artifacts, such as ringing, can arise around sharp features. However, when the number of spherical harmonics is increased, using the raw magnetogram data specified on a grid that is uniform in the sine of the latitude coordinate can result in inaccurate and unreliable results, especially in the polar regions close to the Sun (Tóth et al. 2011). To mitigate this problem, one choice is to remesh the magnetogram onto a uniform latitude–longitude grid, which can limit the highest order of the spherical harmonics to the alias-free limit (Suda & Takami 2002). The other approach is the iterative finite difference method, i.e., the Finite Difference Iterative Potential-field Solver (FDIPS) available at the website http://csem.engin.umich.edu/tools/FDIPS. The FDIPS can accept any magnetogram data represented in the Flexible Image Transport System (FITS) file format as input so that it can optionally remesh the magnetogram onto a grid with uniform resolution in latitude according to the way of spatial discretization. Compared with spherical harmonics, the advantage of finite differences is that the boundary data given by the magnetogram directly affect the solution only locally, while the spherical harmonics are global functions and their amplitudes depend on all of the magnetogram data (Tóth et al. 2011).

In the present study, we embark on developing a hyperbolic cell-centered finite volume method for the Poisson equation. The hyperbolic method, originally proposed by Nishikawa (2007), is mainly used to solve the diffusion equation for fluid flow. Its idea is relaxing the order of spatial derivatives by introducing auxiliary variables for the solution gradient, reformulating the second-order diffusion equation into a first-order hyperbolic system in pseudo time, and solving the reformulated hyperbolic system to get the steady-state solution, which is the true solution of the original diffusion equation. Up to now, many works on the hyperbolic method have been carried out by Nishikawa and his coworkers especially for compressible Navier–Stokes solvers (Mohammed & Ismail 2015; Liu & Nishikawa 2016; Nakashima et al. 2016), for theoretical investigation of the hyperbolic method (Mazaheri et al. 2016; Nishikawa & Roe 2016), for time-dependent problems (Mazaheri & Nishikawa 2014), for dimensional scaling and numerical similarity in hyperbolic method (Nishikawa & Nakashima 2018), for a diffusion equation with discontinuous coefficients (Nishikawa 2018), and for magnetohydrodynamics (Kawashima et al. 2015; Baty & Nishikawa 2016). Additionally, Lee et al. (2018) presented a cell-centered high-order hyperbolic finite volume method for the steady diffusion equation on unstructured grids. Based on the cell-centered finite volume numerical framework, the aim of the present paper is to apply the hyperbolic method to the Poisson equation. Hereafter, the developed algorithm is called a hyperbolic cell-centered finite volume solver (HCCFVS).

With the overarching goal of solving a Poisson equation, HCCFVS could be naturally employed for the Poisson projection method proposed by Brackbill & Barnes (1980) to eliminate the divergence error of the magnetic field during the MHD simulation. The Poisson projection method act as a correction to the magnetic field after the time step is completed by some arbitrary numerical scheme (Tóth 2000). The idea of this Poisson projection method is that the ${{\boldsymbol{B}}}^{* }$ field provided by the base scheme in time step n + 1 is projected to a divergence-free ${{\boldsymbol{B}}}^{n+1}$ field and that this corrected ${{\boldsymbol{B}}}^{n+1}$ is taken as the newly updated magnetic field. It's well known that the magnetic field can be decomposed into the sum of a curl and a gradient ${{\boldsymbol{B}}}^{* }={\rm{\nabla }}\times {\boldsymbol{A}}+{\rm{\nabla }}{\rm{\Phi }}$. Hence, taking the divergence of both sides, we would obtain a Poisson equation

$$\begin{eqnarray*}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}={\rm{\nabla }}\cdot {{\boldsymbol{B}}}^{* },\end{eqnarray*}$$

which can be solved for the scalar function Φ. Then, the magnetic field is corrected by

$$\begin{eqnarray*}&&{{\boldsymbol{B}}}^{n+1}={{\boldsymbol{B}}}^{* }-{\rm{\nabla }}{\rm{\Phi }}.\end{eqnarray*}$$

The Poisson projection method works well, but does come at the cost of solving a Poisson equation at each time step. In the present paper, HCCFVS is used to eliminate the divergence error in solar coronal MHD simulation by solving the Poisson equation involved in the projection method. Besides the Poisson projection method, there exist other built-in divergence-cleaning approaches such as the constrained transport (CT) method (Evans & Hawley 1988; Zhang & Feng 2016), diffusive method (van der Holst & Keppens 2007; Zhang & Feng 2016), hyperbolic divergence-cleaning method (Dedner et al. 2002; Zhang & Feng 2016; Li et al. 2018), constrained-gradient method (Hopkins 2016), and solenoidality-preserving approach (Feng et al. 2019). For the comparative study of the diffusive method, the Poisson projection method, hyperbolic divergence-cleaning method, and CT method in solar coronal simulation, we can refer to Zhang & Feng (2016).

The rest of this paper is organized as follows. In Section 2, we describe the entire procedure of HCCFVS in detail. In Section 3, we test two cases with analytical solutions, a 3D Poisson problem with Dirichlet boundary and a 3D Poisson equation with Neumann boundary, to preliminarily demonstrate the feasibility and the accuracy of HCCFVS. In Section 4, we employ HCCFVS to extrapolate the solar coronal magnetic fields based on the synoptic magnetograms, and compare the resulting large-scale geometry of the solar coronal magnetic fields with that obtained from the spherical harmonic approach and FDIPS. In Section 5, by using the potential magnetic fields from HCCFVS and the SHE as initial inputs, we present a comparative study on the steady-state solar corona during Carrington rotation (CR) 2098. Finally, we draw our conclusions in Section 6.

In this section, the general framework for HCCFVS is described, including the mesh grid system, the hyperbolic formulation for Poisson equation and the discretization in the cell-centered finite volume scheme, the reconstruction for variables, and the implicit solver for the residual equation.

2.1. Mesh Grid System

Following Feng et al. (2010, 2017, 2019), the computational domain is partitioned into the six-component grid system, as shown in Figure 1. Such a six-component grid consists of six identical component meshes with partial overlapping areas. Each component can be treated separately as a low-latitude spherical mesh,

$$\begin{eqnarray*}&&\left(\displaystyle \frac{\pi }{4}-\delta \leqslant \theta \leqslant \displaystyle \frac{3\pi }{4}+\delta \right)\bigcap \left(\displaystyle \frac{3\pi }{4}-\delta \leqslant \phi \leqslant \displaystyle \frac{5\pi }{4}+\delta \right),\end{eqnarray*}$$

where δ is proportionally dependent on the grid spacing entailed for the minimum overlapping area and is taken to be 3Δθ. The six components have the same features, and they can be transformed into each other by coordinate transformation such that all numerical assignments are identical on each component.

Zoom In Zoom Out Reset image size

Each component is divided in the spherical coordinates with the same procedure. In both $\theta$- and $\phi \mbox{-} \mathrm{directions}$, grid points are uniformly distributed and started from $\tfrac{\pi }{4}-\delta$ and $\tfrac{3\pi }{4}-\delta$, respectively: ${\theta }_{i}={\theta }_{\min }+(i-1){\rm{\Delta }}\theta ,i=-2,\cdots ,{N}_{\theta }+2$, ${\phi }_{i}\,={\phi }_{\min }+(i-1){\rm{\Delta }}\phi$, $i=-2,\cdots ,{N}_{\phi }+2$, with Δθ = $({\theta }_{\max }\,-{\theta }_{\min })/({N}_{\theta }-2)$, ${\rm{\Delta }}\phi =({\phi }_{\max }-{\phi }_{\min })/({N}_{\phi }-2)$, where N_θ and N_ϕ are the mesh numbers in latitudinal and longitudinal directions, respectively. ${\theta }_{\min }=\tfrac{\pi }{4},{\theta }_{\max }\,=\tfrac{3\pi }{4},{\phi }_{\min }=\tfrac{3\pi }{4},{\phi }_{\max }=\tfrac{5\pi }{4}$, and N_θ = N_ϕ = 92. In the $r \mbox{-} \mathrm{direction},r(1)=1,r(i+1)=r(i)\,+{\rm{\Delta }}r(i)$, where $i\,=1,\cdots ,{N}_{r},{\rm{\Delta }}r(i)=0.01$ if $r(i)\lt 1.1;$ ${\rm{\Delta }}r(i)\,=\min (A\times \mathrm{lg}(r(i-1)),{\rm{\Delta }}\theta \times r(i-1))$ with $A=0.01/\mathrm{lg}(1.09)$ if $r(i)\lt 3.5$. After the grid mesh partition in the spherical coordinates, the small volumes surrounded by the spherical grid points obtained in the spherical coordinates are interpreted as the corresponding control volume elements in the Cartesian coordinates. In this grid system, every grid cell is a hexahedron like a quadrangular frustum pyramid cell.

2.2. Hyperbolic Formulation for Poisson Equation and Discretization

On the basis of the hyperbolic method (Nishikawa 2007), the solution to the Poisson equation,

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {y}^{2}}+\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}=f(x,y,z),\end{eqnarray} \tag{ 4 }$$

where Φ is a scalar solution variable and f denotes a source term, can be transformed to the steady-state solution to the following pseudo-time-dependent first-order system form:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial \tau } & = & \displaystyle \frac{\partial {{\rm{\Phi }}}_{1}}{\partial x}+\displaystyle \frac{\partial {{\rm{\Phi }}}_{2}}{\partial y}+\displaystyle \frac{\partial {{\rm{\Phi }}}_{3}}{\partial z}-f(x,y,z),\\ \displaystyle \frac{\partial {{\rm{\Phi }}}_{1}}{\partial \tau } & = & \displaystyle \frac{1}{{T}_{r}}\left(\displaystyle \frac{\partial {\rm{\Phi }}}{\partial x}-{{\rm{\Phi }}}_{1}\right),\\ \displaystyle \frac{\partial {{\rm{\Phi }}}_{2}}{\partial \tau } & = & \displaystyle \frac{1}{{T}_{r}}\left(\displaystyle \frac{\partial {\rm{\Phi }}}{\partial y}-{{\rm{\Phi }}}_{2}\right),\\ \displaystyle \frac{\partial {{\rm{\Phi }}}_{3}}{\partial \tau } & = & \displaystyle \frac{1}{{T}_{r}}\left(\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}-{{\rm{\Phi }}}_{3}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

where τ is a pseudo time. T_r, the relaxation time, is a free parameter introduced solely to accelerate the convergence to the pseudo steady state. The hyperbolic method is initially employed to compute the steady-state solution of the diffusion equation, u_t = νΔu, where ν is a positive diffusion coefficient, by solving an equivalent first-order hyperbolic system. As Nishikawa (2007) argued, theoretically speaking, ν does not affect the steady-state solution of the first-order hyperbolic system, but it influences the transient solution. Moreover, the dependency on ν can be eliminated numerically if the entire right-hand side of the first-order hyperbolic system is proportional to ν. Therefore, Nishikawa (2007) set ${T}_{r}={L}_{r}^{2}/\nu$, where L_r is a length scale introduced to keep dimensional consistency, and Nishikawa (2014) set ${L}_{r}=\tfrac{1}{2\pi }$, which could ensure that Fourier modes are propagated to maximize the effect of error propagation. Here, following Nishikawa (2007, 2014), we define T_r to be L_r² and ${L}_{r}=\tfrac{1}{2\pi }$ for simplicity. In Equation (5), $({{\rm{\Phi }}}_{1},{{\rm{\Phi }}}_{2},{{\rm{\Phi }}}_{3})=\left(\tfrac{\partial {\rm{\Phi }}}{\partial x},\tfrac{\partial {\rm{\Phi }}}{\partial y},\tfrac{\partial {\rm{\Phi }}}{\partial z}\right)$ are the auxiliary variables representing the solution gradients introduced to make Equation (4) first order.

Obviously, Equation (5) can be written as follows,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{Q}}}{\partial \tau }+\displaystyle \frac{\partial {\boldsymbol{f}}}{\partial x}+\displaystyle \frac{\partial {\boldsymbol{g}}}{\partial y}+\displaystyle \frac{\partial {\boldsymbol{h}}}{\partial z}={\boldsymbol{S}},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{\boldsymbol{Q}} & = & \left(\begin{array}{c}{\rm{\Phi }}\\ {{\rm{\Phi }}}_{1}\\ {{\rm{\Phi }}}_{2}\\ {{\rm{\Phi }}}_{3}\end{array}\right),\ {\boldsymbol{f}}=\left(\begin{array}{c}-{{\rm{\Phi }}}_{1}\\ -\displaystyle \frac{1}{{T}_{r}}{\rm{\Phi }}\\ 0\\ 0\end{array}\right),\ {\boldsymbol{g}}=\left(\begin{array}{c}-{{\rm{\Phi }}}_{2}\\ 0\\ -\displaystyle \frac{1}{{T}_{r}}{\rm{\Phi }}\\ 0\end{array}\right),\\ {\boldsymbol{h}} & = & \left(\begin{array}{c}-{{\rm{\Phi }}}_{3}\\ 0\\ 0\\ -\displaystyle \frac{1}{{T}_{r}}{\rm{\Phi }}\end{array}\right),\ {\boldsymbol{S}}=\left(\begin{array}{c}-f(x,y,z)\\ -\displaystyle \frac{1}{{T}_{r}}{{\rm{\Phi }}}_{1}\\ -\displaystyle \frac{1}{{T}_{r}}{{\rm{\Phi }}}_{2}\\ -\displaystyle \frac{1}{{T}_{r}}{{\rm{\Phi }}}_{3}\end{array}\right).\end{array}\end{eqnarray*}$$

Numerical solutions of Equation (6) are sought by applying a finite volume spatial discretization procedure in conjunction with a hybrid polynomial reconstruction (Lee et al. 2018) and a Riemann-solver based flux function. Therefore, the set of ordinary differential equations resulting from the application of a finite volume formulation to Equation (6) for cell i surrounded by cell j is given by

$$\begin{eqnarray*}&&\displaystyle \frac{d{\overline{{\boldsymbol{Q}}}}_{i}}{d\tau }{V}_{i}=-\displaystyle \sum _{j=1}^{{N}_{\mathrm{faces}}}{{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R}){\rm{\Delta }}{A}_{{ij}}+{\overline{{\boldsymbol{S}}}}_{i}{V}_{i}={{\boldsymbol{R}}}_{i}({\boldsymbol{Q}}),\end{eqnarray*}$$

where ${\overline{{\boldsymbol{Q}}}}_{i}=\tfrac{1}{{V}_{i}}{\int }_{{{ \mathcal V }}_{i}}{\boldsymbol{Q}}\ {dV}$ and ${\overline{{\boldsymbol{S}}}}_{i}=\tfrac{1}{{V}_{i}}{\int }_{{{ \mathcal V }}_{i}}{\boldsymbol{S}}\ {dV}$ are, respectively, the cell-averaged values of Q and S in cell i, V_i is the volume of cell $i,{\rm{\Delta }}{A}_{{ij}}$ is the area of the interface ij that separates cells i and j, N_faces is the number of faces surrounding cell $i,{{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})$ stands for the numerical flux in the direction normal to the interface ij, the subscripts L and R indicate the values evaluated at the interface ij by the reconstruction polynomials in cells i and j, respectively, and R_i(Q) denotes the solution residual. Here, Roe solver is employed and thus ${{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})$ is computed as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R}) & = & \displaystyle \frac{{\boldsymbol{F}}({{\boldsymbol{Q}}}_{L})\cdot {{\boldsymbol{n}}}_{{ij}}+{\boldsymbol{F}}({{\boldsymbol{Q}}}_{R})\cdot {{\boldsymbol{n}}}_{{ij}}}{2}\\ & & -\displaystyle \frac{| {\boldsymbol{A}}({{\boldsymbol{n}}}_{{ij}})| ({{\boldsymbol{Q}}}_{R}-{{\boldsymbol{Q}}}_{L})}{2},\end{array}\end{eqnarray} \tag{ 7 }$$

where F = (f, g, h), and ${{\boldsymbol{n}}}_{{ij}}={({n}_{x},{n}_{y},{n}_{z})}^{T}$ is the unit normal to the interface pointing from cell i to cell j. By using the eigen-decomposition of the flux Jacobian $\tfrac{\partial ({\boldsymbol{F}}\cdot {\boldsymbol{n}})}{\partial {\boldsymbol{Q}}}$, the absolute flux Jacobian is given by

$$\begin{eqnarray*}| {\boldsymbol{A}}({{\boldsymbol{n}}}_{{ij}})| =\displaystyle \frac{1}{{L}_{r}}\left(\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & {n}_{x}^{2} & {n}_{x}{n}_{y} & {n}_{x}{n}_{z}\\ 0 & {n}_{x}{n}_{y} & {n}_{y}^{2} & {n}_{y}{n}_{z}\\ 0 & {n}_{x}{n}_{z} & {n}_{y}{n}_{z} & {n}_{z}^{2}\end{array}\right).\end{eqnarray*}$$

2.3. State Variable Reconstructions

In this section, we describe the variable reconstructions. Following Lee et al. (2018), we apply the $P2+P1$ hybrid polynomial reconstruction to our current study, that is, the solution variable, Φ, is represented by a quadratic polynomial and the gradient variables, Φ₁, Φ₂, Φ₃, are reconstructed by linear polynomials. In this hybrid polynomial reconstruction, the least-squares procedures are performed only for the gradient variables, and the gradient information for the solution variable is directly recycled from the gradient variables. Obviously, the hybrid reconstruction method demands less work than its routine reconstruction counterpart because it avoids gradient computation for the solution variable.

In the cell-centered finite volume framework, using the least-squares reconstruction method we first obtain the reconstruction polynomial of the gradient variable in cell i as follows,

$$\begin{eqnarray}\begin{array}{rcl}{u}_{i}({\boldsymbol{x}}) & = & {\overline{u}}_{i}+{\left.\displaystyle \frac{\partial u}{\partial x}\right|}_{i}(x-{x}_{i})+{\left.\displaystyle \frac{\partial u}{\partial y}\right|}_{i}(y-{y}_{i})\\ & & +\ {\left.\displaystyle \frac{\partial u}{\partial z}\right|}_{i}(z-{z}_{i}),\end{array}\end{eqnarray} \tag{ 8 }$$

where $u\in \{{{\rm{\Phi }}}_{1},{{\rm{\Phi }}}_{2},{{\rm{\Phi }}}_{3}\}$, ${\overline{u}}_{i}=\tfrac{1}{{V}_{i}}{\int }_{{{ \mathcal V }}_{i}}u({\boldsymbol{x}})\ {dV}$ is the cell-averaged value of u in cell i, x = (x, y, z) is the coordinate vector, (x_i, y_i, z_i) are the coordinates of the centroid of cell i. Here we use ${\left.\tfrac{\partial u}{\partial x}\right|}_{i},{\left.\tfrac{\partial u}{\partial y}\right|}_{i}$ and ${\left.\tfrac{\partial u}{\partial z}\right|}_{i}$ to denote the point value of the first-order partial derivatives at the centroid of cell i. Up to now, we complete the linear reconstruction for the gradient variables by obtaining ${\left.\tfrac{\partial u}{\partial x}\right|}_{i},{\left.\tfrac{\partial u}{\partial y}\right|}_{i}$ and ${\left.\tfrac{\partial u}{\partial z}\right|}_{i}$ with the least-squares technique. Since $u{| }_{i},{\left.\tfrac{\partial u}{\partial x}\right|}_{i},{\left.\tfrac{\partial u}{\partial y}\right|}_{i}$ and ${\left.\tfrac{\partial u}{\partial z}\right|}_{i}$ (u ∈ {Φ₁, Φ₂, Φ₃} ) are available from the reconstruction procedure (8) for the gradient variables, then the reconstructed solution variable, Φ, follows readily

$$\begin{eqnarray*}\begin{array}{rcl}{{\rm{\Phi }}}_{i}({\boldsymbol{x}}) & = & {\overline{{\rm{\Phi }}}}_{i}+{{\rm{\Phi }}}_{1}{| }_{i}(x-{x}_{i})+{{\rm{\Phi }}}_{2}{| }_{i}(y-{y}_{i})+{{\rm{\Phi }}}_{3}{| }_{i}(z-{z}_{i})\\ & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{1}}{\partial x}\right|}_{i}\left[{\left(x-{x}_{i}\right)}^{2}-{\left(\overline{{x}^{2}{y}^{0}{z}^{0}}\right)}_{i}\right]\\ & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{1}}{\partial y}\right|}_{i}\left[(y-{y}_{i})(x-{x}_{i})-{\left(\overline{{x}^{1}{y}^{1}{z}^{0}}\right)}_{i}\right]\\ & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{1}}{\partial z}\right|}_{i}\left[(z-{z}_{i})(x-{x}_{i})-{\left(\overline{{x}^{1}{y}^{0}{z}^{1}}\right)}_{i}\right]\\ & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{2}}{\partial x}\right|}_{i}\left[(x-{x}_{i})(y-{y}_{i})-{\left(\overline{{x}^{1}{y}^{1}{z}^{0}}\right)}_{i}\right]\\ & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{2}}{\partial y}\right|}_{i}\left[{\left(y-{y}_{i}\right)}^{2}-{\left(\overline{{x}^{0}{y}^{2}{z}^{0}}\right)}_{i}\right]\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl} & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{2}}{\partial z}\right|}_{i}\left[(z-{z}_{i})(y-{y}_{i})-{\left(\overline{{x}^{0}{y}^{1}{z}^{1}}\right)}_{i}\right]\\ & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{3}}{\partial x}\right|}_{i}\left[(x-{x}_{i})(z-{z}_{i})-{\left(\overline{{x}^{1}{y}^{0}{z}^{1}}\right)}_{i}\right]\\ & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{3}}{\partial y}\right|}_{i}\left[(y-{y}_{i})(z-{z}_{i})-{\left(\overline{{x}^{0}{y}^{1}{z}^{1}}\right)}_{i}\right]\\ & & +{\left.\displaystyle \frac{1}{2}\displaystyle \frac{\partial {{\rm{\Phi }}}_{3}}{\partial z}\right|}_{i}\left[{\left(z-{z}_{i}\right)}^{2}-{\left(\overline{{x}^{0}{y}^{0}{z}^{2}}\right)}_{i}\right],\end{array}\end{eqnarray*}$$

where ${\left(\overline{{x}^{{p}_{1}}{y}^{{p}_{2}}{z}^{{p}_{3}}}\right)}_{i}=\tfrac{1}{{V}_{i}}{\int }_{{{ \mathcal V }}_{i}}{\left(x-{x}_{i}\right)}^{{p}_{1}}{\left(y-{y}_{i}\right)}^{{p}_{2}}{\left(z-{z}_{i}\right)}^{{p}_{3}}\ {dV}$ with p₁, p₂, p₃ being the integer parameters.

2.4. Implicit Solver

An implicit approach (Ivan et al. 2011) is used to find the steady-state solutions to the discrete form of Equation (6). Applying Newton's method to the residual equation, R(Q) = 0, we obtain the following linear system of equations,

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\partial {\boldsymbol{R}}}{\partial {\boldsymbol{Q}}}\right)}^{m}{\rm{\Delta }}{{\boldsymbol{Q}}}^{m}=-{\boldsymbol{R}}({{\boldsymbol{Q}}}^{m}),\end{eqnarray} \tag{ 9 }$$

for the solution change ${\rm{\Delta }}{{\boldsymbol{Q}}}^{m}={{\boldsymbol{Q}}}^{m+1}-{{\boldsymbol{Q}}}^{m}$ at Newton iteration level m. Using an initial guess, ${{\boldsymbol{Q}}}^{m}={{\boldsymbol{Q}}}^{0}$, successively improved estimates for the solution, ${{\boldsymbol{Q}}}^{m+1}$, are obtained by solving Equation (9) at each step, m, of the Newton method, where $\tfrac{\partial {\boldsymbol{R}}}{\partial {\boldsymbol{Q}}}$ is the residual Jacobian. The iterative procedure is repeated until an appropriate norm of the solution residual is sufficiently small, i.e., ${\parallel {\boldsymbol{R}}({{\boldsymbol{Q}}}^{m+1})\parallel }_{2}\lt \epsilon {\parallel {\boldsymbol{R}}({{\boldsymbol{Q}}}^{0})\parallel }_{2}$ where is a small parameter ( = 10⁻⁸), or the condition ${\parallel {\boldsymbol{R}}({{\boldsymbol{Q}}}^{m+1})\parallel }_{2}\,\gt {\parallel {\boldsymbol{R}}({{\boldsymbol{Q}}}^{m})\parallel }_{2}$ holds.

Below, we present the actual formula of the residual Jacobian in Equation (9). The diagonal block at cell i is given by

$$\begin{eqnarray*}&&{\left(\displaystyle \frac{\partial {\boldsymbol{R}}}{\partial {\boldsymbol{Q}}}\right)}_{{ii}}=\displaystyle \frac{\partial {{\boldsymbol{R}}}_{i}}{\partial {{\boldsymbol{Q}}}_{i}}=\displaystyle \frac{\partial {\overline{{\boldsymbol{S}}}}_{i}}{\partial {{\boldsymbol{Q}}}_{i}}{V}_{i}-\displaystyle \sum _{j=1}^{{N}_{\mathrm{faces}}}\displaystyle \frac{\partial {{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})}{\partial {{\boldsymbol{Q}}}_{i}}{\rm{\Delta }}{A}_{{ij}},\end{eqnarray*}$$

where

$$\begin{eqnarray*}\displaystyle \frac{\partial {\overline{{\boldsymbol{S}}}}_{i}}{\partial {{\boldsymbol{Q}}}_{i}}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & -\displaystyle \frac{1}{{T}_{r}} & 0 & 0\\ 0 & 0 & -\displaystyle \frac{1}{{T}_{r}} & 0\\ 0 & 0 & 0 & -\displaystyle \frac{1}{{T}_{r}}\end{array}\right),\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{\partial {{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})}{\partial {{\boldsymbol{Q}}}_{i}} & = & \displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {{\boldsymbol{Q}}}_{i}}\left({\boldsymbol{F}}({{\boldsymbol{Q}}}_{L})\cdot {{\boldsymbol{n}}}_{{ij}}\right)+\displaystyle \frac{1}{2}| {\boldsymbol{A}}({{\boldsymbol{n}}}_{{ij}})| \\ & = & \displaystyle \frac{1}{2}\left(\begin{array}{cccc}1/{L}_{r} & -{n}_{x} & -{n}_{y} & -{n}_{z}\\ -{n}_{x}/{L}_{r}^{2} & {n}_{x}^{2}/{L}_{r} & {n}_{x}{n}_{y}/{L}_{r} & {n}_{x}{n}_{z}/{L}_{r}\\ -{n}_{y}/{L}_{r}^{2} & {n}_{x}{n}_{y}/{L}_{r} & {n}_{y}^{2}/{L}_{r} & {n}_{y}{n}_{z}/{L}_{r}\\ -{n}_{z}/{L}_{r}^{2} & {n}_{x}{n}_{z}/{L}_{r} & {n}_{y}{n}_{z}/{L}_{r} & {n}_{z}^{2}/{L}_{r}\end{array}\right),\end{array}\end{eqnarray*}$$

and the off-diagonal block is given by

$$\begin{eqnarray*}&&{\left(\displaystyle \frac{\partial {\boldsymbol{R}}}{\partial {\boldsymbol{Q}}}\right)}_{{ij}}=\displaystyle \frac{\partial {{\boldsymbol{R}}}_{i}}{\partial {{\boldsymbol{Q}}}_{j}}=-\displaystyle \frac{\partial {{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})}{\partial {{\boldsymbol{Q}}}_{j}}{\rm{\Delta }}{A}_{{ij}},\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{\partial {{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})}{\partial {{\boldsymbol{Q}}}_{j}}=\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {{\boldsymbol{Q}}}_{j}}\left({\boldsymbol{F}}({{\boldsymbol{Q}}}_{R})\cdot {{\boldsymbol{n}}}_{{ij}}\right)-\displaystyle \frac{1}{2}| {\boldsymbol{A}}({{\boldsymbol{n}}}_{{ij}})| \\ \quad =\,\displaystyle \frac{1}{2}\left(\begin{array}{cccc}-1/{L}_{r} & -{n}_{x} & -{n}_{y} & -{n}_{z}\\ -{n}_{x}/{L}_{r}^{2} & -{n}_{x}^{2}/{L}_{r} & -{n}_{x}{n}_{y}/{L}_{r} & -{n}_{x}{n}_{z}/{L}_{r}\\ -{n}_{y}/{L}_{r}^{2} & -{n}_{x}{n}_{y}/{L}_{r} & -{n}_{y}^{2}/{L}_{r} & -{n}_{y}{n}_{z}/{L}_{r}\\ -{n}_{z}/{L}_{r}^{2} & -{n}_{x}{n}_{z}/{L}_{r} & -{n}_{y}{n}_{z}/{L}_{r} & -{n}_{z}^{2}/{L}_{r}\end{array}\right).\end{array}\end{eqnarray*}$$

Here, both Dirichlet and Neumann boundary conditions are considered, and implemented weakly by the upwind flux with the right (boundary) state specified by the boundary condition, as illustrated in Figure 2. According to Nishikawa & Nakashima (2018), the boundary state in the Dirichlet case is defined as

$$\begin{eqnarray*}&&\left\{\begin{array}{l}{{\rm{\Phi }}}_{R}=2{{\rm{\Phi }}}_{B}-{{\rm{\Phi }}}_{L}\\ {{\rm{\Phi }}}_{1,R}={{\rm{\Phi }}}_{1,L}\\ {{\rm{\Phi }}}_{2,R}={{\rm{\Phi }}}_{2,L}\\ {{\rm{\Phi }}}_{3,R}={{\rm{\Phi }}}_{3,L}\end{array}\right.,\end{eqnarray*}$$

where Φ_B is the solution given at the boundary specified by the Dirichlet condition. In this case, the numerical flux (7) becomes

$$\begin{eqnarray*}&&{{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})=\left(\begin{array}{c}-({{\rm{\Phi }}}_{1,L},{{\rm{\Phi }}}_{2,L},{{\rm{\Phi }}}_{3,L})\cdot {{\boldsymbol{n}}}_{B}-\tfrac{1}{{L}_{r}}({{\rm{\Phi }}}_{B}-{{\rm{\Phi }}}_{L})\\ -\tfrac{1}{{T}_{r}}{{\rm{\Phi }}}_{B}{n}_{x}\\ -\tfrac{1}{{T}_{r}}{{\rm{\Phi }}}_{B}{n}_{y}\\ -\tfrac{1}{{T}_{r}}{{\rm{\Phi }}}_{B}{n}_{z}\end{array}\right),\end{eqnarray*}$$

where ${{\boldsymbol{n}}}_{B}={({n}_{x},{n}_{y},{n}_{z})}^{T}$ is the unit vector normal to the boundary face. The flux Jacobian contributing to the diagonal block is given by

$$\begin{eqnarray*}\displaystyle \frac{\partial {{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})}{\partial {{\boldsymbol{Q}}}_{i}}=\left(\begin{array}{cccc}\displaystyle \frac{1}{{L}_{r}} & -{n}_{x} & -{n}_{y} & -{n}_{z}\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right).\end{eqnarray*}$$

In the Neumann case, the boundary state is defined as

$$\begin{eqnarray*}&&\left\{\begin{array}{l}{{\rm{\Phi }}}_{R}={{\rm{\Phi }}}_{L}\\ {\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{R}=2{\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{B}-{\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{L}\end{array}\right.,\end{eqnarray*}$$

where $\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}={{\rm{\Phi }}}_{1}{n}_{x}+{{\rm{\Phi }}}_{2}{n}_{y}+{{\rm{\Phi }}}_{3}{n}_{z}$ is the normal gradient given at the boundary specified by the Neumann condition. In this case, the numerical flux (7) becomes

$$\begin{eqnarray*}&&{{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})=\left(\begin{array}{c}-{\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{B}\\ -\tfrac{1}{{T}_{r}}{{\rm{\Phi }}}_{L}{n}_{x}-\tfrac{1}{{L}_{r}}{n}_{x}\left({\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{B}-{\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{L}\right)\\ -\tfrac{1}{{T}_{r}}{{\rm{\Phi }}}_{L}{n}_{y}-\tfrac{1}{{L}_{r}}{n}_{y}\left({\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{B}-{\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{L}\right)\\ -\tfrac{1}{{T}_{r}}{{\rm{\Phi }}}_{L}{n}_{z}-\tfrac{1}{{L}_{r}}{n}_{z}\left({\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{B}-{\left.\tfrac{\partial {\rm{\Phi }}}{\partial {\boldsymbol{n}}}\right|}_{L}\right)\end{array}\right),\end{eqnarray*}$$

and the flux Jacobian contributing to the diagonal block is

$$\begin{eqnarray*}\displaystyle \frac{\partial {{ \mathcal F }}_{{ij}}({{\boldsymbol{Q}}}_{L},{{\boldsymbol{Q}}}_{R})}{\partial {{\boldsymbol{Q}}}_{i}}=\displaystyle \frac{1}{{L}_{r}}\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ -\tfrac{{n}_{x}}{{L}_{r}} & {n}_{x}^{2} & {n}_{x}{n}_{y} & {n}_{x}{n}_{z}\\ -\tfrac{{n}_{y}}{{L}_{r}} & {n}_{x}{n}_{y} & {n}_{y}^{2} & {n}_{y}{n}_{z}\\ -\tfrac{{n}_{z}}{{L}_{r}} & {n}_{x}{n}_{z} & {n}_{y}{n}_{z} & {n}_{z}^{2}\end{array}\right).\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

In our computation, the linearized system (9) at each step of Newton's method is solved by the block Jacobi iteration, and the Jacobi iterative procedure is repeated until the difference between two adjacent approximations is less than some small parameter, i.e., $\parallel {\rm{\Delta }}{{\boldsymbol{Q}}}^{m,(\mathrm{Jac}+1)}-{\rm{\Delta }}{{\boldsymbol{Q}}}^{m,(\mathrm{Jac})}{\parallel }_{\infty }\lt \sigma$ where σ is the small parameter (σ = 10⁻⁸) and Jac is the Jacobi iteration counter.

In this section, two 3D test cases with exact solutions are performed to verify the effectiveness and accuracy of the developed code. In numerical simulations, the error for the solution and gradient variables is measured in terms of L₁ and L₂ norms, which are defined as

$$\begin{eqnarray*}\begin{array}{rcl}{\mathrm{Error}}_{{L}_{1}} & = & \displaystyle \frac{{\displaystyle \sum }_{i=1}^{{N}_{\mathrm{total}}}{\displaystyle \int }_{{{ \mathcal V }}_{i}}| {s}_{i}({\boldsymbol{x}})-{s}^{\mathrm{exact}}({\boldsymbol{x}})| \ {dV}}{{\displaystyle \sum }_{i=1}^{{N}_{\mathrm{total}}}{V}_{i}},\\ {\mathrm{Error}}_{{L}_{2}} & = & \sqrt{\displaystyle \frac{{\displaystyle \sum }_{i=1}^{{N}_{\mathrm{total}}}{\displaystyle \int }_{{{ \mathcal V }}_{i}}{[{s}_{i}({\boldsymbol{x}})-{s}^{\mathrm{exact}}({\boldsymbol{x}})]}^{2}\ {dV}}{{\displaystyle \sum }_{i\,=\,1}^{{N}_{\mathrm{total}}}{V}_{i}}},\end{array}\end{eqnarray*}$$

where $s\in \{{\rm{\Phi }},{{\rm{\Phi }}}_{1},{{\rm{\Phi }}}_{2},{{\rm{\Phi }}}_{3}\}$ and N_total is the total number of cells in the computational domain.

The order of accuracy is measured as follows

$$\begin{eqnarray*}&&\mathrm{Order}\ \mathrm{of}\ \mathrm{accuracy}=\displaystyle \frac{{\mathrm{log}}_{10}({E}_{\mathrm{coarse}}/{E}_{\mathrm{fine}})}{{\mathrm{log}}_{10}{\left({N}_{\mathrm{fine}}/{N}_{\mathrm{coarse}}\right)}^{1/3}},\end{eqnarray*}$$

where E_coarse and E_fine are the norm errors described above on the coarse mesh N_coarse and the fine mesh N_fine, respectively. N_coarse and N_fine are the total number of cells under the coarse mesh and the fine mesh, respectively.

3.1. Dirichlet Boundary Value Problem to 3D Poisson Equation with Manufactured Solutions

This case is a Dirichlet boundary value problem to 3D Poisson equation

$$\begin{eqnarray*}&&\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {y}^{2}}+\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}=6,\end{eqnarray*}$$

with the exact solution given by

$$\begin{eqnarray*}&&{\rm{\Phi }}(x,y,z)={x}^{2}+{y}^{2}+{z}^{2}.\end{eqnarray*}$$

The computational domain in this case is defined by inner and outer spheres of radius R_in = 1 and R_out = 3, respectively. In the actual computation, the solution is initialized by zero for all variables, i.e., the solution and gradient variables, and the value for the Dirichlet boundary condition is provided by the exact solution. The L₁ and L₂ norm errors of the numerical solution for each variable are presented in Table 1. From Table 1, it can be seen that second-order accuracy for each variable is obtained as the mesh is refined. For the solution variable, although it is represented by a quadratic reconstruction polynomial, its order of accuracy is second-order for two possible reasons: one is that the overall error in the asymptotic regime is being governed not only by the solution variable but also by the accuracy of the gradient variables (Lee et al. 2018), and the other is that one point quadrature rule is used for the flux integration on the cell interface. Figure 3 shows the comparison of the exact solution and the approximated one on the 180° − 0° plane of 54 × 62 × 62 mesh resolution per component. From this figure, we can see that the numerical solutions are very similar to the exact ones. In order to clearly show the difference between numerical and exact solutions, Figure 4 presents the relative error of them at cell centroid, $\left|\tfrac{{{\rm{\Phi }}}^{\mathrm{num}}-{{\rm{\Phi }}}^{\mathrm{exact}}}{{{\rm{\Phi }}}^{\mathrm{exact}}}\right|$, along the radial direction with different θ = 7575, 9075, and 10575 and the same longitude ϕ = 0°. From this figure we see that the relative errors are very small and the larger errors locate near inner and outer boundaries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Error and Order of Accuracy for Each Variable

Variable	Grid	L1 Error	Order	L2 Error	Order
	(Mesh resolution per component)
Φ	46 × 22 × 22	4.0210E-003		4.7288E-003
	46 × 42 × 42	1.1738E−003	2.8562	1.4289E−003	2.7762
	54 × 62 × 62	5.5548E−004	2.3896	6.7454E−004	2.3975
	64 × 82 × 82	3.1515E−004	2.3322	3.8256E−004	2.3337
Φ1	46 × 22 × 22	3.9827E−003		5.8529E−003
	46 × 42 × 42	1.3301E−003	2.5441	2.1882E−003	2.2823
	54 × 62 × 62	6.3175E−004	2.3780	1.0270E−003	2.4160
	64 × 82 × 82	3.6535E−004	2.2534	5.8897E−004	2.2879
Φ2	46 × 22 × 22	2.3133E−003		3.4966E−003
	46 × 42 × 42	6.2289E−004	3.0436	1.0791E−003	2.7272
	54 × 62 × 62	2.8299E−004	2.5199	4.9827E−004	2.4681
	64 × 82 × 82	1.6168E−004	2.3035	2.8533E−004	2.2940
Φ3	46 × 22 × 22	2.3150E−003		3.5870E−003
	46 × 42 × 42	8.1875E−004	2.4111	1.2712E−003	2.4064
	54 × 62 × 62	3.8973E−004	2.3710	6.0804E−004	2.3555
	64 × 82 × 82	2.2591E−004	2.2439	3.5307E−004	2.2367

Download table as:  ASCIITypeset image

3.2. Neumann Boundary Value Problem to 3D Poisson Equation with Manufactured Solutions

The second case is a Neumann boundary problem to 3D Poisson equation

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {y}^{2}}+\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}=0,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{\partial {\rm{\Phi }}}{\partial r}\right|}_{r={R}_{\mathrm{in}}} & = & -\sin \theta \cos \phi \times \left(\displaystyle \frac{2}{{r}^{3}}+\displaystyle \frac{1}{{R}_{\mathrm{out}}^{3}}\right),\\ {\rm{\Phi }}{| }_{r={R}_{\mathrm{out}}} & = & \sin \theta \cos \phi \times \left(\displaystyle \frac{1}{{r}^{2}}-\displaystyle \frac{r}{{R}_{\mathrm{out}}^{3}}\right),\end{array}\end{eqnarray} \tag{ 11 }$$

with the exact solution given by

$$\begin{eqnarray*}&&{\rm{\Phi }}(x,y,z)=\sin \theta \cos \phi \times \left(\displaystyle \frac{1}{{r}^{2}}-\displaystyle \frac{r}{{R}_{\mathrm{out}}^{3}}\right).\end{eqnarray*}$$

The computational domain is the same as that in the test case of Section 3.1. In the simulation, the solution is also initialized by zero for all variables and the boundary conditions are provided by Equation (11). The errors of each variable in terms of L₁ and L₂ norms and the rates of convergence are listed in Table 2. The results in Table 2 show that all the variables achieve the second-order accuracy as well. Meanwhile we present the comparison of the exact solution and the approximated one on the 180° − 0° plane of 54 × 62 × 62 mesh resolution per component in Figure 5. Besides, Figure 6 shows the relative error between numerical and exact solutions at cell centroid, $\left|\tfrac{{{\rm{\Phi }}}^{\mathrm{num}}-{{\rm{\Phi }}}^{\mathrm{exact}}}{{{\rm{\Phi }}}^{\mathrm{exact}}}\right|$, along the radial direction with different θ = 7575, 9075, 10575 and the same ϕ = 0°. From this figure we can see that the relative error increases gradually near the inner boundary, then begins to decrease, but starts to increase again near the outer boundary. Generally speaking, the value of the relative error is less than 4.5 × 10⁻⁴. Therefore, this case demonstrates the feasibility of the developed code again.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Error and Order of Accuracy for Each Variable

Variable	Grid	L1 Error	Order	L2 Error	Order
	(Mesh Resolution per Component)
Φ	46 × 22 × 22	1.0083E−004		2.0599E−004
	46 × 42 × 42	5.4027E−005	1.4474	1.0750E−004	1.5086
	54 × 62 × 62	2.5263E−005	2.4279	5.7002E−005	2.0263
	64 × 82 × 82	1.4663E−005	2.2385	3.4312E−005	2.0887
Φ1	46 × 22 × 22	3.1304E−004		5.6294E−004
	46 × 42 × 42	1.3080E−004	2.0243	2.1095E−004	2.2769
	54 × 62 × 62	6.7153E−005	2.1294	1.1856E−004	1.8404
	64 × 82 × 82	4.0943E−005	2.0360	7.6004E−005	1.8296
Φ2	46 × 22 × 22	2.0964E−004		3.7166E−004
	46 × 42 × 42	9.5207E−005	1.8311	1.5366E−004	2.0489
	54 × 62 × 62	4.8060E−005	2.1834	8.4215E−005	1.9208
	64 × 82 × 82	2.8887E−005	2.0947	5.2725E−005	1.9269
Φ3	46 × 22 × 22	2.1345E−004		3.6034E−004
	46 × 42 × 42	1.0098E−004	1.7363	1.5808E−004	1.9113
	54 × 62 × 62	5.1800E−005	2.1321	8.9169E−005	1.8288
	64 × 82 × 82	3.1350E−005	2.0664	5.6388E−005	1.8857

Download table as:  ASCIITypeset image

In this section, we apply the developed HCCFVS code to extrapolate the solar coronal magnetic fields, in which the computational domain stretches from 1 to 2.5 R_s as done in previous studies for the FDIPS and PFSS models (Hoeksema et al. 1983; Tóth et al. 2011). As the input to our code, we use the integral Global Oscillation Network Group (GONG) synoptic magnetograms to obtain the observed radial component of the magnetic field. These magnetograms are provided in the FITS format with the 180 × 360 grid, and the grid spacing is uniform in sine of the latitude and in longitude. In FDIPS, the monopole in magnetograms used for the inner boundary is cleaned by removing the average field, $\tfrac{{\sum }_{j,k}{B}_{r,j,k}{({\rm{\Delta }}\cos \theta )}_{j}{({\rm{\Delta }}\phi )}_{k}}{4\pi }$, from the observational field B_r. We adopt the same method as that in FDIPS to remove the monopole in our GONG synoptic magnetograms used for the inner boundary of HCCFVS. The initial conditions are set to zero. For maintaining the divergence- and current-free properties of the potential magnetic field solutions, we incorporate these two properties into the linear reconstruction of gradient variables by the least-squares technique when the magnetic field configuration in the Newton iteration procedure has been basically formed. Actually, this strategy is very similar to the globally solenoidality-preserving (GSP) method proposed by Feng et al. (2019). However, what is different from the GSP method is that three more equations derived from the current-free property are also added into the linear reconstruction of the gradient variables. Its implementation is briefly described in the following.

In the cell-centered second-order finite volume framework, the magnetic divergence in cell $i,{\overline{({\rm{\nabla }}\cdot {\boldsymbol{B}})}}_{i}$, is commonly evaluated by

$$\begin{eqnarray*}\begin{array}{rcl}{\overline{({\rm{\nabla }}\cdot {\boldsymbol{B}})}}_{i} & = & \displaystyle \frac{1}{2{V}_{i}}\displaystyle \sum _{j}\left[{\overline{{\boldsymbol{B}}}}_{i}+({\rm{\nabla }}\otimes {\boldsymbol{B}}){| }_{i}\cdot ({{\boldsymbol{x}}}_{{ij}}-{{\boldsymbol{x}}}_{i})\right.\\ & & \left.+{\overline{{\boldsymbol{B}}}}_{j}+({\rm{\nabla }}\otimes {\boldsymbol{B}}){| }_{j}\cdot ({{\boldsymbol{x}}}_{{ij}}-{{\boldsymbol{x}}}_{j})\right]\cdot {{\boldsymbol{n}}}_{{ij}}{\rm{\Delta }}{A}_{{ij}},\end{array}\end{eqnarray*}$$

with j denoting the face neighbor, and $({\rm{\nabla }}\otimes {\boldsymbol{B}}){| }_{i}$ and $({\rm{\nabla }}\otimes {\boldsymbol{B}}){| }_{j}$ representing the gradient tensors calculated at the centroids of cell i and j, respectively. For ∇ × B, in this study, we assume that

$$\begin{eqnarray*}&&{\overline{({\rm{\nabla }}\times {\boldsymbol{B}})}}_{i}=\displaystyle \frac{1}{{V}_{i}}{\int }_{{{ \mathcal V }}_{i}}{\rm{\nabla }}\times {\boldsymbol{B}}\ {dV}=({\rm{\nabla }}\times {\boldsymbol{B}}){| }_{i}.\end{eqnarray*}$$

In fact, maintaining the divergence- and current-free properties means that

$$\begin{eqnarray*}&&{\overline{({\rm{\nabla }}\cdot {\boldsymbol{B}})}}_{i}=0,\quad {\overline{({\rm{\nabla }}\times {\boldsymbol{B}})}}_{i}=0\end{eqnarray*}$$

holds, that is,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{j}\left[\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{i}\cdot ({{\boldsymbol{x}}}_{{ij}}-{{\boldsymbol{x}}}_{i})\right]\cdot {{\boldsymbol{n}}}_{{ij}}{\rm{\Delta }}{A}_{{ij}}\\ =-\displaystyle \sum _{j}\left[{\overline{{\rm{\nabla }}{\rm{\Phi }}}}_{i}+{\overline{{\rm{\nabla }}{\rm{\Phi }}}}_{j}+\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{j}\cdot ({{\boldsymbol{x}}}_{{ij}}-{{\boldsymbol{x}}}_{j})\right]\cdot {{\boldsymbol{n}}}_{{ij}}{\rm{\Delta }}{A}_{{ij}},\\ \qquad {\left.\left(\displaystyle \frac{\partial {{\rm{\Phi }}}_{3}}{\partial y}-\displaystyle \frac{\partial {{\rm{\Phi }}}_{2}}{\partial z}\right)\right|}_{i}=0,\\ \qquad {\left.\left(\displaystyle \frac{\partial {{\rm{\Phi }}}_{1}}{\partial z}-\displaystyle \frac{\partial {{\rm{\Phi }}}_{3}}{\partial x}\right)\right|}_{i}=0,\\ \qquad {\left.\left(\displaystyle \frac{\partial {{\rm{\Phi }}}_{2}}{\partial x}-\displaystyle \frac{\partial {{\rm{\Phi }}}_{1}}{\partial y}\right)\right|}_{i}=0.\end{array}\end{eqnarray} \tag{ 12 }$$

Note that B in Equation (12) has been replaced by ∇Φ = (Φ₁, Φ₂, Φ₃)^T. Then we add Equation (12) into the usual least-squares problem about Equation (8) related Φ₁, Φ₂ and Φ₃, and get an original overdetermined system searching for $\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{i}$. In order to make all the formulas of Equation (12) be satisfied as accurately as possible, we apply the elimination method to the original overdetermined system. Here, utilizing ${\left.\tfrac{\partial {{\rm{\Phi }}}_{1}}{\partial x}\right|}_{i},{\left.\tfrac{\partial {{\rm{\Phi }}}_{1}}{\partial y}\right|}_{i},{\left.\tfrac{\partial {{\rm{\Phi }}}_{1}}{\partial z}\right|}_{i}$ and ${\left.\tfrac{\partial {{\rm{\Phi }}}_{2}}{\partial z}\right|}_{i}$ for such elimination, we get a reduced overdetermined system. Because $\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{j}$ in the first formula of Equation (12) are unknown when solving the reduced overdetermined system, we need to iteratively solve it. The iterative algorithm for computing gradients is summarized as below: (1) set iteration counter, k = 0, and centroid gradients, $\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{i}^{(0)}$, at all cells. $\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{i}^{(0)}$ is obtained by solving the original overdetermined system; (2) first, solve the reduced overdetermined system in the least-squares sense with the singular value decomposition, and achieve ${\left.\tfrac{\partial {{\rm{\Phi }}}_{2}}{\partial x}\right|}_{i}^{(k+1)},{\left.\tfrac{\partial {{\rm{\Phi }}}_{2}}{\partial y}\right|}_{i}^{(k+1)},{\left.\tfrac{\partial {{\rm{\Phi }}}_{3}}{\partial x}\right|}_{i}^{(k+1)},{\left.\tfrac{\partial {{\rm{\Phi }}}_{3}}{\partial y}\right|}_{i}^{(k+1)}$ and ${\left.\tfrac{\partial {{\rm{\Phi }}}_{3}}{\partial z}\right|}_{i}^{(k+1)}$. Then, we obtain ${\left.\tfrac{\partial {{\rm{\Phi }}}_{1}}{\partial x}\right|}_{i}^{(k+1)},{\left.\tfrac{\partial {{\rm{\Phi }}}_{1}}{\partial y}\right|}_{i}^{(k+1)},{\left.\tfrac{\partial {{\rm{\Phi }}}_{1}}{\partial z}\right|}_{i}^{(k+1)}$ and ${\left.\tfrac{\partial {{\rm{\Phi }}}_{2}}{\partial z}\right|}_{i}^{(k+1)}$ by back substitution via Equation (12). This operation is carried out on all the computational cells; (3) check if convergence is achieved. If yes, stop the iterations, else set k = k + 1 and go to step 2. In this study, the iterations are assumed to be converged when ${\max }_{i}{\left|\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right)\right|}_{i}^{(k+1)}\,-\,\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{i}^{(k)}{| }^{2}\lt \varepsilon$, where ε is a small parameter (ε = 10⁻⁸) and $| \cdot {| }^{2}$ is defined as $| {\boldsymbol{G}}{| }^{2}\,\equiv {\sum }_{a}{\sum }_{b}{\left({G}_{{ab}}\right)}^{2},{\boldsymbol{G}}$ = $\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{i}^{(k+1)}-\left({\rm{\nabla }}\otimes ({\rm{\nabla }}{\rm{\Phi }})\right){| }_{i}^{(k)}$. In the actual computation, besides this convergence criterion, we also set up a maximal number of iterations, 100, to avoid infinite iterations.

Here we simulate the large-scale solar coronal magnetic field structures during CRs 2060, 2112, 2140, and 2180. These four CRs belong to the minimum, rising, maximum, and declining phases of solar activity, respectively. To validate HCCFVS, we compare the potential magnetic fields from it with those obtained by the SHE approach and the FDIPS method. In the SHE approach, the magnetograms are first remeshed into a uniform 181 × 361 latitude–longitude grid. Hence, we can obtain accurate solutions up to N_max = 120° spherical harmonics (Suda & Takami 2002). For the remeshed SHE with N_max = 120, it probably provides a better and much more reasonable solution, which is justified to some extent by Tóth et al. (2011) and Nikolić & Trichtchenko (2012). It is interesting to note that, in HCCFVS and FDIPS, the boundary data given by the magnetogram directly affect the solution only locally, while the spherical harmonics are global functions and their amplitudes depend on all of the magnetogram data. In practice, all the computations are accomplished on the Th-1A supercomputer from the National Supercomputing Center in TianJin, China, in which each computing node is configured with two Intel Xeon X5670 CPUs (2.93 GHz, six-core). For SHE, we only have a serial version. The wall o'clock time of running the serial code of SHE with 120 orders of harmonics is about 9 minutes. For HCCFVS and FDIPS, we use 120 MPI processes to obtain potential field solutions, and the wall o'clock times are 1.5 hr and 22 s, respectively. Obviously, the computational efficiency of HCCFVS is inferior to that of FDIPS. This is mainly due to the fact that there are four unknowns, ${\rm{\Phi }},{{\rm{\Phi }}}_{1},{{\rm{\Phi }}}_{2}$, and Φ₃, to be solved in HCCFVS, and the associated variable reconstructions and data transfer are very time consuming.

Figures 7 and 8 compare the coronal magnetic field topologies obtained with SHE (left), HCCFVS (middle), and FDIPS (right) on two different meridional planes of ϕ = 180° − 0° and ϕ = 270° − 90° for CRs 2060, 2112, 2140, and 2180, respectively. The backgrounds are colored with the base-10 logarithm of the magnetic field strength. The unit of magnetic field is Gauss, G for short. The black lines represent the magnetic field lines, on which the arrowheads denote the directions. Figures 7 and 8 reveal that the large-scale magnetic topologies from SHE, HCCFVS, and FDIPS are almost identical on the two selected meridional planes. By visual inspection, the results in these two figures from three methods achieve almost the same locations and coverages of the corresponding polar holes, the middle- and low-latitude open regions, and the unipolar and bipolar coronal streams from the three methods for the four specified CRs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 presents the relative differences of the magnetic fields along the radial direction in the open- (upper panels) and closed-field (lower panels) regions for CRs 2060 (first column, from the left), 2112 (second column), 2140 (third column), and 2180 (fourth column). The solid, dash–dot–dot, and dashed lines denote $\tfrac{| {{\boldsymbol{B}}}^{\mathrm{HCCFVS}}-{{\boldsymbol{B}}}^{\mathrm{SHE}}| }{| {{\boldsymbol{B}}}^{\mathrm{SHE}}| }$, $\tfrac{| {{\boldsymbol{B}}}^{\mathrm{HCCFVS}}-{{\boldsymbol{B}}}^{\mathrm{FDIPS}}| }{| {{\boldsymbol{B}}}^{\mathrm{FDIPS}}| }$, and $\tfrac{| {{\boldsymbol{B}}}^{\mathrm{FDIPS}}-{{\boldsymbol{B}}}^{\mathrm{SHE}}| }{| {{\boldsymbol{B}}}^{\mathrm{SHE}}| }$, respectively. Here, the superscripts "HCCFVS," "SHE", and "FDIPS" represent the magnetic fields obtained by HCCFVS, SHE, and FDIPS, respectively. From this figure, we can find that beyond 1.3 R_s, the magnetic fields from the three methods are very close in both the open- and closed-field regions for each of the four CRs. The larger differences among them mainly exist near 1 R_s, i.e., the bottom boundary of the computational domain, where the magnetic fields are dependent on high-order harmonic components of the magnetic fields.

Zoom In Zoom Out Reset image size

In order to present a comprehensive picture of the magnetic fields near the bottom boundary from the three methods, in Figure 10 we display the synoptic maps of the radial magnetic field and the superposed quiver representations of (B_θ, B_ϕ) in the "plane" of (θ, ϕ) at 1.005 R_s. From this figure, we can see that the radial magnetic field distributions obtained by HCCFVS agree well with those attained by SHE and FDIPS for all four CRs. The tangential vectors of the magnetic fields from the three methods are roughly consistent with each other except in the regions of extremely strong magnetic fields. The relative lengths of the arrowheads in Figure 10 represent the magnitudes of tangential vectors of the magnetic fields. By calculating the relative differences of the magnitudes of tangential vectors between HCCFVS and SHE, FDIPS and SHE, and HCCFVS and FDIPS, we find that such relative differences less than 40% account for about 80 percent for the four CRs. We further calculate the intersection angles in the selected grid points of the 3D magnetic fields between HCCFVS and SHE, FDIPS and SHE, and HCCFVS and FDIPS, respectively, and list the percentages of the angles less than 20°, 30°, and 35° in Table 3. This table demonstrates that there are about 60 percent of angles less than 20°, over 76 percent less than 30°, and over 82 percent less than 35°. Therefore, there are insignificant deviations among the results from the three methods. In addition, the results from HCCFVS are basically a compromise between those from SHE and FDIPS.

Zoom In Zoom Out Reset image size

In order to further validate the capability of HCCFVS, we, in this section, employ the potential fields from HCCFVS and SHE as the initial inputs to present a preliminary comparison of the simulated steady-state solar coronal solutions for CR 2098 by using the MHD model with a rotated-hybrid scheme developed in Feng et al. (2019). The divergence-free constraint on the magnetic field is kept by solving the Poisson equation in projection method with HCCFVS for both simulations.

Figure 11 shows the modeled magnetic field lines on two different meridional planes at ϕ = 180° − 0° (panels (a), (b), (e), (f)) and ϕ = 270° − 90° (panels (c), (d), (g), (h)) with the background colored with logarithms of plasma number densities to the base 10 (panels (a), (c), (e), (g)) and radial velocity (panels (b), (d), (f), (h)) by using the potential fields from the HCCFVS (upper panels) and SHE (lower panels) as the initial conditions. The arrowheads on the black lines represent the directions of the magnetic fields. The figure reveals that there is no significant difference between the two modeled results. From the figure, we can find that the high- and low-speed plasma flows distribute roughly axisymmetrically in respect to the line with latitudes 28° and −28° in the plane of ϕ = 180° − 0° and the line with latitudes 12° and −12° in the plane of ϕ = 270° − 90° for both simulations. The angular widths of middle- and low-speed plasma flows (v_r < 550 km s⁻¹ at 20 R_s) range between 60° and 70°. Therefore, high-speed plasma flow can be detected in low-latitudinal regions, which is consistent with the fact that the magnetic neutral line (MNL) at the solar rising phases tends to tilt more and the width of low-speed solar flow grows wider when compared with those at the solar minimum phases (e.g., Hoeksema et al. 1983 and Hu et al. 2008). Meticulous inspections find that high-speed flows occupy a slightly larger area in the coronal solutions with the SHE-derived magnetic fields as input.

Zoom In Zoom Out Reset image size

In Figure 12, we present synoptic maps of the observed (a), (b) and modeled results with the initial magnetic fields derived from HCCFVS (c)–(f) and SHE (g)–(j). Panels (a) and (b) show the white-light polarized brightness (pB) at the east and west limbs from the inner coronagraph (COR1) of the instrument suit of the Sun Earth Connection Coronal and Heliospheric aboard the Solar Terrestrial Relations observatory Behind (STEREO B). The red and black dashed lines denote the MNLs of the solutions with HCCFVS- and SHE-derived magnetic fields as input, respectively. Both MNLs are almost coincident and basically follow the trend of the pB bright structures. It should be noted that they are not completely across the brightest structures, which has been pointed out by previous studies (e.g., Sime & McCabe 1990 and Morgan & Habbal 2010). In addition, the bright areas away from the MNLs around the longitudes of 120° and 270° result from pseudo-streamers by inspecting the two-dimensional magnetic field topologies near the Sun (not shown here).

Zoom In Zoom Out Reset image size

Panels (c) and (d) display the modeled proton number density and radial speed at 2.6 R_s and (e) and (f) show the corresponding quantities at 20 R_s. Panels (g)–(j) are the same as (c)–(f), but for the solution with the SHE-derived magnetic field as input. From panels (c)–(j), we can easily find that both simulations achieve the similar distributions of the number density and radial velocity on the surface of 2.6 R_s and 20 R_s. The MNL can be roughly described as a cosine curve with small-amplitude fluctuations. The high-density low-speed coronal plasma flows pervade the zone around the MNLs and their latitudinal widths vary between about 50° and 90°. However, the high-density low-speed plasma flows from the simulation with HCCFVS-derived magnetic field as input span slightly larger domain than those from the simulation initiated by SHE-derived magnetic field. This is probably due to the fact that the solar wind acceleration function used here is dependent on the initial magnetic fields.

To demonstrate the performance of HCCFVS in controlling ∇ · B, Figure 13 presents the temporal evolution of the average relative divergence errors of the magnetic field for the two simulations. Following Zhang & Feng (2016), the average relative divergence error is defined as $\mathrm{Error}{({\boldsymbol{B}})}^{\mathrm{ave}}\,={\sum }_{i=1}^{M}\tfrac{\left|{\int }_{{{ \mathcal V }}_{i}}{\rm{\nabla }}\cdot {\boldsymbol{B}}\ {dV}\right|}{{\int }_{\partial {{ \mathcal V }}_{i}}| {\boldsymbol{B}}| \ {dA}}$/M, where M is the total number of cells in the computational domain. From this figure we can find that the average relative divergence errors increase rapidly in the first 10 hr, then begin to decrease, and tend to remain stable until the end of the simulation (Li et al. 2018). The final average relative divergence errors, Error (B)^ave, are around 10^−4.07 and 10^−4.22 for the simulations with HCCFVS- and SHE-derived magnetic fields as input, respectively. When employing HCCFVS to solve the Poisson equation involved in the projection method at every time step, we make the Newton iteration step of Equation (9) equal to 1, which is sufficient to control ∇ · B in our solar coronal simulation. As far as the relative divergence error is concerned, the performance of solving the Poisson equation in projection method with HCCFVS can be comparable to some commonly used divergence-cleaning approaches (Zhang & Feng 2016).

Zoom In Zoom Out Reset image size

We further measure the average relative differences between the corresponding parameters from the two simulations, where the average relative difference is defined by the formula ${\sum }_{i}\tfrac{| {W}_{i}^{(\mathrm{HCCFVS})}-{W}_{i}^{(\mathrm{SHE})}| /| {W}_{i}^{(\mathrm{SHE})}| }{M}$, with M being the total number of computational cells, and ${W}_{i}^{(\mathrm{HCCFVS})}$ and ${W}_{i}^{(\mathrm{SHE})}$ referring to the simulations initiated by the magnetic fields from HCCFVS and SHE, respectively. After calculation, the average relative difference is about 4% for number density, 3% for pressure, 10% for velocity, and 12% for the magnetic field. The average relative differences of the physical variables are small except for the magnetic field. The slightly larger relative difference for the magnetic field is probably attributed to the relative large difference for the initial magnetic fields from HCCFVS and SHE (13%). In fact, the magnetic fields reverse the direction when crossing the MNL, and thus the relative difference for the magnetic field tends to be larger than the ones for other physical quantities if the MNLs from both methods are not coincident absolutely.

In this paper, we have developed a new HCCFVS code with a six-component grid system to obtain the potential magnetic field solutions based on the solar magnetograms. By introducing solution gradients as additional unknowns and adding a pseudo-time derivative, HCCFVS reformulates the elliptic nature of the Poisson equation into a set of augmented equations that makes the entire system hyperbolic. A $P2+P1$ hybrid polynomial reconstruction is employed to seek the representation of the solution and gradient variables within a cell, which can reduce the computational demands to some extent by avoiding gradient computation for the solution variable. The residual equation is solved by Newton's method and the linear system of equations at each step of Newton's method is worked out by a simple block Jacobi iteration. For the established code, we first test two cases with exact solutions to preliminarily demonstrate its effectiveness and accuracy. The accuracy has been shown to achieve the second-order convergence. Then, we apply the code to the solar potential magnetic field problem. In this problem, for keeping the magnetic field divergence-free and current-free in the finite volume sense, we take into account the divergence- and current-free properties when using the least-squares method to construct the reconstruction polynomials of gradient variables. Comparison between the potential field solutions obtained by HCCFVS, SHE, and FDIPS shows that HCCFVS is capable of achieving the solar coronal magnetic fields. In order to further validate the capability of HCCFVS, we employ the potential fields from HCCFVS and SHE as the initial inputs and present a preliminary comparison of the simulated steady solar coronal solutions for CR 2098. Meanwhile, the divergence-free constraint on the magnetic field is maintained by solving the Poisson equation in projection method with HCCFVS. Both simulations capture large-scale solar coronal structures, and there is no significant difference between the two modeled results. We measure the average relative differences between the corresponding parameters from the two simulations, and the average relative differences are small except for the magnetic field. The slightly larger relative difference for the magnetic field is probably due to the relative large difference for the initial magnetic fields from HCCFVS and SHE. Finally, the average relative divergence errors are around 10⁻⁴ for the simulations with HCCFVS- and SHE-derived magnetic fields as input.

Since HCCFVS is established in the cell-centered finite volume framework, an attractive feature of HCCFVS is that it is applicable to both structured and unstructured grids. As a solver of a hyperbolic system in pseudo time, HCCFVS can employ many numerical techniques suitable for hyperbolic system. In addition, it can be extended to the high-order numerical scheme by increasing the spatial accuracy with high-order reconstructions and high-order flux integration on cell interface. However, the computational cost of HCCFVS is much higher than that of the FDIPS, although the solutions from the two methods are similar. Therefore, in order for HCCFVS to become an efficiently alternative method for the same problem, it is necessary to enhance the applicability and lower the computational cost of HCCFVS in future development. To this end, more efficient approaches, such as the GMRES (generalized minimum residual) algorithm with a LU-SGS (lower-upper symmetric Gauss–Seidel) preconditioner (Luo et al. 1998; Wang et al. 2019a), can be utilized to solve the Newtonian linearized system. Besides, the pseudo-time FAS (full approximation storage) multigrid iterative solver proposed by Kifonidis & Müller (2012) can be a promising way to efficiently solve the nonlinear residual system resulted from HCCFVS. Due to the great computational power of the GPU-acceleration technique (Feng et al. 2013; Wang et al. 2019b), power-efficient computing for compute-intensive GPU applications can be also explored to speed up HCCFVS.

The work is jointly supported by the National Natural Science Foundation of China (grant Nos. 41531073, 41731067, 41874202, 41861164026, 41574171, 41774184, and 41874205) and the Specialized Research Fund for State Key Laboratories. This work utilizes data obtained by the Global Oscillation Network Group (GONG) program, managed by the National Solar Observatory, which is operated by AURA, Inc., under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Inter-American Observatory. The work was carried out at the National Supercomputer Center in Tianjin, China, and the calculations were performed on TianHe-1 (A). The authors thank the reviewer for constructive suggestions for the improvement of the manuscript.