Direct discontinuous Galerkin method for potential magnetic field solutions

Abstract

In this paper, we employ the direct discontinuous Galerkin (DDG) method for the first time to extrapolate the coronal potential magnetic field (PF) with the source surface (SS) and call the developed numerical model as the DDG-PFSS solver. In this solver, the Laplace’s equation is solved by means of the time-dependent method, i.e., introducing a pseudo-time term into the Laplace’s equation and changing the boundary value problem into the initial-boundary value problem. The steady-state solution of the initial-boundary value problem is the solution of the Laplace’s equation to be solved. This formulation facilitates the implementation of the DDG discretization. In order to validate the DDG-PFSS solver, we test a problem with the exact solution, which demonstrates the effectiveness and third-order accuracy of the solver. Then we apply it to the extrapolation for the coronal potential magnetic field. We use the integral GONG synoptic magnetogram of Carrington rotation (CR) 2060 as the boundary condition and achieve the global potential magnetic field solution by the DDG-PFSS solver. The numerical results such as the coronal holes and streamer belts derived from the DDG-PFSS solver are in good agreement with those obtained from the spherical harmonic expansion method. Also, based on the numerical magnetic field and Wang-Sheeley-Arge model, the obtained solar wind speed is found to basically capture the structures of the high- and low-speed streams observed at 1 AU. These results suggest that the DDG-PFSS solver can be seen as a contribution to the numerical methods for obtaining the global potential magnetic field solutions of the solar corona.

1 Introduction

Potential field source surface (PFSS) model, proposed by Schatten, Wilcox, and Ness (1969) and Altschuler and Newkirk (1969), is the most popular model for reconstructing the potential solutions of the solar coronal magnetic field, which has the ability of forecasting the magnetic polarity and solar wind speed near the Earth in operational applications (Arge and Pizzo, 2000; Hakamada et al., 2002; Arge et al., 2003; MacNeice, 2009; Norquist and Meeks, 2010). The PFSS model starts with the synoptic magnetograms of the solar corona and extrapolates the coronal magnetic field to the “source surface”, i.e., a sphere with the particular height R_ss. Usually, R_ss is taken as 2.5R_s (R_s, the radius of the Sun) away from the photosphere, and at this source surface, the magnetic field is assumed to be purely radial. In a potential model, the common key assumption is that there are no currents in the region 1R_s ≤ r ≤ R_ss. Mathematically, the problem of obtaining the current free solution is described as follows: given the magnetogram data B_r (θ, ϕ) at r = 1R_s, find the scalar potential Φ so thatHere, θ ∈ [0, π] and ϕ ∈ [0, 2π] are the colatitude and longitude coordinates, respectively. After the solution of Eq. 1 is found, the potential magnetic field solution is obtained aswhich trivially satisfies both the divergence-free and the current-free properties

In order to obtain the potential magnetic field solutions, one usually uses the spherical harmonic expansion method (referred to as the SHE solver in this work) to solve the Laplace’s Eq. 1 with boundary conditions (2) and (3). However, this method is sensitive to the choice of the number of spherical harmonics, because the mismatches between the number of the harmonics and the resolution of the magnetogram can give rise to the ring-like patterns, especially near the strong magnetic field regions (Tóth, van der Holst, and Huang, 2011; Nikolić and Trichtchenko, 2012; Caplan et al., 2021). To mitigate or completely avoid these ringing effect, a finite difference iterative potential-field solver (FDIPS) has been developed by Tóth, van der Holst, and Huang (2011) and applied to solve PFSS. Encouraged by the first-order hyperbolic system formulation proposed by Nishikawa (2007) in solving the diffusion equation for fluid flow, Liu et al. (2019) developed a hyperbolic cell-centered finite volume solver for the PFSS model. In this paper, along with the purpose of enriching the numerical algorithms for obtaining the potential magnetic field solutions, we attempt to utilize the direct discontinuous Galerkin (DDG) method (Liu and Yan, 2009) to solve Eq. 1.

The DG method has been widely developed for hyperbolic problems since it was initially introduced in 1973 by Reed and Hill (1973), which is a class of finite element methods using a completely discontinuous piecewise polynomial space for the numerical solution and the test functions. This method can be easily extended to high-order approximation. The spatial accuracy of DG method is obtained by means of high-order polynomial approximation within an element rather than by wide stencils used in the finite difference and finite volume methods, which dramatically simplifies the use of high-order methods to some extent. Besides, the DG method has the capacity to handle the domain with complex geometry. Moreover, it is compact in the sense that each cell is treated independently and elements communicate only with direct adjacent elements regardless of the order of accuracy. For the other merits of DG method, we can refer to the literature (e.g., Cockburn, Karniadakis, and Shu, 2000; Hartmann, 2006; Luo et al., 2010; Dolejší and Feistauer, 2015). However, the application of the DG method to diffusion problems has been a challenging task because of the subtle difficulty in defining appropriate numerical fluxes for diffusion terms. To solve this problem, several approaches have been proposed, including the interior penalty method (Douglas and Dupont, 1976; Arnold, 1982), local DG method (Cockburn and Shu, 1998), the first Bassi-Rebay (BR1) scheme (Bassi and Rebay, 1997), the second Bassi-Rebay (BR2) scheme (Bassi et al., 2005), and the DDG method as used in this paper.

