Schedule for: 22w5170 - Outstanding Challenges in Computational Methods for Integral Equations

A unified treatment of boundary integral operators in 2D and 3D using a simple trapezoidal quadrature method is presented. The method is based on generalized Euler-Maclaurin formulas and can be applied to weakly singular as well as hypersingular operators. The construction of such a quadrature rule for a given kernel can be done systematically in simple steps, which we will demonstrate in this talk and show numerical examples.

In the application of the Fast Multipole Method to the computation of potentials for elliptic PDEs and systems thereof, opportunities exist for lowering cost through knowledge of the kernel and the PDE operator. We present two methods that, given various small amounts of user-supplied problem knowledge (e.g. symbolic expression of the kernel, symbolic PDE), automatically exploit these opportunities. The first is devoted to the automatic synthesis of translation operators (e.g. multipole-to-local, point-to-multipole, etc.) for arbitrary kernels. We describe the asymptotic cost of variants of our algorithm available given certain pieces of information, as well as the methods by which they are attained. We present theoretical cost bounds as well as numerical evidence that our algorithms attain them. The second builds on the first and is devoted to the automatic synthesis of execution plans for expressions of potential operators involving multiple inputs and outputs, multiple different kernels, composition, as well as source and target derivatives. Given a symbolic description of such an operator, our system outputs a sequence of operations that realizes cost savings through an algebraic procedure based on syzygies. Finally, we describe the application of the combination of these approaches in a system for the high-order accurate evaluation of layer potentials from unstructured geometries in two and three dimensions.

We introduce a fast direct solver for variable-coefficient elliptic partial differential equations on surfaces based on the hierarchical Poincaré–Steklov method. The method takes as input a high-order quadrilateral mesh of a surface and discretizes surface differential operators on each element using a high-order spectral collocation scheme. Elemental solution operators and Dirichlet-to-Neumann maps tangent to the surface are precomputed and merged in a pairwise fashion to yield a hierarchy of solution operators that may be applied in $\mathcal{O}(N \log N)$ operations for a mesh with $N$ elements. The resulting fast direct solver may be used to accelerate implicit time-stepping schemes, as the precomputed operators can be reused for fast elliptic solves on surfaces. We apply the method to a range of problems on both smooth surfaces and surfaces with sharp corners and edges, including the static Laplace–Beltrami problem, the Hodge decomposition of a tangential vector field, and some time-dependent reaction–diffusion systems.

The talk presents a fast direct solution technique for solving two-dimensional wave scattering problems from quasi-periodic multilayered structures. The computational cost of creating the direct solver scales O(N) where N is the total number of discretization points on all interfaces. The bulk of the precomputation can be re-used for any choice of the incident wave. As a result, the direct solver can solve over 200 scattering problems involving an eleven-layer geometry with complex interfaces 100 times faster than building a new fast direct solver from scratch for each new set of boundary data. An added benefit of the presented solver is that building an updated solver for a new geometry involving a replaced interface or a change in material property in one layer is inexpensive compared to building a new fast direct solver from scratch. Numerical results illustrate the improved performance of the new solver over some previous approaches.

Bandlimited functions arise in a wide variety of applications in scientific computing and signal processing. In this talk, we present a fast algorithm for computing quadratures for bandlimited functions, based on recent advances in the numerical treatment of prolate spheroidal wave functions. The resulting quadrature rules are capable of integrating functions with a given bandlimited to high accuracy and can be computed rapidly, with only O(n) operations required to compute an n-point rule with fixed bandlimit.

In this paper we present the transpose method for quasi-Newton optimization problems in acoustics and electromagnetism. In particular we use a fast algorithm for optimizing the location of p.e.c. spheres in free space such that a certain goal function is optimized. The goal function can be used to maximize the total field in a certain point in space (focal point) or any other function that depends on the scattered field produced by the spheres. The method can be used to create new materials (metamaterials) or goal-oriented materials with uncommon electromagnetic/acoustic properties.

