Domain Decomposition for the Closest Point Method

TL;DR

This paper extends domain decomposition methods to solve elliptic PDEs on surfaces using the Closest Point Method, enabling parallel computation for applications like reaction-diffusion systems and shape analysis.

Contribution

It introduces novel domain decomposition techniques tailored for surface PDEs discretized by the Closest Point Method, enhancing parallel solver efficiency.

Findings

Extended two DDMs to surface PDEs with CPM

Demonstrated potential for parallel computation in surface PDEs

Facilitated faster solutions for applications like shape classification

Abstract

The discretization of elliptic PDEs leads to large coupled systems of equations. Domain decomposition methods (DDMs) are one approach to the solution of these systems, and can split the problem in a way that allows for parallel computing. Herein, we extend two DDMs to elliptic PDEs posed intrinsic to surfaces as discretized by the Closest Point Method (CPM) \cite{SJR:CPM,CBM:ICPM}. We consider the positive Helmholtz equation $\left(c-\Delta_\mathcal{S}\right)u = f$, where $c\in\mathbb{R}^+$ is a constant and $\Delta_\mathcal{S}$ is the Laplace-Beltrami operator associated with the surface $\mathcal{S}\subset\mathbb{R}^d$. The evolution of diffusion equations by implicit time-stepping schemes and Laplace-Beltrami eigenvalue problems \cite{CBM:Eig} both give rise to equations of this form. The creation of efficient, parallel, solvers for this equation would ease the investigation of reaction-diffusion equations on surfaces \cite{CBM:RDonPC}, and speed up shape classification \cite{Reuter:ShapeDNA}, to name a couple applications.

Figures

Tables

Table 1. Table 1: Here the iterations to convergence of the ORAS solver are gathered for various parameters. Convergence was declared when the 2 − limit-from 2 2- norm of the residual was reduced by a factor of 10 6 superscript 10 6 10^{6} .

$N_O=4$, $\alpha =16$				$N_S=64$, $\alpha =16$				$N_S=64$, $N_O=4$
$N_S$	$64$	$96$	$128$	$N_O$	$4$	$6$	$8$	$\alpha$	$16$	$32$	$64$
RAS	$992$	$1237$	$1533$		$992$	$747$	$610$		$992$	$992$	$992$
ORAS	$526$	$672$	$833$		$526$	$418$	$393$		$526$	$707$	$868$

Equations

$$\left(c-{\Delta }_{\mathcal{S}}\right)u=f,$$

