A fast multipole method for the three-dimensional Stokes equations
Introduction
During the last two decades, a substantial body of work has emerged on fast algorithms for potential theory. The bulk of these have considered the Coulombic or gravitational N-body problem or the solution of related boundary integral equations. The N-body problem requires evaluations of sums on the formwhere , is an arbitrary orientation vector, and . A typical boundary integral equation takes the formwhere Γ is a surface in , , refers to the normal derivative at the point , and are constants that depend on the particular boundary condition being imposed. In general, F(x) is given as data and σ is an unknown surface density. In physical terms, the first sum in (1) is the field due to a collection of charges and the second sum is the field due to a collection of dipoles. Likewise, the first integral in (2) is the field due to a surface charge distribution (a single layer potential) and the second integral is the field due to a continuous distribution of normally oriented dipoles (a double layer potential). Discretizing the integrals in (2) by some quadrature rule, we obtain discrete sums over the quadrature points of the same form as (1).
In short, N-body calculations arise in two main contexts: the evaluation of particle interactions and the evaluation of layer potentials (often as part of the iterative solution of boundary integral equations). Since the governing equation is the Laplace equation, we will refer to these as harmonic N-body calculations.
The straightforward evaluation of in (1) for clearly requires work. There are, however, a variety of fast algorithms that are capable of reducing that cost to O(N) or . These include fast multipole methods, tree codes, the method of local corrections, multigrid methods, panel clustering methods, particle-in-cell methods, particle-mesh Ewald methods, and pre-corrected FFT methods. We will not attempt to review the literature on fast algorithms, but refer the reader to a few selected papers on the fast multipole method (FMM) [1], [2], [6], [8], [15], since it is the algorithm we will rely on here. It is worth noting, in the present context, that both particle-mesh-Ewald methods and pre-corrected FFT methods have been extended to the case of Stokes flow and used with great effect [10], [13]. The principal advantage of the FMM is that it is fully adaptive, and handles highly inhomogeneous source distributions as easily as it does homogeneous ones. It has the disadvantage that it is much more complex to implement efficiently.
In Stokes flow and linear elasticity, computations analogous to (1) arise naturally. Both involve vector versions of the N-body problem. In this short paper, we concentrate on the Stokes problem.
One fundamental solution for the Stokes equations is the Stokesletwhere is the Kronecker delta, and . The corresponding summation problem is the calculation of the vector at each source location from the formulawhere are the vector source strengths. A second fundamental solution is the Stressletwith an orientation vector. The corresponding summation problem is to evaluate as given bywhere is an orientation vector at the th source location and are the vector source strengths. and can be thought of as velocity fields induced by point forces or surface stresses, respectively [9]. (Integrals involving Sij and Dij are referred to as single and double layer potentials, by analogy with the electrostatic case).
Several fast multipole methods for these kernels (or the closely related ones of linear elasticity) have been developed based essentially on expansion of the Green’s function for the biharmonic equation [3], [7], [8], [12], [16]. A clever rearrangement of terms allows for the use of a total of four sets of multipole/local expansion coefficients. The cost of the method is, therefore, approximately four times that of a harmonic FMM. Recently, Wang et al. [11] have presented an efficient, parallel implementation for Stokeslet and Stresslet summations along these lines. However, none of these schemes make direct use of existing “black-box” harmonic FMMs.
Another useful approach, developed by Ying et al. [15], is the kernel independent fast multipole algorithm that can be applied to essentially any non-oscillatory kernel. This allows for straightforward “black-box” application to the Stokes equations, but at a cost of approximately six to nine scalar FMM calls (three calls with 3N sources each).
Finally, one can rewrite the Stokeslet and Stresslet summation formulas in such a way that the harmonic FMM (or any other fast electrostatic method) can be used in “black-box” fashion. This approach was taken by Wang and LeSar for dislocation dynamics [14], by Fu et al. for linear elasticity [4], and by Fu and Rodin for Stokes flow [5]. The latter formulation requires four harmonic FMM calls for Stokeslets and twelve harmonic FMM calls for Stresslets.
Here, we follow the last approach and present a new, simple and efficient decomposition of the Stokeslet and Stresslet summations into four harmonic N-body problems each, although the Stresslet case requires a little care (see below). It is, perhaps, worth citing three reasons for choosing to revisit this problem. First, the reformulation of the Stresslet in terms of harmonic interactions appears to be new. Second, using this approach, improvements to any fast algorithm for harmonic interactions, such as low-level code optimization, hardware acceleration, or parallel implementation, become available for fluid dynamic and elasticity applications. Third, the storage cost of the harmonic FMM is about half that of the biharmonic case, providing some advantage for large-scale simulations.
Section snippets
Decomposition of the stokeslet and stresslet sums
We begin by observing that the Stokeslet (3) can be written asFrom this, we show in Appendix A that the summation in (4) can be written aswhere and .
That is, to compute , , we need to evaluate four harmonic sums over the locations with source strengths , and , respectively. For this, we use the FMM as described in [2], which returns both the potential and
Numerical results
In this section, we present some timing results for the FMM-based Stokeslet and Stresslet sums described by (4), (6). All calculations were carried out on a laptop computer with a 1.2 GHz Pentium M processor and 500 Mb of RAM. We also present the time for an efficient implementation of direct summation, which uses only one square root evaluation to generate all nine matrix entries in (3), (5). For the sake of comparison, we also present timings of the FMM and the direct method for simple
Conclusions
In this paper, we have presented a fast multipole method (FMM) for Stokeslet and Stresslet calculations, based on the use of the harmonic FMM. Each of the Stresslet and Stokeslet summations requires four calls to such an FMM, which has been trivially modified so that the net cost scales like three harmonic interactions. A feature of our approach is that additional improvements to any fast algorithm for harmonic interactions, from code optimization to specialized hardware, become immediately
Acknowledgement
L.G. acknowledges support from DOE Grant DE-FG02-88ER25053. A.K.T. acknowledges support from NSF grant DMS-0412203 and an Alfred P. Sloan Research Fellowship. A.K.T. is also a Royal Swedish Academy of Sciences Research Fellow supported by a grant from the Knut and Alice Wallenberg Foundation.
References (16)
- et al.
A fast adaptive multipole algorithm in three-dimensions
J. Comput. Phys.
(1999) A fast multipole implementation of the qualocation mixed-velocity-traction approach for exterior Stokes flow
Eng. Anal. Bound. Elem.
(2005)- et al.
A fast algorithm for particle simulations
J. Comput. Phys.
(1987) - et al.
Fast multipole method for the biharmonic equation in three dimensions
J. Comput. Phys.
(2006) - et al.
A kernel independent adaptive fast multipole algorithm in two and three-dimensions
J. Comput. Phys.
(2004) - R. Beatson, L. Greengard, A short course on fast multipole methods, In: Wavelets, Multilevel methods and Elliptic PDEs,...
- et al.
A fast solution method for three-dimensional many-particle problems of linear elasticity
Int. J. Num. Meth. Eng.
(1998) - et al.
Fast solution methods for three-dimensional Stokesian many-particle problems
Commun. Numer. Meth. En.
(2000)
Cited by (94)
Efficient convergent boundary integral methods for slender bodies
2024, Journal of Computational PhysicsAdjoint-based control of three dimensional Stokes droplets
2023, Journal of Computational PhysicsFast Ewald summation for Stokes flow with arbitrary periodicity
2023, Journal of Computational PhysicsA generalised drift-correcting time integration scheme for Brownian suspensions of rigid particles with arbitrary shape
2022, Journal of Computational PhysicsA numerical method for suspensions of articulated bodies in viscous flows
2022, Journal of Computational PhysicsA domain decomposition solution of the Stokes-Darcy system in 3D based on boundary integrals
2022, Journal of Computational Physics