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) Ruuth and Merriman (2008); Macdonald and Ruuth (2010). We consider the positive Helmholtz equation
where is a constant and is the Laplace-Beltrami operator associated with the surface . The evolution of diffusion equations by implicit time-stepping schemes and Laplace-Beltrami eigenvalue problems Macdonald et al. (2011) 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 Macdonald et al. (2013), and speed up shape classification Reuter et al. (2006), 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 Floater and Hormann (2005). 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 Dziuk and Elliott (2007). This approach leads to a sparse and symmetric system but is sensitive to the quality of the triangulation. Level set methods for surface PDEs Bertalmio et al. (2001) 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 and Macdonald (2015). 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 Gander (2006); St-Cyr et al. (2007), 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 Balay et al. (2018, 1997) 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 (2018).
2 The closest point method
The CPM was introduced in Ruuth and Merriman (2008) 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, for , 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 consisting of all points within a distance of the surface, where is an upper bound on the principal curvatures of the surface Chu and Tsai (2018).
From this mapping, an extension operator can be defined that sends functions defined on the surface, , to functions defined on the embedding space via composition with the closet point mapping, . 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 Ruuth and Merriman (2008).
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 interpolationBerrut and Trefethen (2004) 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 interpolation and a grid spacing of requires that the tube contain all points within a distance of from the surface, thus limiting the acceptable grid spacings in relation to .
The grid points within the computational tube form the set of active nodes, . For (1), we need only discretize the regular Laplacian on . Here we consider the second order accurate centered difference approximation requiring points. Around and lying outside the tube, a set of ghost nodes, , is formed from any incomplete differencing stencils. With a total of active nodes and ghost nodes, we define the discrete Laplacian and extension operators, and , where applies the centered difference Laplacian over all active nodes, and is the discretization of . 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 , which was used successfully for parabolic equations with explicit time-stepping in Ruuth and Merriman (2008). However, for implicit time-stepping Macdonald and Ruuth (2010) and eigenvalue problems Macdonald et al. (2011) a modified form is needed. In Macdonald and Ruuth (2010) it was recognized that there was a redundant interpolation being performed, and that its removal could stabilize the discretization. The stabilized form 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 St-Cyr et al. (2007)). 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 is decomposed into disjoint subdomains, , for . These disjoint subdomains are then grown to form overlapping subdomains , whose boundaries are labelled depending on where they lie in the disjoint partitioning. Taking gives and allows the definition of the local problems
where are generally linear boundary operators defining the TCs. RAS is achieved by choosing as identity operators, corresponding to Dirichlet TCs, while ORAS uses Robin TCs, , where is the outward pointing boundary normal on and is a constant weight on the Dirichlet contribution.
The subproblems in equation (2) are initialized with a guess for the global solution (defined at least over the boundaries ), which is usually just taken as . After all of the subproblems have been solved a new global solution is constructed with respect to the disjoint partitioning, , 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 . (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 , overlapping partitions by , and the boundaries of the overlapping partitions by .
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 Karypis and Kumar (1998) is a graph partitioner that is frequently used within the DD community to partition meshes Dolean et al. (2015). The stencils of and 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 . Fig. 1 shows a portion of one such disjoint partition, in black circles, for a circular surface.
With obtained from METIS, overlapping subdomains can be formed. This construction proceeds in the following steps:
All nodes in are added to .
layers of overlap nodes are added around . Layers are added one at a time from globally active nodes neighboring .
A subset of the ghost nodes, , are placed in which consists of nodes that neighbor a member of .
The shapes of the disjoint and overlapping subdomains are not known in advance. The boundary is approximated discretely by the closest points of the final layer of overlap nodes, and held in the set .
Nodes needed to complete stencils from the ambient Laplacian or extension operator, including extension from the points , are placed in the set .
For ORAS a layer of ghost nodes around are also placed in .
The active nodes in the subdomain consist of and the active portion of . 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 at one of its boundaries.
The Robin TCs, to be defined in Section 5, need some final information about the subdomain. Every node in is identified with the point in 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 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
makes this construction difficult. Instead we define the conormal vectors from the point of view of the boundary nodes. Takeas a node whose associated conormal direction, , is sought. Let be its closest point in , and be the unit surface normal there. Connecting the boundary location to the boundary node via , we obtain a usable approximation to the conormal by computing the component of that is orthogonal to and normalizing, i.e., In the unlikely event that lies perfectly in the surface normal direction, we set 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 beyond the surface boundary Macdonald et al. (2011). As such, the local operators will take the form
where is the extension operator for the nodes in as inherited from the global operator and is the modified extension operator for the nodes in . When solving for the local correction to the solution the right hand side of the local problem, , will be the restriction of the residual to . The final rows of the right hand side, those lying over , become zeros corresponding to the homogenous TCs.
Homogeneous Dirichlet TCs can be enforced to first order accuracy by extending zeros over all of . With the right hand side already set to zero there, the modified extension reduces to the identity mapping,
We discretize the Robin condition
using a forward difference in the direction for each node in and the first order accurate Dirichlet condition from above. Taking the partial derivative , and applying the change of variables , allows one to write the Neumann term in equation (4) in terms of the displacement vector from Section 4
. Assuming for the moment thatand are not perpendicular, the derivative in the conormal direction can be approximated by where denotes the modified closest point function identifying points in with points in . Combining this with (4), and applying the identity extension for the Dirichlet component, , we find that must enforce the extension with replaced by the same interpolation used in the global scheme discussed in Section 2.
As approaches the surface normal direction, will tend to zero. In this event, the extension reduces to , which is just the standard extension corresponding to the interior. Fortuitously, this case arises when the point lies adjacent to the interior points where this condition would be applied anyway, and in our experience this ensures that the method remains robust.
The solvers described in the previous sections were implemented in C++, with PETSc Balay et al. (2018, 1997) providing the linear algebra data structures and MPI parallelization, and Umfpack Davis (2004) 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 Turk and Levoy (1994), 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 , which paired with tri-quadratic interpolation gives 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 was used after extending it to be constant along the surface normals.
Table 1 shows the effects of subdomain count , overlap width , and Robin parameter . For comparison, GMRES preconditioned with the standard block-Jacobi method with , , and blocks requires more than 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 . The partitioning, the initial error, and the error in the approximate solution after and iterations are visible in Fig. 2 for one run of the solver.
Choosing an optimal value for is non-trivial as it depends on the value of , 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 , but determining precise values a priori is limited to simple splittings Gander (2006); Gander and Kwok (2012). An upcoming work from the same authors explores this in much greater detail.
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 (2018), 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.
- PETSc users manual. Technical report Technical Report ANL-95/11 - Revision 3.9, Argonne National Laboratory. Cited by: §1, §6.
- Efficient management of parallelism in object oriented numerical software libraries. In Modern Software Tools in Scientific Computing, E. Arge, A. M. Bruaset, and H. P. Langtangen (Eds.), pp. 163–202. Cited by: §1, §6.
- Barycentric Lagrange interpolation. SIAM Review 46 (3), pp. 501–517 (eng). External Links: Cited by: §2.
Variational problems and partial differential equations on implicit surfaces. Journal of Computational Physics 174 (2), pp. 759–780 (eng). External Links: Cited by: §1.
- The closest point method and multigrid solvers for elliptic equations on surfaces. SIAM Journal on Scientific Computing 37 (1), pp. A134–A155. Cited by: §1.
- Volumetric variational principles for a class of partial differential equations defined on surfaces and curves. Research in the Mathematical Sciences 5 (2), pp. 19. External Links: Cited by: §2.
- Algorithm 832: UMFPACK v4.3—an unsymmetric-pattern multifrontal method. ACM Transactions on Mathematical Software (TOMS) 30 (2), pp. 196–199. Cited by: §6.
- An introduction to domain decomposition methods: algorithms, theory, and parallel implementation. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. External Links: Cited by: §4.
- SURFACE finite elements for parabolic equations. Journal of Computational Mathematics 25 (4), pp. 385–407 (eng). External Links: Cited by: §1.
- Surface parameterization: a tutorial and survey. Mathematics and Visualization Advances in Multiresolution for Geometric Modelling, pp. 157–186. Cited by: §1.
- Best Robin parameters for optimized Schwarz methods at cross points. SIAM Journal on Scientific Computing 34 (4), pp. 1849–1879 (eng). External Links: Cited by: §6.
- Optimized Schwarz methods. SIAM Journal on Numerical Analysis 44 (2), pp. 699–731 (eng). External Links: Cited by: §1, §6.
- A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM Journal on Scientific Computing 20 (1), pp. 359–392 (eng). External Links: Cited by: §4.
- Solving eigenvalue problems on curved surfaces using the closest point method. Journal of Computational Physics 230 (22), pp. 7944–7956 (eng). External Links: Cited by: §1, §2, §5.
- Simple computation of reaction diffusion processes on point clouds. Proceedings of the National Academy of Sciences 110 (23) (eng). External Links: Cited by: §1.
- The implicit closest point method for the numerical solution of partial differential equations on surfaces. SIAM Journal on Scientific Computing 31 (6), pp. 4330–4350 (eng). External Links: Cited by: §1, §2.
- Domain decomposition solvers and preconditioners for the implicit closest point method. Master’s Thesis, Simon Fraser University. Cited by: §1, §7.
- Laplace–Beltrami spectra as ‘Shape-DNA’ of surfaces and solids. Computer-Aided Design 38 (4), pp. 342–366 (eng). External Links: Cited by: §1.
- A simple embedding method for solving partial differential equations on surfaces. Journal of Computational Physics 227 (3), pp. 1943–1961 (eng). External Links: Cited by: §1, §2, §2, §2.
- Optimized multiplicative, additive, and restricted additive Schwarz preconditioning. SIAM Journal on Scientific Computing 29 (6), pp. 2402–2425 (eng). External Links: Cited by: §1, §3.
- Zippered polygon meshes from range images. In Proceedings of the 21st annual conference on computer graphics and interactive techniques, SIGGRAPH ’94, pp. 311–318 (eng). External Links: Cited by: §6.