Anna‐Karin Tornberg

David Krantz, Anna-Karin Tornberg

When numerically evaluating layer potentials at target points close to the domain boundary, specialized quadrature techniques are required for accuracy because of rapid variations in the integrand. To efficiently achieve a prescribed error tolerance, we introduce an adaptive quadrature method for smooth axisymmetric surfaces in which all algorithmic choices are determined automatically from the requested error tolerance. Standard quadrature is used wherever it is sufficient, while a specialized near-quadrature correction is applied only for those target points where additional accuracy is required. This correction combines singularity swap quadrature in the azimuthal direction with adaptive refinement in the polar direction; on the resulting refined polar grid, either standard quadrature or singularity swap quadrature is used depending on the predicted quadrature error. The method is coupled to a standard quadrature based on the trapezoidal rule in the azimuthal direction and Gauss--Legendre quadrature in the polar direction, and is activated only when that rule is predicted to be insufficient. Quadrature and interpolation error predictors are derived using complex analysis and are used to control both activation and refinement. While each surface is assumed to be axisymmetric, the layer density and the overall geometry need not be, allowing applications to configurations with multiple smooth axisymmetric bodies and patchwise discretizations. Numerical examples for Laplace, Helmholtz, and Stokes layer potentials demonstrate reliable error control across a range of geometries, including multi-body configurations.

Ludvig af Klinteberg, Anna-Karin Tornberg

In boundary integral methods it is often necessary to evaluate layer potentials on or close to the boundary, where the underlying integral is difficult to evaluate numerically. Quadrature by expansion (QBX) is a new method for dealing with such integrals, and it is based on forming a local expansion of the layer potential close to the boundary. In doing so, one introduces a new quadrature error due to nearly singular integration in the evaluation of expansion coefficients. Using a method based on contour integration and calculus of residues, the quadrature error of nearly singular integrals can be accurately estimated. This makes it possible to derive accurate estimates for the quadrature errors related to QBX, when applied to layer potentials in two and three dimensions. As examples we derive estimates for the Laplace and Helmholtz single layer potentials. These results can be used for parameter selection in practical applications.

Anna Broms, Anna-Karin Tornberg, Alex H. Barnett

We tackle two key difficulties in the simulation of the viscous hydrodynamics of a large dense collection of rigid particles: (i) the poor convergence rate of an iterative solution of the discretized linear system as particle gaps shrink, and (ii) the large number of unknowns needed to accurately discretize the resulting lubrication-driven flows. Our focus is the 2D Stokes resistance and mobility boundary value problems for nearly-touching disks. To address both challenges, we introduce a general two-body preconditioning strategy, and implement it with the method of fundamental solutions. For each close particle pair, the hard-to-resolve interaction is represented in a basis precomputed by solving a local boundary value problem on a fine grid. In an iterative solve, the resulting flow field corrects that obtained from a coarse representation of all particles. The local fine-grid correction can even be compressed so that all particles except the pair itself are affected by an equivalent set of coarse sources. Numerical experiments demonstrate rapid GMRES convergence in challenging multi-particle settings, with iteration counts remaining low even in densely packed suspensions. For example, the mobility problem is solved for a random close packing with area fraction $φ= 0.65$, $P = 10000$ monodisperse disks, and minimum separation $10^{-3}$, in just 47 GMRES iterations, achieving five digits of accuracy with 72 vector unknowns per body.

Davood Saffar Shamshirgar, Berk Hess, Anna-Karin Tornberg

The smooth particle mesh Ewald (SPME) method is an FFT based method for the fast evaluation of electrostatic interactions under periodic boundary conditions. A highly optimized implementation of this method is available in GROMACS, a widely used software for molecular dynamics simulations. In this article, we compare a more recent method from the same family of methods, the spectral Ewald (SE) method, to the SPME method in terms of performance and efficiency. We consider serial and parallel implementations of both methods for single and multiple core computations on a desktop machine as well as the Beskow supercomputer at KTH Royal Institute of Technology. The implementation of the SE method has been well optimized, however not yet comparable to the level of the SPME implementation that has been improved upon for many years. We show that the SE method is very efficient whenever used to achieve high accuracy and that it already at this level of optimization can be competitive for low accuracy demands.

Dag Lindbo, Anna-Karin Tornberg

