Volume Integrals in Complex Geometry: A Case Study of Poisson's Equation Travis Askham, Department of Mathematics Courant Institute of Mathematical Sciences CBMS-NSF Conference: Fast Direct Solvers for Elliptic PDEs Abstract Notes on Steps 2 and 3 The solution of many problems in scientic computing requires the evaluation of a volume integral over some domain. For free-space problems, the use of an adaptive tree structure with tensor product grids dened on the leaves allows for a simple, lightweight discretization of the integral. When coupled with an FMM, such a discretization can be integrated rapidly. Volume integrals dened on domains with complex geometry cannot be immediately approached in this manner. The main issue is that the boundary of the domain will intersect the interior of many of the boxes in the discretization, much like the cut-cells encountered by nite volume schemes. We discuss global function extension as a solution to this problem and present a scheme for computing such an extension. We present the method in the context of Poisson's equation. Our focus is step 1, so we will briey discuss steps 2 and 3 rst. Applications With a Volume Component The solution to the following problems can be obtained, in part, by computing a volume integral: • Inhomogeneous problems, e.g. Poisson's equation. • Time-stepping for linear and nonlinear problems, e.g. IMEX schemes for advection diusion problems. • Variable coecient elliptic problems, e.g. of the form ∇ · ((x)∇φ(x)) = f (x). The solution of step 2 is perhaps the best understood of the three steps and we will say little about it. Because we consider multiply connected domains, we choose a combination of single and double layer potentials to represent the solution uH . In our implementation, we employ Alpert quadrature (Alpert, 1999) to approximate the integrals and use the fast-direct solver package HODLR (Ambikasaran and Darve, 2013) to solve the system. Step 3, setting u = v + uH , may appear trivial. What is implied here, however, is the evaluation of the layer H potential describing u at the discretization points of the domain Ω. Given a quadrature scheme, this step can be done quickly with the fast multipole method. For points near to the boundary, designing the quadrature is dicult because the integral kernels for the layer potentials are nearly singular. There are a few good options for this in 2D and for our implementation we use the quadrature by expansion (QBX) scheme of (Klöckner, Barnett, Greengard, O'Neil, 2013). Volume Integral for Poisson's Equation It is well known that v as in step 1 can be calculated as Z v(x) = G(x, y)˜ g (y) dy R2 Prototypical Approach to Poisson's Equation For the sake of context, we outline the typical integral equations-based approach to computing the solution of Poisson's equation. To solve ∆u = g(x) in Ω u = h(x) on ∂Ω, (1) (2) where g˜ = g in Ω and G is the free-space Green's function for Poisson's equation. To discretize the integral, there are two basic approaches. One could use a bodytted mesh (e.g. a nite element style mesh) or one could use an embedded boundary mesh, in which the volume mesh is aligned with the coordinate axes without regard for the boundary. We opt for the latter approach, choosing an adaptive quad-tree based mesh with tensor product grids, for the following reasons: • There are specialized FMMs (Ethridge and Green- follow these steps: 1. Compute v (typically via a volume integral) such that ∆v = g in Ω. 2. Solve an appropriate boundary integral equation for uH such that ∆uH = 0 in Ω and uH = h − v on ∂Ω. 3. Evaluate u = u + v . H gard, 2001) which are competitive with FFTs in terms of work per grid point for such a mesh. • Because the local interactions are always of the same form (up to scale), precomputed tables are used and it is simple to accurately compute the gradient and Hessian of v . • The creation of the mesh itself is cheap. · June 28, 2014 · Dartmouth College · Hanover, New Hampshire The Cut-Cell Problem Solver Demonstration and Discussion Because the boundary is embedded in the mesh, there are cells (leaf boxes of the quad-tree) of the discretization which are intersected by the boundary and have grid points both inside and outside the domain. For the exterior points, we must invent values for the right hand side g . If we simply set g to zero for exterior points, then it is unlikely that within any box the right hand side will be resolved. Further, discontinuities across boxes cause the volume integral v to be nonsmooth (at most C 1 ) in this region. Therefore, setting g to zero for exterior points requires heavy renement near the boundary. Further, local extension of the function g based on nearby points does not alleviate the issue of discontinuities across boxes. In light of this, we seek a global extension of g which can be computed rapidly. To demonstrate the ecacy of global function extension, we implemented a solver for Poisson's equation using the tools described above. For steps 2 and 3, 10,000 boundary points were used with 16th order Alpert quadrature. The QBX and HODLR subroutines were run with a tolerance of 1.0E-14. The volume integral was computed with a modied version of (Ethridge and Greengard), with 4th order grids. The experiments were run on an Intel Core i5 1.7GHz processor on a laptop with 4Gb of RAM. The right hand side g is as in the previous section. The adaptive discretization hasn't been perfected for complex geometries yet, so we report accuracy for both the unextended right hand side (setting g = 0 in the exterior) and the continuously extended right hand side (via a Laplace solve in the exterior) on the same grid. Proposed Global Function Extension Method Take Ω2 to be a larger set containing Ω and g a function ¯ . Let w solve the following Laplace problem dened on Ω ∆w = 0 in Ω2 \ Ω w = g on ∂Ω, (3) (4) setting either w ≡ 0 on ∂Ω2 or choosing growth conditions at innity. Then g˜ = g on Ω and w on Ω2 is a continuous extension of g to Ω2 . To compute w, we simply use the same technology as steps 2 and 3. Above: The original function g is shown on the left. We computed the extension g˜ on the right for the exterior Dirichlet boundary value problem (Ω2 = R2 ) with appropriate growth conditions. It is clear that by looking at polyharmonic Dirichlet problems on Ω2 \ Ω one can nd smoother extensions of g . However, the normal derivatives of the data g would have to be supplied or computed. Above: plot of log10 (error) in the solution using the unextended g (left) and the extended g (right). Let Nv be the number of points for the volume integral and Nin be the number of points inside the domain; let Eext and E0 be the error with and without extension; and let TEG , TQ , and TQE be the time for the volume integral, QBX, and QBX for the extension, respectively. Nv Nin Eext E0 TEG TQ TQE 65,152 31,517 1.02E-05 1.21E-03 0.232 6.924 6.26 359,392 180.989 1.57E-08 3.27E-06 1.08 11.01 9.06 3,056,512 1,529,310 2.84E-11 4.45E-08 8.93 13.9 11.2 Based on these early experiments, we believe that the framework presented here, with an embedded boundary mesh and global function extension, will allow for fast, accurate computations in complex geometries. Future work will revolve around applying these tools to an adaptive Navier Stokes scheme in complex geometry in 2D. Acknowledgments This work was supported by the Applied Mathematical Sciences program of the U.S. Department of Energy under Contract DEFG0288ER25053 and by the Oce of the Assistant Secretary of Defense for Research and Engineering and AFOSR under NSSEF Program Award FA9550-10-1-0180.