$$\{\begin{array}{ll}\left(c-{\Delta }_{\mathcal{S}}\right)u_j^{(n+1)}=f,& \mathrm{in}{\mathcal{S}}_j,\\ {\mathcal{T}}_{jk}{u}_j^{(n+1)}={\mathcal{T}}_{jk}{u}^{(n)},& \mathrm{on}{\Gamma }_{jk},k=1,\dots ,N_S,k\ne j,\end{array}$$

$${\mathbf{A}}_j=\left(c+\frac{2d}{h^2}\right)-\left(\frac{2d}{h^2}+{\Delta }_j^h\right)\left[\begin{array}{c}{\mathbf{E}}_j\\ {\mathbf{T}}_j\end{array}\right],$$

$${\frac{\partial u}{\partial {\hat{q}}_i}|}_{C{P}_{{\mathcal{S}}_j}(x_i)}+\alpha u\left(C{P}_{{\mathcal{S}}_j}(x_i)\right)=0,$$

Full text

11institutetext: Ian May 22institutetext: Simon Fraser University, 8888 University Dr, Burnaby, BC V5A 1S6, 22email: mayianm@sfu.ca 33institutetext: Ronald Haynes 44institutetext: Memorial University of Newfoundland, 230 Elizabeth Ave, St. Johns, NL A1C 5S7 44email: rhaynes@mun.ca 55institutetext: Steven Ruuth 66institutetext: Simon Fraser University, 8888 University Dr, Burnaby, BC V5A 1S6, 66email: sruuth@sfu.ca

Domain Decomposition for the Closest Point Method

Ian May and Ronald D. Haynes and Steven J. Ruuth

1 Introduction

The discretization of elliptic PDEs leads to large coupled systems of equations. Domain decomposition methods (DDMs) are one approach to the solution of these systems, and can split the problem in a way that allows for parallel computing. Herein, we extend two DDMs to elliptic PDEs posed intrinsic to surfaces as discretized by the Closest Point Method (CPM) SJR:CPM ; CBM:ICPM . We consider the positive Helmholtz equation

[TABLE]

where $c\in\mathbb{R}^{+}$ is a constant and $\Delta_{\mathcal{S}}$ is the Laplace-Beltrami operator associated with the surface $\mathcal{S}\subset\mathbb{R}^{d}$. The evolution of diffusion equations by implicit time-stepping schemes and Laplace-Beltrami eigenvalue problems CBM:Eig both give rise to equations of this form. The creation of efficient, parallel, solvers for this equation would ease the investigation of reaction-diffusion equations on surfaces CBM:RDonPC , and speed up shape classification Reuter:ShapeDNA , to name a couple applications.

Several methods exist for the discretization of surface intrinsic PDEs. The surface may be parametrized to allow the use of standard methods in the parameter space FloaterHormann:Para . Unfortunately, many surfaces of interest do not have simple, or even known, parametrizations. Given a triangulation of the surface, a finite element discretization can be formed DziukElliot . This approach leads to a sparse and symmetric system but is sensitive to the quality of the triangulation. Level set methods for surface PDEs Cheng:LSM solve the problem in a higher dimensional embedding space over a narrow band containing the surface. The solution of model equation (1) by this method requires using gradient descent, as the approach was formulated only for parabolic problems. The CPM is also discretized over a narrow band in the embedding space, but has the advantage of using a direct discretization of equation (1).

The solution of the linear system arising from the CPM discretization of the model equation (1) has relied primarily on direct methods, although a multigrid method was discussed in Chen:MG . Herein we formulate restricted additive Schwarz (RAS) and optimized restricted additive Schwarz (ORAS) solvers compatible with the CPM to step towards efficient iterative solvers and to allow for parallelism. The optimized variant of the classical RAS solver uses Robin transmission conditions (TCs) to pass additional information between the subdomains Gand:OptPar ; Cyr:OMSORAS , and can accelerate convergence dramatically. This formulation is described in Sections 4 and 5 after reviewing the CPM in Section 2 and (O)RAS solvers in Section 3. Then, we discuss a PETSc petsc-user-ref ; petsc-efficient implementation and show some numerical examples in Section 6. A more thorough exploration of these solvers, and an initial look at their use as preconditioners, can be found in May’s thesis May:Thesis .

2 The closest point method

The CPM was introduced in SJR:CPM as an embedding method for surface intrinsic PDEs. It allows the reuse of standard flat space discretizations of differential operators and provides a surface agnostic implementation. At the core of this method is the closest point mapping, $CP_{\mathcal{S}}(x)=\underset{y\in\mathcal{S}}{\operatorname*{arg\,min}}|x-y|$ for $x\in\mathbb{R}^{d}$, which identifies the closest point on the surface for (almost) any point in the embedding space. This mapping exists and is continuous in the subset of $\mathbb{R}^{d}$ consisting of all points within a distance $\kappa_{\infty}^{-1}$ of the surface, where $\kappa_{\infty}$ is an upper bound on the principal curvatures of the surface Chu:Vari .

From this mapping, an extension operator $E$ can be defined that sends functions defined on the surface, $f:\mathcal{S}\rightarrow\mathbb{R}$, to functions defined on the embedding space via composition with the closet point mapping, $Ef=f\circ CP_{\mathcal{S}}$. The extended functions are constant in the surface normal direction and retain their original values on the surface. This extension operator can be used to define surface intrinsic differential operators from their flat space analogs SJR:CPM .

Discretization typically requires a Cartesian grid on the embedding space within a narrow tube surrounding the surface. The extension operator can be defined by any suitable interpolation scheme, with tensor product barycentric Lagrangian interpolation Tref:Bary being used here. As such, the computational tube must be wide enough to contain the interpolation stencil for any point on the surface. Using degree $p$ interpolation and a grid spacing of $\Delta x$ requires that the tube contain all points within a distance of $\gamma=\Delta x(p+2)\sqrt{d}/2$ from the surface, thus limiting the acceptable grid spacings in relation to $\kappa_{\infty}$.

The grid points within the computational tube form the set of active nodes, $\Sigma_{A}$. For (1), we need only discretize the regular Laplacian on $\mathbb{R}^{d}$. Here we consider the second order accurate centered difference approximation requiring $2d+1$ points. Around $\Sigma_{A}$ and lying outside the tube, a set of ghost nodes, $\Sigma_{G}$, is formed from any incomplete differencing stencils. With a total of $N_{A}$ active nodes and $N_{G}$ ghost nodes, we define the discrete Laplacian and extension operators, $\Delta^{h}:\mathbb{R}^{N_{A}+N_{G}}\rightarrow\mathbb{R}^{N_{A}}$ and ${\rm\bf E}:\mathbb{R}^{N_{A}}\rightarrow\mathbb{R}^{N_{A}+N_{G}}$, where $\Delta^{h}$ applies the centered difference Laplacian over all active nodes, and ${\rm\bf E}$ is the discretization of $E$. ${\rm\bf E}$ extends data on the active nodes to both the active and ghost nodes, and has entries consisting of the interpolation weights for each node’s closest point.

The Laplace-Beltrami operator can be directly discretized as $\Delta^{h}_{\mathcal{S},dir}=\Delta^{h}{\rm\bf E}$, which was used successfully for parabolic equations with explicit time-stepping in SJR:CPM . However, for implicit time-stepping CBM:ICPM and eigenvalue problems CBM:Eig a modified form is needed. In CBM:ICPM it was recognized that there was a redundant interpolation being performed, and that its removal could stabilize the discretization. The stabilized form $\Delta^{h}_{\mathcal{S}}=-\frac{2d}{\Delta x^{2}}{\rm\bf I}+\left(\frac{2d}{\Delta x^{2}}{\rm\bf I}+\Delta^{h}\right){\rm\bf E},$ will be used in the remainder of this work.

3 (Optimized) Restricted additive Schwarz

Both RAS and ORAS are overlapping DDMs, and can work on the same set of subdomains (given an additional overlap condition for ORAS Cyr:OMSORAS ). We define these solvers from the continuous point of view and subsequently discretize, rather than defining them purely algebraically. This will ease the discussion of TCs within the context of the CPM later in Section 5.

First, the whole surface $\mathcal{S}$ is decomposed into $N_{S}$ disjoint subdomains, $\widetilde{\mathcal{S}}_{j}$, for $j=1,\ldots,N_{S}$. These disjoint subdomains are then grown to form overlapping subdomains $\mathcal{S}_{j}$, whose boundaries are labelled depending on where they lie in the disjoint partitioning. Taking $\Gamma_{jk}=\partial\mathcal{S}_{j}\cap\widetilde{\mathcal{S}}_{k}$ gives $\partial\mathcal{S}_{j}=\bigcup\limits_{k}\Gamma_{jk}$ and allows the definition of the local problems

[TABLE]

where $\mathcal{T}_{jk}$ are generally linear boundary operators defining the TCs. RAS is achieved by choosing $\mathcal{T}_{jk}$ as identity operators, corresponding to Dirichlet TCs, while ORAS uses Robin TCs, $\mathcal{T}_{jk}=\left(\frac{\partial}{\partial\hat{{\rm\bf n}}_{jk}}+\alpha\right)$, where $\hat{{\rm\bf n}}_{jk}$ is the outward pointing boundary normal on $\Gamma_{jk}$ and $\alpha\in\mathbb{R}^{+}$ is a constant weight on the Dirichlet contribution.

The subproblems in equation (2) are initialized with a guess for the global solution $u^{(0)}$ (defined at least over the boundaries $\Gamma_{jk},\forall j,k$), which is usually just taken as $u^{(0)}=0$. After all of the subproblems have been solved a new global solution is constructed with respect to the disjoint partitioning, $u^{(n+1)}=\sum\limits_{j}\left.u_{j}^{(n+1)}\right|_{\widetilde{\mathcal{S}}_{j}}$, where the use of the term restricted indicates that the portion of the local solutions in the overlap regions are discarded. From this new approximation for the global solution the local problems may be solved again with new boundary data, and the process repeats until the global solution is satisfactory.

4 Subdomain construction

To solve problems arising from the CPM we first need to decompose the global set of active nodes $\Sigma_{A}$. (O)RAS solvers rely on both a disjoint partitioning of the active nodes and an induced overlapping partitioning. Following the notation in Section 3, disjoint partitions will be denoted by $\widetilde{\Sigma}_{j}$, overlapping partitions by $\Sigma_{j}$, and the boundaries of the overlapping partitions by $\Lambda_{j}$.

To ensure the solvers work on a variety of surfaces, we seek an automated and surface agnostic partitioning scheme to generate the disjoint partitions. METIS metis is a graph partitioner that is frequently used within the DD community to partition meshes DoleanNataf . The stencils of $\Delta^{h}$ and ${\rm\bf E}$ may be used to induce connectivity between the active nodes and define a graph. Here we only consider nearest neighbor coupling through the stencil for $\Delta^{h}$. Fig. 1 shows a portion of one such disjoint partition, in black circles, for a circular surface.

With $\widetilde{\Sigma}_{j}$ obtained from METIS, overlapping subdomains $\Sigma_{j}$ can be formed. This construction proceeds in the following steps:

All nodes in $\widetilde{\Sigma}_{j}$ are added to $\Sigma_{j}$. 2. 2.

$N_{O}$ layers of overlap nodes are added around $\Sigma_{j}$. Layers are added one at a time from globally active nodes neighboring $\Sigma_{j}$. 3. 3.

A subset of the ghost nodes, $\Sigma_{G}$, are placed in $\Sigma_{j}^{G}$ which consists of nodes that neighbor a member of $\Sigma_{j}$. 4. 4.

The shapes of the disjoint and overlapping subdomains are not known in advance. The boundary $\partial\mathcal{S}_{j}$ is approximated discretely by the closest points of the final layer of overlap nodes, and held in the set $\Lambda_{j}$. 5. 5.

Nodes needed to complete stencils from the ambient Laplacian or extension operator, including extension from the points $x_{i}\in\Lambda_{j}$, are placed in the set $\Sigma_{j}^{BC}$. 6. 6.

For ORAS a layer of ghost nodes around $\Sigma_{j}^{BC}$ are also placed in $\Sigma_{j}^{BC}$.

The active nodes in the $j^{\rm th}$ subdomain consist of $\Sigma_{j}$ and the active portion of $\Sigma_{j}^{BC}$. $\Sigma_{j}^{BC}$ is kept separate as that is where the TCs in Section 5 are defined. Each of these sets are shown in Fig. 1, which shows a portion of one subdomain on a circle in the vicinity of the points in $\Lambda_{j}$ at one of its boundaries.

The Robin TCs, to be defined in Section 5, need some final information about the subdomain. Every node in $\Sigma_{j}^{BC}$ is identified with the point in $\Lambda_{j}$ that is closest to it. This identification will be used to override the global closest point function in the following section. For each point in $\Lambda_{j}$ we also need to know the direction that is simultaneously orthogonal to the boundary and the surface normal direction. We call this the conormal direction. It is in this direction that the Neumann component of the Robin condition will be enforced. However, the discrete nature of $\Lambda_{j}$ makes this construction difficult. Instead we define the conormal vectors from the point of view of the boundary nodes. Take $x_{i}\in\Sigma_{j}^{BC}$ as a node whose associated conormal direction, $\hat{q}_{i}$, is sought. Let $y_{i}$ be its closest point in $\Lambda_{j}$, and $\hat{n}_{i}$ be the unit surface normal there. Connecting the boundary location to the boundary node via $d_{i}=x_{i}-y_{i}$, we obtain a usable approximation to the conormal by computing the component of $d_{i}$ that is orthogonal to $\hat{n}_{i}$ and normalizing, i.e., $\hat{q}_{i}=(d_{i}-\left(d_{i}\cdot\hat{n}_{i}\right)\hat{n}_{i})/(\left|d_{i}-\left(d_{i}\cdot\hat{n}_{i}\right)\hat{n}_{i}\right|).$ In the unlikely event that $d_{i}$ lies perfectly in the surface normal direction, we set $\hat{q}_{i}=0$ which recovers the natural boundary condition on the computational tube as discussed at the end of Section 5.

5 Transmission conditions

Boundary conditions in the CPM are imposed by modifying the extension operator over the nodes $\Sigma_{j}^{BC}$ beyond the surface boundary CBM:Eig . As such, the local operators will take the form

[TABLE]

where ${\rm\bf E}_{j}$ is the extension operator for the nodes in $\Sigma_{j}$ as inherited from the global operator and ${\rm\bf T}_{j}$ is the modified extension operator for the nodes in $\Sigma_{j}^{BC}$. When solving for the local correction to the solution the right hand side of the local problem, ${\rm\bf A}_{j}v_{j}=r_{j}$, will be the restriction of the residual to $\Sigma_{j}$. The final rows of the right hand side, those lying over $\Sigma_{j}^{BC}$, become zeros corresponding to the homogenous TCs.

Homogeneous Dirichlet TCs can be enforced to first order accuracy by extending zeros over all of $\Sigma_{j}^{BC}$. With the right hand side already set to zero there, the modified extension reduces to the identity mapping, ${\rm\bf T}_{j}=\begin{bmatrix}{\bf 0}&~{}&{\bf I}\end{bmatrix}$, with the zero matrix padding the columns corresponding to the interior nodes.

We discretize the Robin condition

[TABLE]

using a forward difference in the $\hat{q}_{i}$ direction for each node in $\Sigma_{j}^{BC}$ and the first order accurate Dirichlet condition from above. Taking the partial derivative $\frac{\partial u}{\partial d_{i}}$, and applying the change of variables $d_{i}=\hat{q}_{i}+\hat{n}_{i}$, allows one to write the Neumann term in equation (4) in terms of the displacement vector $d_{i}$ from Section 4. Assuming for the moment that $d_{i}$ and $\hat{q}_{i}$ are not perpendicular, the derivative in the conormal direction can be approximated by $\left.\frac{\partial u}{\partial\hat{q}_{i}}\right|_{CP_{\mathcal{S}_{j}}(x_{i})}\approx\frac{u(x_{i})-u\left(CP_{\mathcal{S}_{j}}(x_{i})\right)}{d_{i}\cdot\hat{q}_{i}}$ where $CP_{\mathcal{S}_{j}}$ denotes the modified closest point function identifying points in $\Sigma_{j}^{BC}$ with points in $\Lambda_{j}$. Combining this with (4), and applying the identity extension for the Dirichlet component, $u(x_{i})=u(CP_{\mathcal{S}_{j}}(x_{i}))$, we find that ${\rm\bf T}_{j}$ must enforce the extension $u(x_{i})=\frac{u\left(CP_{\mathcal{S}_{j}}(x_{i})\right)}{1+\alpha d_{i}\cdot\hat{q}_{i}},$ with $u\left(CP_{\mathcal{S}_{j}}(x_{i})\right)$ replaced by the same interpolation used in the global scheme discussed in Section 2.

As $d_{i}$ approaches the surface normal direction, $d_{i}\cdot\hat{q}_{i}$ will tend to zero. In this event, the extension reduces to $u(x_{i})=u\left(CP_{\mathcal{S}_{j}}(x_{i})\right)$, which is just the standard extension corresponding to the interior. Fortuitously, this case arises when the point $x_{i}$ lies adjacent to the interior points where this condition would be applied anyway, and in our experience this ensures that the method remains robust.

6 Results

The solvers described in the previous sections were implemented in C**++**, with PETSc petsc-user-ref ; petsc-efficient providing the linear algebra data structures and MPI parallelization, and Umfpack Davis:Umfpack providing the local solutions. Here we focus on evaluating the solver, though in practice one should accelerate the solver with a Krylov method. The (O)RAS solver was placed into a PETSc PCSHELL preconditioner, allowing it to be embedded in any of their Krylov methods, and we have found coupling with GMRES to be a favorable pair.

Equation (1) was solved over the Stanford Bunny bunny , which has been scaled to be two units tall. The original triangulation has not been modified in any way beyond this scaling. This surface has several holes and is complicated enough to stress the solvers, making it a good test case. Our chosen grid spacing was $\Delta x=1/120$, which paired with tri-quadratic interpolation gives $N_{A}=947,964$ active nodes in the global problem. The origin was placed at the center of the bounding box containing the bunny and the right hand side $f=\phi(\pi-\phi)\sin(3\phi)(\sin\theta+\cos(10\theta))/2$ was used after extending it to be constant along the surface normals.

Table 1 shows the effects of subdomain count $N_{S}$, overlap width $N_{O}$, and Robin parameter $\alpha$. For comparison, GMRES preconditioned with the standard block-Jacobi method with $64$, $96$, and $128$ blocks requires more than $10000$ iterations. The solvers display the expected behavior with the iteration count increasing for larger subdomain counts and decreasing with larger overlap widths. ORAS consistently requires fewer iterations than RAS, though the final sub-table shows the dependence of this performance on the appropriate choice of Robin weight $\alpha$. The partitioning, the initial error, and the error in the approximate solution after $10$ and $550$ iterations are visible in Fig. 2 for one run of the solver.

Choosing an optimal value for $\alpha$ is non-trivial as it depends on the value of $c$, the mesh width, and the geometry. Additionally, the presence of cross points in decomposition, where more than two subdomains meet, complicate the matter. From the planar case, it is known that $\alpha\sim\mathcal{O}\left(\Delta x^{-1/2}\right)$, but determining precise values a priori is limited to simple splittings Gand:OptPar ; Gand:XPT . An upcoming work from the same authors explores this in much greater detail.

7 Conclusion

Restricted additive Schwarz and optimized restricted additive Schwarz solvers were formulated for the closest point method applied to (1). These solvers provide a solution mechanism for larger problem sizes and will allow users of the CPM to leverage large scale parallelism. Table 1 shows the dramatic reduction in iteration count when Robin TCs are used. These solvers were more completely evaluated in May:Thesis , which includes an exploration of their utility as preconditioners. The optimized conditions come at the cost of some additional complexity in the implementation, and even the standard RAS solver brings parallel capabilities to the user. Interesting extensions to this work include multiplicative methods, non-overlapping Robin schemes, two-level solvers, and inclusion of advective terms in the model equation.

Acknowledgements The authors gratefully acknowledge the financial support of NSERC Canada (RGPIN 2016-04361 and RGPIN 2018-04881), and the preliminary work of Nathan King that inspired this project.