A new method for interface tracking is presented. The interface representation, based on domain decomposition, provides the interface location explicitly, yet is Eulerian. This allows for well established finite difference methods on uniform grids to be used for the numerics of advecting the interface and other computations. CFL-stable and second order accurate explicit and implicit time-stepping methods are derived. Numerical results are given to substantiate stated convergence properties, as well as convergence in interface curvature and mass conservation. Our method is applied to a boundary integral formulation for Stokes flow, and the resulting integrals are analyzed and treated numerically to second order accuracy. Finally, we embed our method in the familiar immersed boundary- and immersed interface methods for two-phase Navier-Stokes flow.

Erik Boström, Anna-Karin Tornberg, Ludvig af Klinteberg

Fast Ewald summation efficiently evaluates Coulomb interactions and is widely used in molecular dynamics simulations. It is based on a split into a short-range and a long-range part, where evaluation of the latter is accelerated using the fast Fourier transform (FFT). The accuracy and computational cost depend critically on the mollifier in the kernel split and the window function used in the spreading and interpolation steps that enable the use of the FFT. The first prolate spheroidal wavefunction (PSWF) has optimal concentration in real and Fourier space simultaneously, and is used when defining both a mollifier and a window function. We provide a complete description of the method and derive rigorous error estimates. In addition, we obtain closed-form approximations of the Fourier truncation and aliasing errors, yielding explicit parameter choices for the achieved error to closely match the prescribed tolerance. Numerical experiments confirm the analysis: PSWF-based Ewald summation achieves a given accuracy with significantly fewer Fourier modes and smaller window supports than Gaussian- and B-spline-based approaches, providing a superior alternative to existing Ewald methods for particle simulations.

David Krantz, Ludvig af Klinteberg, Anna-Karin Tornberg

Dual-space multilevel kernel-splitting (DMK) is a fast summation framework that combines ideas from the fast multipole method, Ewald summation, and multilevel summation. Originally formulated for free-space problems, and later extended to fully periodic problems on a cube, it decomposes the kernel interaction into a smooth global contribution and a hierarchy of localized interactions evaluated on an octree. We extend DMK to problems on rectangular cuboids with periodic boundary conditions in one, two, or three coordinate directions. The periodization leverages the fact that interactions on all tree levels below the root are localized, allowing for their evaluation with minimal modification on a cubical tiling of the domain. The remaining smooth root-level far-field contribution is evaluated in Fourier space, with Fourier series in the periodic directions and Fourier integrals in the free directions. For reduced periodicity, truncated kernels are used to regularize singular and near-singular Fourier kernels, yielding rapidly convergent trapezoidal discretizations and a unified treatment of all periodicities. For large-aspect-ratio cuboids, the root-level sum can be accelerated using the fast Fourier transform. We validate the method for the electrostatic potential and Stokeslet, stresslet and rotlet potentials, for all periodicities and a wide range of aspect ratios. Numerical experiments show that the periodization adds only a small overhead to the original free-space DMK algorithm, also for high-aspect-ratio cuboids. The resulting method provides a framework for applying DMK to problems with mixed periodicity on rectangular cuboids, and extends naturally to other non-oscillatory kernels for which a kernel split is available.

Edith Frisk Gärtner, Fredrik Fryklund, Anna-Karin Tornberg

We present a fast and accurate potential theory-based method for the two-dimensional modified Helmholtz equation, treating the involved singular and nearly singular layer evaluations together with volume potentials within a single computational framework. The method is based on a decomposition of the free-space Green's function into a short-range local part and a smooth long-range part. The long-range contribution is evaluated efficiently using the non-uniform fast Fourier transform (NUFFT), while the local contribution is treated by asymptotic expansions. For the layer potentials, an intermediate telescoping sum over dyadic refinement levels is added, where the resulting difference kernels are smooth and rapidly decaying, allowing the dyadic levels to be evaluated without specialized quadrature rules. The volume potential is evaluated on triangular cut-cell meshes, where the mesh only enters the scheme as quadrature rule for smooth data. This makes the method robust with respect to small and distorted mesh cells, without the need for stabilization or cell-merging techniques. Numerical experiments demonstrate the expected convergence rates, high throughput of the potential evaluations, and robustness with respect to mesh quality.