A kernel-independent treecode based on barycentric Lagrange interpolation

TL;DR

This paper introduces a kernel-independent treecode that uses barycentric Lagrange interpolation to efficiently compute particle interactions, reducing complexity from quadratic to near-linear time, suitable for non-oscillatory kernels.

Contribution

The paper presents a novel kernel-independent treecode based on barycentric Lagrange interpolation, enabling efficient and scalable particle interaction computations for non-oscillatory kernels.

Findings

Reduces interaction computation from O(N^2) to O(N log N)

Demonstrates effectiveness for 3D Stokeslet and rotlet systems

Low memory usage allows easy parallelization

Abstract

A kernel-independent treecode (KITC) is presented for fast summation of particle interactions. The method employs barycentric Lagrange interpolation at Chebyshev points to approximate well-separated particle-cluster interactions. The KITC requires only kernel evaluations, is suitable for non-oscillatory kernels, and it utilizes a scale-invariance property of barycentric Lagrange interpolation. For a given level of accuracy, the treecode reduces the operation count for pairwise interactions from $O(N^2)$ to $O(N\log N)$, where $N$ is the number of particles in the system. The algorithm is demonstrated for systems of regularized Stokeslets and rotlets in 3D, and numerical results show the treecode performance in terms of error, CPU time, and memory consumption. The KITC is a relatively simple algorithm with low memory consumption, and this enables a straightforward OpenMP parallelization.

Tables

Table 1. Table 1: Example 1, regularized Stokeslets in a cube ( 20 ), MRS parameter ϵ = 0.02 italic-ϵ 0.02 \epsilon=0.02 , system size N 𝑁 N , CPU time (s) for direct sum and treecode ( d , t 𝑑 𝑡 d,t ), speedup ( d / t 𝑑 𝑡 d/t ), treecode error ( E 𝐸 E ), results are from Fig.  3 with MAC parameter θ = 0.7 𝜃 0.7 \theta=0.7 , degree n = 7 𝑛 7 n=7 , maximum leaf size N 0 = 2000 subscript 𝑁 0 2000 N_{0}=2000 .

$N$	direct sum (s), $d$	treecode (s), $t$	speedup, $d/t$	error, $E$
10K	0.54	0.51	1.06	2.25e-6
80K	34.01	9.78	3.48	1.44e-5
640K	2248.93	149.14	15.08	3.17e-5
5.12M	157,542.13	1892.57	83.24	1.62e-5
40.96M	10,082,696.32	20,966.04	480.91	7.06e-6

Table 2. Table 2: Example 1, regularized Stokeslets in a cube ( 20 ), MRS parameter ϵ = 0.02 italic-ϵ 0.02 \epsilon=0.02 , direct sum and KITC memory consumption ( MB ) MB \rm{(MB)} , system size N 𝑁 N , treecode MAC parameter θ = 0.7 𝜃 0.7 \theta=0.7 , degree n = 1 , 3 , 5 , 7 , 9 𝑛 1 3 5 7 9 n=1,3,5,7,9 , maximum leaf size N 0 = 2000 subscript 𝑁 0 2000 N_{0}=2000 .

system size $N$	10K	80K	640K	5.12M	40.96M
direct sum memory (MB)	0.88	7.04	56.36	450.60	3604.80
treecode memory (MB)
$n=1$	0.89	7.09	56.74	453.69	3629.24
$n=3$	0.91	7.27	58.17	465.09	3720.46
$n=5$	0.98	7.74	61.96	495.37	3962.69
$n=7$	1.16	8̇.66	69.28	553.96	4431.40
$n=9$	1.44	10.16	81.32	650.30	5202.10

Table 3. Table 3: Example 1, regularized Stokelets in a cube ( 20 ), system size N = 𝑁 absent N\,{=}\, 640K, MRS parameter ϵ italic-ϵ \epsilon , microorganism length l = ϵ 𝑙 italic-ϵ l=\epsilon , error E 𝐸 E , treecode parameters θ = 0.7 , n = 7 formulae-sequence 𝜃 0.7 𝑛 7 \theta=0.7,n=7 , maximum leaf size N 0 = 2000 subscript 𝑁 0 2000 N_{0}=2000 .

MRS parameter, $ϵ$	0.005	0.01	0.02	0.04	0.08
error, $E$	7.68e-6	2.08e-5	3.17e-5	2.78e-5	1.93e-5

Table 4. Table 4: Example 1, regularized Stokeslets in a cube ( 20 ), parallel performance, system size N = 𝑁 absent N\,{=}\, 640K, MRS parameter ϵ = 0.02 italic-ϵ 0.02 \epsilon=0.02 , KITC parameters θ = 0.7 , n = 7 , N 0 = 2000 formulae-sequence 𝜃 0.7 formulae-sequence 𝑛 7 subscript 𝑁 0 2000 \theta=0.7,n=7,N_{0}=2000 , error E 𝐸 E = 3.17e-5, number of cores ( n ​ c 𝑛 𝑐 nc ), CPU time (s) is shown for direct sum and KITC ( d , t ) 𝑑 𝑡 (d,t) , parallel speedup = d 1 / d n ​ c , t 1 / t n ​ c absent subscript 𝑑 1 subscript 𝑑 𝑛 𝑐 subscript 𝑡 1 subscript 𝑡 𝑛 𝑐 =d_{1}/d_{nc},t_{1}/t_{nc} , parallel efficiency (PE %) = ( d 1 / d n ​ c ) / n ​ c , ( t 1 / t n ​ c ) / n ​ c subscript 𝑑 1 subscript 𝑑 𝑛 𝑐 𝑛 𝑐 subscript 𝑡 1 subscript 𝑡 𝑛 𝑐 𝑛 𝑐 (d_{1}/d_{nc})/nc,(t_{1}/t_{nc})/nc , parallel KITC speedup ( d / t ) 𝑑 𝑡 (d/t) .

