Current work with my postdoc advisor, Adrianna Gillman.

The Hierarchical Poincaré-Steklov (HPS) method is a discretization technique based on domain decomposition and classical spectral collocation methods. It is accurate and robust, even for highly oscillatory solutions. An associated fast direct solver is built from hierarchically merging Poincaré-Steklov operators on box boundaries and storing the necessary information; the solve stage then involves applying the operators in a series of small matrix-vector multiplications.

We are working on parallelizing, optimizing, and accelerating the algorithm, especially the build stage. This will be crucial in 3D. Our 2D shared memory parallelization technique has already reduced the build stage times from approximately 10 minutes to approximately 30 seconds for over 3 million unknowns.

Speedup achieved in the build and solve portions of the algorithm with 56 available threads.

We developed and implemented a new, high order finite element-integral equation coupling method for elliptic interface problems. The method leverages the strengths of FE and IE methods to handle different aspects of the problem. It can handle general jump conditions at the interface, and the jump conditions appear only in the right hand side of the system to be solved. Additionally, method does not suffer from loss of convergence order near the boundary when the boundary is smooth and the problem data can be extended as necessary across the interface. It also doesn't require the construction or use of any special basis functions; in fact, much of the standard machinery for both the FE and IE solvers can be used, simplifying implementation from a software standpoint.

An interface problem on a starfish-shaped domain.
An example interface problem with jump conditions along a starfish-shaped interface.

The goal of this project was to create an FE-based version of the particle-particle–particle-mesh (P3M) methods for N-body problems. P3M methods separate particle interactions into short-range (only required for near neighbors) and far-field (smooth, and solved easily on a mesh with your favorite numerical method). A common choice is to achieve the P-P/P-M splitting through Gaussian screen functions (as in the classic Ewald sum) and solve the mesh problem with FFTs.

We developed a method which constructs special polynomial screen functions out of finite element basis functions. Thus the screens are represented exactly on the finite element mesh, removing several possible sources of error in a finite element solution to the mesh problem. Indeed, we aimed to make the mesh problem especially suited to being solved with finite elements, avoiding the geometry restrictions and parallel communication burdens of the FFT.

Example polynomial screens.
Polynomial screens formed from increasing order of basis functions for a charge located at the starred location.