The main idea of the DDG method (Liu and Yan, 2009) for diffusion or elliptic equation is to use the direct weak formulation for solution of equations in each computational cell and let cells communicate through a numerical flux only. The main feature in the DDG scheme lies in the numerical flux choices for the solution gradient, which involve higher order derivatives calculated across cell interfaces (Liu, 2021). The DDG method has the advantage of easier formulation and implementation, as well as efficient computation of the solution. The DDG method holds the compactness in DG formulations, which is conducive to high efficiency massively parallel computation (He et al., 2022). Based on the finite difference scheme, there have been three publicly available solvers, including the FDIPS Fortran code publicly available at http://csem.engin.umich.edu/tools/FDIPS Web site, the PFSSPY Python package (Stansby, Yeates, and Badman, 2020) and the POT3D Fortran code (Caplan et al., 2021). Here we devote to establishing a numerical PFSS model based on the DDG method, which is termed as the DDG-PFSS solver. And, the magnetic fields derived by the DDG-PFSS solver are compared with that obtained by one of the three public methods. In this study, we choose the FDIPS code for comparision.

The remainder of this paper is organized as follows. We review the grid system of our code in Section 2. The discontinuous Galerkin discretization process is described in Section 3. In Section 4, we present the implicit time integration method used in our model. The numerical results are provided in Section 5. Finally, we state our conclusions in Section 6.

2 Grid system

The computational domain is set as a spherical shell geometrical region ranging from 1R_s to 2.5R_s. Following Feng et al. (2010), Feng et al. (2019), such spherical shell geometry is partitioned into the six-component grid system shown in Figure 1. This system consists of six identical component meshes with partial overlapping areas. Each component is a low-latitude spherical meshwhere δ is proportionally dependent on the grid spacing entailed for the minimum overlapping area. These 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.

FIGURE 1

In both θ and ϕ directions, grid points are uniformly distributed: θ_i = θ_min + (i − 1)Δθ, i = 1, …, N_θ − 1, ϕ_i = ϕ_min + (i − 1)Δϕ, i = 1, …, N_ϕ − 1, with Δθ = (θ_max − θ_min)/(N_θ − 2), Δϕ = (ϕ_max − ϕ_min)/(N_ϕ − 2), where N_θ and N_ϕ are the mesh numbers in latitudinal and longitudinal directions, respectively. In this paper, , , , , δ = 3Δθ. In the r direction, r (1) = 1R_s, r (i + 1) = r(1) + Δr(i), where i = 1, …, N_r. Δr(1) = 0.01R_s if r(1) < 1.1R_s; Δr(1) = min (A × log (r (i − 1)), Δθ × r (i − 1)) with A = 0.01/log (1.09) if r(1) < 3.5R_s.

There are two types of cells in our grid system after the computation area is partitioned. One is the hexahedron with six planes, shown in Figure 2A, and the other is the hexahedron with five planes and a spherical surface, shown in Figure 2B. The former denotes the interior computational cell and the latter the boundary cell near the inner boundary of the computational domain. The difference between them is that the boundary cell has one face lying on a spherical surface, and the associated area is labeled by the blue color.

FIGURE 2

3 Discontinuous Galerkin discretization

Introducing a pseudo-time derivative term to the Laplace’s Eq. 1 makes this equation into the following form:where τ is the pseudo time. Obviously, the steady-state solution of Eq. 4 is the solution of the Laplace’s Eq. 1.

In DG method, the numerical solution of Eq. 4 in each cell of computational grid is represented as a linear combination of the local polynomial basis functions ψ_l(x):Coefficients of this expansion, w_l(τ), are the main unknown values in the DG methods. Usually, the DG method is termed as the DG (P_k) method if the variable Φ is represented using a piecewise polynomials of degree k. The number of basis functions, K, depends on the degree k of the expanded polynomials and the spatial dimensions d,

For DG (P₂) case in three-dimensional space, K = 10, where P₂ indicates that a piecewise polynomial of degree of two is used to compute the fluxes. The basis functions we choose are Taylor basis functions and are listed as follows. For cell Ω_i,where x_i, y_i and z_i are the coordinates of the centroid of cell Ω_i. Δx_i, Δy_i, and Δz_i are the length scales in x −, y −, and z − directions, Δx_i = 0.5 (x_max − x_min), Δy_i = 0.5 (y_max − y_min), Δz_i = 0.5 (z_max − z_min), and x_max, x_min, y_max, y_min, z_max, and z_min are the maximum and minimum coordinates in the local cell in three directions, respectively. And, the coefficients in Eq. 5, i.e., w_l(τ), are the cell-average value and the scaled derivatives at the centroid of cell:

The weak formulation can be given through multiplying Eq. 4 by each basis function ψ_l(x), integrating it over cell Ω_i with boundary Γ_i = ∂Ω_i, and performing integration by parts:where n_ij denotes the unit outward normal vector to the boundary Γ_i. Substituting Eq. 5 into Eq. 9 results in

The integral in Eq. 10 is computed exactly or approximately by using suitable numerical quadratures (Xia, 2013; Liu et al., 2017). To complete the DG space discretization, we need to define the numerical flux , for which we adopt the DDG method in this work.

In order to describe DDG method in a more convenient way, we introduce the following jump and average operators:where the subscripts ⋅|_L and ⋅|_R denote the left and right states of a common face, respectively. With this definition, the x-direction, y-direction and z-direction components of flux for the DDG method can be given aswhere β₀ and β₁ are regarded as two constant penalty parameters, and Δ is the characteristic face length. According to Cheng et al. (2016), we choosewhere d_ij = (x_j − x_i, y_j − y_i, z_j − z_i) is the displacement vector from the centroid of Ω_i to the centroid of Ω_j. Hence, the numerical flux is calculated asAt this moment, the flux Jacobian matrices for Φ and ∇Φ contributing to the diagonal block are respectivelyand the flux Jacobian flux for Φ and ∇Φ contributing to the off-diagonal block are respectively

As for boundary face, the numerical flux is constructed by using only the first derivatives and by assuming the jump of the second derivatives to be zero, which is as follows:withTherefore, at the outer boundary face,where Φ_B is given at the boundary specified by the Dirichlet condition. In this case, the numerical flux given by Eq. 19 becomesand the flux Jacobian matrices for Φ and ∇Φ contributing to the diagonal block are respectivelyAt the inner boundary face,where is the normal gradient given at the boundary specified by the Neumann condition. In this case, the numerical flux given by Eq. 19 becomesand the flux Jacobian matrices for Φ and ∇Φ contributing to the diagonal block areAt an inner boundary face as a patch of the spherical surface, the normal direction of the boundary coincides with the radial direction. And, the outer normal direction of the inner boundary is negative to the radial direction, i.e., pointing to the interior of the Sun. As a result, . At boundary faces, Δ is chosen as the distance from cell centroid to the boundary faces.

The spatial discretization in Eq. 10 leads to a system of ordinary differential equations:where M_i denotes the mass matrix whose entries are , W is the solution vector made up of the coefficients in Eq. 5, and R_i is the residual vector over cell Ω_i, consisting of K components shown in the right hand side of Eq. 10. In order to speed up the convergence in steady problems, implicit temporal discretization is applied to the pseudo-time system Eq. 29 in this work.

4 Implicit time integration

Using Euler implicit time integration, the semi-discrete system of ordinary differential Eq. 29 is rewritten aswhere Δτ is the pseudo-time increment, and with CFL = 1000 in this work. Δh_min denotes the minimum value of the sizes of all of the computational cells. is the difference of unknown vector between levels n and n + 1. Each component of the residual vector R_i (Φⁿ⁺¹) can be linearized in pseudo-time aswhere denotes the gradient of Φ, andwith the superscripts i and j representing cell Ω_i and Ω_j, respectively.

Substituting Eq. 31 into Eq. 30, we obtain the following linear system over cell Ω_iwhere

For boundary cells, note that there are only five matrices on the left hand side of Eq. 36, corresponding to the five face-neighboring cells respectively. Considering all cells in the computational domain, we obtain the large linear system at each pseudo-time step. In terms of solving large linear systems, iterative methods are vital tools. In this study, we employ the block successive over-relaxation (SOR) scheme which can be written at cell Ω_i aswhere m denotes the iteration step, and ω the over-relaxation factor which is taken to be 1.3 in this work. This relaxation is performed to a specified tolerance, i.e., , or to a specified maximum number 20 of relaxations. The implicit solver is taken as converged when the L₂ norm of the Φ residual satisfies , where ϵ is a small positive number and considered in the range 10^–6 ∼ 10^–10.

5 Numerical results

In this section, a 3D test case with the exact solutions is first performed to verify the effectiveness and accuracy of the developed DDG-PFSS code. Then we apply it to the extrapolation for solar potential magnetic field. All the computations in this study 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).

5.1 Validation of DDG-PFSS

To validate the solutions generated with the DDG-PFSS solver, we use an analytic solution, which is constructed according to Caplan et al. (2021), given by

The potential solution is given by

