A Two-Level Direct Solver for the Hierarchical Poincaré-Steklov Method

In this subsection we will describe a global spectral method to approximate an elliptic PDE on a box domain. We will apply this later to individual subdomains within a larger boundary value problem. Suppose we have a PDE such as in (1), where $\Omega$ is a rectangular domain such as $[0,1]^{d}$, $d=2$ or $3$. We discretize $\Omega$ into a product grid of $p^{d}$ Chebyshev nodes $\{x_{k}\}_{k=1}^{p^{d}}$, where $p$ is the chosen polynomial order (other works have used polynomial coefficients as degrees of freedom instead [18]). Let $\mathbf{u}\in\mathbb{R}^{p^{d}}$ be the numerical approximation of $u$ on this discretization, such that $\mathbf{u}_{k}\approx u(x_{k})$. There exist spectral differentiation matrices $\mathbf{D}^{(1)}$, $\mathbf{D}^{(2)}$, and (in 3D) $\mathbf{D}^{(3)}$ such that $[\mathbf{D}^{(i)}\mathbf{u}]\approx[\partial_{i}u](x_{k})$, and more generally $\mathbf{D}^{(i)}\mathbf{u}$ is a discretization of $\partial_{i}u$ [19]. Thus we can approximate the partial differential operator with a spectral differentiation matrix. For example, in the 3D Helmholtz equation shown in (2), we can discretize the partial differential operator with

where $\mathbf{I}$ is the $p^{3}\times p^{3}$ identity matrix. Since Chebyshev nodes include the endpoints of an interval, there will always be points on $\partial\Omega$: $4p$ points on the edges for $d=2$, and $6p^{2}$ points on the faces for $d=3$. Let $B$ denote the indices of $\mathbf{u}$ that lie on $\partial\Omega$, and $I$ denote the indices of $\mathbf{u}$ on the interior of $\Omega$. We can then show

$\mathbf{S}$ is sometimes called the solution operator [20], and is of size $\#\text{ of interior points }\times\#\text{ of boundary points}=(p-2)^{d}\times 2dp^{d-1}$. Thus given a boundary value problem such as (1) with partial differential operator $A$ and a discretization $\mathbf{A}$ based on Chebyshev nodes $\{x\}_{k=1}^{p^{d}}$, we can derive $\mathbf{S}$ that maps the problem’s Dirichlet data to its solution on the interior of $\Omega$ (with the addition of a body load, if inhomogeneous). This is a global spectral method, and it excels at approximating highly-oscillatory problems such as the Helmholtz equation. However, it is poorly conditioned for larger $p$, and it can struggle to capture function details near the center of $\Omega$ since there are fewer Chebyshev nodes in that region. In addition, both $\mathbf{A}$ and $\mathbf{S}$ are dense matrices, making use of this method for large $p$ computationally expensive.

To mitigate the shortcomings of a global spectral method in solving an elliptic PDE such as in (1), we will partition $\Omega$ into two rectangular non-overlapping subdomains $\Omega^{(1)}$ and $\Omega^{(2)}$ (we will alternatively call these subdomains “leaves” or “boxes”). We will discretize each subdomain using $p^{d}$ Chebyshev nodes, as we did for $\Omega$ in Section 2.1. This leads to a total of $2p^{d}-p^{d-1}$ points, since there is one shared face of points between $\Omega^{(1)}$ and $\Omega^{(2)}$. A diagram of this discretization in 2D is presented in Figure 1, with $\mathbf{u}_{5}$ denoting the shared face. We already know the data of the other faces ($\mathbf{u}_{3},\mathbf{u}_{4}$ in the figure) on both boxes as part of the Dirichlet boundary condition. If we can solve for the points on the shared face first, then we can solve for the interior of each subdomain separately using a spectral method as in Section 2.1, enforcing the PDE locally. To help solve for the shared face, we will add an additional condition to our problem: the Neumann derivative of $u$ on $\partial\Omega^{(1)}$ and $\partial\Omega^{(2)}$ must be equal on the shared face. In other words, we enforce continuity of the Neumann derivative.

