Fast Ewald summation for Stokes flow with arbitrary periodicity

TL;DR

The paper addresses the computation of three Stokes-flow potentials under arbitrary periodicity by extending the Spectral Ewald method to $D\in\{0,1,2,3\}$, achieving $O(N\log N)$ scaling. It develops unified, modified-kernel formulations, efficient discretization via gridding, PKB windowing, and adaptive Fourier transforms, along with automated parameter selection and rigorous error estimates. Key contributions include improved truncation-error bounds for $S$, $T$, and $\Omega$, analytical validation formulas for reduced periodicity, and a precomputation strategy for the fully aperiodic case, enabling accurate, fast simulations of Stokes flow in complex geometries. The practical impact is substantial for boundary-integral and potential-method simulations, where fast, spectrally accurate evaluation of Stokes potentials with arbitrary periodicity is essential.

Abstract

A fast and spectrally accurate Ewald summation method for the evaluation of stokeslet, stresslet and rotlet potentials of three-dimensional Stokes flow is presented. This work extends the previously developed Spectral Ewald method for Stokes flow to periodic boundary conditions in any number (three, two, one, or none) of the spatial directions, in a unified framework. The periodic potential is split into a short-range and a long-range part, where the latter is treated in Fourier space using the fast Fourier transform. A crucial component of the method is the modified kernels used to treat singular integration. We derive new modified kernels, and new improved truncation error estimates for the stokeslet and stresslet. An automated procedure for selecting parameters based on a given error tolerance is designed and tested. Analytical formulas for validation in the doubly and singly periodic cases are presented. We show that the computational time of the method scales like O(N log N) for N sources and targets, and investigate how the time depends on the error tolerance and window function, i.e. the function used to smoothly spread irregular point data to a uniform grid. The method is fastest in the fully periodic case, while the run time in the free-space case is around three times as large. Furthermore, the highest efficiency is reached when applying the method to a uniform source distribution in a primary cell with low aspect ratio. The work presented in this paper enables efficient and accurate simulations of three-dimensional Stokes flow with arbitrary periodicity using e.g. boundary integral and potential methods.

Figures (15)

• Figure 1: Illustration of how the truncation radius $R$ is set, two-dimensional case (corresponding to $D=1$). The rectangle is the extended box $\mathcal{B}'_{D\mathcal{P}}$ projected on the $x_2x_3$-plane, and $R$ is its diagonal, given by \ref{['eq:truncation-radius']}. In (a), the box is a square, while in (b), it has a higher aspect ratio. Since the modified kernels are based on radial kernels, $R$ must be the same in all directions (it must define a circle, not an ellipse). This means that the upsampling factor $s_0$, given by \ref{['eq:upsampling-factor-s0']}, will tend to be larger for boxes with high aspect ratio. For this reason, the SE method will be most efficient when applied to a box with low aspect ratio.
• Figure 2: Plot of window function $w_0$ with support $[-3h, 3h]$ and evaluation points as in \ref{['eq:1d-window']} indicated by points. The offset of $-x_{n,i}+p_i$ with respect to the uniform grid determines the offset of the evaluation points in this figure. The number of evaluation points (window size) is here $P=2 a_w/h=6$. The window function shown here is \ref{['eq:KB-window']} with $\beta=15$.
• Figure 3: Schematic illustration of adaptive Fourier transform (AFT), two-dimensional example. The upsampling factor $s$ in the free direction is different for different wavenumbers $\boldsymbol{k}^\mathcal{P}$ in the periodic direction, as given by \ref{['eq:adaptive-upsampling-factor']}. The areas of the shaded rectangles represent the number of discrete modes that are stored; for example, $\lvert \mathcal{K}_* \rvert s_* M'$ modes are stored for the set $\mathcal{K}_*$.
• Figure 4: Fourier-space part truncation errors and estimates. Actual errors are shown as colored symbols, estimates \ref{['eq:fs-trunc-est-rotlet']}--\ref{['eq:fs-trunc-est-stresslet']} are shown as solid black curves. For comparison, black crosses ($\times$) mark stokeslet and stresslet errors for $D=0$ when using the modified biharmonic \ref{['eq:general-biharmonic-0p']} with $a_B=b_B=0$, which was the case in afKlinteberg2017. Plots (a)--(c) show the three kernels for $\xi = 10$, $D=0,1,2,3$; (d) shows the stokeslet for $\xi = 5, 10, 15, 20$ (left to right), $D=3$. In all plots, there are $N=10^4$ random (see explanation in text) sources which also serve as evaluation points, and $L=2$, $Q=100$. Parameters other than $\xi$ and $h=\pi/k_\infty$ are selected such that errors other than the truncation error are negligible. The reference potential is given by the SE method with $k_\infty = 4\pi \xi$.
• Figure 5: Real-space part truncation errors and estimates. (a)--(c) show the three kernels for $\xi = 5, 10, 15, 20$ (right to left), $D=3$ (errors for other periodicities are almost identical and not shown). In all plots, there are $N=10^4$ random sources which also serve as evaluation points, and $L=2$, $Q=100$. The reference potential is given by the same method with $r_\mathrm{c} = L$. Actual errors are shown as symbols. Estimates \ref{['eq:rs-trunc-est-stokeslet']}--\ref{['eq:rs-trunc-est-rotlet']} are shown as solid curves.