Adrianna Gillman

Adrianna Gillman, Patrick Young, Per-Gunnar Martinsson

An algorithm for the direct inversion of the linear systems arising from Nystrom discretization of integral equations on one-dimensional domains is described. The method typically has O(N) complexity when applied to boundary integral equations (BIEs) in the plane with non-oscillatory kernels such as those associated with the Laplace and Stokes' equations. The scaling coefficient suppressed by the "big-O" notation depends logarithmically on the requested accuracy. The method can also be applied to BIEs with oscillatory kernels such as those associated with the Helmholtz and Maxwell equations; it is efficient at long and intermediate wave-lengths, but will eventually become prohibitively slow as the wave-length decreases. To achieve linear complexity, rank deficiencies in the off-diagonal blocks of the coefficient matrix are exploited. The technique is conceptually related to the H and H^2 matrix arithmetic of Hackbusch and co-workers, and is closely related to previous work on Hierarchically Semi-Separable matrices.

Adrianna Gillman, Alex Barnett

We consider the numerical solution of the scattering of time-harmonic plane waves from an infinite periodic array of reflection or transmission obstacles in a homogeneous background medium, in two dimensions. Boundary integral formulations are ideal since they reduce the problem to $N$ unknowns on the obstacle boundary. However, for complex geometries and/or higher frequencies the resulting dense linear system becomes large, ruling out dense direct methods, and often ill-conditioned (despite being 2nd-kind), rendering fast multipole-based iterative schemes also inefficient. We present an integral equation based solver with O(N) complexity, which handles such ill-conditioning, using recent advances in "fast" direct linear algebra to invert hierarchically the isolated obstacle matrix. This is combined with a recent periodizing scheme that is robust for all incident angles, including Wood's anomalies, based upon the free space Green's function kernel. The resulting solver is extremely efficient when multiple incident angles are needed, as occurs in many applications. Our numerical tests include a complicated obstacle several wavelengths in size, with $N=10^5$ and solution error of $10^{-10}$, where the solver is 66 times faster per incident angle than a fast multipole based iterative solution, and 600 times faster when incident angles are chosen to share Bloch phases.

Adrianna Gillman, Per-Gunnar Martinsson

