An alternative extended linear system for boundary value problems on locally perturbed geometries
AAn alternative extended linear system for boundary value problems onlocally perturbed geometries
Y. Zhang and A. Gillman
Abstract
This manuscript presents a new extended linear system for integral equation based techniques forsolving boundary value problems on locally perturbed geometries. The new extended system is similarto the one presented by Greengard et al. (2009) [4] for which a fast direct solver was constructed byZhang and Gillman (2018) [9] but the corresponding fast direct solver is more efficient than the onefor the original extended system. Another important feature of the new extended system is that itcan easily handle problems that require specialized quadrature for weakly singular kernels. Numericalresults illustrate the improved performance of the fast direct solver for the new extended system whencompared to the fast direct solver for the original extended system.
This manuscript presents an integral equation based solution technique for elliptic boundary value prob-lems on locally-perturbed geometries. Such problems arise in applications such as optimal shape design.In each iteration or optimization cycle the changes to the object shape often stay local to certain partsof the object. The proposed approach formulates an extended linear system that allows for the boundaryvalue problem on the new geometry to be expressed in terms of a linear system on the original geometryplus a correction to account for the local perturbation. This idea was first proposed in [4] and a fastdirect solver was constructed for the resulting formulation in [9]. Unfortunately, the fast direct solverfor the original extended system required inverting a matrix the size of the number of points removedfrom the original geometry which is expensive if the removed portion is large. Another difficulty of theoriginal extended system is that care is required when the technique is applied to systems discretizedusing quadrature for weakly singular kernels. The extended linear system proposed in this manuscriptovercomes these two difficulties. Additionally, a fast direct solver for the new extended system can beconstructed from the tools presented in [9] but is more efficient than the original fast direct solver.This manuscript briefly reviews a boundary integral formulation for a Laplace problem with Dirichletboundary data and the linear system that results from the discretization in section 2. Next, the originalextended system and the new extended system are presented in section 3. Finally numerical resultsillustrate the efficiency of the fast direct solver for the new extended system in section 4.
Consider the interior Laplace problem with Dirichlet boundary condition − ∆ u ( x ) = 0 for x ∈ Ω ,u ( x ) = g ( x ) for x ∈ Γ (1)where Ω denotes the interior of the geometry, and Γ denotes the boundary of Ω, as illustrated in Figure1(a). Let G ( x , y ) = − π log | x − y | denote the Green’s function for the Laplace operator and D ( x , y ) =1 a r X i v : . [ m a t h . NA ] J u l n ( y ) G ( x , y ) denote the double layer kernel where n ( x ) denotes the outward facing normal vector at thepoint x ∈ Γ. The solution to (1) can be expressed as u ( x ) = (cid:90) Γ D ( x , y ) σ ( y ) ds ( y ) for x ∈ Ω , (2)where σ ( x ) is some unknown density defined only on the boundary Γ. Enforcing that u ( x ) satisfies theboundary condition results in the following integral equation for σ ( x ); − σ ( x ) + (cid:90) Γ D ( x , y ) σ ( y ) ds ( y ) = g ( x ) . (3)Upon discretization via a Nystr¨om or boundary element method, one is left with solving a dense linearsystem Aσ = g , (4)where A is the discretized boundary integral operator and σ is the vector approximating σ at the dis-cretization points. ΩΓ x n ( x )(a) Γ k Γ p Γ c (b) Figure 1: (a) A sample geometry Ω with boundary Γ and outward facingnormal vector n ( x ) at the point x ∈ Γ. (b) A sample locally perturbedgeometry where the original boundary is Γ o = Γ k ∪ Γ c , the portion of theboundary being removed is Γ c , the portion of the original boundary remain-ing is Γ k and the newly added boundary is Γ p . Consider a boundary value problem on a geometry with a local perturbation as illustrated in Figure 1(b).Let Γ o denote the boundary of the original geometry, Γ k denote the portion of the boundary that is notchanging and Γ c denote the portion that is cut or removed. So Γ o = Γ c ∪ Γ k . Let Γ p denote the newportion of the boundary. Then the new geometry has a boundary Γ n defined by Γ n = Γ k ∪ Γ p .The discretized linear systems can be partitioned according to this notation. In other words, the originalsystem can be expressed as A oo σ o = (cid:20) A kk A kc A ck A cc (cid:21) (cid:18) σ k σ c (cid:19) = (cid:18) g k g c (cid:19) = g o , (5)and the linear system for the perturbed geometry can be expressed as A nn σ n = (cid:20) A kk A kp A pk A pp (cid:21) (cid:18) σ k σ p (cid:19) = (cid:18) g k g p (cid:19) = g n (6)where σ k denotes the vector whose entries are the approximate solution at the discretization points onΓ k , σ c denotes the vector whose entries are the approximate solution at the discretization points on Γ c ,etc. Likewise A kk is the submatrix of the discretized integral equation corresponding to the interaction ofΓ k with itself, A kc is the submatrix of the discretized integral equation corresponding to the interactionof Γ k with Γ c , etc. 2 .1 Original extended linear system The discretized problem on Γ n can be expressed as an extended linear system [4] by (cid:20) A oo A pp (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ˜ A + (cid:18) − A kc − B cc (cid:19) A op A pk (cid:124) (cid:123)(cid:122) (cid:125) Q orig σ k σ c σ p (cid:124) (cid:123)(cid:122) (cid:125) σ ext = g k g p (cid:124) (cid:123)(cid:122) (cid:125) g ext (7)where A kc denotes the submatrix of A oo corresponding to the interaction between Γ k and Γ c , A op denotesthe discretization of the double layer integral operator on Γ p evaluated on Γ o , A pk denotes the discretiza-tion of the double layer integral operator on Γ k evaluated on Γ p , and B cc denotes the sub-matrix of A oo corresponding to the interaction of Γ c with itself but the diagonal entries are set to zero. The matrix Q orig is called the update matrix . The extended system (7) is obtained by subtracting the contributionsfrom Γ c in A oo and adding the contributions from Γ p . Upon solving (7), only σ k and σ p are used toevaluate the solution inside of Γ n . Effectively σ c is a dummy vector. Details of the derivation of (7) areprovided in [4, 9]. We exploit the fact that the contribution from Γ c is not used to find the solution inside of Γ n in thenew extended linear system. Specifically, we introduce the vector σ dum c fully knowing a priori that it willcontain useless information. Then solving (6) is equivalent to solving the following A kk A kp A ck A cc A pk A pp σ k σ dumc σ p = g k g p . (8)The expanded form of (8) is (cid:20) A oo A pp (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ˜ A + − A kc A kp A pk (cid:124) (cid:123)(cid:122) (cid:125) Q new σ k σ dumc σ p (cid:124) (cid:123)(cid:122) (cid:125) σ ext = g k g p (cid:124) (cid:123)(cid:122) (cid:125) g ext . (9)Here Q new is the new update matrix. Notice that Q new has a zero row. When constructing the fast direct solver for the locally perturbed boundary value problem, there areadvantages to writing the system in the form of (7) and (9). Since the matrix A oo is the system resultingfrom the discretization of the integral equation on the original geometry, we assume that a fast directsolver was already computed for it. Any fast direct solver such as Hierarchically Block Separable (HBS) [3],
Hierarchically Semi-Separable (HSS) [8, 1],
Hierarchical Interpolative Factorization (HIF) [6] and H and H - matrix methods [5] can be used. Additionally, the update matrices Q orig and Q new are low rank. Thisallows for the inverse of the extended systems to be applied rapidly via a Sherman-Morrison-Woodburyformula σ ext = (cid:16) ˜ A + Q (cid:17) − g ext ≈ (cid:16) ˜ A + LR (cid:17) − g ext ≈ ˜ A − g ext − ˜ A − L (cid:16) I + R ˜ A − L (cid:17) − R ˜ A − g ext , (10)3here I is an identity matrix, and LR denotes the low rank factorization of the update matrix Q .The low rank property of the update matrices Q orig and Q new can be observed by noting that the matrices A kc , A kp , A pk and A op are low rank. The only full rank matrix in the update matrices is B cc . Since Q new does not contain that matrix, it’s rank is smaller than Q orig and the Sherman-Morrison-Woodburyformula can be applied more rapidly. The details for efficiently creating the low rank factorizations canbe found in [9]. This section illustrates the performance of the fast direct solver for the proposed extended linear systemfor a collection of problems. The integral equations are discretized via the Nystr¨om method with a 16 th order composite Gaussian quadrature. For all problems, the original geometry is discretized with enoughpoints in order for the boundary value problem to be solved to 10 digits of accuracy. The HBS directsolver [3] was used in the examples in this section. For all tests, the tolerance for HBS compression andlow-rank approximation is set to (cid:15) = 10 − .Roughly speaking the cost of building fast direct solvers is split into two parts: precomputation andsolve. The time for precomputation is the time for constructing all the parts of the fast direct solver.For the extended systems, this includes constructing the low rank factorizations of the update matrices Q orig or Q new and inverting the small matrix in the Sherman-Morrison-Woodbury formula (10). For theHBS solver, the precomputation includes creating a compressed approximation of the discretized systemon the new geometry and inverting that approximation. The time for the solve is the time for applyingthe resulting solver to one vector (or right-hand-side). For the extended systems, this is the time forapplying the Sherman-Morrison-Woodbury formula (10). For the HBS solver, it is the time for applyingthe approximate inverse.To illustrate the efficiency of the proposed technique, we compare the performance of the new solutiontechnique with the fast solver developed for the original extended system and building an HBS solverfrom scratch for the new geometry. We report the following: • N o : the number of discretization points on the original geometry; • N c : the number of discretization points cut from the original geometry; • N p : the number of discretization points added; • T new ,p : the time in seconds for the precomputation of the proposed solver; • T orig ,p : the time in seconds for the precomputation of the original fast solver; • T hbs ,p : the time in seconds for the precomputation of HBS from scratch for the new geometry; • r p = T hbs ,p T new ,p ; • T new ,s : the time in seconds for applying the proposed solver to one right-hand-side; • T orig ,s : the time in seconds for applying the original solver to one right-hand-side; • T hbs ,s : the time in seconds for applying the HBS inverse to one right-hand-side; • r s = T hbs ,s T new ,s .The ratios r p and r s are measures for the speed-up (or slow-down) by using the proposed solver versusbuilding a new fast direct solver from scratch for the new geometry. If r p is greater than 1, the precom-putation of the proposed solver is faster than building a fast direct solver from scratch. If r p is less than1, the precomputation of the proposed solver is slower than building a fast direct solver from scratch, etc.4ll experiments were run on a dual 2.3 GHz Intel Xeon Processor E5-2695 v3 desktop workstation with256 GB of RAM. The code is implemented in MATLAB, apart from the randomized linear algebra utilizedin creating low rank factorization rapidly which is implemented in Fortran. Consider the Laplace boundary value problem (1) on the geometry illustrated in Figure 2. The cornersare smoothed via the scheme in [2]. A detailed description of this geometry is given in [9].In the first experiment, the number of points cut remains fixed, N c = 16, while the number of discretizationpoints on Γ k grows. In Figure 2, this corresponds to the nose height d decreasing as N k grows. The resultsare reported in Table 1. All three solution techniques are linear with respect to N o and the precomputationtime for the new solution technique is about the same as the original extended system solver. It is roughly3.5 times faster than building a new direct solver from scratch for the new geometry. The cost of applyingthe proposed solver is almost as fast applying the HBS approximate inverse.In the next example, N c grows by the same factor as N k . The nose height d in Figure 2 remains fixed.Table 2 reports on the performance of all three solvers for this geometry. The proposed solution techniqueis the fastest for the precomputation step. It is much faster than the solver based on the original extendedsystem formulation, especially for the case where N c is large. A factor of roughly 2.9 speed up in theprecomputation is observed. Again applying the proposed solver is slightly slower than applying the HBSapproximate inverse. ΩΓ k Γ c Γ c Γ p Γ p d (cid:101) (cid:82) (cid:38)(cid:37)(cid:39)(cid:36) Figure 2:
The square with nose geometry. A nose of height d is smoothlyattached to the a square. N o T orig ,p T new ,p T hbs ,p r p T orig ,s T new ,s T hbs ,s r s Table 1:
Times for applying the solution technique to (1) on the square withthinning nose geometry.
Next we consider applying the proposed solution technique to the Laplace boundary value problem (1)where the local perturbation is a refinement in a portion of the geometry. Figure 3(a) illustrates the5 o N c T orig ,p T new ,p T hbs ,p r p T orig ,s T new ,s T hbs ,s r s Table 2:
Times for applying the solution techniques to (1) on the squarewith fixed nose geometry. geometry under consideration. It is given by the following parameterization: x ( t ) = (cid:18) r ( t ) cos( t ) r ( t ) sin( t ) (cid:19) , with r ( t ) = 1 + 0 . t ) for t ∈ [0 , π ] . The portion of the boundary being refined is highlighted in red. Figure 3(b) is a zoomed in illustrationof that region. Figure 3(c) illustrates the local refinement. Three Gaussian panels ( N c = 48) are replacedwith N p discretization points ( N p /
16 Gaussian panels). The number of discretization points on Γ k remainsfixed; N k = 6352.Table 3 reports on the performance of all three solution techniques for this problem. The proposedsolution technique is 13 to 21 times faster than building a new solver from scratch while applying thesolver is less than a factor two slower than applying the HBS approximate inverse.Γ k (a) (b) (c) Γ c Γ p Figure 3: (a) The sunflower geometry with the portion of the boundary tobe refined in red. (b) The three Gaussian panels in the boxed region from theoriginal discretization. (c) Six Gaussian panels replacing the original threepanels. N p N p N o T orig ,p T new ,p T hbs ,p r p T orig ,s T new ,s T hbs ,s r s
96 0.015 6.06e-01 5.03e-01 7.52e+00 14.9 1.10e-02 1.30e-02 1.32e-02 1.02192 0.03 6.16e-01 3.62e-01 7.77e+00 21.4 1.17e-02 1.25e-02 9.30e-03 0.74384 0.06 6.83e-01 3.90e-01 7.72e+00 19.8 1.36e-02 1.42e-02 9.13e-03 0.64768 0.12 7.60e-01 4.11e-01 7.78e+00 18.9 2.01e-02 1.20e-02 9.06e-03 0.761536 0.24 1.01e+00 6.09e-01 8.03e+00 13.2 4.72e-02 1.66e-02 1.00e-02 0.60
Table 3:
Times for applying the solution techniques to (1) on the geometryin Figure 3 with local refinement. .3 A Helmholtz problem with a locally refined discretization Besides being faster than the solver for the original extended system, the proposed solver has the advantagethat it can easily handle problems that are using specialized quadrature for weakly singular kernels. Theissue that arises for the original extended system is that it would be cumbersome to evaluate the entriesof the matrix A op corresponding to the interaction of Γ c with Γ p . This matrix does not arise in the newextended system.To illustrate the efficiency of the solver for systems that involve specialized quadrature we consider thefollowing exterior Dirichlet Helmholtz boundary value problem − ∆ u ( x ) + ω u = 0 for x ∈ Ω c ,u ( x ) = g ( x ) for x ∈ Γ (11)with Sommerfeld radiation condition on the sunflower geometry illustrated in Figure 3 where ω denotesthe wave number. We chose to represent the solution with the following combined field u ( x ) = (cid:90) Γ D ω ( x , y ) σ ( y ) ds ( y ) − iω (cid:90) Γ S ω ( x , y ) σ ( y ) ds ( y ) , (12)where D ω and S ω denote the double and single layer Helmholtz kernel and σ ( x ) is the unknown boundarycharge distribution.The integral equation that results from enforcing the Dirichlet boundary condition is12 σ ( x ) + (cid:90) Γ D ω ( x , y ) σ ( y ) ds ( y ) − iω (cid:90) Γ S ω ( x , y ) σ ( y ) ds ( y ) = g ( x ) . (13)We discretize the operator via Nystr¨om with a composite generalized Gaussian quadrature [7]. The wavenumber is set to ω = 20 which corresponds to the geometry being approximately 8.3 wavelengths in size.Again, we consider the local refinement problem. Table 4 reports on the performance of the proposedsolution technique and building a fast direct solver from scratch. For this problem, the proposed solveris anywhere from 15 to 35 times faster than building the fast direct solver from scratch. This speed up isthe result of the increased ranks associated with Helmholtz problems. Applying the proposed solver to aright-hand-side is roughly 1.5 times slower than applying the HBS solver. N p N p N o T new ,p T hbs ,p r p T new ,s T hbs ,s r s
96 0.015 1.13e+00 3.97e+01 35.2 4.06e-02 2.86e-02 0.71192 0.03 1.36e+00 4.08e+01 29.9 4.64e-02 2.93e-02 0.63384 0.06 1.44e+00 4.08e+01 28.4 3.91e-02 2.54e-02 0.65768 0.12 1.64e+00 4.17e+01 25.4 3.69e-02 2.70e-02 0.731536 0.24 2.64e+00 4.08e+01 15.4 4.39e-02 3.33e-02 0.76
Table 4:
Times for applying the solution techniques to (11) on the geometryin Figure 3 with local refinement.
This manuscript presented a new extended linear system for integral equation based solution techniquesfor boundary value problems on locally perturbed geometries. A fast direct solver for the new extendedsystem is significantly faster than a fast direct solver for the original system and can be easily applied todiscretizations utilizing quadrature for weakly singular kernels.7 eferences [1] S. Chandrasekaran and M. Gu A divide-and-conquer algorithm for the eigendecomposition of sym-metric block-diagonal plus semiseparable matrices.
Numerische Mathematik , 96(4):723-731, 2004[2] C. Epstein and M. O’Neil. Smoothed Corners and Scattered Waves.
SIAM Journal on ScientificComputing , 38, 2015[3] A. Gillman, P. Young and P.G. Martinsson A direct solver O ( N ) complexity for integral equationson one-dimensional domains. Frontiers of Mathematics in China , 7:217-247, 2012[4] L. Greengard, D. Gueyffier, P.G. Martinsson and V. Rokhlin Fast direct solvers for integral equationsin complex three-dimensional domains.
Acta Numerica , 18:243-275, 2009[5] W. Hackbusch A Sparse Matrix Arithmetic Based on H-Matrices; Part I: Introduction to H-Matrices.
Computing , 62:89-108, 1999[6] K. Ho and L. Ying Hierarchical Interpolative Factorization for Elliptic Operators: Integral Equations.
Communications on Pure and Applied Mathematics , 69, 2013[7] P. Kolm and V. Rokhlin Numerical quadratures for singular and hypersingular integrals.
Computersand Mathematics with Applications , 41:327-352, 2001[8] Z. Sheng, P. Dewilde and S. Chandrasekaran Algorithms to solve hierarchically semi-separablesystems.
Operator Theory: Advances and Applications , 176: 255-294, 2007[9] Y. Zhang and A. Gillman. A fast direct solver for boundary value problems on locally perturbedgeometries.