$nc$	$d$ time (s)	$d_1/d_{nc}$	$d$ PE (%)	$t$ time (s)	$t_1/t_{nc}$	$t$ PE (%)	$d/t$
1	2249.71	1.00	100.00	149.75	1.00	100.00	15.02
2	1123.72	2.00	100.00	75.01	2.00	99.82	14.98
4	569.58	3.95	98.74	38.83	3.86	96.41	14.67
8	311.50	7.22	90.28	21.44	6.98	87.31	14.53
12	215.49	10.43	87.00	16.34	9.16	76.37	13.19
24	108.22	20.79	86.62	10.19	14.70	61.23	10.62

Table 5. Table 5: Example 2, regularized Stokeslets and rotlets on an array of helical rods ( 21a )-( 21b ) following  [ 47 ] , parallel performance, N r = 80 2 subscript 𝑁 𝑟 superscript 80 2 N_{r}=80^{2} rods, system size N = 966400 𝑁 966400 N=966400 , MRS parameter ϵ = 0.3 italic-ϵ 0.3 \epsilon=0.3 , KITC parameters θ = 0.7 , n = 7 formulae-sequence 𝜃 0.7 𝑛 7 \theta=0.7,n=7 , N 0 = 1000 subscript 𝑁 0 1000 N_{0}=1000 , error E 𝐸 E = 2.24e-5, number of cores ( n ​ c 𝑛 𝑐 nc ), CPU time (s) is shown for direct sum and KITC ( d , t ) 𝑑 𝑡 (d,t) , parallel speedup = d 1 / d n ​ c , t 1 / t n ​ c subscript 𝑑 1 subscript 𝑑 𝑛 𝑐 subscript 𝑡 1 subscript 𝑡 𝑛 𝑐 d_{1}/d_{nc},t_{1}/t_{nc} , parallel efficiency (PE %) = ( d 1 / d n ​ c ) / n ​ c , ( t 1 / t n ​ c ) / n ​ c subscript 𝑑 1 subscript 𝑑 𝑛 𝑐 𝑛 𝑐 subscript 𝑡 1 subscript 𝑡 𝑛 𝑐 𝑛 𝑐 (d_{1}/d_{nc})/nc,(t_{1}/t_{nc})/nc , parallel KITC speedup ( d / t ) 𝑑 𝑡 (d/t) .

$nc$	$d$ time (s)	$d_1/d_{np}$	$d$ PE (%)	$t$ time (s)	$t_1/t_{np}$	$t$ PE (%)	$d/t$
1	15651.61	1.00	100.00	720.87	1.00	100.0	21.71
2	7907.10	1.98	98.97	371.91	1.94	96.91	21.26
4	3962.36	3.95	98.75	191.89	3.76	93.91	20.65
8	2120.36	7.38	92.27	110.88	6.50	81.27	19.12
12	1424.27	10.99	91.58	80.35	8.97	74.76	17.73
24	742.85	21.07	87.79	43.63	16.52	68.84	17.03

Equations

$$u({\mathbf{x}}_i)=\sum _{j=1}^{N}k({\mathbf{x}}_i,{\mathbf{x}}_j)f_j,i=1,\dots ,N,$$

$$k(\mathbf{x},\mathbf{y})\approx \sum _{k=0}^n{\varphi }_k(\mathbf{x}){\psi }_k(\mathbf{y}).$$

$$p_n(t)=\sum _{k=0}^n{f}_k{L}_k(t),f_k=f(s_k),k=0,\dots ,n,$$

$$L_k(t)=\frac{{\Pi }_{j=0,j\ne k}^n(t-s_j)}{{\Pi }_{j=0,j\ne k}^n(s_k-s_j)},k=0,\dots ,n,$$

$$k(x,y)\approx \sum _{k=0}^{n}k(x,s_k)L_k(y).$$

$$k(\mathbf{x},\mathbf{y})\approx \sum _{k_1=0}^n\sum _{k_2=0}^n\sum _{k_3=0}^{n}k(\mathbf{x},{\mathbf{s}}_{\mathbf{k}})L_{k_1}(y_1)L_{k_2}(y_2)L_{k_3}(y_3).$$

$$L_k(t)=\frac{\frac{w_k}{t-s_k}}{\sum _{k=0}^n\frac{w_k}{t-s_k}},w_k=\frac{1}{{\Pi }_{j=0,j\ne k}^n(s_k-s_j)},k=0,\dots ,n,$$

$$s_k=\cos {\theta }_k,{\theta }_k=\frac{k\pi }{n},k=0,\dots ,n.$$