To enforce this continuity, we need to find a numerical approximation to the Neumann derivative. This can be done using spectral differentiation matrices. In particular we form a Neumann operator matrix $\mathbf{N}$ on each subdomain such that

Here $l,r,d,u,b,f$ are indices for each face of the subdomain boundary, with $l,r$ pointing in the $x_{1}$ direction, $d,u$ pointing in the $x_{2}$ direction, and $b,f$ pointing in the $x_{3}$ direction (respectively “left”, “right”, “down”, “up”, “back”, and “front” - in 2D $b$ and $f$ are omitted). $\mathbf{N}$ obtains the Neumann derivative for each face, taking in both the boundary points $\mathbf{u}_{B}$ and interior points $\mathbf{u}_{I}$ as inputs. To obtain the Neumann derivatives using only $\mathbf{u}_{B}$, we combine $\mathbf{N}$ with the solution operator $\mathbf{S}$ to get

$\mathbf{Id}$ is the $\text{size}(\mathbf{u}_{B})\times\text{size}(\mathbf{u}_{B})$ identity matrix and $\mathbf{T}$ denotes the Dirichlet-to-Neumann (DtN) map for $\Omega^{(i)}$ (also called the Poincaré-Steklov operator). For homogeneous PDEs we approximate the Neumann derivative entirely with $\mathbf{T}\mathbf{u}_{B}$, whereas the inhomogeneous case requires some contribution from the body load as seen in (2.2). Since $\mathbf{T}$ maps from the box boundary to itself, it is a square matrix of size $2dp^{d-1}\times 2dp^{d-1}$.

Borrowing terminology from Figure 1, let $\mathbf{u}_{1},\mathbf{u}_{2},\mathbf{u}_{3},\mathbf{u}_{4},\mathbf{u}_{5}$ denote the various subsets of $\mathbf{u}$ on discretized $\Omega$. Let $\mathbf{v}_{3},\mathbf{v}_{4},\mathbf{v}_{5}$ denote the Neumann derivatives on the same points as $\mathbf{u}_{3},\mathbf{u}_{4},\mathbf{u}_{5}$, and let $\mathbf{T}^{(1)},\mathbf{T}^{(2)}$ denote the DtNs on $\Omega^{(1)},\Omega^{(2)}$. Then, with a reindexing of some rows and columns of the DtNs,

By setting $\mathbf{v}_{5}$ (Neumann derivative of the shared face) equal to itself, we get

Thus we can compute the shared face $\mathbf{u}_{5}$ $\mathbf{u}_{3}$ and $\mathbf{u}_{4}$, given by the Dirichlet boundary condition, and a new solution operator $\mathbf{S}^{(3)}$ computed entirely from submatrices of the DtNs. This assumes a homogeneous PDE - if the problem is inhomogeneous, additional terms with the body load are present such as in (2.2). Once $\mathbf{u}_{5}$ is solved, we can treat it as Dirichlet data for $\Omega_{1}$ and $\Omega_{2}$, along with $\mathbf{u}_{3}$ and $\mathbf{u}_{4}$, to get $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$.

This divide-and-conquer approach has a few advantages over a global spectral method. The solution operators and DtN maps for each subdomain can be computed separately, providing an easy source of parallelism. While the merged solution operator $\mathbf{S}^{(3)}$ is still somewhat ill-conditioned, it is better conditioned than a global spectral matrix with the same domain size and degrees of freedom. In addition, nodes are more evenly spaced across the entire domain than in a global spectral method.

One interpretation of this merge operation is as an analogy to $LU$ factorization [21]. Consider Figure 1, simplifying it by considering the domain Dirichlet boundary terms $\mathbf{u}_{3},\mathbf{u}_{4}$ parts of $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ respectively. Then we can represent our numerical partial differential operator (PDO) $\mathbf{A}$ with

Here $\mathbf{u}$ is the solution and $\mathbf{f}$ is the body load on the interior of $\Omega$, split into the two subdomains and shared boundary. Because an elliptic PDO is a local operator there is no interaction through $\mathbf{A}$ on $\mathbf{f}_{1}$ by $\mathbf{u}_{2}$ or $\mathbf{f}_{2}$ by $\mathbf{u}_{1}$, so $\mathbf{A}=0$ in those submatrices. However, both $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$, are local to the shared boundary where $\mathbf{u}_{5}$ lies, so they both interact with the body load there, $\mathbf{f}_{5}$, through $\mathbf{A}$. To interpret the first two block rows, note that