To test the code, we set the initial guess Φ = 0 and solve the system using the analytic B_r of Eq. 40 as the inner boundary condition in the manner of Eq. 2 and using Eq. 3 as the outer boundary condition. The implicit solver is taken as converged when the L₂ norm of the Φ residual is reduced by six orders of magnitude, i.e., ϵ = 1.0 × 10^–6. For any solution variable s, the L₂ norm error is defined aswhere s(x) and s(x)^exact are respectively the numerical and exact solutions evaluated at point x. N_global denotes the total number of the interior computational cells. Starting with a resolution of 40 × 12 × 12 × 6, we test the code with successively increased resolution. All runs are performed using 96 MPI processes. The results are presented in Table 1 and Figure 3. In Table 1, we present the L₂ norm error of Φ on each level of meshes, as well as the order of accuracy for Φ. Between the coarser and finer meshes, following Lee, Ahn, and Luo (2018), the order of accuracy is measured as followswith the superscripts “coarse” and “fine” denote the coarse mesh and fine mesh, respectively. Figure 3 displays the three-dimensional topologies of the gradient B = ∇Φ of Φ on the level of mesh with the resolution of 40 × 42 × 42 × 6. The left panel in this figure represents the analytic solution and the right panel the numerical solution. In these two panels, the starting points of the black lines representing ∇Φ are the same in order to make a visual comparison between the numerical solution and the exact solution. We see that the numerical solution is almost identical with the analytic solution, and the code exhibits the third-order accuracy as the mesh is gradually refined.

FIGURE 3

5.2 Coronal potential magnetic field extrapolation

In this section, we extrapolate the coronal magnetic fields for CR 2060 (14 August To 10 September 2007) using the developed DDG-PFSS solver. The integral GONG synoptic magnetograms (in “mrmqs” format) are used as the inner boundary condition. These magnetograms are provided in the FITS format with the uniform 180 × 360, sin θ − ϕ mesh. With the help of MATLAB for reading the FITS files, we remesh the magnetograms to a uniform 181 × 361 latitude-longitude grid. In our DDG-PFSS solver, the magnetogram needs to be pre-processed, which is to remove the monopole. Therefore, the observed radial field B_r (R_s, θ_j, ϕ_k) at the solar surface is corrected into . The initial condition is set to zero, and the parameter ϵ of the convergence criterion is set to 1.0 × 10^–10. When extrapolating the coronal magnetic fields with the DDG-PFSS solver, we adopt the uniform radial grid, which is similar to the finite difference solver FDIPS. The mesh resolution we used in the DDG-PFSS solver is (Δr, Δθ, Δϕ) = (0.025R_s, 2.25°, 2.25°). In this solver, we use 168 MPI processes to obtain the potential field solutions, and the wall o’clock times are about 4.7 h. The divergence-free condition is basically satisfied to the level of the solver tolerance. The maximum value of the divergence errors over all cells in the simulation domain, , is about 4.2 × 10^–7 and the global divergence L₂ norm of the magnetic field, , is about 5.05 × 10^–10, where . The magnetic field B at the cell interfaces is calculated in the same way as the first term of the right residual in Eq. 10 when we calculate this surface integral. In order to validate the numerical magnetic fields generated by the DDG-PFSS solver, we compare these results with those derived by the FDIPS code. The magnetic fields obtained by the SHE method with a number n_max = 120 of harmonics on the mesh resolution of (Δr, Δθ, Δϕ) = (0.025R_s, 1°, 1°) serve as the reference solution in this study. The results from the FDIPS code are obtained by the version of (UseCosTheta, nR, nTheta, nPhi) = (T, 60, 180, 360). For the SHE solver and the FDIPS code, we adopt the serial versions, and their wall o’clock times are about 10.5 min and 7.5 min, respectively. Obviously, the computational efficiency of the DDG-PFSS solver is inferior to that of the FDIPS code. This is the main weakness of the DG method since it requires solutions of systems of equations with more unknowns.

According to Schrijver et al. (2006) and Liu et al. (2011), we calculate some metrics that are used for checking the performance of the extrapolation solvers. They are C_vec, C_cs, E_m and E_n, respectively, which are briefly defined below. C_vec is used to quantify the vector correlation,where and are the magnetic field derived by the SHE solver and the extrapolated field at the grid point i, respectively. If the field vectors are identical, then C_vec = 1; if , then C_vec = 0. C_cs is based on the Cauchy-Schwarz inequality and is mostly applied for measuring the differences of the vector fields,where M_global is the total number of vectors in the domain to be calculated. C_cs = 1 means that and are parallel; C_cs = −1 means that and are anti-parallel; C_cs = 0 means that . E_n is a normalized vector error,E_m is a mean of the normalized vector errors,The agreement is perfect when E_n and E_m are equal to zero, which is unlike the first two metrics.

In Table 2, we show the metrics calculated using the extrapolated magnetic fields of the FDIPS and DDG-PFSS solvers. From the metrics, we see that the potential magnetic fields derived by the DDG-PFSS solver are slightly better than that obtained by the FDIPS method. Figure 4 presents the coronal magnetic field topologies obtained with (A) the SHE solver, (B) the FDIPS code and (C) the DDG-PFSS solver on the meridional plane of ϕ = 180°–0°. The backgrounds are colored with the radial component of magnetic field (Unit: G). The black lines represent the magnetic field lines, on which the arrowheads denote the directions. In order to make direct comparisons between these two numerical methods, the starting points at 1R_s of the magnetic fields are the same. By visual inspection, the large-scale magnetic topologies obtained from the three solvers are almost identical, such as the open-field regions, bipolar streamers and pseudo-streamers. In addition, we also present the comparison of the radial magnetic field between the two numerical solvers in Figure 5. The black, green and magenta lines represent the radial fields obtained by the SHE, FDIPS and DDG-PFSS solvers, respectively. The backgrounds are colored with the radial field of the SHE solver. From the low-latitude regions to the high-latitude regions, we can see that the radial field derived by the DDG-PFSS solver is closer to that of the SHE solver when compared to that derived by FDIPS. These results suggest the credibility of our established DDG-PFSS solver. In the following section, we continue to obtain other large-scale structures using our DDG-PFSS solver.