The large sparse linear systems arising from the finite element or finite difference discretization of elliptic PDEs can be solved directly via, e.g., nested dissection or multifrontal methods. Such techniques reorder the nodes in the grid to reduce the asymptotic complexity of Gaussian elimination from $O(N^{2})$ to $O(N^{1.5})$ for typical problems in two dimensions. It has recently been demonstrated that the complexity can be further reduced to O(N) by exploiting structure in the dense matrices that arise in such computations (using, e.g., $\mathcal{H}$-matrix arithmetic). This paper demonstrates that such \textit{accelerated} nested dissection techniques become particularly effective for boundary value problems without body loads when the solution is sought for several different sets of boundary data, and the solution is required only near the boundary (as happens, e.g., in the computational modeling of scattering problems, or in engineering design of linearly elastic solids.

Tracy Babb, Adrianna Gillman, Sijia Hao et al.

This paper presents a numerical method for variable coefficient elliptic PDEs with mostly smooth solutions on two dimensional domains. The PDE is discretized via a multi-domain spectral collocation method of high local order (order 30 and higher have been tested and work well). Local mesh refinement results in highly accurate solutions even in the presence of local irregular behavior due to corner singularities, localized loads, etc. The system of linear equations attained upon discretization is solved using a direct (as opposed to iterative) solver with $O(N^{1.5})$ complexity for the factorization stage and $O(N \log N)$ complexity for the solve. The scheme is ideally suited for executing the elliptic solve required when parabolic problems are discretized via time-implicit techniques. In situations where the geometry remains unchanged between time-steps, very fast execution speeds are obtained since the solution operator for each implicit solve can be pre-computed.

Peter Geldermans, Adrianna Gillman

This manuscript presents an adaptive high order discretization technique for elliptic boundary value problems. The technique is applied to an updated version of the Hierarchical Poincaré-Steklov (HPS) method. Roughly speaking, the HPS method is based on local pseudospectral discretizations glued together with Poincaré-Steklov operators. The new version uses a modified tensor product basis which is more efficient and stable than previous versions. The adaptive technique exploits the tensor product nature of the basis functions to create a criterion for determining which parts of the domain require additional refinement. The resulting discretization achieves the user prescribed accuracy and comes with an efficient direct solver. The direct solver increases the range of applicability to time dependent problems where the cost of solving elliptic problems previously limited the use of implicit time stepping schemes.

Adrianna Gillman

This paper presents an integral formulation for Helmholtz problems with mixed boundary conditions. Unlike most integral equation techniques for mixed boundary value problems, the proposed method uses a global boundary charge density. As a result, Calderón identities can be utilized to avoid the use of hypersingular integral operators. More importantly, the formulation avoids spurious resonances. Numerical results illustrate the performance of the proposed solution technique.

Natalie Beams, Adrianna Gillman, Russell J. Hewett

The recently developed Hierarchical Poincaré-Steklov (HPS) method is a high-order discretization technique that comes with a direct solver. Results from previous papers demonstrate the method's ability to solve Helmholtz problems to high accuracy without the so-called pollution effect. While the asymptotic scaling of the direct solver's computational cost is the same as the nested dissection method, serial implementations of the solution technique are not practical for large scale numerical simulations. This manuscript presents the first parallel implementation of the HPS method. Specifically, we introduce an approach for a shared memory implementation of the solution technique utilizing parallel linear algebra. This approach is the foundation for future large scale simulations on supercomputers and clusters with large memory nodes. Performance results on a desktop computer (resembling a large memory node) are presented.

Yabin Zhang, Adrianna Gillman

Many applications involve solving several boundary value problems on geometries that are local perturbations of an original geometry. The boundary integral equation for a problem on a locally perturbed geometry can be expressed as a low rank update to the original system. A fast direct solver for the new linear system is presented in this paper. The solution technique utilizes a precomputed fast direct solver for the original geometry to efficiently create the low rank factorization of the update matrix and to accelerate the application of the Sherman-Morrison formula. The method is ideally suited for problems where the local perturbation is the same but its placement on the boundary changes and problems where the local perturbation is a refined discretization on the same geometry. Numerical results illustrate that for fixed local perturbation the method is three times faster than building a new fast direct solver from scratch.

Carlos Borges, Adrianna Gillman, Leslie Greengard

We describe a fast, stable algorithm for the solution of the inverse acoustic scattering problem in two dimensions. Given full aperture far field measurements of the scattered field for multiple angles of incidence, we use Chen's method of recursive linearization to reconstruct an unknown sound speed at resolutions of thousands of square wavelengths in a fully nonlinear regime. Despite the fact that the underlying optimization problem is formally ill-posed and non-convex, recursive linearization requires only the solution of a sequence of linear least squares problems at successively higher frequencies. By seeking a suitably band-limited approximation of the sound speed profile, each least squares calculation is well-conditioned and involves the solution of a large number of forward scattering problems, for which we employ a recently developed, spectrally accurate, fast direct solver. For the largest problems considered, involving 19,600 unknowns, approximately one million partial differential equations were solved, requiring approximately two days to compute using a parallel MATLAB implementation on a multi-core workstation.

Gary Marple, Alex Barnett, Adrianna Gillman et al.

This paper presents a new boundary integral equation (BIE) method for simulating particulate and multiphase flows through periodic channels of arbitrary smooth shape in two dimensions. The authors consider a particular system---multiple vesicles suspended in a periodic channel of arbitrary shape---to describe the numerical method and test its performance. Rather than relying on the periodic Green's function as classical BIE methods do, the method combines the free-space Green's function with a small auxiliary basis, and imposes periodicity as an extra linear condition. As a result, we can exploit existing free-space solver libraries, quadratures, and fast algorithms, and handle a large number of vesicles in a geometrically complex channel. Spectral accuracy in space is achieved using the periodic trapezoid rule and product quadratures, while a first-order semi-implicit scheme evolves particles by treating the vesicle-channel interactions explicitly. New constraint-correction formulas are introduced that preserve reduced areas of vesicles, independent of the number of time steps taken. By using two types of fast algorithms, (i) the fast multipole method (FMM) for the computation of the vesicle-vesicle and the vesicle-channel hydrodynamic interaction, and (ii) a fast direct solver for the BIE on the fixed channel geometry, the computational cost is reduced to $O(N)$ per time step where $N$ is the spatial discretization size. Moreover, the direct solver inverts the wall BIE operator at $t = 0$, stores its compressed representation and applies it at every time step to evolve the vesicle positions, leading to dramatic cost savings compared to classical approaches. Numerical experiments illustrate that a simulation with $N=128, 000$ can be evolved in less than a minute per time step on a laptop.