The Poisson equation is central in numerous physics and engineering applications, such as computational fluid dynamics and acoustic wave propagation, where efficient and accurate solutions are essential. This study focuses on the numerical solution of the 2D Poisson equation with Dirichlet boundary conditions using a fourth-order compact Implicit Finite Difference scheme. Finite difference methods, particularly high-order schemes, are advantageous for solving the Poisson equation due to their efficiency and suitability for structured grids. To address the computational demands of large-scale problems, we incorporate domain decomposition and the Multicolor Successive Over Relaxation method, facilitating parallel computation. Through numerical experiments, we demonstrate that our approach significantly enhances both accuracy and computational efficiency when compared to traditional second-order methods.
keywords
Fourth-order finite difference, Domain decomposition, MPI, Poisson equation, Speedup, Elapsed time
The Poisson equation plays a crucial role in a wide range of physics and engineering applications, including computational fluid dynamics, acoustic wave propagation, and theoretical physics [1]. However, the implementation of simple, accurate, and efficient numerical solver remains fundamental to the successful resolution of these problems. In this work, we study the numerical solution of the 2D Poisson equation, which is formulated as follows:
|
(1) |
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega = (\mathfrak{a},\mathfrak{b})\times (\mathfrak{c},\mathfrak{d})}
and Dirichlet boundary conditions are applied on Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \partial \Omega }
. The right-hand side Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle f}
is a know function defined over Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega }
. The study of large-scale physical phenomena, which often involves solving Poisson's equation (1), requires significant computational effort. As a result, these computations demand parallel computing systems with high-speed, massively parallel processors and large memory capacities to efficiently implement numerical approaches within this framework [2]. Domain decomposition is a popular strategy for dividing the overall problem into smaller sub-domain problems, where parallelization can be effectively applied [3]. This technique offers advantages such as simpler construction and greater speed-up. Many numerical methods have been developed to exploit these benefits. However, finite difference (FD) methods are preferred for solving the Poisson equation over rectangular domains, as they typically lead to a large, block-structured, and sparse system of equations [4]. In addition, for smooth problems, high-order FD methods can compute numerical solutions hundreds of times faster than second-order FD schemes for a fixed level of accuracy. Implicit finite difference (IFD) methods provide high-order schemes without the drawback of requiring large stencils; in fact, some of them are associated with compact schemes [5].
In this paper, we present a compact implicit finite difference (IFD) scheme for numerically solving the 2D Poisson problem, which is primarily based on the domain decomposition strategy. Furthermore, the Successive Over-Relaxation (SOR) method is chosen for its efficiency and simplicity as a solver. However, since SOR is inherently sequential, it cannot be directly parallelized. To address this limitation, a coloring method is employed, as suggested in previous works [6,7,8,9].
This paper is organized as follows: Section 2 provides the high-order finite difference discretization. Section 3 is devoted to explaining the domain decomposition technique. In Section 4, we test our algorithm with several numerical examples and evaluate its speedup and efficiency. Finally, Section 5 presents the conclusions and suggestions for future work.
The classic one-dimensional finite-difference approximation for a second derivative of a real-valued function Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \phi }
at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle x}
with a small value Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h > 0}
is given by
|
(2) |
Using Taylor series, it follows that the second-order derivative at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle x}
can be approximated by (2) with the following local truncation error
|
(3) |
Note that the accuracy of centered finite difference approximation (3) can be fourth-order re-arranging the terms as follows
|
(4) |
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle b=1/12} . Furthermore, it can be shown that the second-order derivative operator Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\cdot )_{xx}}
on the left-hand side can be replaced by a finite difference operator Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta ^2}
without compromising the overall order of accuracy [5], as follows
|
(5) |
This formulation is called implicit because approximations to the second-order derivative appear in both sides of equation (5). Finally, we remark that the application of the implicit operator for a general real-valued function Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u}
(in particular Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u=\phi _{xx}}
) is given by
|
(6) |
In this section, we present the numerical scheme for solving the 2D Poisson problem. The computational domain, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega } , is divided using a uniform mesh with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N}
sub-divisions along the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle x}
- and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle y} -axes, where the grid points are defined as follows:
|
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_x = (\mathfrak{b}-\mathfrak{a})/N_x}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_y=(\mathfrak{d}-\mathfrak{c})/N_y}
. To maintain a simple explanation, we assume that the computational domain is a square, so Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=N_x=N_y}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h = h_x = h_y}
. We use the following notation: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_{i,j} = u(x_i,y_j)}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle f_{i,j} = f(x_i,y_j)}
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (x_i,y_j)}
represents the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ij}
-th grid point.
Let us consider the finite difference operators
|
(9) |
The subscripts at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathfrak{D}}
are not derivatives, they only indicate the dependence of the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle x}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle y}
direction respectively. Applying simultaneously these operators at equation (1), it yields
|
(10) |
Next, we use formula (5) to the second derivatives in Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle x}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle y}
to get
|
(11) |
Thus the two-dimensional fourth-order implicit finite difference scheme (IFD) results from applying (6) in 11. The following theorem states the final numerical discretization for 2D Poisson problem using the IFD method.
Theorem 1: Let us assume that the solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u}
has sufficient regularity and the grid size Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h}
is positive. The full discretization of the problem (1) is given by the fourth-order scheme
|
(12) |
where
|
(13) |
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle U_{i,j}}
is an approximation of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_{i,j}}
. More details about the proof of the Theorem 1 can be found in [5]. Additionally, if we set Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle b=0}
in (12) the resulting scheme becomes the classical second-order method. This parameter Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle b} allows us directly compare both numerical schemes by setting it either Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 0} or Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 1/12}
. Furthermore scheme (12) is compact scheme because its stencil requires only the surrounding nine points. The stencil of the classical second-order and the fourth-order implicit finite-difference schemes are illustrated in Figure 1.
| Figure 1: Second- and fourth-order finite-difference stencil for the two-dimensional Poisson equation. The numbers correspond to the weights next to unknown variables. |
The discretization given in (12) and (13) results in large linear system Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{AU} = \mathbf{F}}
as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N}
becomes larger. Therefore, to efficiently obtain the solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle U_{i,j}}
, this system should be solved using an appropiate high-performance numerical solver. It is not difficult to demonstrate that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{A}}
is a diagonally dominant matrix, thus a wide range of iterative linear solvers are available in this case.
Due to their robustness, efficiency and ease of implementation, we consider the classical stationary methods: the Jacobi method, the Gauss-Seidel (GS) method, and the Successive Over-Relaxation (SOR) method [10]. The Jacobi method is a well-known iterative approach that is straightforward to implement and parallelize. However, its slow convergence rate makes it unsuitable for practical, large-scale applications. The GS and SOR methods offer more efficient solvers, with SOR converging faster than GS when an appropriate relaxation factor is chosen. The main challenge with the SOR method is that, in its original form, it cannot be easily parallelized due to its inherently sequential nature. Thus, alternative approaches must be considered for parallel implementation. In this work, we utilize the Multicolor SOR method due to its simplicity in parallelization and its minimal modifications compared to the Jacobi method.
We begin by considering the Jacobi method. While its convergence is slow, it provides the advantages of high execution speed and straightforward implementation, particularly when it comes to vectorization and parallelization. Let the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k} -th equation of the resulting linear system
|
(14) |
then, an iteration of the Jacobi method is given by
|
(15) |
where the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (n+1)}
values are from the current iteration, and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (n)}
values are from the previous iteration. Here the vectors Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \phi }
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \ell }
are constructed using a natural indexation of matrices Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{U}}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{F}}
, respectively. The iteration process will be stopped once a convergence tolerance Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon }
is achieved. Note that the Jacobi method can be easily implemented in parallel since all updated values depend solely on the previous iteration Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (n)}
. The convergence rate of the Jacobi method can be improved by employing the SOR method. In these iterative solvers, each component at the new iteration depends on some previously computed components within the same iteration. The SOR method is expressed as:
|
(16) |
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 0<\omega{<2}}
is the relaxation factor. The Gauss-Seidel method is recovered by setting Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega = 1}
in the SOR method. It is important to note that, unlike the Jacobi method, the computations in the SOR method are inherently sequential, meaning the updates cannot be performed simultaneously.
To apply this solver to our compact fourth-order implicit scheme, we can rewrite (16) but now using a residual vector as follows:
|
(17) |
where
|
(18) |
The stopping criterion for (17) with a given tolerance Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon }
is that the residual (18) is sufficiently small, i.e., Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle |R_{i,j}|<\varepsilon }
, for all indices Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i,j} . In addition to the tolerance, we include a maximum number of iterations as a complementary stopping criterion.
Several algorithms exist for parallelizing the SOR method on structured grids, and due to their simplicity, coloring methods have often been a popular choice [6]. However, there are limited works on the implementation and performance of these techniques when applied to compact schemes using a nine-point stencil [7,8,9].
The classical two-color SOR method (RBSOR) divides the domain into a chessboard pattern of red and black points. In this configuration, each unknown is updated based on values of the opposite color, allowing the points of the same color to be computed in parallel. However, this algorithm is only applicable to stencils like the classical second-order finite-difference scheme, see Figure 2. Each iteration is carried out in two steps. Updating the first set of points (red) relies solely on previous black values. In contrast, the remaining points (black) require only the already updated red points for their calculations. The challenge with a nine-point stencil lies in the use of corner points, which complicates the direct application of the RBSOR method.
| Figure 2: Second-order finite-difference stencil for the two-dimensional Poisson equation using RBSOR iterative method. |
In the case of our nine-point compact scheme, four colors are sufficient to decouple the grid points. Thus, the four-color SOR (MCSOR) method is employed. For simplicity, we assign the colors red (R), black (B), green (G), and orange (O). For each grid point of a given color, the nearest grid points along both coordinate directions must be assigned different colors to ensure independence during the update process, as shown in Figure 3. This figure also illustrates the four sub-steps employed to update all grid points for each iteration Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (n)} .
| Figure 3: Fourth-order finite-difference stencil for the two-dimensional Poisson equation. The numbers correspond to the weights next to unknown variables. |
It is important to remark that we use the same variable to store all color values at each iteration. This approach allows the grid points of each color to be updated using the most recently computed values from the other colors. For example, when updating a reference grid point colored green, the neighboring grid points colored red and black are taken from the current iteration, while the orange values are taken from the previous iteration, as presented below:
|
(19) |
where
|
(20) |
and the superscripts in Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle U}
indicate the corresponding color. This configuration is expected to lead to faster convergence than the original SOR method, as demonstrated by the numerical results presented in the next section. However, further theoretical investigations are needed in this direction to better understand and confirm the potential improvements in convergence. Finally, we emphasize that, similar to the Jacobi method, the MCSOR method is well-suited to straightforward parallel implementation because all new values only depends on calculations previous sub-steps.
Instead of computing the solution over the entire domain Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega } , we subdivide the domain into smaller parts and compute the solution over each subdomain Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega _k} , allowing for parallel computation. Let us define a domain decomposition for the problem using Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle M}
strips, as follows
|
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega _k}
is an vertical strip. Then we find a solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u}
of the 2D Poison problem such Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u = \sum _{k=1}^{M} u_k}
, where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_k}
is the solution of the problem over the strip Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega _k}
. In practice, we solve the discrete problem (12) over the computational domain Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega ^h}
with step size Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h}
, which is decomposed into Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle M}
sub-grids Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega ^h_k}
, see Figure 4(a).
| Figure 4: Sketch of rectangular domain decomposition into smaller sub-domains with overlapping grid points. | |
The compact IFD numerical scheme using the MCSOR method, combined with the domain decomposition method, requires that each element of the sub-domain obtains information from its nine neighbors in each iteration. Consequently, the boundary (or artificial boundary) approximations of each sub-domain need information from the elements in the columns of neighboring sub-domains. This presents a challenge when parallelizing the scheme, as the processes computing the solution in one part of the domain do not have direct access to the memory of other processes. To address this issue, we propose creating sub-domains that overlap, ensuring they include the adjacent neighboring elements, see Figure 4(b). The adjacent information is exchanged using the Message Passing Interface (MPI), meaning that neighboring processors communicate with each other. These communications occur after each sub-step of the MCSOR solver within each sub-domain.
In this section, we perform several tests to analyze the capabilities of the numerical scheme and describe the differences between standard finite difference methods and the implicit finite difference method within the parallel paradigm.
The error is analyzed by comparing the exact and numerical solutions, and it is measured using the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l^\infty } - and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l^2} -norms at all grid points. The order of accuracy is calculated as follows:
|
(21) |
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Vert \cdot \Vert _{N_1}}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Vert \cdot \Vert _{N_2}}
represent the error norms with a resolution relative to grid resolution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N_1}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N_2}
respectively.
Next, we conduct several examples to test the robustness of the proposed methodology. We present different examples using iterative solvers such as Jacobi, SOR, RBSOR, and MCSOR methods to find the solutions of the resulting linear system. However, only Jacobi, RBSOR, and MCSOR are suitable for parallelization.
In this example, we show that the numerical method can recover the numerical solution with great accuracy. We consider the 2D Poisson problem over the unit square Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega = [0,1]\times [0,1]}
given by
|
Notice that we impose the Dirichlet boundary conditions at the boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \partial \Omega } . The exact solution of problem (22) and (23) is given by the function
|
Figure 5 shows the numerical and the absolute error for example using a grid resolution with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=80}
. Notice that even when we use a coarse grid resolution, we obtain an accurate numerical solution. Table 1 shows the convergence analysis for Example 1. As expected, the numerical order for the FD and IFD methods corresponds to second- and fourth-order accuracy, respectively. This confirms that the proposed IFD scheme provides an accurate approximation.
| FD (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): b=0
) |
IFD (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): b=1/2
) | |||||||
| N | Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l^\infty
-norm |
Order | Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l^2
-norm |
Order | Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l^\infty
-norm |
Order | Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l^2
-norm |
Order |
| 20 | 8.27E-03 | –- | 4.13E-02 | –- | 2.69E-05 | –- | 1.34E-05 | –- |
| 40 | 2.06E-03 | 2.01 | 1.03E-02 | 2.01 | 1.69E-06 | 3.99 | 8.44E-07 | 3.99 |
| 80 | 5.14E-04 | 2.00 | 2.57E-03 | 2.00 | 1.06E-07 | 4.00 | 5.28E-08 | 4.00 |
| 160 | 1.29E-04 | 2.00 | 6.43E-04 | 2.00 | 6.64E-09 | 4.00 | 3.30E-09 | 4.00 |
| 320 | 3.21E-05 | 2.00 | 1.61E-04 | 2.00 | 4.44E-10 | 3.98 | 2.06E-010 | 4.00 |
The SOR iterative solver has a parameter Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 0<\omega{<2}}
that can be tuned to improve performance. Figure 6 illustrates how the number of iterations depends on the parameter Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega } for a given tolerance of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon = 10^{-9}} and a grid resolution of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N = 320}
. Figure 6(a) compares the FD and IFD methods, showing that the implicit method not only provides better accuracy but also achieves greater efficiency in terms of the number of iterations required for convergence when the using SOR method. Figure 6(b) shows that both multicolor methods exhibit similar behavior to the standard SOR method. Here, MCSOR requires fewer iterations for convergence compared to RBSOR.
Figure 7(a) shows residual errors for the SOR and Multicolor SOR methods with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega = 1.8}
and a grid resolution of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=320}
. As expected, in all methods, the residual errors decrease as the number of iterations increases, with the IFD methods offering better efficiency in terms of requiring fewer iterations. On the other hand, Figure 7(b) shows the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l^\infty } -norm of the absolute error between the numerical and exact solutions of the problem. Note that the IFD method not only improves efficiency but also offers greater accuracy, a distinct advantage inherent to this implicit technique. Although the results are not shown here, the Jacobi iteration solver demonstrates similar behavior between the IFD and FD methods.
In the second example, we test our numerical scheme on a more challenging problem characterized by oscillatory behavior, a key aspect often encountered in computational fluid dynamics simulations. Here, we demonstrate not only the performance of the parallel implementation of the implicit finite-difference method but also a comparison with the standard central finite-difference method (the five-point stencil scheme).
Let us consider the 2D Poisson problem over the unit square Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega }
defined by
|
Notice that we impose the Dirichlet boundary conditions at the boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \partial \Omega } . Here, the exact solution of problem (24) and (25) is given by the function
|
Problems with large oscillations (high frequency) bring challenges for numerical approximation, as capturing these oscillations accurately requires a large number of grid points. While this may seem like a minor issue for static problems, it becomes significant when solving large Poisson problems as part of evolutionary simulations, such as those involving the Navier-Stokes equations over long time periods.
Figure 8 shows the exact solution and the absolute errors obtained from both the standard and implicit finite difference schemes using a resolution of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=80} . As expected, the implicit method yields an accurate approximation. It is also noteworthy that, in both cases, the errors are distributed throughout the entire domain.
Table 2 shows the norm errors and simulation time in seconds of the proposed FD methods using a stopping criterion of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon = 10^{-9}}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega=1.85}
. Similar to the previous example, IFD shows higher accuracy as the standard FD method. For instance, to obtain an approximation with an absolute error less than Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 10^{-3}} , results shows that the FD method requires Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=640}
grid points, whereas the IFD method only needs Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=160}
grid points. Furthermore, the computational time required to achieve this absolute error using the Jacobi method with the FD method (47.24 seconds) is approximately 337 times longer than with the IFD method (0.14 seconds). On the other hand, using more advanced and faster solvers such as RBSOR and MCSOR provides the same precision, but both are significantly faster than the Jacobi method. For instance, MCSOR is approximately 30 times faster than Jacobi when using Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=640}
.
| FD (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): b=0
) |
IFD (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): b=1/2
) | |||||||
| N | Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l^\infty
-norm |
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l^2
-norm |
Jacobi (sec.) | RBSOR(sec.) | Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l^\infty
-norm |
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l^2
-norm |
Jacobi(sec.) | MCSOR(sec.) |
| 40 | 2.34E-01 | 1.17E-01 | 0.000125 | 0.00424 | 1.30E-02 | 6.52E-03 | 0.000131 | 0.00539 |
| 80 | 5.30E-02 | 2.65E-02 | 0.01 | 0.0173 | 1.01E-03 | 5.07E-04 | 0.00878 | 0.0224 |
| 160 | 1.30E-02 | 6.48E-03 | 0.17 | 0.0683 | 6.54E-05 | 3.27E-05 | 0.14 | 0.0956 |
| 320 | 3.22E-03 | 1.61E-03 | 2.81 | 0.28 | 4.13E-06 | 2.06E-06 | 2.42 | 0.39 |
| 640 | 8.04E-04 | 4.02E-04 | 47.24 | 1.19 | 2.58E-07 | 1.29E-07 | 39.91 | 1.58 |
To analyze the performance of the proposed method, the resolution of the numerical problem will be set into Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N = 641}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=1281}
(more than 400 thousand and 1.6 million of grid points, respectively). The stopping criterion is set to Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon =10^{-9}}
. This stringent tolerance is selected to ensure that both FD and IFD methods reach their theoretical absolute error.
Figure 9 shows the speedup and elapsed time for solving the 2D Poisson problem using different schemes and resolutions with the Jacobi solver. Note that the speedup and elapsed time are nearly identical in both cases; however, the IFD method yields more accurate results. In this case, the maximum error for the FD method is approximately Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 10^{-4}} , while for the same setting with the IFD method, the maximum error is close to Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 10^{-6}} . Finally, it is worth noting that the computational time with 32 processors is 24.6 seconds, which is 26 times faster than the sequential solver (641.46 seconds) when using the IFD with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=1281} . This represents a reduction of nearly 96% in the original runtime.
Figure 10 illustrates the performance for the RBSOR and MCSOR solvers. It is important to note that RBSOR and MCSOR are quite different solvers, as they consider the stencil structure for both FD and IFD techniques. Therefore, they recover the solution by accessing the global matrix in different ways and perform different arithmetic operations. However, both methods achieve the same accuracy as described for the Jacobi solver. In terms of parallel performance, the MCSOR method achieves a computational time of 2.37 seconds with 32 processors, making it 16 times faster than the sequential solver (37.99 seconds) for the IFD with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N=1281} . This represents a nearly 94% reduction compared to MCSOR on a single processor. Moreover, this parallel execution is 270 times faster than the sequential Jacobi solver, yielding a 99.6% reduction in runtime.
We present an implicit finite-difference scheme combined with domain decomposition and the Multicolor SOR method for solving the 2D Poisson equation. This scheme provides a fourth-order approximation, consistent with its sequential counterpart. Our experiments demonstrate that, under the same conditions, the IFD method produces more accurate solutions than the standard finite-difference method. Moreover, the speedup and elapsed time are comparable for both numerical schemes, making the IFD method an attractive option since it delivers better results with the same computational effort. In future work, we propose developing a sixth-order implicit finite-difference scheme to solve 2D and 3D Poisson problems.
This work was partially supported by the Consejo Nacional de Humanidades, Ciencias y Tecnologías (CONAHCYT) under the program Investigadores e Investigadoras por México, and CONAHCYT México under Ciencia de Frontera Project Number: CF-2023-I-2639. We also acknowledge to the supercomputer facilities provided by CIMAT, specifically Cluster Merida TOLOK, for their invaluable support in conducting our research.
[10] Uh Zapata, M., Pham Van Bang, D., & Nguyen, K. D. (2016). Parallel SOR methods with a parabolic-diffusion acceleration technique for solving an unstructured-grid Poisson equation on 3D arbitrary geometries. International Journal of Computational Fluid Dynamics, 30(5), 370-385.
Published on 06/11/24
Submitted on 09/10/24
Licence: CC BY-NC-SA license
Are you one of the authors of this document?