By tracing the magnetic field lines from the source surface to the photosphere, we obtain the coronal holes (CHs) regions, i.e., the footpoints of open magnetic field lines, which are shown in the dark regions in Figure 6. In Figure 6A, we present the 195 Å EUVI observation obtained by Solar Terrestrial Relations Observatory (STEREO)-A/Extreme Ultraviolet Imager (EUVI) instrument. In Figures 6B,C, we display the distributions of CHs derived from the SHE and DDG-PFSS solvers, respectively. From Figures 6B,C, we find no apparent differences between them. In addition, we calculate the CHs areas, shown as fraction of the Sun’s total surface area, which is about 13.77% for the SHE solver and 14.28% for the DDG-PFSS solver, respectively. From Figure 6A, we see that large CHs are mainly located in polar regions (polar CHs) during the solar minimum. The observed polar CHs and the extension of the southern polar CHs at about the longitudinal range of 90°–180° to lower latitude are well caught by the SHE and DDG-PFSS solvers. In observation, we also see some small and isolated CHs in the low-latitude regions, which we call as low-latitude CHs. Comparing the SHE- and DDG-derived results with observation shows that the configuration of low-latitude CHs is basically reproduced. However, some observed low-latitude CHs marked by red rectangles in Figure 6A are not captured by the SHE and DDG-PFSS solvers, which is also evidenced by Nikolić (2019). The red rectangles appearing in Figures 6B,C are located at the same positions as that shown in Figure 6A, with the purpose of making the visual comparison more convenient.

In Figure 7A, we show the synoptic map of polarized brightness (pB) observation at 2.5R_s. This synoptic map is created using data from the east limb of the SOHO/LASCO C2 instrument. The bright structures in the pB synoptic map may include the bipolar streamer and the pseudo-streamer (Wang, Sheeley, and Rich, 2007; Li, Feng, and Wei, 2021). The red line superimposed on this pB synoptic map denotes the magnetic neutral line (B_r = 0) from the SHE solver, while the blue line represents the magnetic neutral line from the DDG-PFSS solver. From Figure 7A, we see that the blue line almost coincides with the red line and they are surrounded by the bright structures. Usually, the position of the magnetic neutral line is considered as the bipolar streamer belt. In comparison with the observation, the positions of SHE-derived and DDG-derived bipolar streamer belts agree reasonably with the observed pB bright structures. In addition, we use a streamer identification algorithm (Owens, Crooker, and Lockwood, 2014; Li, Feng, and Wei, 2021) to further identify the bipolar streamer and pseudo-streamer from the SHE- and DDG-derived results. That is, the regions where the parameter log (ds) > 1 are identified as bipolar streamer and pseudo-streamer belts. The parameter log (ds) is defined as with d_ss denoting the distance between two adjacent grid points on the source surface and d_ph the distance between the photospheric footpoints of their corresponding magnetic field lines. These results from the SHE solver and the DDG-PFSS solver are presented in Figures 7B,C, respectively. Apart from the bipolar streamer belt, i.e., the position of the magnetic neutral line, we find that the pseudo-streamer belt derived from the DDG-based results is not only similar to that obtained by the SHE-based result, but also agrees well with the observation.

The PFSS model plays an important role in forecasting solar wind. Here, we take advantage of the Wang-Sheeley-Arge (WSA) model (Arge et al., 2003; Riley, Linker, and Arge, 2015) and the Arge-Pizzo kinematic evolution model (Arge and Pizzo, 2000; Riley and Lionello, 2011) to predict the solar wind speed v_r at 1 AU. More specifically, by using the SHE- and DDG-derived magnetic fields we first configure the WSA model to obtain the solar wind speed at 30R_s, and then map the solar wind streams from 30R_s to 1 AU using the Arge-Pizzo kinematic evolution model. The general formulation of WSA model for solar wind speed is as followswhere a₁ ∼ a₈ are the empirical and adjustable parameters. f_s represents the magnetic flux tube expansion factor and θ_b the minimum angular separation between an open magnetic field footpoint and its nearest coronal hole boundary. In Figure 8, we display the synoptic maps of two parameters f_s and θ_b at 2.5R_s, in which f_s and θ_b in Figure 8(A1, A2) are obtained by the SHE-derived magnetic fields, and those in Figure 8(B1, B2) from the DDG-derived magnetic fields. According to Riley and Lionello (2011), the Arge-Pizzo kinematic evolution model is expressed aswhere V_i+1,j denotes the speed at r_i+1 for longitude j. The iterative number in radial direction, which is applied between 30R_s and 1 AU, is set here to be 120. The mapped results are exhibited in Figure 9, which demonstrates the comparison between the OMNI data and the numerical results. The red lines in Figures 9A,B represent the mapped results based on the SHE solver and the DDG-PFSS solver, respectively. The black line denotes the OMNI data. The parameters in Eq. 48 which are chosen to produce the best match to the OMNI data are listed in Table 3.

