3
boundary there is a relatively simple image construction for Gdue to Blake [27]. This image
construction has also enabled the development of fast methods for evaluating (2) for a single wall,
based on either (flexible periodicity) Fast Multipole Method [28,29] or the Fast Fourier Transform
(FFT) [30,31]. For triply-periodic (TP) domains, the Positively Split Ewald (PSE) method [26]
for evaluating (2) (and also generating Brownian velocities) for the RPY kernel provides the basis
for Fast Stokesian Dynamics [15]. However, it remains a challenge to construct a similarly-efficient
method for confined suspensions since PSE uses Fourier representations in a key way for each
component of the method (Stokes solver, generating Brownian increments, and Ewald splitting).
In principle, a greater flexibility with respect to boundary conditions can be achieved by re-
placing analytical Green’s functions with a grid-based Stokes solver [32]. This requires replacing
the singular surface delta function form of δa(r) with a smooth function that can be resolved on a
grid. In the FCM, a Gaussian envelope function δais used, which allows for analytical calculation
of Oin an unbounded domain [6,7]; in numerical methods the Gaussian is truncated. By con-
trast, in the IBM, δais a discrete grid function specifically constructed to maximize grid invariance
[5,12,17], and the double convolution in (1) is discrete. This makes all IBM results grid- and
solver-dependent, without a direct continuum limit.
We develop a method that combines favorable features of FCM and IBM. Namely, we maintain
the continuum representation (1) from FCM, however, we do not use a Gaussian envelope but
rather use the “exponential of a semi-circle” (ES) kernel proposed by Barnett for its efficiency
in the context of non-uniform FFTs [11]. Not only does this kernel allow for greater flexibility in
tuning the hydrodynamic radius of the particles [33], but also allows us to use many fewer grid cells
per particle than for a Gaussian kernel, comparable to the IBM, while still solving the continuum
equations to several digits of accuracy. Specifically, for about two-three digits of accuracy we need
four grid cells per particle (in each dimension) if only translational velocity is required, and five
or six if rotational velocities are also required. The grid independence of the results makes them
transferable, and allows us to separate the Stokes solver from the physics (i.e., the problem has a
solution without a discretization). We find that Ocomputed with the ES and Gaussian kernels
are indistinguishable in practice. We employ an image construction following Yeo and Maxey [8]
(also used in the IBM [17,34]) to generalize (1) to the relevant-in-practice case when some of the
blob kernels overlap the boundaries; this is considerably harder to do for the RPY kernel, and has
required ad hoc fixes in past work by some of us [35].
For TP domains one can easily solve the Stokes equations spectrally in Fourier space, but this
is considerably harder when boundaries are present. In Sec. III we use some of the ideas applied to
the Poisson equation in [13] to develop a Stokes solver for DP geometries that are unbounded or
confined by one or two boundaries in the aperiodic direction. Specifically, we focus on bottom wall
(BW) DP geometries, with a single wall at z= 0, and slit channel (SC) DP geometry, confined by
two walls at z= 0 and z=H, but the Stokes solver can handle more flexible boundary conditions
in the aperiodic direction. Our novel fluid solver makes it possible to handle geometries that are
partially unbounded in one direction, unlike in existing IBM or FCM implementations based on
more traditional grid solvers like finite differences or finite elements [36]. Our fluid solver has the
additional advantage that its implementation requires only calls to the three-dimensional (3D) FFT
with an oversampling factor of 2 in the aperiodic direction, combined with trivially parallelizable
one-dimensional pentadiagonal boundary value solvers in the zdirection. We use this to implement
the method on Graphical Processing Units, achieving linear scaling up to as many as one million
particles, with scaling constants much better than existing fast methods [29] for the types of
problems we study here. It is important to note that Srinivasan and Tornberg have recently
developed a sophisticated method based on 3D FFTs for the singular Stokes Green’s function with
a bottom wall geometry (using images) and flexible periodicity [30,31]. Our approach is different
and specialized to doubly-periodic geometries, is conceptually simpler, and affords flexibility in the