Biomolecular electrostatics is key in protein function and the chemical processes affecting it. Implicit-solvent models expressed by the Poisson-Boltzmann (PB) equation can provide insights with less computational power than full atomistic models, making large-system studies -- at the scale of viruses, for example -- accessible to more researchers. In this talk we present a Galerkin BEM approach to compute solvation for large protein models based on the black-box coupling of the Bempp software with the Exafmm kernel-independent FMM library. We discuss the implementation of the coupling and scalability for large problems. We conclude the talk with remarks on ongoing work to develop FEM/BEM coupled solvers for solvation problems with inhomogeneities, and an outlook to software development for extreme scale solvation models.

Oscillatory systems are ubiquitous in physics: they arise in celestial and quantum mechanics, electrical circuits, molecular dynamics, and beyond. Yet even in the simplest case, when the frequency of oscillations changes slowly but is large, the vast majority of numerical methods struggle to solve such equations. Methods based on approximating the solution with polynomials are forced to take $\mathcal{O}(k)$ timesteps, where $k$ is the characteristic frequency of oscillations. This scaling can generate unacceptable computational costs when the ODE in question needs to be solved billions of times, e.g.\ as the forward modelling step of Bayesian parameter estimation. In this talk I will introduce an efficient method for solving 2nd order, linear ODEs with highly oscillatory solutions. The solver employs two methods: in regions where the solution varies slowly, it uses a spectral method based on Chebyshev nodes and with an adaptive stepsize, but in the highly oscillatory phase it automatically switches over to an asymptotic method. The asymptotic method constructs a nonoscillatory phase function solution of the Riccati equation associated with the ODE. In the talk I will present how the method fits in the landscape of oscillatory solvers, the theoretical underpinnings of the asymptotic solver, a summary of the switching and stepsize-update algorithms, some examples, and a brief error analysis.

Time-harmonic wave phenomena in inhomogenous media are governed by Helmholtz and Maxwell equations with variable coefficients, and are typically simulated with finite-difference (FD)/finite-element-based differential equation solvers or volume integral equation based solvers (VIE). Fast, accurate and stable algorithms for solving these problems in the high-frequency regime are computationally very challenging. Although a few recent works have leveraged the so-called butterfly compression techniques to construct fast FD and VIE-based direct solvers, they suffer from a few other computational issues. The FD-based solver is plagued with numerical dispersion, PML truncation error and zero-pivoting during sparse matrix inversion, hence cannot handle large systems with high-order accuracy. The VIE solver, on the other hand, requires inverting a large dense linear system and turns out to be still expensive even with butterfly acceleration. In this work, we consider another approach called Hadamard-Babich integrator, which represents a high-frequency ansatz of the Green's function for inhomogenous media. We first construct low-rank products of the phase and amplitude ingredients of the Babich integrator and use the results to construct butterfly compression of the discretized integral operator for the entire computation domain. The resulting Babich integrator-based solver is very accurate for smoothly varying media and computationally very efficient compared to FD or VIE solvers. When further combined with surface integral equation (SIE) formulations, the proposed solver also applies to large 2D and 3D domains with surface inclusions or multiple regions. This is a joint work with Jianliang Qian, Jian Song from MSU and Robert Burridge from UNM.

Many-body Green's function methods are of central importance in modern approaches to computational quantum physics which have attempted to reach higher accuracies than those provided by effective one-body approximations like density functional theory. These Green's functions satisfy nonlinear integral equations, called Dyson equations, and present an interesting set of algorithmic challenges which may be possible to address, in part, using ideas developed in the computational integral equations community. This talk will discuss a few of these challenges, compare and contrast them with those arising in more standard computational integral equations problems, and introduce new fast algorithms for Dyson equations.

Moderator: A. Gillman Panelists: Andreas Kloeckner, Timo Betcke, and Manas Rachh