FIGURE 4

FIGURE 5

FIGURE 6

FIGURE 7

FIGURE 8

FIGURE 9

The comparison between the modeled and observational results suggests that the simulated results basically capture the structures of the high- and low-speed streams occurring in the observation, apart from the small deviation of the high-speed stream at about longitudinal range of 150°–190° shown in Figure 9A.

6 Conclusion

In this paper, we develop a DDG-PFSS solver for obtaining the potential magnetic field solutions. By introducing the pseudo-time into the Laplace’s Eq. 1, the Laplace’s equation is changed into a diffusion-like Eq. 4. The initial-boundary value problem for the modified equation is obtained to replace the boundary value problem for the Laplace’s equation. Then we apply the DG method to the space discretization of the modified equation. The steady-state solution of the initial-boundary value problem is the sought solution of the Laplace’s equation. As for the space discretization, we choose the DDG method to calculate the numerical flux across the interfaces because of its simplicity in implementation and efficiency in computation. In order to speed up the convergence, implicit backward Euler temporal discretization is applied to the pseudo-time system, and the resulting large linear system at each pseudo-time step is solved by the block SOR scheme.

To validate the DDG-PFSS solver, we test a problem with the exact solution, which demonstrates the effectiveness and third-order accuracy of the solver. Then we utilize this solver to calculate the global coronal magnetic field of CR 2060. The numerical magnetic fields generated by the DDG-PFSS solver are compared with those derived by the FDIPS code. The magnetic fields obtained by the SHE method serve as a reference solution. By comparison, we see that the magnetic fields derived by the DDG-PFSS solver are closer to that of the SHE solver. Finally, we also obtain the large-scale structures including CHs, bipolar streamer belts and pseudo-streamer belts using the DDG-PFSS solver. These structures have a good agreement with solar observations except for a few small and isolated low-latitude CHs. In addition, we obtain the solar wind speed at 1 AU using the WSA model and the Arge-Pizzo kinematic evolution model. The derived solar wind speed basically captures the structures of the high- and low-speed streams as occurred in the OMNI data. These results strongly suggest the credibility of our established DDG-PFSS solver. Compared to other numerical methods including the finite difference and finite volume methods, DDG discrete method possesses some attractive features, such as the compactness and the easy extension to high-order approximation. The numerical viscous flux defined by the DDG method is consistent and conservative (Liu and Yan, 2009). The most prominent feature of the DDG method is its simplicity in implementation and its efficiency in computation (Cheng et al., 2017). In a word, our DDG-PFSS model can be seen as an alternative contribution to computational space weather modelling in space physics.

The work is jointly supported by the National Natural Science Foundation of China (Grant Nos. 42030204, 41874202, 41731067, 41861164026, and 41531073), the strategic priority research program of CAS with Grant No. XDB09000000 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 Interamerican Observatory. The STEREO/SECCHI data are produced by a consortium of the NRL (United States), LMSAL (United States), NASA/GSFC (United States), RAL (United Kingdom), UBHAM (United Kingdom), MPS (Germany), CSL (Belgium), IOTA (France), and IAS (France). This work also utilizes data obtained by LASCO C2/SOHO. The OMNI data are obtained from the GSFC/SPDF CDAWeb interface at https://cdaweb.gsfc.nasa.gov/index.html. The work was carried out at the National Supercomputer Center in Tianjin, China, and the calculations were performed on TianHe-1 (A).

References

• 1

  AltschulerM. D.NewkirkG. (1969). Magnetic fields and the structure of the solar corona. Sol. Phys.9 (1), 131–149. 10.1007/BF00145734

  • CrossRef
  • Google Scholar

• 2

  ArgeC. N.OdstrcilD.PizzoV. J.MayerL. R. (2003). “Improved method for specifying solar wind speed near the Sun,”. Editors VelliM.BrunoR.MalaraF.BucciB. (American Institute of Physics Conference Series), Sol. Wind Ten. 679. 190–193. 10.1063/1.1618574

  • CrossRef
  • Google Scholar

• 3

  ArgeC. N.PizzoV. J. (2000). Improvement in the prediction of solar wind conditions using near-real time solar magnetic field updates. J. Geophys. Res.105 (A5), 10465–10479. 10.1029/1999JA000262

  • CrossRef
  • Google Scholar

• 4

  ArnoldD. N. (1982). An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal.19 (4), 742–760. 10.1137/0719052

  • CrossRef
  • Google Scholar