This is similar to our solution operator formulation shown in (2.1), with the part of the subdomain boundary on $\partial\Omega$ now folded into $\mathbf{u}_{1}$. Next, let $\mathbf{\tilde{u}}_{1}=\mathbf{A}_{11}^{-1}\mathbf{f}_{1}$ and $\mathbf{\tilde{u}}_{2}=\mathbf{A}_{22}^{-1}\mathbf{f}_{2}$. If we supply these vectors in place of $\mathbf{u}_{1},\mathbf{u}_{2}$ in (8) and set $\mathbf{u}_{5}=\mathbf{0}$, then the equation still holds with the same body loads $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$. Thus $\mathbf{\tilde{u}}_{1}$ and $\mathbf{\tilde{u}}_{2}$ are the particular solutions to our PDE on the interiors of $\Omega^{(1)}$, and $\Omega^{(2)}$ with $\mathbf{u}_{5}=\mathbf{0}$. Now consider an upper triangular matrix $\mathbf{U}$ that satisfies

This linear system can be confirmed with (9). $\mathbf{U}$ decouples our solutions $\mathbf{u}_{1},\mathbf{u}_{2}$ into the particular solutions $\mathbf{\tilde{u}}_{1},\mathbf{\tilde{u}}_{2}$ which are derived from our global Dirichlet BC without the subdomain-only (shared face) Dirichlet BC, and $\mathbf{u}_{5}$ which is that subdomain-only boundary condition. Thus $\mathbf{U}^{-1}$ collects both components of our solutions - in effect it represents the end of solving the merged system, where we use $\mathbf{u}_{3}$ and $\mathbf{u}_{5}$ to get $\mathbf{u}_{1}$, and $\mathbf{u}_{4}$ and $\mathbf{u}_{5}$ to get $\mathbf{u}_{2}$. Next, let $\mathbf{L}$ be a lower triangular matrix defined such that

To interpret this, note that

This has some resemblance to (7) when noting that $\mathbf{\tilde{u}}_{1}$ encodes the impact of Dirichlet data $\mathbf{u}_{3}$ without $\mathbf{u}_{5}$ (and likewise $\mathbf{\tilde{u}}_{2}$ is $\mathbf{u}_{4}$ without $\mathbf{u}_{5}$). If an additional operator was composed onto $[\mathbf{A}_{51}|\mathbf{A}_{52}]$ that resembled a Neumann derivative approximation following an inverse of the PDO, then it would closely resemble $[-\mathbf{T}_{53}^{(1)}|\mathbf{T}_{54}^{(2)}]$ from (7). Another way to view (12) is that it breaks loads from decoupled $\mathbf{\tilde{u}}_{1},\mathbf{\tilde{u}}_{2}$ on the two subdomains from $\mathbf{f}_{5}$, leaving another decoupled term in $\mathbf{\tilde{f}}_{5}$. By giving us $\mathbf{\tilde{f}}_{5}$ from $\mathbf{f}$, $\mathbf{L}^{-1}$ can enable us to formulate a reduced problem that will solve for $\mathbf{u}_{5}$ given $\mathbf{\tilde{f}}_{5}$ - it performs the first part of the merge operation. To define this reduced problem, consider a block diagonal matrix $\mathbf{D}$ where

This now represents a fully decoupled system for the particular solutions $\mathbf{\tilde{u}}_{1},\mathbf{\tilde{u}}_{2}$ and shared boundary $\mathbf{u}_{5}$. However, we need to find a submatrix $\mathbf{T}_{55}$ that satisfies this system. Such a $\mathbf{T}_{55}$ can be found as

This formulation for $\mathbf{T}_{55}$ includes the effect of $\mathbf{u}_{5}$ on its body load ($\mathbf{A}_{55}$) as well as terms that resemble the solution operators in (2.1) acting on $\mathbf{u}_{5}$ only ($\mathbf{A}_{11}^{-1}\mathbf{A}_{15},\mathbf{A}_{22}^{-1}\mathbf{A}_{25}$). If we combine it with (12) then we see