$$w_k=(-1{)}^k{\delta }_k,{\delta }_k=\{\begin{array}{ll}1/2,& k=0\text{or}k=n,\\ 1,& k=1,\dots ,n-1.\end{array}$$

$$p_n(t)=\sum _{k=0}^n\frac{\frac{w_k}{t-s_k}}{\sum _{k=0}^n\frac{w_k}{t-s_k}}{f}_k,$$

$${\bar{s}}_k=\cos {\bar{\theta }}_k,{\bar{\theta }}_k=\frac{(2k-1)\pi }{2n},k=1,\dots ,n,$$

$$p_{n-1}(t)=\sum _{k=1}^n{f}_k{\bar{L}}_k(t),{\bar{L}}_k(t)=\frac{1}{n}+\frac{2}{n}\sum _{j=1}^{n-1}{T}_j({\bar{s}}_k)T_j(t),k=1,\dots ,n,$$

$$u({\mathbf{x}}_i)=\sum _{j=1}^{N}k({\mathbf{x}}_i,{\mathbf{x}}_j)f_j=\sum _{C}u({\mathbf{x}}_i,C),i=1,\dots ,N,$$

$$u({\mathbf{x}}_i,C)=\sum _{{\mathbf{y}}_j\in C}k({\mathbf{x}}_i,{\mathbf{y}}_j)f_j,$$

$$u({\mathbf{x}}_i,C)\approx \sum _{k_1=0}^n\sum _{k_2=0}^n\sum _{k_3=0}^{n}k({\mathbf{x}}_i,{\mathbf{s}}_{\mathbf{k}}){\hat{f}}_{\mathbf{k}},$$

$${\hat{f}}_{\mathbf{k}}=\sum _{{\mathbf{y}}_j\in C}{L}_{k_1}(y_{j1})L_{k_2}(y_{j2})L_{k_3}(y_{j3})f_j,k_1,k_2,k_3=0,\dots ,n,$$

$$L_{k_{\ell }}(y_{j\ell })=\frac{\frac{w_{k_{\ell }}}{y_{j\ell }-s_{k_{\ell }}}}{\sum _{k=0}^n\frac{w_k}{y_{j\ell }-s_k}}=\frac{a_{j,k_{\ell },\ell }}{\sum _{k=0}^n{a}_{j,k,\ell }},a_{j,k,l}=\frac{w_k}{y_{jl}-s_k},$$

$$\frac{r}{R}\le \theta ,$$

$$H_1(r)$$

$$D_1(r)$$

$$\mathbf{u}(\mathbf{x})=\sum _{j=1}^N\left({\mathbf{f}}_j{H}_1(r_j)+\left[{\mathbf{f}}_j\cdot (\mathbf{x}-{\mathbf{y}}_j)\right](\mathbf{x}-{\mathbf{y}}_j)H_2(r_j)\right),$$

$$\mathbf{u}(\mathbf{x})$$

$$\mathbf{w}(\mathbf{x})$$

$$E=\left(\sum\limits_{i=1}^{N}|{\bf u}^{d}({\bf x}_{i})-{\bf u}^{t}({\bf x}_{i})|^{2}\bigg{/}\sum\limits_{i=1}^{N}|{\bf u}^{d}({\bf x}_{i})|^{2}\right)^{\!\!1/2},$$

Full text

1 Introduction

Consider the problem of evaluating the sum

[TABLE]

where $u({\bf x}_{i})$ is a velocity (or potential or force) and $\{{\bf x}_{i}\}\subset\mathbb{R}^{d}$ is a set of particles with weights $\{f_{i}\}$. Depending on the application, the velocity and weights may be scalars or vectors, and the kernel may be a tensor. The kernel $k({\bf x},{\bf y})$ describes the interaction between a target particle ${\bf x}$ and a source particle ${\bf y}$, and we are interested in non-oscillatory kernels that are smooth for ${\bf x}\neq{\bf y}$ and decay slowly for $|{\bf x}-{\bf y}|\to\infty$. It is understood that if the kernel is singular for ${\bf x}={\bf y}$, then the sum omits the $i=j$ term.

These types of sums arise in particle simulations involving point masses, point charges, and point vortices, as well as in boundary element methods where the particles are quadrature points. Evaluating (1) by direct summation requires $O(N^{2})$ operations which is prohibitively expensive when $N$ is large, and several fast methods have been developed to reduce the cost. One can distinguish between two types of methods, particle-mesh methods in which the particles are projected onto a uniform mesh where the FFT or multigrid can be used (e.g. P3M [29], particle-mesh Ewald [16], spectral Ewald [2], multilevel summation [9, 27]), and tree-based methods in which the particles are partitioned into a hierarchy of clusters with a tree structure and the particle-particle interactions are replaced by particle-cluster or cluster-cluster approximations (e.g. treecode [4], fast multipole method (FMM) [24], panel clustering [26]).

Tree-based methods. The present work is concerned with tree-based methods that rely on degenerate kernel approximations of the form,

[TABLE]

Such approximations can be classified as near-field/local or far-field/multipole depending on their domain of validity in the variables ${\bf x},{\bf y}$. The treecode originally used a far-field monopole approximation for the Newtonian potential [4], while the FMM improved on this by employing higher-order multipole and local approximations, in particular using Laurent series for the 2D Laplace kernel and spherical harmonics for the 3D Laplace kernel [24, 25]. Later versions of the FMM used plane wave expansions for the 3D Laplace kernel [10] and spherical Bessel function expansions for the Yukawa potential [23]. Methods based on Cartesian Taylor expansions were also developed for some common kernels [13, 15, 37, 51, 35, 30, 57].

Kernel-independent methods. The tree-based methods cited above rely on analytic series expansions specific to each kernel and alternative approximation methods have been investigated. An early example in this direction was an FMM for Laplace kernels based on discretizing the Poisson integral formula [3], and this was followed by a pseudoparticle method that reproduces the multipole moments for these kernels [38]. Later work developed approximations suitable for a wide class of non-oscillatory kernels. One approach based on polynomial interpolation ([32], section 11.4) has been applied in the context of multilevel approximation [21], hierarchical matrices [7], and the black-box FMM (bbFMM) [19]. An alternative method employed in the kernel-independent FMM (KIFMM) uses equivalent densities [60, 61], and other kernel-independent FMMs use Legendre expansions [22], matrix compression based on skeletonization [41], and truncated Fourier series [62]. Recently an FMM based on the Cauchy integral formula and Laplace transform was proposed for general analytic functions [34], and a kernel-independent treecode was developed using approximate skeletonization for particle systems in high dimensions [40].

Present work. There is ongoing interest in exploring different strategies for fast summation of particle interactions, and the present work contributes a kernel-independent treecode (KITC) with operation count $O(N\log N)$ in which the far-field approximation uses barycentric Lagrange interpolation at Chebyshev points [6, 55]. The barycentric Lagrange interpolant can be efficiently implemented and has good stability properties [46, 28, 42]; the 1D case is reviewed in [6, 55] and here we apply it in 3D using a tensor product to compute well-separated particle-cluster approximations.

It should be noted that the bbFMM [19] and KITC both use polynomial interpolation, but they differ in two ways. The first difference concerns the interpolating polynomial; the bbFMM uses Chebyshev points of the 1st kind (roots of Chebyshev polynomials) and expresses the Lagrange polynomials in terms of Chebyshev polynomials, while the KITC uses Chebyshev points of the 2nd kind (extrema of Chebyshev polynomials) and expresses the Lagrange polynomials in barycentric form; as explained in Section 3, this enables the KITC to take advantage of the scale-invariance property of barycentric Lagrange interpolation. The second difference concerns the algorithm structure; the KITC uses only far-field approximations and avoids the multipole-to-local translations and SVD compression steps in the bbFMM [19]. With these choices the KITC is a relatively simple algorithm with low memory consumption, and this enables a straightforward OpenMP parallelization.

We present numerical results motivated by the method of regularized Stokeslets (MRS) for slow viscous flow [11, 12]. The MRS has been applied to simulate cilia- and flagella-driven flow [18, 49], helical swimming [12], slender body flow [8], coupled Stokes-Darcy flow [52], and flow around elastic rods [44]. Due to the complexity of the MRS kernels, they are prime candidates for kernel-independent fast summation methods, but as far as we know only recently has the KIFMM been applied to MRS simulations [47]. Here we apply the KITC to systems of regularized Stokeslets and rotlets from [47]. The results demonstrate the method’s good performance in terms of accuracy, efficiency, and memory consumption in serial and parallel simulations.

The paper is organized as follows. Section 2 discusses polynomial interpolation and its application to kernel approximation. Section 3 reviews barycentric Lagrange interpolation following [6, 55]. Section 4 explains how the interpolant is used to approximate particle-cluster interactions. Section 5 describes how some quantities called modified weights are computed. Section 6 presents the KITC algorithm. Section 7 reviews the MRS kernels (regularized Stokeslet and rotlet). Section 8 presents numerical results for two examples motivated by recent MRS simulations [47]. A summary is given in Section 9.

2 Polynomial interpolation and kernel approximation

We begin by recalling some basic facts about polynomial interpolation in 1D [55]. Given a function $f(t)$ and $n+1$ distinct points $s_{k}\in[-1,1]$ for $k=0,\ldots,n$, there is a unique polynomial $p_{n}(t)$ of degree at most $n$ that interpolates the function at these points, $p_{n}(s_{k})=f(s_{k}),k=0,\ldots,n$. The Lagrange form of the interpolating polynomial is

[TABLE]

where the Lagrange polynomials,

[TABLE]

have degree $n$ and satisfy $L_{k}(s_{j})=\delta_{jk}$. We view the interpolating polynomial $p_{n}(t)$ as an approximation to $f(t)$, and applying this idea to a kernel $k(x,y)$ in 1D, we hold $x$ fixed and interpolate with respect to $y$ to obtain

[TABLE]

The approximation on the right is a polynomial of degree $n$ in the variable $y$ and it interpolates the kernel at $y=s_{k}$; this idea is well known (e.g. [32], section 11.4).

Now consider a kernel $k({\bf x},{\bf y})$ in 3D and a tensor product set of grid points ${\bf s}_{\bf k}=(s_{k_{1}},s_{k_{2}},s_{k_{3}})\in[-1,1]^{3}$, where ${\bf k}=(k_{1},k_{2},k_{3})$ is a multi-index with $k_{\ell}=0,\ldots,n$ for $\ell=1,2,3$. As above we hold ${\bf x}$ fixed and interpolate with respect to ${\bf y}=(y_{1},y_{2},y_{3})$ to obtain

[TABLE]

In this case the approximation on the right is a polynomial of degree $n$ in each variable $(y_{1},y_{2},y_{3})$ and it interpolates the kernel at the grid points ${\bf y}={\bf s}_{\bf k}$. Moreover, the Lagrange polynomial expressions (5) and (6) are degenerate kernel approximations of the form (2).

3 Barycentric Lagrange interpolation

The expression for the interpolating polynomial described in Section 2 is not well-suited for practical computing due to cost and stability issues; the problem though is not with the Lagrange form (3) for $p_{n}(t)$, but rather with the expression (4) for the Lagrange polynomials $L_{k}(t)$. Berrut and Trefethen [6] advocated using instead the 2nd barycentric form of the Lagrange polynomials,

[TABLE]

where $w_{k}$ are the barycentric weights. This form is mathematically equivalent to (4), with the understanding that the removable singularity at $t=s_{j}$ is resolved by setting $L_{k}(s_{j})=\delta_{jk}$; to enforce this condition, following [6] the code is written so that if the argument $t$ is closer to an interpolation point $s_{j}$ than some tolerance, this is flagged and the correct value $L_{k}(s_{j})=\delta_{jk}$ is provided. In our implementation the tolerance is the minimum positive IEEE double precision floating point number, DBL_MIN = 2.22507e-308; further details will be given in Section 5.

We work with Chebyshev points of the 2nd kind,

[TABLE]

In this case the interpolating polynomial $p_{n}(t)$ converges rapidly and uniformly on $[-1,1]$ to $f(t)$ as $n$ increases, under mild smoothness assumptions on the given function [55]. Note that computing the barycentric weights $w_{k}$ by the definition (7) requires $O(n^{2})$ operations, but this expense disappears for the $s_{k}$ chosen in (8) because in that case the following simple weights can be used instead [48, 55],

[TABLE]

This relies on a scale-invariance property of the barycentric form of $L_{k}(t)$ in (7); namely, if the weights $w_{k}$ have a common constant factor $\alpha\neq 0$, then $\alpha$ can be cancelled from the numerator and denominator, and (7) stays the same.

Hence the barycentric Lagrange form of the interpolating polynomial is

[TABLE]

and with the simple weights (9), evaluating $p_{n}(t)$ requires $O(n)$ operations [6, 55].

The scale-invariance property is also important when working on different intervals; if $[-1,1]$ is linearly mapped to $[a,b]$ by $t\to\frac{1}{2}(a+b+t(b-a))$, then the weights $w_{k}$ defined in (7) gain a factor of $2^{n}/(b-a)^{n}$ which could lead to overflow or underflow, but this factor can be safely omitted due to scale-invariance. This means that the simple weights (9) can be used for any interval $[a,b]$, along with the linearly mapped Chebyshev points; this is important in the present work because the treecode uses intervals of different sizes. Note also that barycentric Lagrange interpolation is stable in finite precision arithmetic [46, 28, 42], and the Chebfun software package uses this form of polynomial interpolation [14].

For comparison, the bbFMM [19] uses Chebyshev points of the 1st kind,

[TABLE]

and a Chebyshev Lagrange form of the interpolating polynomial,

[TABLE]

where $T_{j}(t)$ is the $j$th degree Chebyshev polynomial and $\bar{L}_{k}(\bar{s}_{j})=\delta_{jk}$. The cost of evaluating $p_{n-1}(t)$ by directly summing (12) is $O(n^{2})$, or $O(n\log n)$ if a fast transform is used. We are not aware of an analog of the barycentric scale-invariance property.

4 Particle-cluster interactions

In a treecode the particles ${\{\bf x}_{i}\}$ are partitioned into a hierarchy of clusters $\{C\}$, and the sum (1) is written as

[TABLE]

where

[TABLE]

is the interaction between a target particle ${\bf x}_{i}$ and a source cluster $C=\{{\bf y}_{j}\}$. The sum over $C$ in (13) denotes a suitable subset of clusters depending on the target particle. Figure 1a depicts a particle-cluster interaction with target particle ${\bf x}_{i}$, cluster center ${\bf y}_{c}$, cluster radius $r$, and particle-cluster distance $R=|{\bf x}_{i}-{\bf y}_{c}|$. In this work the clusters are rectangular boxes whose sides are aligned with the coordinate axes. When the particle and cluster are well-separated (the criterion is given in Section 6), the interaction (14) is computed using the kernel approximation (6); this is depicted in Fig. 1b where the Chebyshev grid points ${\bf s}_{\bf k}$ are mapped to the cluster. Hence the target particle ${\bf x}_{i}$ interacts with the interpolation points ${\bf s}_{\bf k}$ rather than the source particless ${\bf y}_{j}$.

In more detail, we substitute the kernel approximation (6) into the particle-cluster interaction (14) and switch the order of summation to obtain the far-field particle-cluster approximation,

[TABLE]

where the modified weights are

[TABLE]

and ${\bf y}_{j}=(y_{j1},y_{j2},y_{j3})$. Note that when the target ${\bf x}_{i}$ is well-separated from the sources ${\bf y}_{j}$, the kernel $k({\bf x},{\bf y})$ is interpolated on a subdomain where it is smooth and this ensures the accuracy of the approximation.

We see that (15) has the form of a particle-cluster interaction (14), where the source particles and weights $\{{\bf y}_{j},f_{j}\}$ are replaced by the Chebyshev grid points and modified weights $\{{\bf s}_{k},\widehat{f}_{\bf k}\}$. This is advantageous for two reasons. First, the modified weights $\widehat{f}_{\bf k}$ for a given cluster are independent of the target particle ${\bf x}_{i}$, so they can be precomputed and re-used for different targets; details are given below. Second, the cost of the direct sum (14) is $O(N_{c})$, where $N_{c}$ is the number of particles in cluster $C$, while the cost of the approximation (15) is $O(n^{3})$ assuming the modified weights are known, so there is a cost reduction when $N_{c}\,{>>}\,n^{3}$. Finally note that the approximation (15) depends only on kernel evaluations and hence qualifies as kernel-independent.

5 Computation of modified weights

Algorithm 1 describes how the modified weights are computed. Before proceeding, note that the Lagrange polynomial terms in (16) in barycentric form are

[TABLE]

where the variables $a_{j,k,l}$ are introduced for clarity. There are three loops, the outer loop over source particles $j=1,\ldots,N_{c}$, the 1st inner loop over Chebyshev points $k=0,\ldots,n$, and the 2nd inner loop over coordinate indices $\ell=1,2,3$. To handle the removable singularity, we adapt the procedure in [6] (page 510).

In Algorithm 1, the source particles and weights ${\bf y}_{j},f_{j}$ for a given cluster are input, and the modified weights $\widehat{f}_{\bf k}$ are output. The modified weights are initialized to zero on line 4, the loop over source particles starts on line 6, and the flags and sums are initialized on line 7. Lines 9-13 loop over the Chebyshev points and coordinate indices to compute the terms $a_{j,k,\ell}$; this requires $O(nN_{c})$ operations and a temporary storage array of size $O(nN_{c})$ which is reused for different clusters. Line 10 checks whether a source particle coordinate is close to a Chebyshev point; if so, then a flag records the index. Then on lines 15-17 the code checks to see whether a flag was set; if so, then the temporary variables are adjusted to enforce the condition $L_{k}(s_{j})=\delta_{jk}$. Then on lines 20-22 the code loops over the tensor product of Chebyshev point indices and increments the modified weights $\widehat{f}_{\bf k}$ as indicated in (16); this requires $O(n^{3}N_{c})$ operations and $O(n^{3})$ storage for each cluster. The modified weights are computed this way for each cluster when the tree is constructed and in practice this requires only a small fraction of the total KITC CPU time; for example in one case with $N$ = 640 K particles, the total CPU time was 149 s, while computing the modified weights took less than 2 s.

The quantities $L_{k_{\ell}}(y_{j\ell})$ could be computed instead using the Chebyshev form (12), but the argument of the Chebyshev polynomials appearing there lies in the unit interval $[-1,1]$, so the particles ${\bf y}_{j}$ must be mapped to the unit cube $[-1,1]^{3}$, with cost $O(N_{c})$. The barycentric form (7) does not require this mapping; the particles stay in place in the cluster; hence the barycentric form has an advantage over the Chebyshev form in that it avoids the cost of mapping the particles ${\bf y}_{j}$.

6 Kernel-independent treecode algorithm

Aside from using barycentric Lagrange interpolation for the far-field approximation, the present algorithm is similar to previous treecodes based on analytic series expansions [4, 13, 35]. The procedure is outlined in Algorithm 2. After inputting the particle data and treecode parameters, a hierarchical tree of particle clusters is built and the modified weights for each cluster are computed. In our implementation the root cluster of the tree is the smallest rectangular box enclosing the particles. The root is bisected in any coordinate direction for which the side length is greater than $\ell_{max}/\sqrt{2}$, where $\ell_{max}$ is the maximum side length of the root. Hence the root is bisected in its long directions, but not its short directions, and the child clusters are bisected in the same way. The process continues until a cluster has fewer than $N_{0}$ particles, a user-specified parameter. The clusters obtained this way are rectangular boxes and a cluster may have 8, 4, or 2 children depending on which sides were bisected. The flexibility to use rectangular boxes instead of cubes is important in cases like Example 2 below where the particles lie in a rectangular slab.

Turning to Algorithm 2, after building the tree of clusters (line 5), the code cycles through the particles (line 7), and each particle interacts with clusters, starting at the root and proceeding to the child clusters. The multipole acceptance criterion (MAC),

[TABLE]

determines whether a particle ${\bf x}_{i}$ and cluster $C$ are well-separated, where (recall Fig. 1) $r$ is the cluster radius, $R$ is the particle-cluster distance, and $\theta$ is a user-specified parameter. If the MAC (18) is satisfied, the particle-cluster interaction is computed by barycentric Lagrange interpolation (15) with user-specified degree $n$; otherwise the code checks the child clusters, or if the cluster is a leaf (no children), then the interaction is computed directly by (14).

This is essentially the Barnes-Hut algorithm [4], extended to higher-order particle-cluster approximations computed by barycentric Lagrange interpolation. There are two stages; stage 1 is building the tree and computing the modified weights for each cluster, and stage 2 is computing the particle-cluster interactions. In principle both stages scale like $O(N\log N)$, although in practice stage 1 is much faster than stage 2. There are three user-specified parameters ($\theta,n,N_{0}$) that control the accuracy and efficiency of the treecode; the present work uses representative values with no claim that they are optimal.

7 Regularized Stokes kernels

The method of regularized Stokeslets (MRS) uses regularized kernels to represent point forces and torques in slow viscous flow [11, 12, 44]. In particular, consider the following functions defined on ${\mathbb{R}}^{3}$,

[TABLE]

where $r=|{\bf x}-{\bf y}|$, ${\bf x}$ is a target, ${\bf y}$ is a source, and $\epsilon$ is the MRS regularization parameter. Given a set of source particles $\{{\bf y}_{j}\}$ with forces $\{{\bf f}_{j}\}$ and torques $\{{\bf n}_{j}\}$, for $j=1,\ldots,N$, the velocity induced by regularized Stokeslets is

[TABLE]

where $r_{j}=|{\bf x}-{\bf y}_{j}|$, while the linear velocity and angular velocity induced by regularized Stokeslets and rotlets are

[TABLE]

These sums have the pairwise interaction form (1) suitable for a fast summation method; in this case the kernels are tensors, and the particle weights and output velocities are vectors. However the MRS functions (19) are somewhat complicated and we know of only one application [47] using a version of the KIFMM in which the equivalent densities are defined on coronas (or shells) around each cluster [61]. Note also that when $\epsilon\,{=}\,0$, the MRS kernels reduce to the usual Stokes kernels which are homogeneous (i.e. $k(\alpha{\bf x},\alpha{\bf y})=\alpha^{\lambda}k({\bf x},{\bf y})$ for some constant $\lambda$ and all $\alpha>0$); some versions of the FMM use this property to improve performance (e.g. [39, 50]), but this optimization is not available for the MRS kernels because they are non-homogeneous when $\epsilon\neq 0$. In the next section we demonstrate the capability of the KITC in evaluating the MRS sums (20)-(21) with $\epsilon\neq 0$.

8 Numerical results

We present results for two examples motivated by recent MRS simulations [47]. In these examples the targets and sources coincide, but this is not an essential restriction. To quantify the accuracy of the KITC we define the relative error,

[TABLE]

where ${\bf u}^{d}({\bf x}_{i})$ is the exact velocity computed by direct summation, ${\bf u}^{t}({\bf x}_{i})$ is the treecode approximation, and $|{\bf u}|$ is the Euclidean norm. In Example 2 where the linear velocity (21a) and angular velocity (21b) are computed, they are combined into a single vector and the error is computed as in (22). All lengths are nondimensional. References to “CPU time” are the total wall-clock run time in seconds.

The algorithm was coded in C++ in double precision and compiled using the Intel compiler icpc with $-$O2 optimization. The computations were performed on the University of Wisconsin-Milwaukee Mortimer Faculty Research Cluster; each node is a Dell PowerEdge R430 server with two 12-core Intel Xeon E5-2680 v3 processors at 2.50 GHz and 64 GB RAM; all compute and I/O nodes are linked by Mellanox FDR Infiniband (56Gb/s) and gigabit Ethernet networks. Parallel runs used OpenMP with up to 24 cores on a single node. The source code is available for download (github.com/Treecodes/stokes-treecode).

8.1 Example 1

The first example simulates microorganisms randomly located in a cube of side length $L=10$, where each microorganism is a pair of particles representing its body and flagella, and the particles exert unit forces in opposite directions along the organism length [47]. The total number of particles is $N$ and we consider five systems with $N\,{=}\,$10K, 80K, 640K, 5.12M, 40.96M. The maximum number of particles in a leaf was set to $N_{0}=2000$; note however that in this case each cluster has 8 children and the actual leaf size is approximately $N/8^{d}=1250$, where $d$ is the depth of the tree (some leaves may have a few more or less since the particle locations are random). The microorganism length is $\ell\,{=}\,0.02$ and the MRS parameter is $\epsilon=0.02$. In this example we compute the Stokeslet velocity (20).

8.1.1 CPU time versus error

Figure 2 focuses on intermediate system size $N\,{=}\,$640K, showing the treecode CPU time versus the error $E$, for MAC parameters in the range $0.4\leq\theta\leq 0.8$ and interpolation degrees $n=1,\ldots,10$ (increasing from right to left). For a given MAC parameter $\theta$, the error $E$ decreases as the degree $n$ increases, but the CPU time increases. The lower envelope of the data gives the most efficient treecode parameters. For example to achieve error $E\leq$ 1e-4, we can choose MAC $\theta=0.7$ and degree $n=7$. Figure 2 also shows the direct sum CPU time (2249 s) as a horizontal line (the value is independent of the error); the treecode is faster than direct summation for this range of parameters.

8.1.2 Error and CPU time versus system size

Figure 3 plots the treecode error $E$ (a), and CPU time (b) versus system size $N$, for MAC parameter $\theta=0.7$ and interpolation degree $n=1,3,5,7,9$. Since it is impractical to carry out the direct sum for the largest system with $N\,{=}\,$40.96M, in that case the error was computed using a sample of 2000 particles and the CPU time was extrapolated from smaller systems. In Fig. 3a for a given degree $n$, the error first increases with increasing system size $N$ and then it decreases slightly for the two largest systems; nonetheless for a given system size $N$, the error decreases as the interpolation degree $n$ increases.

Figure 3b shows that the treecode is faster than direct summation except for the smallest system with $N\,{=}\,$10K (those runs all take less than 1 s). Two reference lines are shown; the direct sum CPU time is parallel to the slope 2 line consistent with $O(N^{2})$ scaling, while the KITC time is slightly steeper than the slope 1 line consistent with $O(N\log N)$ scaling. As expected, for a given degree $n$ the speedup provided by the treecode increases with the system size $N$. To quantify this observation, Table 1 records the direct sum and treecode CPU time and the treecode error for the five systems using MAC parameter $\theta=0.7$ and degree $n\,{=}\,7$. For the largest system with $N\,{=}\,$40.96M particles, the treecode is more than 480 times faster than direct summation while achieving error $E\,{=}\,$7.06e-6.

8.1.3 Memory consumption

The treecode memory consumption has two main components which are estimated as follows. First, the memory due to the particle coordinates and weights is $O(N)$, where $N$ is the number of particles in the system; this is the same as in direct summation. Second, the memory due to the modified weights is $O(n^{3}N)$, where $n$ is the degree of polynomial interpolation and in this case $N$ represents the number of clusters in the tree (this assumes the number of levels in the tree is $O(\log N)$). Hence the theoretical memory consumption of the treecode is $O(n^{3}N)$, although in practice the $O(N)$ particle data dominates when $n$ is small.

To see the actual values, Table 2 presents the direct sum and treecode memory consumption (MB) for Example 1; the results were obtained using the Valgrind Massif tool (www.valgrind.org). For the largest system with $N\,{=}\,$40.96M particles, the direct sum memory consumption (3604.80MB) was obtained by extrapolating from the smaller systems. It can be verified that the memory consumption in Table 2 is consistent with the estimates in the previous paragraph; for a given degree $n$, the memory consumption scales like $O(N)$, while for a given system size $N$, beyond the baseline particle data, the additional memory consumption scales like $O(n^{3})$. Over this parameter range, the treecode uses less than 1.65 times as much memory as direct summation, and a KITC computation with $N\,{=}\,$40.96M particles and degree $n=9$ uses about 5.2GB of memory. For comparison, a bbFMM computation with 1M Stokeslets ran out of memory for degree $n>7$ ([19], Fig. 8), although the precise memory consumption value was not given.

8.1.4 Effect of MRS parameter

Table 3 shows the error $E$ for values of the MRS parameter in the range $0.005\leq\epsilon\leq 0.08$, with microorganism length $l=\epsilon$ and system size $N$ = 640K [47]. The KITC used MAC $\theta=0.7$ and degree $n=7$, and the error is on the order of 1e-5. As $\epsilon$ increases, the error first increases and then decreases; some error variation with $\epsilon$ is expected, but the reason for this particular variation is not clear. It can be noted that error variation with $\epsilon$ also occurred in KIFMM simulations ([47], Table 3).

8.1.5 Parallel simulations

We parallelized the direct sum and KITC using OpenMP with up to 24 cores on a single node. In both methods, the loop over target particles was parallelized, and the computation of the modified weights in the KITC was also parallelized.

We consider a system of $N\,{=}\,$640K regularized Stokeslets; the treecode parameters are $\theta=0.7,n=7,N_{0}=2000$, yielding $E\,{=}\,$3.17e-5. Table 4 shows the CPU time for direct sum and KITC $(d,t)$, parallel speedup and efficiency, and ratio of direct sum and treecode CPU time. Parallelizing the direct sum reduces the CPU time from 2249 s with 1 core to 108 s with 24 cores, yielding parallel efficiency $86\%$. Parallelizing the KITC reduces the CPU time from 149 s with 1 core to 10 s with 24 cores, yielding parallel efficiency $61\%$. The KITC is 15 times faster than direct summation with 1 core and 10 times faster with 24 cores.

For comparison, the KIFMM results in [47] employed the Matlab Parallel Computing Toolbox with up to 8 cores. For a system of size $N\,{=}\,$640K using a $4^{3}$ grid to solve for the equivalent densities, the error was $E\,{=}\,$7.33e-4, and the CPU time was reduced from 537 s with 1 core to 100 s with 8 cores, yielding parallel efficiency $67$% ([47], Tables 4 and 5). In comparing the KITC results and KIFMM results in [47], the point is not to claim that one method is better than the other; that would require more extensive testing beyond the scope of this work, but we believe the present results indicate that the KITC is a competitive option for fast summation of MRS kernels.

8.2 Example 2

The second example models an array of rods representing cilia or free-swimming flagella [44, 47]. The number of rods is $N_{r}$ and each rod is a helical curve with $M$ segments, so the total number of particles is $N=N_{r}(M+1)$. Each rod is parametrized by the $z$-coordinate and has the form $(x_{0}+0.3\cos 2z,y_{0}+0.3\sin 2z,z)$ for $0\leq z\leq 9$, where the base point $(x_{0},y_{0})$ lies at a regular Cartesian grid point in the $xy$-plane and the rods extend in the $z$-direction. Figure 4 shows an example with $N_{r}=16$ rods and $M=25$ segments on each rod.

In this example each particle is a superposition of a regularized Stokeslet and rotlet, and the KITC is applied to compute the linear velocity (21a) and angular velocity (21b) for a given rod configuration. Starting with $N_{r}=15^{2}$ rods with base points in the domain $[-8,8]\times[-8,8]$ as in [47], we increase the number of rods ($N_{r}=15^{2},20^{2},30^{2},40^{2},60^{2},80^{2}$) while expanding the horizontal dimensions of the domain to maintain constant rod density. Each rod has $M=150$ segments, the MRS parameter is $\epsilon=5L/M=0.3$, and each component of the force ${\bf f}_{j}$ and torque ${\bf n}_{j}$ is a random number in $[-1,1]$.

Unlike Example 1 where the particles lie in a cube, in Example 2 the particles lie in a rectangular slab (the $z$-direction is shorter than the $xy$-directions). In particular, the slab dimensions vary from $16^{2}\times 9$ for the smallest system $(N_{r}=15^{2},N=33975)$ to approximately $85^{2}\times 9$ for the largest system $(N_{r}=80^{2},N=966400)$. The tree construction scheme described above yields clusters that are well adapted to the slab geometry of this example.

Also note that the regularized Stokeslet/rotlet kernel (21) in Example 2 has more terms and is more expensive to evaluate than the regularized Stokeslet kernel (20) in Example 1. Hence to better balance the cost of direct particle-particle interactions (14) and particle-cluster approximations (15), the maximum leaf size was reduced to $N_{0}=1000$ in Example 2.

8.2.1 Error and CPU time versus system size

Figure 5 presents the treecode error (a) and the direct sum and treecode CPU time (b) versus system size $N$. The treecode MAC parameter is $\theta=0.7$ and the degree is $n=1,3,5,7,9$. The error increases slightly with system size $N$, but for a given degree $n$, the error amplitude is comparable to the results in Example 1 for this range of system size. As before the direct sum CPU time is parallel to the slope 2 line, while the KITC time is slightly steeper than the slope 1 line and is consistent with $O(N\log N)$ scaling.

8.2.2 Parallel simulations

Table 5 presents results for $N_{r}=80^{2}$ rods, system size $N=966400$, and with treecode parameters $\theta=0.7,n=7,N_{0}=1000$, yielding error $E$ = 2.24e-5. Parallelizing the direct sum reduces the CPU time from 15651 s with 1 core to 742 s with 24 cores, yielding parallel efficiency $87\%$, while parallelizing the KITC reduces the CPU time from 720 s with 1 core to 43 s with 24 cores, yielding parallel efficiency $68\%$. The KITC is 21 times faster than direct summation with 1 core and 17 times faster with 24 cores.

9 Summary

We presented a kernel-independent treecode (KITC) for fast summation of particle interactions in 3D. The method employs barycentric Lagrange interpolation at Chebyshev points to approximate well-separated particle-cluster interactions. The KITC requires only kernel evaluations, is suitable for non-oscillatory kernels, and it utilizes on a scale-invariance property of barycentric Lagrange interpolation. Numerical results were presented for the non-homogeneous kernels arising in the Method of Regularized Stokeslets (MRS) [47]. For a given level of accuracy, the treecode CPU time scales like $O(N\log N)$, where $N$ is the number of particles, and a substantial speedup over direct summation is achieved for large systems. The KITC is a relatively simple algorithm with low memory consumption, and this enables a straightforward parallelization; here we employed OpenMP with up to 24 cores on a single node; alternative approaches including distributed memory parallelization for larger systems will be considered in the future [17, 43, 58].

Other kernels for which the KITC may be suitable are the regularized Green’s functions used in computing nearly singular integrals (e.g. [5, 53, 54]), the RPY (Rotne-Prager-Yamakawa) tensor for hydrodynamic interactions [36], and the generalized Born potential for implicit solvent modeling [59]. We expect that the KITC can take advantage of several techniques employed in advanced implementations of the KIFMM and bbFMM such as blocking operations, utilizing BLAS routines, AVX and SSE vectorization, and GPU and Phi coprocessing (e.g. [50, 1, 33, 39]). It should be noted that an FMM could be implemented using barycentric Lagrange interpolation for both the target and source variables, and this is also an interesting direction for future study. Finally we mention that the extension to barycentric Hermite interpolation has been carried out for scalar electrostatic kernels [31, 56].

Acknowledgments

We are grateful for a computer allocation on the University of Wisconsin-Milwaukee Mortimer Faculty Research Cluster. We thank Nathan Vaughn and Leighton Wilson for help with programming. Partial support is acknowledged from start-up funds provided by the University of Wisconsin-Milwaukee, NSF grant DMS-1819094, the Michigan Institute for Computational Discovery and Engineering (MICDE), and the Mcubed program at the University of Michigan.