• 5

  BassiF.CrivelliniA.RebayS.SaviniM. (2005). Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k − ω turbulence model equations. Comput. Fluids34 (4), 507–540. Residual Distribution Schemes, Discontinuous Galerkin Schemes and Adaptation. 10.1016/j.compfluid.2003.08.004

  • CrossRef
  • Google Scholar

• 6

  BassiF.RebayS. (1997). A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations. J. Comput. Phys.131 (2), 267–279. 10.1006/jcph.1996.5572

  • CrossRef
  • Google Scholar

• 7

  CaplanR. M.DownsC.LinkerJ. A.MikicZ. (2021). Variations in finite-difference potential fields. Astrophys. J.915 (1), 44. 10.3847/1538-4357/abfd2f

  • CrossRef
  • Google Scholar

• 8

  ChengJ.LiuX.LiuT.LuoH. (2017). A parallel, high-order direct discontinuous Galerkin method for the Navier-Stokes equations on 3D hybrid grids. Commun. Comput. Phys.21 (5), 1231–1257. 10.4208/cicp.OA-2016-0090

  • CrossRef
  • Google Scholar

• 9

  ChengJ.YangX.LiuX.LiuT.LuoH. (2016). A direct discontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids. J. Comput. Phys.327, 484–502. 10.1016/j.jcp.2016.09.049

  • CrossRef
  • Google Scholar

• 10

  CockburnB.KarniadakisG. E.ShuC. W. (Editors) (2000). Discontinuous Galerkin methods (Berlin, Heidelberg: Springer). 978-3-642-64098-8. 10.1007/978-3-642-59721-3

  • CrossRef
  • Google Scholar

• 11

  CockburnB.ShuC. W. (1998). The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J. Numer. Anal.35 (6), 2440–2463. 10.1137/S0036142997316712

  • CrossRef
  • Google Scholar

• 12

  DolejšíV.FeistauerM. (2015). Discontinuous Galerkin method. Cham: Springer. 978-3-319-19267-3. 10.1007/978-3-319-19267-3

  • CrossRef
  • Google Scholar

• 13

  DouglasJ.DupontT. (1976). Interior penalty procedures for elliptic and parabolic Galerkin methods. Editors GlowinskiR.LionsJ. L., 58, 207. 10.1007/BFb0120591

  • CrossRef
  • Google Scholar

• 14

  FengX. S.YangL. P.XiangC. Q.WuS. T.ZhouY. F.ZhongD. K. (2010). Three-dimensional solar wind modeling from the Sun to Earth by a SIP-CESE MHD model with a six-component grid. Astrophys. J.723, 300–319. 10.1088/0004-637x/723/1/300

  • CrossRef
  • Google Scholar

• 15

  FengX.LiX.XiangC.LiH.WeiF. (2019). A new MHD model with a rotated-hybrid scheme and solenoidality-preserving approach. Astrophys. J.871 (2), 226. 10.3847/1538-4357/aafacf

  • CrossRef
  • Google Scholar

• 16

  HakamadaK.KojimaM.TokumaruM.OhmiT.YokobeA.FujikiK. (2002). Solar wind speed and expansion rate of the coronal magnetic field in solar maximum and minimum phases. Sol. Phys.207, 173–185. 10.1023/A:1015511805863

  • CrossRef
  • Google Scholar

• 17

  HartmannR. (2006). Adaptive discontinuous Galerkin methods with shock-capturing for the compressible Navier-Stokes equations. Int. J. Numer. Methods Fluids51 (9-10), 1131–1156. 10.1002/fld.1134

  • CrossRef
  • Google Scholar

• 18

  HeX.WangK.FengY.LvL.LiuT. (2022). An implementation of mpi and hybrid openmp/mpi parallelization strategies for an implicit 3d ddg solver. Comput. Fluids241, 105455. 10.1016/j.compfluid.2022.105455

  • CrossRef
  • Google Scholar

• 19

  LeeE.AhnH. T.LuoH. (2018). Cell-centered high-order hyperbolic finite volume method for diffusion equation on unstructured grids. J. Comput. Phys.355, 464–491. 10.1016/j.jcp.2017.10.051

  • CrossRef
  • Google Scholar

• 20

  LiH.FengX.WeiF. (2021). Comparison of synoptic maps and PFSS solutions for the declining phase of solar cycle 24. JGR. Space Phys.126 (3). 10.1029/2020JA028870

  • CrossRef
  • Google Scholar

• 21

  LiuH. (2021). Analysis of direct discontinuous Galerkin methods for multi-dimensional convection–diffusion equations. Numer. Math. (Heidelb).147, 839–867. 10.1007/s00211-021-01183-x

  • CrossRef
  • Google Scholar

• 22

  LiuH.YanJ. (2009). The direct discontinuous Galerkin (DDG) methods for diffusion problems. SIAM J. Numer. Anal.47 (1), 675–698. 10.1137/080720255

  • CrossRef
  • Google Scholar