This more closely resembles (7). Lastly, by combining the matrices together,

Thus $\mathbf{L}\mathbf{D}\mathbf{U}$ is an LDU factorization of $\mathbf{A}$, with $(\mathbf{L}\mathbf{D})^{-1}$ resembling the merge operations we derive via DtN maps, and $\mathbf{U}^{-1}$ resembling the interior solve for $\mathbf{u}_{1},\mathbf{u}_{2}$ once we obtain the shared boundary data. This conceptualization will prove useful in justifying our algorithm for HPS.

Given an elliptic PDE such as in (1) with a rectangular domain $\Omega$, we can hierarchically partition $\Omega$ into a series of rectangular “box” subdomains, illustrated in Figure 2. Note that the boxes in each level are children of boxes in the previous level, with each box in level $\ell$ having 2 children in level $\ell+1$, until we reach the final level of “leaf” boxes.

The HPS algorithm is split into two stages: the “build” stage, and the “solve” stage. In the build stage, we first construct solution and DtN maps for every leaf box in the final level of our discretization ($\ell=4$ in Fig. 2). We then proceed in an upward pass, constructing solution maps for each box in the next level up using the leaf DtNs, following (7). These solution maps solve for the boundaries between paired leaf boxes (boxes 8 and 9, 10 and 11, 12 and 13, and 14 and 15 in Fig. 2). We then recursively construct solution and DtN maps for each level of the discretization, until we get a solution operator for $\ell=1$ which solves for one large face in the middle of the domain (between boxes 2 and 3 in Fig. 2). Once all solution operators are constructed, we proceed with the solve stage that follows a downward pass: using the Dirichlet boundary condition and the $\ell=1$ solution operator to solve for the middle boundary, then proceeding down and using each level’s solution operators to solve for subsequent box faces, until every leaf box has all faces solved. Then we solve for the leaf interiors using their local solution operators.

The method’s upward and downward passes recursively extend the formulation and use of solution operators we derived in Section 2.2. Similarly, they also resemble the LDU factorization described in Section 2.3. We can imagine the upward pass as similar to $\mathbf{L}^{-1}$ described there, with both encoding merge operations that reduce the body loads of the problem onto shared boundaries. $\mathbf{D}^{-1}$ then resembles formulating the final components of the merged solution operators and solving for the shared face points - it combines elements of the upward and downward pass. Lastly, $\mathbf{U}^{-1}$ is the end of the downward pass, obtaining the total solution. Given the similarities, it is not surprising that the classical algorithm for obtaining the DtN and solution operators of traditional HPS resembles multifrontal methods for sparse LU factorization [22, 23].

The original HPS method has a few significant drawbacks. The solution and DtN maps can be formed in parallel within each level, and prior work has used shared-memory models or MPI to accelerate this component [24, 25]. However, across levels these operators must be constructed recursively. This forces a large serial component to the algorithm. In addition, the method on a discretization with $L$ levels requires the storage of $2^{L}-1$ dense solution operators, which can be large especially for 3D problems. One can choose to recompute certain levels during the downward pass to reduce storage cost, but this greatly increases the computation time as a trade-off. In addition, higher level DtN maps tend to get increasingly ill-conditioned in this method, which can cause numerical challenges in the method’s stability (formulations using other merge operators such as impedance maps may be better conditioned [26]). Since each solution operator is applied separately in the downward pass, techniques such as pivoting are limited to improving the conditioning of one solution operator at a time - this is one way HPS differs significantly from sparse LU algorithms, despite the other similarities. Our contributions mitigate these problems with the original HPS method.

In this section we will discuss the contribution of our work: designing a direct solver for elliptic PDEs in 2D or 3D that utilizes the HPS method in a 2 level framework. Instead of constructing higher level solution and DtN operators, we use the leaf DtN maps to formulate a system that solves for all leaf boundaries concurrently. We use this system to acquire the leaf boundaries, then apply the leaf solution operators in parallel to solve for leaf interiors. With this approach we use the discretization for HPS to reduce the problem in a way analogous to static condensation, visualized in Figure 3. Similar to traditional HPS, our method can naturally be divided into a build stage and a solve stage.

The build stage has an asymptotic cost of

We formulate the DtNs with embarrassingly parallel dense matrix operations that are well-suited for GPUs using batched linear algebra routines. As $p$ increases, its effect on the DtN size and run time is substantial - this is because of the need to invert a large block of the spectral differentiation submatrix $\mathbf{A}^{(i)}$, as shown in Algorithm 3. Nevertheless, we can still easily formulate DtNs in batch for $p\leq 15$. Even for larger $p$ where batched operations are impractical, formulating the DtNs in serial is relatively fast as shown in our numerical results. In practice the steepest cost in the build stage comes from the factorization of the sparse system, which for a 3D problem using a multilevel solver is $O(N^{2})$. However, the sparse factorization is unaffected by the choice of $p$ and this cost is only required once per differential operator and domain discretization. Since the sparse system is relatively small (it has $<\frac{6N}{p}\times\frac{6N}{p}$ total entries with $\approx 11p^{2}$ nonzero entries per row), it can be stored in memory and reused for multiple solves. For the solve stage, our asymptotic cost is

If we store the solution operators in memory then the interior box solves are only $O(p^{2}N)$, but in practice it is typically preferred to reformulate the solution operators and avoid having to store and potentially transfer additional memory.

We investigate the accuracy of our HPS solver on the homogeneous Poisson equation,

Here the Dirichlet boundary condition is defined so that it has a solution in its Green’s function. We compute numerical solutions to this equation across a range of block sizes $h$ and polynomial orders $p$, comparing their relative $\ell^{2}$-norm error to the Green’s function. The accuracy results are in Figure 7. We notice convergence rates of $\approx h^{7}$ for $p=6$ and $h^{12}$ for $p=8$ until the relative error reaches $\approx 10^{-13}$, at which point it flattens (for larger $p$ the error reaches $10^{-13}$ before a convergence rate can be fairly estimated). Subsequent refinement does not generally improve accuracy, but it does not lead to an increase in error either. This is likely due to the HPS discretization being fairly well-conditioned, especially compared to global spectral methods.

To assess our method’s accuracy on oscillatory problems we also evaluate it on the homogeneous Helmholtz equation with its Green’s function used as a manufactured solution with wavenumber $\kappa$,

In Figure (8) we show relative $\ell^{2}$-norm errors for our solver on Equation (20) with fixed wavenumbers $\kappa=16$ and $\kappa=30$, which correspond to roughly 2.5 and 4.5 wavelengths on the domain, respectively. We see comparable convergence results to the Poisson equation, though the lower bound on relative error is somewhat higher for each problem ($\approx 10^{-11}$ for $\kappa=16$ and $\approx 10^{-10}$ for $\kappa=30$). However, increasing the resolution further does not cause the error to diverge. Along with the results for Equation (19), this suggests our method remains stable when problems are over-resolved. Figure 7 shows the relative error for the Helmholtz equation with a varying $\kappa$, set to provide $\approx 10$ discretization points per wavelength. We see higher $p$ corresponds to a better numerical accuracy, while increasing $h$ does not improve results (since the wavenumber is also increased), but errors remain relatively stable. Since our method maintains accuracy while increasing degrees of freedom linearly with $\kappa$, it appears to not suffer from the pollution effect.

We then consider the gravity Helmholtz equation,

This equation adds a spatially-varying component to the wavenumber variation based on $x_{3}$ analogous to the effect of a gravitational field, hence the name [36]. Since $1<1-x_{3}<2.2$, the variation is greater through $\Omega$ than in the constant Helmholtz equation with the same $\kappa$. Unlike the previous Poisson and Helmholtz equations, equation (21) does not have a manufactured solution. It also produces relatively sharp gradients near the domain boundaries since there is a zero Dirichlet boundary condition, and it has a variable coefficient in spatial coordinates (specifically $x_{3}$). We investigate our solver’s accuracy by computing the relative $\ell^{2}$-norm errors of our results to an over-resolved solution at select points. Figure 9 shows the solution of this problem using select $p$ and $h$ values, while Figure 10 shows these errors in the case of $p$-refinement, for a range of values of $h$. We can see our solver shows convergence for each $h$.