• 23

  LiuS.ZhangH. Q.SuJ. T.SongM. T. (2011). Study on two methods for nonlinear force-free extrapolation based on semi-analytical field. Sol. Phys.269, 41–57. 10.1007/s11207-010-9691-4

  • CrossRef
  • Google Scholar

• 24

  LiuX.FengX.XiangC.ShenF. (2019). Hyperbolic cell-centered finite volume method for obtaining potential magnetic field solutions. Astrophys. J.887 (1), 33. 10.3847/1538-4357/ab4b53

  • CrossRef
  • Google Scholar

• 25

  LiuX.XuanL.XiaY.LuoH. (2017). A reconstructed discontinuous Galerkin method for the compressible Navier-Stokes equations on three-dimensional hybrid grids. Comput. Fluids152, 217–230. 10.1016/j.compfluid.2017.04.027

  • CrossRef
  • Google Scholar

• 26

  LuoH.LuoL.NourgalievR.MousseauV. A.DinhN. (2010). A reconstructed discontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids. J. Comput. Phys.229 (19), 6961–6978. 10.1016/j.jcp.2010.05.033

  • CrossRef
  • Google Scholar

• 27

  MacNeiceP.JianL. K.AntiochosS. K.ArgeC. N.Bussy-ViratC. D.DeRosaM. L.et al (2009). Validation of community models: 2. Development of a baseline using the wang-sheeley-arge model. Space weather.7 (12), 1644–1667. 10.1029/2018SW002040

  • CrossRef
  • Google Scholar

• 28

  NikolićL. (2019). On solutions of the PFSS model with GONG synoptic maps for 2006-2018. Space weather.17 (8), 1293–1311. 10.1029/2019SW002205

  • CrossRef
  • Google Scholar

• 29

  NikolićL.TrichtchenkoL. (2012). Development of coronal field and solar wind components for MHD interplanetary simulations. Int. J. Environ. Chem. Ecol. Geol. Geophys. Eng.6 (11), 698–702.

  • Google Scholar

• 30

  NishikawaH. (2007). A first-order system approach for diffusion equation. I: Second-order residual-distribution schemes. J. Comput. Phys.227, 315–352. 10.1016/j.jcp.2007.07.029

  • CrossRef
  • Google Scholar

• 31

  NorquistD. C.MeeksW. C. (2010). A comparative verification of forecasts from two operational solar wind models. Space weather.8 (12). 10.1029/2010SW000598

  • CrossRef
  • Google Scholar

• 32

  OwensM. J.CrookerN. U.LockwoodM. (2014). Solar cycle evolution of dipolar and pseudostreamer belts and their relation to the slow solar wind. J. Geophys. Res. Space Phys.119 (1), 36–46. 10.1002/2013JA019412

  • CrossRef
  • Google Scholar

• 33

  ReedW. H.HillT. R. (1973). Triangular mesh methods for the neutron transport equation.

  • Google Scholar

• 34

  RileyP.LinkerJ. A.ArgeC. N. (2015). On the role played by magnetic expansion factor in the prediction of solar wind speed. Space weather.13 (3), 154–169. 10.1002/2014SW001144

  • CrossRef
  • Google Scholar

• 35

  RileyP.LionelloR. (2011). Mapping solar wind streams from the sun to 1 au: A comparison of techniques. Sol. Phys.270, 575–592. 10.1007/s11207-011-9766-x

  • CrossRef
  • Google Scholar

• 36

  SchattenK. H.WilcoxJ. M.NessN. F. (1969). A model of interplanetary and coronal magnetic fields. Sol. Phys.6 (3), 442–455. 10.1007/BF00146478

  • CrossRef
  • Google Scholar

• 37

  SchrijverC. J.De RosaM. L.MetcalfT. R.LiuY.McTiernanJ.RégnierS.et al (2006). Nonlinear force-free modeling of coronal magnetic fields Part I: A quantitative comparison of methods. Sol. Phys.235 (1-2), 161–190. 10.1007/s11207-006-0068-7

  • CrossRef
  • Google Scholar

• 38

  StansbyD.YeatesA.BadmanS. T. (2020). pfsspy: A python package for potential field source surface modelling. J. Open Source Softw.5 (54), 2732. 10.21105/joss.02732

  • CrossRef
  • Google Scholar

• 39

  TóthG.van der HolstB.HuangZ. (2011). Obtaining potential field solutions with spherical harmonics and finite differences. Astrophys. J.732, 102. 10.1088/0004-637X/732/2/102

  • CrossRef
  • Google Scholar

• 40

  WangY. M.SheeleyJ. N. R.RichN. B. (2007). Coronal pseudostreamers. Astrophys. J.658 (2), 1340–1348. 10.1086/511416

  • CrossRef
  • Google Scholar

• 41

  XiaY. (2013). A parallel implicit reconstructed discontinuous Galerkin method for compressible flows on hybrid grids. PhD thesis. Raleigh, NC: North Carolina State University.

  • Google Scholar