To evaluate the computational complexity of our HPS solver, we measure the execution times of its three components: the construction of the DtN operators used for the build stage, the factorization time for our sparse matrix constructed by the DtNs, and the solve time. We know the asymptotic cost of constructing the DtN operators is $O\left(p^{6}N\right)$ as shown in Equation (3.3). Figure (11) illustrates the run times for the build stage in solving the Poisson and Helmholtz equations detailed in Section 4.1, using batched linear algebra routines on an Nvidia V100 GPU. As expected, we observe a linear increase in runtime with $h$-refinement (i.e. a higher $N$ with fixed $p$). The constant for the linear scaling in $h$, though, increases sharply with higher values of $p$.

A challenge with utilizing GPUs for DtN formulation is memory storage. In 3D each DtN map is size $\approx 6p^{2}\times 6p^{2}$. Computing them requires the inversion of the spectral differentiation operator on interior points which is size $\approx(p-2)^{3}\times(p-2)^{3}$. These matrices grow quickly as $p$ increases (at $p=20$ they are $\approx 5.7$ million and $34$ million entries, respectively). For larger problem sizes it is impossible to store all DtN maps on a GPU concurrently. Interior differentiation submatrices do not have to be stored beyond formulating their box’s DtN map, but they are much larger for high $p$ and present considerable overhead. To handle this memory concern we have implemented an aggressive two-level scheduler that stores a limited number of DtNs on the GPU at once, and formulates a smaller number of DtNs in batch at a time. For $p\geq 16$ on a V100 GPU the best results occur when only one DtN is formulated at a time. However, we may still store multiple DtNs on the GPU before transferring them collectively to the CPU.

In Figure (12) we see the time to factorize the resulting sparse matrix for the build stage, using MUMPS. The expected asymptotic cost is $O(N^{2})$. However, our results more closely resemble a lower cost of $O(N^{3/2})$ for $N$ up to $\approx$ 8 million. Overall the factorization time dominates the build stage of our method given its more costly scaling in $N$ and the reliance on CPUs through MUMPS. Lastly in Figure (13) we see the runtime for the leaf interior solves. These follow a linear scaling similar to the batched DtN construction.

The HPS method has been extended to other problems such as surface PDEs [37], inverse scattering problems [38], and more generally non-rectangular domains. We can apply our HPS solver to non-rectangular domain geometries through the use of an analytic parameterization between the domain we wish to model, such as a sinusoidal curve along the $y$-axis shown in Figure (14(a)), and a rectangular reference domain. These parameter maps extend the versatility of the solver, but feature multiple challenges in the implementation: they require variable coefficients for our differential operators; they may use mixed second order differential terms, i.e. $\frac{\partial^{2}u}{\partial x_{i}\partial x_{j}}$ where $i\neq j$; and (depending on the domain) may have areas that are near-singular, such as around the sharp point $(x,y)\approx(2,0)$ in Figure (14(a)). We investigate the accuracy of our solver on this sinusoidal domain for the Helmholtz equation as shown in Eq. (20), but now on the domain

We map (20) on $\Psi$ to a cubic reference domain $\Omega=[1.1,2.1]\times[-1,0]\times[-1.2,0.2]$ using the parameter map

We test this problem with a manufactured solution in the Green’s function of the Helmholtz equation, with corresponding Dirichlet data matching our test in Section 4.1. A plot of this problem and the computed solution for wavenumber $\kappa=26$ is shown in Figure 14(a), while convergence studies for this problem with $\kappa=16$ and $\kappa=30$ are shown in Figure 15. We observe a similar convergence pattern to the Helmholtz test on a cubic domain, although the speed of convergence is somewhat slower.

We also investigate our solver’s accuracy on the Helmholtz equation when the domain is a three-dimensional annulus, shown in Figure 14(b). We use an analytic parametrization to map this problem to a rectangular reference domain with higher resolution in the $y$-axis and a periodic boundary condition. Relative errors in $p$ and $N$ are comparable to those for the problem on a curved domain.