2.1 Saddle point preconditioner

The saddle point matrices described in § 2 are efficiently solved using iterative techniques such as the generalized minimum residuals method (GMRES) (Saad & Schultz Reference Saad and Schultz1986). Preconditioners are commonly applied to improve the convergence of iterative matrix solvers, particularly in applications where the matrix is ill-conditioned. The main idea behind preconditioning is to apply a linear transformation ${\mathcal{P}}^{-1}\approx {\mathcal{A}}^{-1}$ to the matrix equation ${\mathcal{A}}\boldsymbol{\cdot }\boldsymbol{x}=\boldsymbol{b}$ , so that ${\mathcal{P}}^{-1}\boldsymbol{\cdot }{\mathcal{A}}$ has a smaller condition number than  ${\mathcal{A}}$ , and whose solution therefore converges more quickly. Saddle point matrices commonly arise in the formulation of engineering and mathematical problems with constraints, and their preconditioning and solution are well characterized (Benzi et al. Reference Benzi, Golub and Liesen2005). To develop a preconditioner, we start with the block triangular factorization of the saddle point matrix

(2.10) $$\begin{eqnarray}{\mathcal{A}}=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644} & \mathbf{0}\\ {\mathcal{B}}^{\text{T}}\boldsymbol{\cdot }{\mathcal{M}}^{-1} & \unicode[STIX]{x1D644}\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}cc@{}}{\mathcal{M}} & \mathbf{0}\\ \mathbf{0} & {\mathcal{S}}\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644} & {\mathcal{M}}^{-1}\boldsymbol{\cdot }{\mathcal{B}}\\ \mathbf{0} & \unicode[STIX]{x1D644}\end{array}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity matrix. A good preconditioner provides a reasonably accurate approximation of the actual matrix while remaining inexpensive to compute. In the preconditioner, we set ${\mathcal{M}}=\unicode[STIX]{x1D701}^{-1}\unicode[STIX]{x1D644}$ , where $\unicode[STIX]{x1D701}$ is an effective drag coefficient, because ${\mathcal{M}}$ is dense and ill-conditioned, so its inverse is challenging to compute. The projector ${\mathcal{B}}$ is generally straightforward to compute so it is used in the preconditioner without approximation. To apply the preconditioner, it is necessary to invert the approximation for  ${\mathcal{A}}$ , so ${\mathcal{S}}$ is approximated by ${\mathcal{L}}\boldsymbol{\cdot }{\mathcal{L}}^{\text{T}}$ , where ${\mathcal{L}}$ is the zero fill-in incomplete lower Cholesky factor of $\unicode[STIX]{x1D64B}\boldsymbol{\cdot }{\mathcal{S}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}^{\text{T}}$ or a reasonable approximation to it, and $\unicode[STIX]{x1D64B}$ is the unitary permutation matrix representing the reverse Cuthill–McKee (RCM) reordering (Cuthill & McKee Reference Cuthill and McKee1969; George & Liu Reference George and Liu1981) of $\unicode[STIX]{x1D64D}_{FU}^{nf}$ . The RCM reordering reduces the sparse matrix bandwidth and is applied in this case to reduce the number of dropped entries in the incomplete Cholesky factor. This improves the stability of the method because the incomplete Cholesky factorization is not guaranteed to exist even for symmetric positive-definite matrices. The approximation for ${\mathcal{S}}$ used in the preconditioner is

(2.11) $$\begin{eqnarray}\tilde{{\mathcal{S}}}=\unicode[STIX]{x1D701}\unicode[STIX]{x1D644}+\tilde{\unicode[STIX]{x1D64D}}_{FU}^{nf},\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D64D}}_{FU}^{nf}$ is an approximation to $\unicode[STIX]{x1D64D}_{FU}^{nf}$ built from particle pairs whose separation is less than $r_{p}<r_{cut}$ , and $\unicode[STIX]{x1D701}$ is an approximation for the drag coefficient. For monodisperse spheres one might choose the Stokes drag $\unicode[STIX]{x1D701}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}a$ , but Viera (Reference Viera2002) notes that larger values can improve performance in various applications. Here we will simply use the Stokes drag coefficient. This approximation increases the sparsity of $\unicode[STIX]{x1D64D}_{FU}^{nf}$ , reducing memory requirements and accelerating computation, while simultaneously remaining an effective preconditioner. We will show in § 3 that $\tilde{\unicode[STIX]{x1D64D}}_{FU}^{nf}$ is an effective preconditioner because it retains the largest values from $\unicode[STIX]{x1D64D}_{FU}^{nf}$ , which originate with the most narrowly separated particles, and that rather small values of $r_{p}$ are sufficient to accelerate convergence. The preconditioner is

(2.12) $$\begin{eqnarray}{\mathcal{P}}=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644} & \mathbf{0}\\ {\mathcal{B}}^{\text{T}} & \unicode[STIX]{x1D644}\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D701}^{-1}\unicode[STIX]{x1D644} & \mathbf{0}\\ \boldsymbol{ 0} & {\mathcal{L}}\boldsymbol{\cdot }{\mathcal{L}}^{\text{T}}\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644} & {\mathcal{B}}\\ \mathbf{0} & \unicode[STIX]{x1D644}\end{array}\right],\end{eqnarray}$$

and its inverse

(2.13) $$\begin{eqnarray}{\mathcal{P}}^{-1}=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644} & -{\mathcal{B}}\\ \mathbf{0} & \unicode[STIX]{x1D644}\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D701}\unicode[STIX]{x1D644} & \mathbf{0}\\ \mathbf{0} & {\mathcal{L}}^{-\text{T}}\boldsymbol{\cdot }{\mathcal{L}}^{-1}\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644} & \mathbf{0}\\ -{\mathcal{B}}^{\text{T}} & \unicode[STIX]{x1D644}\end{array}\right],\end{eqnarray}$$

can be quickly and inexpensively applied to a vector using the triangular solves of ${\mathcal{L}}$ and the simple projector  ${\mathcal{B}}$ .

2.2 Brownian dynamics in the saddle point formulation

In the SD framework, the stochastic force $\boldsymbol{F}^{B}$ is also decomposed into near-field and far-field parts (Banchio & Brady Reference Banchio and Brady2003)

(2.14) $$\begin{eqnarray}\boldsymbol{F}^{B}=\boldsymbol{F}_{nf}^{B}+\boldsymbol{F}_{ff}^{B},\end{eqnarray}$$

where

(2.15) $$\begin{eqnarray}\boldsymbol{F}_{nf}^{B}=\sqrt{\frac{2k_{B}T}{\unicode[STIX]{x0394}t}}(\unicode[STIX]{x1D64D}_{FU}^{nf})^{1/2}\boldsymbol{\cdot }\unicode[STIX]{x1D733}_{nf}\end{eqnarray}$$

and

(2.16) $$\begin{eqnarray}\boldsymbol{F}_{ff}^{B}=\sqrt{\frac{2k_{B}T}{\unicode[STIX]{x0394}t}}({\mathcal{B}}^{\text{T}}\boldsymbol{\cdot }{\mathcal{M}}^{-1}\boldsymbol{\cdot }{\mathcal{B}})^{1/2}\boldsymbol{\cdot }\unicode[STIX]{x1D733}_{ff},\end{eqnarray}$$

where $\unicode[STIX]{x1D733}_{nf}$ and $\unicode[STIX]{x1D733}_{ff}$ are vectors of standard independent Gaussian random variables. The stochastic quantities can also be represented as stochastic variables with zero mean and a covariance set by the fluctuation dissipation theorem, namely

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{F}_{nf}^{B}(\text{i}\unicode[STIX]{x0394}t)\boldsymbol{F}_{nf}^{B}(j\unicode[STIX]{x0394}t)\rangle =\frac{2k_{B}T}{\unicode[STIX]{x0394}t}\unicode[STIX]{x1D64D}_{FU}^{nf}\unicode[STIX]{x1D6FF}_{ij}, & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{F}_{ff}^{B}(\text{i}\unicode[STIX]{x0394}t)\boldsymbol{F}_{ff}^{B}(j\unicode[STIX]{x0394}t)\rangle =\frac{2k_{B}T}{\unicode[STIX]{x0394}t}{\mathcal{B}}^{\text{T}}\boldsymbol{\cdot }{\mathcal{M}}^{-1}\boldsymbol{\cdot }{\mathcal{B}}\unicode[STIX]{x1D6FF}_{ij}, & \displaystyle\end{eqnarray}$$

where the angle brackets denote an ensemble average over realizations of the stochastic variables and $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta. In this way, computing $\boldsymbol{F}_{nf}^{B}$ and $\boldsymbol{F}_{ff}^{B}$ can be thought of as generating samples of an instantaneously correlated random vector whose covariance is proportional to a specific matrix. The near-field Brownian forces in (2.15) are reasonably straightforward to compute because the near-field resistance tensor is known and can be directly computed. The far-field Brownian forces, however, are more challenging because the far-field interactions are described in terms of the mobility, whose inverse is therefore required to compute the forces, as shown in (2.16). Banchio & Brady (Reference Banchio and Brady2003) compute these forces directly using the Chebyshev polynomial approximation developed by Fixman (Reference Fixman1986), which is computationally expensive because ${\mathcal{M}}$ is ill-conditioned. The saddle point formulation enables a more efficient approach, namely the computation of a far-field slip velocity

(2.19) $$\begin{eqnarray}{\mathcal{U}}^{B}=\sqrt{\frac{2k_{B}T}{\unicode[STIX]{x0394}t}}{\mathcal{M}}^{1/2}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D733}}_{ff},\end{eqnarray}$$

that emerges by rearranging (2.8) to give

(2.20) $$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}{\mathcal{M}} & {\mathcal{B}}\\ {\mathcal{B}}^{\text{T}} & -\unicode[STIX]{x1D64D}_{FU}^{nf}\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}c@{}}\hat{{\mathcal{F}}}\\ \boldsymbol{U}\end{array}\right]=\left[\begin{array}{@{}c@{}}{\mathcal{U}}^{B}\\ -(\boldsymbol{F}^{P}+\boldsymbol{F}_{nf}^{B})\end{array}\right],\end{eqnarray}$$

where $\hat{{\mathcal{F}}}={\mathcal{F}}+{\mathcal{M}}^{-1}\boldsymbol{\cdot }{\mathcal{U}}^{B}$ is the set of far-field hydrodynamic forces plus an additional force required to satisfy the no-slip boundary condition on the particle surfaces and $\hat{\unicode[STIX]{x1D733}}_{ff}$ is another white noise vector. Introducing ${\mathcal{U}}^{B}$ in this way is physically equivalent to imposing a fluctuating stress in the far-field model of the fluid, whose associated particle displacements satisfy the fluctuation–dissipation theorem. For many representations of  ${\mathcal{M}}$ , ${\mathcal{U}}^{B}$ is far easier to compute than $\boldsymbol{F}_{ff}^{B}$ , and accelerated techniques with $O(N)$ or nearly $O(N)$ scaling have been developed for various approximations to  ${\mathcal{M}}$ . Examples include the fluctuating immersed boundary method published by Delong et al. (Reference Delong, Usabiaga, Delgado-Buscalioni, Griffith and Donev2014), the fluctuating force coupling method of Delmotte & Keaveny (Reference Delmotte and Keaveny2015) and our own positively split Ewald approach (Fiore et al. Reference Fiore, Usabiaga, Donev and Swan2017; Fiore & Swan Reference Fiore and Swan2018).

Others have also used the concept of a far-field slip velocity to efficiently model the thermal motion of dispersed particles. For example, Delmotte & Keaveny (Reference Delmotte and Keaveny2015) considered spherical particles with force, torques and stresslets, and introduced a slip velocity that included a non-zero rate of strain on each particle. Sprinkle et al. (Reference Sprinkle, Balboa Usabiaga, Patankar and Donev2017) considered rigid multi-blob particles and introduced a fluid slip velocity in an approach that is a special case of the one presented here. In particular, equation (6) in their paper is the special case of (2.20) in which lubrication interactions are neglected and ${\mathcal{B}}$ represents their geometric block matrix, $-{\mathcal{K}}$ , describing rigid body motions. While the saddle point matrix equation has been previously used by Sprinkle et al. (Reference Sprinkle, Balboa Usabiaga, Patankar and Donev2017) and others in the context of hydrodynamic mobility problems, equation (2.20) is the first time these ideas have been applied to SD or resistance-type problems, in which fluctuating and deterministic near-field forces also need be considered. Equation (2.20) represents a completely general framework for SD simulations of arbitrary particle shapes with relatively few restrictions on the model imposed for near-field HI. No prescription on the particle sizes or shapes has been made, and the method is valid for any symmetric positive semi-definite approximations for ${\mathcal{M}}$ and $\unicode[STIX]{x1D64D}_{FU}^{nf}$ . In §§ 3 and 4 we discuss techniques specific to spherical particles in periodic geometries, both for computing the deterministic HI and the stochastic forces and slip velocities.

2.3 Stochastic sampling preconditioning

As mentioned in the previous section, accelerated sampling methods for the effective far-field Brownian force can be developed by leveraging the algebraic structure of various approximations for ${\mathcal{M}}$ to generate a far-field slip velocity which is then projected onto the set of hydrodynamic constraints. Further discussion of such methods is deferred to § 4. In this section, we describe an efficient and completely general approach to sample the near-field Brownian forces. This approach relies on the iterative Krylov subspace approximations of Chow & Saad (Reference Chow and Saad2014) and leverages the preconditioner that was already developed in § 2.1 for the saddle point problem.

The preconditioning is more simply described using an alternative, but equivalent, notation for the stochastic variables, where an arbitrary random vector $\boldsymbol{x}$ drawn from a normal distribution with mean $\langle \boldsymbol{x}(\text{i}\unicode[STIX]{x0394}t)\rangle =\unicode[STIX]{x1D741}$ and covariance $\langle \boldsymbol{x}(\text{i}\unicode[STIX]{x0394}t)\boldsymbol{x}^{\text{T}}(j\unicode[STIX]{x0394}t)\rangle =\boldsymbol{V}\unicode[STIX]{x1D6FF}_{ij}$ is denoted as

(2.21) $$\begin{eqnarray}\boldsymbol{x}\sim {\mathcal{N}}(\unicode[STIX]{x1D741},\boldsymbol{V}).\end{eqnarray}$$

Standard random number generators are only able to draw random samples with a covariance proportional to  $\unicode[STIX]{x1D644}$ ,

(2.22) $$\begin{eqnarray}\boldsymbol{n}\sim {\mathcal{N}}(\unicode[STIX]{x1D741},\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D644}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ is a positive scalar. Sampling a distribution with arbitrary covariance $\boldsymbol{V}$ is achieved by computing the product of the square root of the covariance matrix with the independent and identically distributed vector $\boldsymbol{n}$ computed the random number generator

(2.23) $$\begin{eqnarray}\boldsymbol{x}=\unicode[STIX]{x1D741}+\boldsymbol{V}^{1/2}\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$

If $\boldsymbol{V}$ is known explicitly, direct calculation of $\boldsymbol{V}^{1/2}$ can be evaluated with eigenvalue or singular value decomposition, incurring a computational cost of $O(N^{3})$ . Iterative Krylov subspace approximations for $\boldsymbol{V}^{1/2}\boldsymbol{\cdot }\boldsymbol{n}$ , such as the commonly used Chebyshev polynomial approach of Fixman (Reference Fixman1986) and the more recently developed Lanczos method of Chow & Saad (Reference Chow and Saad2014), can be applied for an arbitrary  $\boldsymbol{V}$ . While these iterative schemes may converge even for ill-conditioned matrices, they are especially efficient when $\boldsymbol{V}$ is well-conditioned and $\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{x}$ can be rapidly computed. Preconditioning can be used to improve the condition number and accelerate convergence. We will show how this can be achieved to sample  $\unicode[STIX]{x1D64D}_{FU}^{nf}$ .

Important to any sampling approach is that $\unicode[STIX]{x1D64D}_{FU}^{nf}$ is positive semi-definite, which guarantees that the resulting forces are real-valued. Preconditioning is necessary because $\unicode[STIX]{x1D64D}_{FU}^{nf}$ is ill-conditioned. The large condition number arises because it has both very large eigenvalues originating in the divergent forces between nearly touching particles and very small eigenvalues for particles with few neighbours whose separations are close to the near-field cutoff radius, and whose entries in the resistance tensor are nearly zero. The large and small eigenvalues are most effectively preconditioned independently.

The preconditioner developed in § 2.1 is primarily designed to condition the large eigenvalues present in $\unicode[STIX]{x1D64D}_{FU}^{nf}$ , so it can be recycled here with some careful modifications that ensure the resulting variables have the proper covariance. Here, the relevant piece of the saddle point preconditioner is the incomplete Cholesky factor ${\mathcal{L}}$ of the matrix $\unicode[STIX]{x1D64B}\boldsymbol{\cdot }(\unicode[STIX]{x1D64D}_{FU}^{nf}+\unicode[STIX]{x1D701}\unicode[STIX]{x1D644})\boldsymbol{\cdot }\unicode[STIX]{x1D64B}^{\text{T}}$ , which is applied to $\unicode[STIX]{x1D64D}_{FU}^{nf}$ symmetrically in order to reduce the magnitude of the large eigenvalues,

(2.24) $$\begin{eqnarray}{\mathcal{L}}^{-1}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}\boldsymbol{\cdot }\unicode[STIX]{x1D64D}_{FU}^{nf}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}^{\text{T}}\boldsymbol{\cdot }{\mathcal{L}}^{-\text{T}}.\end{eqnarray}$$

To precondition the small eigenvalues, we use a diagonal matrix that is designed to scale each row and column of $\unicode[STIX]{x1D64D}_{FU}^{nf}$ with suitably small non-zero diagonals, $0<(\unicode[STIX]{x1D64D}_{FU}^{nf})_{ii}<\unicode[STIX]{x1D701}$ , by the diagonal element, so that the relative magnitude of the entries within each row and column is preserved but the absolute magnitudes cluster around  $\unicode[STIX]{x1D701}$ . This preconditioner does not alter the rows and columns of $\unicode[STIX]{x1D64D}_{FU}^{nf}$ with diagonals that are either zero or greater than or equal to a nominal drag coefficient, say $\unicode[STIX]{x1D701}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}a$ for mono-disperse spherical particles, i.e.

(2.25) $$\begin{eqnarray}\unicode[STIX]{x1D60B}_{ii}=\left\{\begin{array}{@{}ll@{}}(\unicode[STIX]{x1D64D}_{FU}^{nf})_{ii}^{1/2}, & 0<(\unicode[STIX]{x1D64D}_{FU}^{nf})_{ii}<6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}a,\\ 1, & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{ii}=0,\\ 1, & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{ii}\geqslant 6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}a.\end{array}\right.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D63F}$ is primarily expected to be important at low particle volume fractions, for which particles are more likely to be well separated and have only a few neighbours within the near-field cutoff radius. Combining this preconditioner with the incomplete Cholesky preconditioner for large eigenvalues yields the preconditioned matrix

(2.26) $$\begin{eqnarray}{\mathcal{L}}^{-1}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}^{-1}\boldsymbol{\cdot }\unicode[STIX]{x1D64D}_{FU}^{nf}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}^{-\text{T}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}^{\text{T}}\boldsymbol{\cdot }{\mathcal{L}}^{-\text{T}}.\end{eqnarray}$$

Unfortunately, this preconditioner is not a great approximation for $\unicode[STIX]{x1D64D}_{FU}^{nf}$ because $\unicode[STIX]{x1D64D}_{FU}^{nf}$ only has non-zero rows for particles that are closely separated, whereas the addition of $\unicode[STIX]{x1D701}\boldsymbol{I}$ in the computation of ${\mathcal{L}}$ guarantees that ${\mathcal{L}}$ has no rows or columns that are purely zero. The effect of the preconditioner can be improved by instead sampling forces from the distribution

(2.27) $$\begin{eqnarray}\hat{\boldsymbol{F}}_{nf}^{B}\sim {\mathcal{N}}\left(\mathbf{0},\frac{2k_{B}T}{\unicode[STIX]{x0394}t}(\unicode[STIX]{x1D64D}_{FU}^{nf}+\unicode[STIX]{x1D701}\unicode[STIX]{x1D644}_{nn})\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D644}_{nn}$ is a modified identity tensor with ones where the diagonal elements of $\unicode[STIX]{x1D64D}_{FU}^{nf}$ are zero, i.e. for the rows/columns associated with neighbourless particles, and zero otherwise, i.e.

(2.28) $$\begin{eqnarray}(\unicode[STIX]{x1D644}_{nn})_{ii}=\left\{\begin{array}{@{}ll@{}}1, & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{ii}=0,\\ 0, & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{ii}\neq 0.\end{array}\right.\end{eqnarray}$$

Neighbourless particles experience no near-field Brownian force, so the proper force distribution can be obtained by projecting the corresponding elements of $\boldsymbol{F}_{nf}^{B}$ to zero

(2.29) $$\begin{eqnarray}\boldsymbol{F}_{nf}^{B}=(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D644}_{nn})\boldsymbol{\cdot }\hat{\boldsymbol{F}}_{nf}^{B}\sim {\mathcal{N}}\left(\mathbf{0},\frac{2k_{B}T}{\unicode[STIX]{x0394}t}\unicode[STIX]{x1D64D}_{FU}^{nf}\right).\end{eqnarray}$$

The operator $(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D644}_{nn})$ is a projector that removes the additional noise introduced by adding $\unicode[STIX]{x1D701}\unicode[STIX]{x1D644}_{nn}$ to $\unicode[STIX]{x1D64D}_{FU}^{nf}$ . This operation works because we project out a subset of the Brownian forces. The preconditioned sampling produces random vectors

(2.30) $$\begin{eqnarray}\boldsymbol{z}\sim {\mathcal{N}}\left(\mathbf{0},\frac{2k_{B}T}{\unicode[STIX]{x0394}t}{\mathcal{L}}^{-1}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}^{-1}\boldsymbol{\cdot }(\unicode[STIX]{x1D64D}_{FU}^{nf}+\unicode[STIX]{x1D701}\unicode[STIX]{x1D644}_{nn})\boldsymbol{\cdot }\unicode[STIX]{x1D63F}^{-\text{T}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}^{\text{T}}\boldsymbol{\cdot }{\mathcal{L}}^{-\text{T}}\right),\end{eqnarray}$$

where

(2.31) $$\begin{eqnarray}\boldsymbol{F}_{nf}^{B}=(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D644}_{nn})\boldsymbol{\cdot }\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}^{\text{T}}\boldsymbol{\cdot }{\mathcal{L}}\boldsymbol{\cdot }\boldsymbol{z}\sim {\mathcal{N}}\left(\mathbf{0},\frac{2k_{B}T}{\unicode[STIX]{x0394}t}\unicode[STIX]{x1D64D}_{FU}^{nf}\right).\end{eqnarray}$$

The use of $\unicode[STIX]{x1D644}_{nn}$ and $\unicode[STIX]{x1D644}-\unicode[STIX]{x1D644}_{nn}$ bears some additional discussion. As mentioned, adding $\unicode[STIX]{x1D644}_{nn}$ to $\unicode[STIX]{x1D64D}_{FU}^{nf}$ improves the efficacy of the incomplete Cholesky preconditioner but increases the variance of the stochastic forces. Fortunately, particles that have no neighbours in the near field cannot experience a near-field Brownian force; the operator $\unicode[STIX]{x1D644}-\unicode[STIX]{x1D644}_{nn}$ projects these forces out of $\boldsymbol{F}_{nf}^{B}$ and restores the expected covariance. This approach only works because those forces should be exactly zero. Additional preconditioning may be possible for specific forms of  $\unicode[STIX]{x1D64D}_{FU}^{nf}$ , where knowledge of the specific interactions can be used to further compress the eigenspectrum and improve convergence. The efficacy of the preconditioner presented in this section is demonstrated in § 4 for the Stokesian dynamics lubrication approximation for the near-field resistance.

It is important to note that the approach advocated in this section is not the same as simply omitting the particles that have no neighbours. Rather, Brownian forces are computed for all particles then projected onto the set of particles with neighbours. The motivation behind this approach is to leverage the preconditioner already developed for the saddle point problem in § 2.1, which effectively preconditions many of the features that contribute to large and small eigenvalues in the covariance matrix of the stochastic sampling problem, so as to avoid repeating those calculations. In order to use the saddle point preconditioner, which considers all particles and not solely the ones with neighbours, the Brownian sampling method must at some stage generate noise for these particles without neighbours.

3.1 Far-field hydrodynamics for monodisperse spheres

The far-field mobility tensor ${\mathcal{M}}$ is constructed by multipole expansion. We have previously developed a compact representation of ${\mathcal{M}}$ (Fiore et al. Reference Fiore, Usabiaga, Donev and Swan2017; Fiore & Swan Reference Fiore and Swan2018), which will be summarized here for completeness. Equivalent representations have also been developed by Ladd (Reference Ladd1990) and Mizerski et al. (Reference Mizerski, Wajnryb, Zuk and Szymczak2014). This representation is derived from the integral formulation of the Stokes equations for rigid suspended bodies (Ladyzhenskaya & Silverman Reference Ladyzhenskaya and Silverman1969)

(3.1) $$\begin{eqnarray}u_{i}(\boldsymbol{x})=\int _{S_{\unicode[STIX]{x1D6FD}}}\text{d}S_{\unicode[STIX]{x1D6FD}}(\,\boldsymbol{y})\unicode[STIX]{x1D611}_{im}(\boldsymbol{x},\boldsymbol{y})f_{m}(\,\boldsymbol{y}),\end{eqnarray}$$

where $u_{i}(\boldsymbol{x})$ is the $i$ th component of the fluid velocity at a point $\boldsymbol{x}$ in the fluid propagated by the force density $\boldsymbol{f}(\,\boldsymbol{y})$ at point $\boldsymbol{y}$ on the particle surface $S_{\unicode[STIX]{x1D6FD}}$ through the Green’s function $\unicode[STIX]{x1D645}$ . In this and following equations, Einstein summation notation is assumed for repeated indices. All particles are assumed to be spherical and the first few multipole moments are used to approximate the force density on the surface of each particle

(3.2) $$\begin{eqnarray}f_{j}=\frac{1}{4\unicode[STIX]{x03C0}a^{2}}F_{j}+\frac{3}{4\unicode[STIX]{x03C0}a^{3}}\hat{n}_{n}C_{nj}+\cdots ,\end{eqnarray}$$

where $\hat{\boldsymbol{n}}$ is the normal to the particle surface, $a$ is the particle radius, $\boldsymbol{F}$ is the total force on the particle and $\boldsymbol{C}$ is the traceless couplet, whose symmetric and antisymmetric components are the stresslet and torque, respectively. The multipole expansion is truncated after the couplet for computational efficiency, but can in principle be extended to arbitrary order. The force and couplet are defined by

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle F_{m}=\int _{S_{\unicode[STIX]{x1D6FD}}}\,\text{d}S_{\unicode[STIX]{x1D6FD}}\,f_{m}, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle C_{mn}=\int _{S_{\unicode[STIX]{x1D6FD}}}\,\text{d}S_{\unicode[STIX]{x1D6FD}}\left[(y_{m}-x_{\unicode[STIX]{x1D6FD},m})f_{n}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{mn}(y_{p}-x_{\unicode[STIX]{x1D6FD},p})f_{p}\right], & \displaystyle\end{eqnarray}$$

where $\boldsymbol{x}_{\unicode[STIX]{x1D6FD}}$ is the location of the centre of particle $\unicode[STIX]{x1D6FD}$ , and higher-order moments are defined by higher-order outer products of the distance vector $(\,\boldsymbol{y}-\boldsymbol{x})$ with $\boldsymbol{f}$ . Faxeń laws infer particle motion in response to the stresses exerted on its surface by a flow field, and are effectively integral statements of the no-slip boundary condition,

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle U_{i}=\frac{1}{4\unicode[STIX]{x03C0}a^{2}}\int _{S_{\unicode[STIX]{x1D6FC}}}\,\text{d}S_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x})u_{i}(\boldsymbol{x}), & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D60B}_{ij}=\frac{3}{4\unicode[STIX]{x03C0}a^{3}}\int _{S_{\unicode[STIX]{x1D6FC}}}\,\text{d}S_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x})u_{i}\hat{n}_{j}, & \displaystyle\end{eqnarray}$$

where $S_{\unicode[STIX]{x1D6FC}}$ is the surface of particle $\unicode[STIX]{x1D6FC}$ , and $V_{\unicode[STIX]{x1D6FC}}$ is its volume. The particle velocity gradient, $\unicode[STIX]{x1D63F}$ , is composed of a symmetric part, the rate of strain  $\boldsymbol{E}$ , and an anti-symmetric part, the angular velocity  $\unicode[STIX]{x1D734}$ .

The linearity of the Stokes equations allows independent coupling and superposition of the total force or couplet to the velocity or velocity gradient. Substituting the fluid velocity (3.1) and the multipole expansion (3.2) into the Faxén laws (3.5), (3.6) gives the following relations between the force moments on particle $\unicode[STIX]{x1D6FD}$ and the velocity moments induced on particle  $\unicode[STIX]{x1D6FC}$

(3.7) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}U_{i}^{\unicode[STIX]{x1D6FC}}\\ \unicode[STIX]{x1D60B}_{ij}^{\unicode[STIX]{x1D6FC}}\end{array}\right]=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D614}_{UF,im}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D614}_{UC,imn}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\\ \unicode[STIX]{x1D614}_{DF,ijm}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D614}_{DC,ijmn}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}c@{}}F_{m}^{\unicode[STIX]{x1D6FD}}\\ C_{mn}^{\unicode[STIX]{x1D6FD}}\end{array}\right],\end{eqnarray}$$

where the elements of the pair mobility tensors are given in terms of integrals of the Green’s function over the particle surfaces

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D614}_{UF,im}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\frac{1}{4\unicode[STIX]{x03C0}a^{2}}\int _{S_{\unicode[STIX]{x1D6FC}}}\,\text{d}S_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x})\frac{1}{4\unicode[STIX]{x03C0}a^{2}}\int _{S_{\unicode[STIX]{x1D6FD}}}\,\text{d}S_{\unicode[STIX]{x1D6FD}}(\,\boldsymbol{y})\unicode[STIX]{x1D611}_{im}(\boldsymbol{x},\boldsymbol{y}), & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D614}_{UC,imn}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\frac{1}{4\unicode[STIX]{x03C0}a^{2}}\int _{S_{\unicode[STIX]{x1D6FC}}}\,\text{d}S_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x})\frac{3}{4\unicode[STIX]{x03C0}a^{3}}\int _{S_{\unicode[STIX]{x1D6FD}}}\,\text{d}S_{\unicode[STIX]{x1D6FD}}(\,\boldsymbol{y})\unicode[STIX]{x1D611}_{im}(\boldsymbol{x},\boldsymbol{y})\hat{n}_{n}^{\unicode[STIX]{x1D6FD}}, & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D614}_{DF,ijm}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\frac{3}{4\unicode[STIX]{x03C0}a^{3}}\int _{S_{\unicode[STIX]{x1D6FC}}}\,\text{d}S_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x})\hat{n}_{j}^{\unicode[STIX]{x1D6FC}}\frac{1}{4\unicode[STIX]{x03C0}a^{2}}\int _{S_{\unicode[STIX]{x1D6FD}}}\,\text{d}S_{\unicode[STIX]{x1D6FD}}(\,\boldsymbol{y})\unicode[STIX]{x1D611}_{im}(\boldsymbol{x},\boldsymbol{y}), & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D614}_{DC,ijmn}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\frac{3}{4\unicode[STIX]{x03C0}a^{3}}\int _{S_{\unicode[STIX]{x1D6FC}}}\,\text{d}S_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x})\hat{n}_{j}^{\unicode[STIX]{x1D6FC}}\frac{3}{4\unicode[STIX]{x03C0}a^{3}}\int _{S_{\unicode[STIX]{x1D6FD}}}\,\text{d}S_{\unicode[STIX]{x1D6FD}}(\,\boldsymbol{y})\unicode[STIX]{x1D611}_{im}(\boldsymbol{x},\boldsymbol{y})\hat{n}_{n}^{\unicode[STIX]{x1D6FD}}. & \displaystyle\end{eqnarray}$$

For periodic systems, which are commonly used in simulation, the appropriate Green’s function is given by Hasimoto (Reference Hasimoto1959),

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D645}(\boldsymbol{x},\boldsymbol{y})=\frac{1}{\unicode[STIX]{x1D702}V}\mathop{\sum }_{\boldsymbol{k}\neq \mathbf{0}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{x}-\boldsymbol{y})}\frac{1}{k^{2}}(\unicode[STIX]{x1D644}-\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}),\end{eqnarray}$$

where $\unicode[STIX]{x1D645}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{F}$ is the fluid velocity at point $\boldsymbol{x}$ propagated by a force $\boldsymbol{F}$ located at point $\boldsymbol{y}$ and its periodic images on a lattice, $\unicode[STIX]{x1D644}$ is the identity tensor, $\boldsymbol{k}$ are drawn from the set of reciprocal (Fourier) space lattice vectors excluding the zero vector, $\unicode[STIX]{x1D702}$ is fluid viscosity and $V$ is the periodic cell volume. With Hasimoto’s solution, the integrals in (3.8)–(3.11) can be evaluated in principle. The integrals do not necessarily have analytically tractable solutions for arbitrary particle shapes, but are straightforward to evaluate for spherical particles, which are the focus of this work. Substituting (3.12) into (3.8)–(3.11) and evaluating the surface integrals in Fourier space gives algebraic expressions for the pair mobilities of spherical particles with equal, finite size

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D648}_{pq}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\frac{1}{\unicode[STIX]{x1D702}V}\mathop{\sum }_{\boldsymbol{k}\neq \mathbf{0}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{x}_{\unicode[STIX]{x1D6FC}}-\boldsymbol{x}_{\unicode[STIX]{x1D6FD}})}\frac{1}{k^{2}}\boldsymbol{s}^{p}(ka)(\unicode[STIX]{x1D644}-\hat{\boldsymbol{k}}\hat{\boldsymbol{k}})\boldsymbol{s}^{q}(ka).\end{eqnarray}$$

The shape factors are proportional to spherical Bessel functions of the first kind, and can be related to the order of the moment in the multipole expansion $p$ as,

(3.14) $$\begin{eqnarray}\boldsymbol{s}^{p}(ka)=(-1)^{p/2}(2p+1)!!\text{j}_{p}(ka)\,\hat{\boldsymbol{k}}^{p},\end{eqnarray}$$

where $\hat{\boldsymbol{k}}^{p}$ is an irreducible tensor of rank  $p$ (Ladd Reference Ladd1990) and $\text{j}_{m}(x)$ is the spherical Bessel function of the first kind of order $m$ (Abramowitz & Stegun Reference Abramowitz and Stegun1965). The shape factors for the force and linear velocity ( $p=0$ ) are $\boldsymbol{s}^{0}(ka)=\text{j}_{0}(ka)$ and the shape factors for the couplet and velocity gradient ( $p=1$ ) are $\boldsymbol{s}^{1}(ka)=3\sqrt{-1}\boldsymbol{k}\text{j}_{1}(ka)$ .

The sum in (3.13) is slowly converging, but the convergence can be accelerated using the Ewald sum splitting (Ewald Reference Ewald1921), which filters the slowly decaying, large wavenumber terms out of the sum and explicitly reconstructs them in real space, where they converge quickly because they are associated with high-frequency (short wavelength) motions. Mathematically, the mobility sum $\unicode[STIX]{x1D648}_{pq}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\sum _{\boldsymbol{k}\neq \mathbf{0}}\unicode[STIX]{x1D648}_{k}$ is written in two parts $\unicode[STIX]{x1D648}_{pq}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\sum _{\boldsymbol{k}\neq \mathbf{0}}(1-H(k,\unicode[STIX]{x1D709}))\unicode[STIX]{x1D648}_{k}+\sum _{\boldsymbol{k}\neq \mathbf{0}}H(k,\unicode[STIX]{x1D709})\unicode[STIX]{x1D648}_{k}$ , where the first term is analytically reconstructed in real space by applying the inverse Fourier transform and the second sum is left in wave space. Hasimoto (Reference Hasimoto1959) derived the Ewald sum for periodic Stokes flow using the filter

(3.15) $$\begin{eqnarray}H(k,\unicode[STIX]{x1D709})=(1+k^{2}/4\unicode[STIX]{x1D709}^{2})\text{e}^{-k^{2}/4\unicode[STIX]{x1D709}^{2}},\end{eqnarray}$$

and showed that the terms in the resulting real-space sum are exponentially decaying functions of distance, while the terms in the wave-space sum are exponentially decaying functions of wavenumber. The variable $\unicode[STIX]{x1D709}$ is an arbitrary parameter that controls how many terms are filtered out of the wave-space sum; in the limit that $\unicode[STIX]{x1D709}\rightarrow 0$ , the entire sum is evaluated in real space. Using this filter, the real-space and wave-space contributions to the pair mobility tensor, $\unicode[STIX]{x1D648}_{pq}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=(\unicode[STIX]{x1D648}_{pq}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}})_{r}+(\unicode[STIX]{x1D648}_{pq}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}})_{w}$ , are respectively

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D648}_{pq}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}})_{r}=\frac{1}{\unicode[STIX]{x1D702}V}\mathop{\sum }_{\boldsymbol{k}\neq \mathbf{0}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{x}_{\unicode[STIX]{x1D6FC}}-\boldsymbol{x}_{\unicode[STIX]{x1D6FD}})}\frac{1}{k^{2}}(1-H(k,\unicode[STIX]{x1D709}))\boldsymbol{s}^{p}(ka)(\unicode[STIX]{x1D644}-\hat{\boldsymbol{k}}\hat{\boldsymbol{k}})\boldsymbol{s}^{q}(ka), & \displaystyle\end{eqnarray}$$

(3.17) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D648}_{pq}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}})_{w}=\frac{1}{\unicode[STIX]{x1D702}V}\mathop{\sum }_{\boldsymbol{k}\neq \mathbf{0}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{x}_{\unicode[STIX]{x1D6FC}}-\boldsymbol{x}_{\unicode[STIX]{x1D6FD}})}\frac{1}{k^{2}}H(k,\unicode[STIX]{x1D709})\boldsymbol{s}^{p}(ka)(\unicode[STIX]{x1D644}-\hat{\boldsymbol{k}}\hat{\boldsymbol{k}})\boldsymbol{s}^{q}(ka). & \displaystyle\end{eqnarray}$$

The real-space representation of the pair mobility tensors are computed by taking the inverse Fourier transform of (3.16) for each set of shape factors. The compact representations of these tensors are

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D614}_{UF,im}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}})_{r}=c_{UF}^{(1)}(r)(\unicode[STIX]{x1D6FF}_{im}-\hat{r}_{i}\hat{r}_{m})+c_{UF}^{(2)}(r)\hat{r}_{i}\hat{r}_{m}, & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D614}_{UC,imn}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}})_{r}=c_{UC}^{(1)}(r)(\unicode[STIX]{x1D6FF}_{im}\hat{r}_{n}-\hat{r}_{i}\hat{r}_{m}\hat{r}_{m})+c_{UC}^{(2)}(r)(\unicode[STIX]{x1D6FF}_{in}\hat{r}_{m}+\unicode[STIX]{x1D6FF}_{mn}\hat{r}_{i}-4\hat{r}_{i}\hat{r}_{m}\hat{r}_{n}), & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D614}_{DC,ijmn}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}})_{r} & = & \displaystyle c_{DC}^{(1)}(r)(\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D6FF}_{mn}+\unicode[STIX]{x1D6FF}_{im}\unicode[STIX]{x1D6FF}_{jn}-4\unicode[STIX]{x1D6FF}_{in}\unicode[STIX]{x1D6FF}_{jm})+c_{DC}^{(2)}(r)(\unicode[STIX]{x1D6FF}_{jm}\hat{r}_{i}\hat{r}_{n}-\hat{r}_{i}\hat{r}_{j}\hat{r}_{m}\hat{r}_{n})\nonumber\\ \displaystyle & & \displaystyle +\,c_{DC}^{(3)}(r)(\unicode[STIX]{x1D6FF}_{ij}\hat{r}_{m}\hat{r}_{n}+\unicode[STIX]{x1D6FF}_{im}\hat{r}_{j}\hat{r}_{n}+\unicode[STIX]{x1D6FF}_{jn}\hat{r}_{i}\hat{r}_{m}+\unicode[STIX]{x1D6FF}_{mn}\hat{r}_{i}\hat{r}_{j}+\unicode[STIX]{x1D6FF}_{in}\hat{r}_{j}\hat{r}_{m}-6\hat{r}_{i}\hat{r}_{j}\hat{r}_{m}\hat{r}_{n}-\unicode[STIX]{x1D6FF}_{in}\unicode[STIX]{x1D6FF}_{jm}),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\hat{\boldsymbol{r}}$ is the unit distance vector between the particle centres, $r$ is its magnitude, $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta function and the scalar coefficients $c_{UF}^{(1)}$ , $c_{UF}^{(2)}$ , $c_{UC}^{(1)}$ , $c_{UC}^{(2)}$ , $c_{DC}^{(1)}$ , $c_{DC}^{(2)}$ , $c_{DC}^{(3)}$ are exponentially decaying functions of the interparticle distance derived by analytically evaluating the inverse Fourier transform and projecting onto a minimal tensor basis. These functions are detailed in our previous work (Fiore & Swan Reference Fiore and Swan2018). It is important to note that the pair mobility tensor presented here, combined with the Hasimoto filter for the Ewald sum applied in Fourier space produces what we term a positive splitting, i.e. a splitting for which the real-space and wave-space contributions to the operator are independently positive semi-definite, as discussed previously (Fiore et al. Reference Fiore, Usabiaga, Donev and Swan2017; Fiore & Swan Reference Fiore and Swan2018).

The far-field grand mobility tensor, which describes the set of interactions between all particle pairs

(3.21) $$\begin{eqnarray}{\mathcal{M}}=\left[\begin{array}{@{}cc@{}}{\mathcal{M}}_{UF} & {\mathcal{M}}_{UC}\\ {\mathcal{M}}_{DF} & {\mathcal{M}}_{DC}\end{array}\right],\end{eqnarray}$$

whose sub-blocks are built from the pair mobility tensors, e.g.

(3.22) $$\begin{eqnarray}{\mathcal{M}}_{UF}=\left[\begin{array}{@{}cccc@{}}\unicode[STIX]{x1D648}_{UF}^{11} & \unicode[STIX]{x1D648}_{UF}^{12} & \cdots \, & \unicode[STIX]{x1D648}_{UF}^{1N}\\ \unicode[STIX]{x1D648}_{UF}^{21} & \unicode[STIX]{x1D648}_{UF}^{22} & & \unicode[STIX]{x1D648}_{UF}^{2N}\\ \vdots & & \ddots & \vdots \\ \unicode[STIX]{x1D648}_{UF}^{N1} & \unicode[STIX]{x1D648}_{UF}^{N2} & \cdots \, & \unicode[STIX]{x1D648}_{UF}^{NN}\end{array}\right],\end{eqnarray}$$

is also positive semi-definite for all particle configurations and its Ewald sum

(3.23) $$\begin{eqnarray}{\mathcal{M}}={\mathcal{M}}_{r}+{\mathcal{M}}_{w},\end{eqnarray}$$

also represents a positive splitting. This is an important feature and is crucial for fast sampling of the far-field Brownian noise, as discussed in § 4.1.

The wave-space component of the mobility can be written compactly using linear operators

(3.24) $$\begin{eqnarray}{\mathcal{M}}_{w}=\mathbb{D}^{\dagger }\boldsymbol{\cdot }\mathbb{P}^{\dagger }\boldsymbol{\cdot }\mathbb{B}\boldsymbol{\cdot }\mathbb{P}\boldsymbol{\cdot }\mathbb{D},\end{eqnarray}$$

where $\mathbb{D}$ is the unitary non-uniform discrete Fourier transform (NUDFT) operator, and $\mathbb{D}^{\dagger }$ is its adjoint/inverse. $\mathbb{D}$ transforms the force moments applied on particles from real space to Fourier space. $\mathbb{B}$ is a symmetric block diagonal matrix whose $i$ th diagonal block is

(3.25) $$\begin{eqnarray}\mathbb{B}_{i}=\frac{H(k_{i},\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D702}Vk_{i}^{2}}(\unicode[STIX]{x1D644}-\hat{\boldsymbol{k}}_{i}\hat{\boldsymbol{k}}_{i}),\end{eqnarray}$$

where $\hat{\boldsymbol{k}}_{i}$ is $i$ th wavevector; $\mathbb{B}$ multiplies the Fourier components of the force by the Green’s function $1/\unicode[STIX]{x1D702}Vk_{i}^{2}$ and Hasimoto $H(k_{i},\unicode[STIX]{x1D709})$ scaling factors and projects them onto the space of divergence-free velocity fields via the orthogonal projector $(\unicode[STIX]{x1D644}-\hat{\boldsymbol{k}}\hat{\boldsymbol{k}})$ . The projector $\mathbb{P}$ contains the non-negative shape factors that map the Fourier-space force moments onto the force density, and its adjoint $\mathbb{P}^{\dagger }$ maps the Fourier-space velocity field onto the Fourier-space velocity moments; $\mathbb{P}$ is a rectangular matrix whose only non-zero entries are rectangular blocks along the main diagonal,

(3.26) $$\begin{eqnarray}\mathbb{P}_{i}=\left[\begin{array}{@{}cc@{}}\,\text{j}_{0}(k_{i}a)\unicode[STIX]{x1D644} & 3\sqrt{-1}\text{j}_{1}(k_{i}a)\boldsymbol{k}_{i}\end{array}\right].\end{eqnarray}$$

The matrix $\mathbb{P}_{i}^{\dagger }\boldsymbol{\cdot }\mathbb{B}_{i}\boldsymbol{\cdot }\mathbb{P}_{i}$ is the wave-space representation of the mobility tensor at the $i$ th wave vector, i.e. the summand in (3.17) except for the exponential factor, where the exponential factors and the sum over $\boldsymbol{k}\neq \mathbf{0}$ come from the application of the NUDFT and it inverse, $\mathbb{D}$ and  $\mathbb{D}^{\dagger }$ .

The results presented in this section hold for monodisperse suspensions of spherical particles. Polydispersity can be incorporated in the expressions presented in this section through subtle modifications to the shape factors and scalar functions. This application is discussed in appendix A. Non-spherical particles can be built from rigid assemblies of spherical beads represented using the shape factors presented here. A detailed discussion of this method is deferred to § 6.

3.2 Near-field hydrodynamics

We construct the near-field lubrication tensor $\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{nf}$ separately for each pair of particles $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ whose distance is shorter than a cutoff radius of $r_{lub}$ . The exact series expansion for the hydrodynamic resistance between two spheres (Jeffrey & Onishi Reference Jeffrey and Onishi1984; Jeffrey Reference Jeffrey1992b ) provides the basis for the SD lubrication tensor. The series expansion is equivalent to carrying the multipole expansion for the isolated particle pair out to an infinite number of terms, so the corresponding resistance tensor $\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{2B}$ contains within it all the multipole moments associated with each particle pair. The monopole, dipole and isotropic quadrupole and octupole moments are already included in the far-field approximation to the hydrodynamic mobility for each pair, ${\mathcal{M}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{2B}$ , so their contribution to the resistance tensor for each particle pair, ${\mathcal{R}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{2B}=({\mathcal{M}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{2B})^{-1}$ , is subtracted from the exact result to avoid double counting them. Therefore, the pair resistance tensor used is

(3.27) $$\begin{eqnarray}\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{nf}=\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{2B}-{\mathcal{R}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{2B}.\end{eqnarray}$$

This pair resistance tensor is comprised of four sub-blocks describing the $FU$ , $FE$ , $SU$ and $SE$ couplings for the particle pair,

(3.28) $$\begin{eqnarray}\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{nf}=\left[\begin{array}{@{}cc@{}}(\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & (\unicode[STIX]{x1D64D}_{FE}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\\ (\unicode[STIX]{x1D64D}_{SU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & (\unicode[STIX]{x1D64D}_{SE}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\end{array}\right],\end{eqnarray}$$

where

(3.29) $$\begin{eqnarray}(\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\left[\begin{array}{@{}cc@{}}(\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{11} & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{12}\\ (\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{21} & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{22}\end{array}\right]\end{eqnarray}$$

is the $FU$ component of the resistance tensor for the pair of particles $\unicode[STIX]{x1D6FC}$ and  $\unicode[STIX]{x1D6FD}$ . In this expression, the indices $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are the global indices of the two particles constituting the pair, i.e. with $N$ particles $1\leqslant \unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}\leqslant N$ , while within the pair particles $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are denoted by the indices 1 and 2, respectively. Similar definitions hold for the other three sub-blocks.

The near-field lubrication tensor for all particles is constructed by summing the interactions over all particle pairs

(3.30) $$\begin{eqnarray}\unicode[STIX]{x1D64D}_{FU}^{nf}=\mathop{\sum }_{\unicode[STIX]{x1D6FC}}\mathop{\sum }_{\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FC}}\boldsymbol{S}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\boldsymbol{\cdot }(\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\boldsymbol{\cdot }\boldsymbol{S}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{\text{T}},\end{eqnarray}$$

where $\boldsymbol{S}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ maps the forces for a given particle pair into the list of all particle forces, its transpose projects the velocities for the pair out of the set of all particle velocities and $(\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\mathbf{0}$ if $r_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}>4a$ . Written explicitly in matrix notation,

(3.31) $$\begin{eqnarray}\unicode[STIX]{x1D64D}_{FU}^{nf}=\left[\begin{array}{@{}cccc@{}}\displaystyle \mathop{\sum }_{\unicode[STIX]{x1D6FD}\neq 1}(\unicode[STIX]{x1D64D}_{FU}^{nf})_{1\unicode[STIX]{x1D6FD}}^{11} & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{12}^{12} & \cdots \, & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{1N}^{12}\\ (\unicode[STIX]{x1D64D}_{FU}^{nf})_{21}^{12} & \displaystyle \mathop{\sum }_{\unicode[STIX]{x1D6FD}\neq 2}(\unicode[STIX]{x1D64D}_{FU}^{nf})_{2\unicode[STIX]{x1D6FD}}^{11} & & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{2N}^{12}\\ \vdots & & \ddots & \vdots \\ (\unicode[STIX]{x1D64D}_{FU}^{nf})_{N1}^{12} & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{N2}^{12} & \cdots \, & \displaystyle \mathop{\sum }_{\unicode[STIX]{x1D6FD}\neq N}(\unicode[STIX]{x1D64D}_{FU}^{nf})_{N\unicode[STIX]{x1D6FD}}^{11}\end{array}\right].\end{eqnarray}$$

The SU, FE and SE components are constructed similarly. In these expressions, and throughout the rest of this paper, we use $\boldsymbol{F}$ to denote the generalized force, which is comprised of the force and torque, and $\boldsymbol{U}$ to denote the generalized velocity, which is comprised of the linear and angular velocity. These can also be thought of as the forces and velocities associated with rigid body motion. In the notation of Jeffrey & Onishi (Reference Jeffrey and Onishi1984),

(3.32) $$\begin{eqnarray}(\unicode[STIX]{x1D64D}_{FU})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \tilde{\boldsymbol{B}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\\ \boldsymbol{B}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \boldsymbol{ C}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\end{array}\right],\end{eqnarray}$$

and in the notation of Jeffrey (Reference Jeffrey1992b ),

(3.33) $$\begin{eqnarray}(\unicode[STIX]{x1D64D}_{SU})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\left[\begin{array}{@{}cc@{}}\boldsymbol{G}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \boldsymbol{H}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\end{array}\right].\end{eqnarray}$$

The far-field two-body contribution to the near-field resistance tensor is computed from the inverse of the mobility tensor for a pair of particles

(3.34) $$\begin{eqnarray}{\mathcal{R}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{2B}=\left[\begin{array}{@{}cccc@{}}\unicode[STIX]{x1D64D}_{FU}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D64D}_{FU}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D64D}_{SU}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D64D}_{SU}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\\ \unicode[STIX]{x1D64D}_{FU}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D64D}_{FU}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D64D}_{SU}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D64D}_{SU}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}}\\ \unicode[STIX]{x1D64D}_{FE}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D64D}_{FE}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D64D}_{SE}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D64D}_{SE}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\\ \unicode[STIX]{x1D64D}_{FE}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D64D}_{FE}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D64D}_{SE}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D64D}_{SE}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}}\end{array}\right]=\left[\begin{array}{@{}cccc@{}}\unicode[STIX]{x1D648}_{UF}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D648}_{UF}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D648}_{US}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D648}_{US}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\\ \unicode[STIX]{x1D648}_{UF}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D648}_{UF}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D648}_{US}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D648}_{US}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}}\\ \unicode[STIX]{x1D648}_{EF}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D648}_{EF}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D648}_{ES}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D648}_{ES}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\\ \unicode[STIX]{x1D648}_{EF}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D648}_{EF}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}} & \unicode[STIX]{x1D648}_{ES}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & \unicode[STIX]{x1D648}_{ES}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}}\end{array}\right]^{-1}.\end{eqnarray}$$

As both the mobility and resistance tensors can be written in terms of scalar functions and orthogonal tensors, the inversion in (3.34) can be used to relate the scalar mobility and resistance functions, and therefore does not have to be performed during the simulation. Because the terms of $\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{nf}=\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}-{\mathcal{R}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ decay very quickly (for $r\gg 2a$ , $\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\approx {\mathcal{R}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ ) when constructing the near-field resistance tensor for all particles, only the entries associated with particle pairs whose distance is within a cutoff radius of $r_{lub}=4a$ are retained, consistent with previous implementations of Stokesian dynamics (Sierou & Brady Reference Sierou and Brady2001).

3.3 Hydrodynamic saddle point solution

As described in § 2, the far-field and near-field contributions to the hydrodynamic interactions can be combined in a saddle point matrix to solve for the unknown particle velocities and far-field hydrodynamic force moments. Iterative solutions of this matrix equation can be effectively preconditioned with a block preconditioner based on the Schur factorization. The efficacy of this block saddle preconditioner, equation (2.12), for iterative GMRES inversion of the saddle point matrix, equation (2.8), is shown in figure 1. Specifically, the convergence rate of the relative residual as a function of  $r_{p}$ , normalized by the particle radius  $a$ , is given in figure 1(a) for a random dispersion comprised of 1000 hard spheres at a volume fraction of $\unicode[STIX]{x1D719}=0.2$ . A relatively short cutoff radius, $r_{p}=2.1$ , is sufficient to drastically improve the convergence rate, and is almost as effective as larger cutoffs, $r_{p}=2.5$ or $r_{p}=3.0$ . We choose $r_{p}=2.1$ because it is very effective and requires minimal memory and operation counts to evaluate. It effectively only considers particles within the first neighbour shell. Within a distance of $r_{p}=2.1$ , monodisperse spheres can have at most 12 neighbours, so the memory required to construct the preconditioner is limited. For fixed  $r_{p}$ , the preconditioner performance is expected to deteriorate with increasing $\unicode[STIX]{x1D719}$ as particles gain more neighbours within $r_{lub}$ that are not also within  $r_{p}$ . Figure 1(b) shows the convergence rate of the preconditioner with $r_{p}=2.1$ across a range of volume fractions. The convergence rate decreases modestly with increasing volume fraction but the number of iterations remains acceptably small. For large scale simulations, the most important metric is the convergence rate as a function of system size. This is shown in figure 1(c) for $r_{p}=2.1$ and $\unicode[STIX]{x1D719}=0.2$ . The convergence rate is effectively independent of system size, given that the number of iterations required to converge to a given tolerance increases by less than a factor of two when the number of particles is increased by three orders of magnitude.

Figure 1.

Convergence rate of the relative error in the iterative GMRES solution of the saddle point formulation for the hydrodynamic interactions. (a) Convergence rate dependence on the preconditioner cutoff radius (in units of particle radius) $r_{p}=\{2.0,2.01,2.1,2.5\}$ for a random hard sphere dispersion with $N=1000$ and $\unicode[STIX]{x1D719}=0.2$ . Although not shown on the plot, the convergence rate with $r_{p}=3.0$ is identical to that for $r_{p}=2.5$ . (b) Convergence rate dependence on the volume fraction $\unicode[STIX]{x1D719}=\{0.1,0.2,0.3,0.4\}$ for a random hard sphere dispersion with $N=1000$ and $r_{p}=2.1$ . (c) Convergence rate dependence on the system size $N=\{10^{2},10^{3},10^{4},10^{5}\}$ for random hard sphere dispersions with $\unicode[STIX]{x1D719}=0.2$ and $r_{p}=2.1$ .

4.1 Far-field Brownian slip velocity

The positively split Ewald (PSE) approach to sampling the far-field stochastic velocities has been described in great detail previously (Fiore et al. Reference Fiore, Usabiaga, Donev and Swan2017; Fiore & Swan Reference Fiore and Swan2018), but will be summarized here for clarity. In the PSE approach, the Brownian velocity ${\mathcal{U}}^{B}\sim {\mathcal{N}}(\mathbf{0},{\mathcal{M}})$ is decomposed into two components, corresponding to the real-space and wave-space contributions to the Ewald sum

(4.1) $$\begin{eqnarray}{\mathcal{U}}^{B}={\mathcal{U}}_{r}^{B}+{\mathcal{U}}_{w}^{B},\end{eqnarray}$$

where

(4.2a,b ) $$\begin{eqnarray}{\mathcal{U}}_{r}^{B}\sim {\mathcal{N}}\left(\mathbf{0},\frac{2k_{B}T}{\unicode[STIX]{x0394}t}{\mathcal{M}}_{r}\right),\quad {\mathcal{U}}_{w}^{B}\sim {\mathcal{N}}\left(\mathbf{0},\frac{2k_{B}T}{\unicode[STIX]{x0394}t}{\mathcal{M}}_{w}\right).\end{eqnarray}$$

Samples from these two distributions are drawn using separate techniques tailored to their different operator structures.

The real-space contribution to the displacements is sampled from  ${\mathcal{M}}_{r}$ ,

(4.3) $$\begin{eqnarray}{\mathcal{U}}_{r}^{B}={\mathcal{M}}_{r}^{1/2}\boldsymbol{\cdot }\unicode[STIX]{x1D733}_{r},\end{eqnarray}$$

using the iterative Lanczos sampling technique developed by Chow & Saad (Reference Chow and Saad2014). The rate of convergence of the Lanczos iteration is controlled by the condition number of the sampled matrix. Because ${\mathcal{M}}_{r}$ contains exponentially decaying interactions evaluated only for particles within a small neighbourhood of each other, its condition number is both small and independent of the total number of particles modelled (Fiore et al. Reference Fiore, Usabiaga, Donev and Swan2017; Fiore & Swan Reference Fiore and Swan2018). Therefore, the Lanczos iterates converge quickly for any number of particles.

The wave-space contribution to the displacements is sampled from ${\mathcal{M}}_{w}$ whose compact operator representation is,

(4.4) $$\begin{eqnarray}{\mathcal{M}}_{w}=\mathbb{D}^{\dagger }\boldsymbol{\cdot }\mathbb{P}^{\dagger }\boldsymbol{\cdot }\mathbb{B}\boldsymbol{\cdot }\mathbb{P}\boldsymbol{\cdot }\mathbb{D},\end{eqnarray}$$

which can be written as the product of a matrix and its adjoint

(4.5) $$\begin{eqnarray}{\mathcal{M}}_{w}=(\mathbb{D}^{\dagger }\boldsymbol{\cdot }\mathbb{P}^{\dagger }\boldsymbol{\cdot }\mathbb{B}^{1/2})\boldsymbol{\cdot }(\mathbb{D}^{\dagger }\boldsymbol{\cdot }\mathbb{P}^{\dagger }\boldsymbol{\cdot }\mathbb{B}^{1/2})^{\dagger }.\end{eqnarray}$$

The square root of a matrix $\unicode[STIX]{x1D653}$ useful for stochastic sampling is defined as any matrix $\unicode[STIX]{x1D653}^{1/2}$ such that $\unicode[STIX]{x1D653}^{1/2}\boldsymbol{\cdot }(\unicode[STIX]{x1D653}^{1/2})^{\dagger }=\unicode[STIX]{x1D653}$ . Therefore, one possible (and very useful) square root of ${\mathcal{M}}_{w}$ is given analytically by

(4.6) $$\begin{eqnarray}{\mathcal{M}}_{w}^{1/2}=\mathbb{D}^{\dagger }\boldsymbol{\cdot }\mathbb{P}^{\dagger }\boldsymbol{\cdot }\mathbb{B}^{1/2},\end{eqnarray}$$

so that

(4.7) $$\begin{eqnarray}{\mathcal{U}}_{w}^{B}=\mathbb{D}^{\dagger }\boldsymbol{\cdot }\mathbb{P}^{\dagger }\boldsymbol{\cdot }\mathbb{B}^{1/2}\boldsymbol{\cdot }\unicode[STIX]{x1D733}_{w},\end{eqnarray}$$

which is no more challenging to compute than the product of the mobility tensor with a vector.

Asserting that $\unicode[STIX]{x1D733}_{r}$ and $\unicode[STIX]{x1D733}_{w}$ are statistically independent guarantees that ${\mathcal{U}}_{r}^{B}$ and ${\mathcal{U}}_{w}^{B}$ are also statistically independent. Furthermore, if $\unicode[STIX]{x1D733}_{r}$ and $\unicode[STIX]{x1D733}_{w}$ are computed such that $\unicode[STIX]{x1D733}_{r}\sim \unicode[STIX]{x1D733}_{w}\sim {\mathcal{N}}(\mathbf{0},(2k_{B}T/\unicode[STIX]{x0394}t)\,\unicode[STIX]{x1D644})$ , then

(4.8) $$\begin{eqnarray}{\mathcal{U}}_{r}^{B}+{\mathcal{U}}_{w}^{B}\sim {\mathcal{N}}\left(\mathbf{0},\frac{2k_{B}T}{\unicode[STIX]{x0394}t}({\mathcal{M}}_{r}+{\mathcal{M}}_{w})\right),\end{eqnarray}$$

as required by the fluctuation–dissipation theorem, where ${\mathcal{M}}={\mathcal{M}}_{r}+{\mathcal{M}}_{w}$ . It is important to note that this split sampling approach only works when ${\mathcal{M}}_{r}$ and ${\mathcal{M}}_{w}$ are independently positive semi-definite, which guarantees that they each have a well-defined square root. To our knowledge, the only representation of ${\mathcal{M}}_{r}$ and ${\mathcal{M}}_{w}$ for which this is true is the formulation we have presented previously (Fiore & Swan Reference Fiore and Swan2018) and reviewed in § 3.1.

4.2 Near-field Brownian force

Banchio & Brady (Reference Banchio and Brady2003) first developed a method to sample the near-field Brownian force in the context of ASD, avoiding the $O(N^{3})$ cost for direct evaluation of the square root by individually computing the force for each set of neighbouring particles $\unicode[STIX]{x1D6FC}$ and  $\unicode[STIX]{x1D6FD}$ ,

(4.9) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}(\boldsymbol{F}_{nf}^{B})^{\unicode[STIX]{x1D6FC}}\\ (\boldsymbol{F}_{nf}^{B})^{\unicode[STIX]{x1D6FD}}\end{array}\right]\sim {\mathcal{N}}\left(\mathbf{0},\frac{2k_{B}T}{\unicode[STIX]{x0394}t}\left[\begin{array}{@{}cc@{}}(\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}} & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\\ (\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}} & (\unicode[STIX]{x1D64D}_{FU}^{nf})_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}}\end{array}\right]\right),\end{eqnarray}$$

and showed that the total force vector is appropriately distributed, i.e. $\boldsymbol{F}_{nf}^{B}\sim {\mathcal{N}}(\mathbf{0},2k_{B}T\,\unicode[STIX]{x1D64D}_{FU}^{nf})$ . The sampling for each pair of particles within $r_{lub}$ is evaluated using the complete lower Cholesky decomposition $\boldsymbol{L}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ of the pair resistance tensor

(4.10) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}(\boldsymbol{F}_{nf}^{B})^{\unicode[STIX]{x1D6FC}}\\ (\boldsymbol{F}_{nf}^{B})^{\unicode[STIX]{x1D6FD}}\end{array}\right]=\boldsymbol{L}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\boldsymbol{\cdot }\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D733}_{nf}^{\unicode[STIX]{x1D6FC}}\\ \unicode[STIX]{x1D733}_{nf}^{\unicode[STIX]{x1D6FD}}\end{array}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D733}_{nf}^{\unicode[STIX]{x1D6FC}}$ and $\unicode[STIX]{x1D733}_{nf}^{\unicode[STIX]{x1D6FD}}$ are uncorrelated and $\unicode[STIX]{x1D733}_{nf}^{\unicode[STIX]{x1D6FC}}\sim \unicode[STIX]{x1D733}_{nf}^{\unicode[STIX]{x1D6FD}}\sim {\mathcal{N}}(\mathbf{0},(2k_{B}T/\unicode[STIX]{x0394}t)\,\unicode[STIX]{x1D644})$ . This approach is efficient from a global memory perspective, because only the pairwise resistance tensors need to be stored, but is inefficient from a performance perspective because each particle pair must be individually computed, requiring $O(N)$ Cholesky factorizations. Some parallelization of this scheme is possible at the cost of increased memory usage, but large scale parallelization and acceleration of the pair Cholesky calculation is infeasible on current computer hardware.

Instead, we choose to directly sample $\boldsymbol{F}_{nf}^{B}$ for all particles simultaneously from $\unicode[STIX]{x1D64D}_{FU}^{nf}$ using the iterative Lanczos scheme of Chow & Saad (Reference Chow and Saad2014). $\unicode[STIX]{x1D64D}_{FU}^{nf}$ is positive semi-definite, which guarantees that the square root exists and the resulting forces are real valued. Although the Lanczos scheme can be applied directly to sample Brownian forces, $\unicode[STIX]{x1D64D}_{FU}^{nf}$ is ill-conditioned, so effective preconditioning schemes are required to accelerate convergence. The preconditioner described in § 2.3 is particularly well suited to this problem.

Figure 2.

Effectiveness of the preconditioner for iterative Krylov subspace approximation to the square root of the near-field resistance tensor. (a) Eigenspectrum of the near-field lubrication tensor for a random suspension of $N=1000$ hard spheres at $\unicode[STIX]{x1D719}=0.2$ (solid black line), modified eigenspectrum after applying the incomplete Cholesky preconditioner (dashed light grey line, overlaps the black line for much of the range), modified eigenspectrum after applying both the incomplete Cholesky and diagonal preconditioners (dashed-dot medium grey line). The dotted black line denotes a uniform eigenspectrum with $\unicode[STIX]{x1D706}_{i}=1$ , which would result from ‘perfect’ preconditioning. (b) Condition number of resistance tensor before and after preconditioning for random hard sphere suspension comprised of $N=1000$ particles at various volume fractions. Same line definitions as in panel (a). (c) Number of iterations required for convergence of the method to a relative step-norm error of $10^{-6}$ for the hard sphere suspensions in panel (b). (d) Number of iterations required for convergence to a relative step-norm error of $10^{-3}$ , typical of many simulations, for the fully preconditioned sampling scheme applied to hard sphere dispersions with $\unicode[STIX]{x1D719}=0.2$ at varying system sizes.

Figure 2(a) shows how these preconditioners compress the spectral radius of the near-field resistance tensor. Specifically, as intended, the incomplete Cholesky factor removes the largest magnitude eigenvalues, while the diagonal preconditioner shifts all eigenvalues that are less than one to be clustered around one. The effect of this eigenspectrum compression on the condition number is shown in figure 2(b) for hard sphere suspensions of 1000 particles across a wide range of volume fractions. The preconditioners produce a well-conditioned operator for all volume fraction and the condition number of the preconditioned operator is independent of  $\unicode[STIX]{x1D719}$ . Figure 2(c) shows the number of iterations required for the Krylov subspace sampling technique to converge to a step-norm relative error of  $10^{-6}$ . With our preconditioning scheme, the iteration count is independent of the volume fraction, and is 3–15 times smaller than without preconditioning. Most impressively, figure 2(d) shows that the iteration count is effectively independent of the system size; increasing the number of particles of particles by a factor of 1000 only increases the number of iterations by 50 %.

4.3 Brownian drift

We use the random finite differencing (RFD) scheme of Delong et al. (Reference Delong, Usabiaga, Delgado-Buscalioni, Griffith and Donev2014) to compute the Brownian drift, $k_{B}T\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64D}_{FU}$ . This method was already clearly described by Delong et al., but will be summarized here for completeness. The idea of RFD is to approximate the divergence of a matrix $\unicode[STIX]{x1D63C}(\boldsymbol{q})$ , which is a function of the set of particle coordinates  $\boldsymbol{q}$ , as the difference between the matrix evaluated at two nearby points. The central random finite difference is

(4.11) $$\begin{eqnarray}\frac{\left\langle \unicode[STIX]{x1D63C}\left(\boldsymbol{q}+{\displaystyle \frac{\unicode[STIX]{x1D716}}{2}}\unicode[STIX]{x0394}\boldsymbol{q}\right)\boldsymbol{\cdot }\unicode[STIX]{x0394}\boldsymbol{q}-\unicode[STIX]{x1D63C}\left(\boldsymbol{q}-{\displaystyle \frac{\unicode[STIX]{x1D716}}{2}}\unicode[STIX]{x0394}\boldsymbol{q}\right)\boldsymbol{\cdot }\unicode[STIX]{x0394}\boldsymbol{q}\right\rangle }{\unicode[STIX]{x1D716}}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63C}+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

where $\unicode[STIX]{x0394}\boldsymbol{q}\sim {\mathcal{N}}(\mathbf{0},\unicode[STIX]{x1D644})$ , $\unicode[STIX]{x1D716}$ is a very small discretization parameter, and the finite difference approximates the divergence of $\unicode[STIX]{x1D63C}$ in expectation. Methods for choosing $\unicode[STIX]{x1D716}$ that maximize the accuracy of the RFD calculation are detailed by Sprinkle et al. (Reference Sprinkle, Balboa Usabiaga, Patankar and Donev2017). This scheme applied to FSD nominally requires two saddle point solves, one at each of the displaced points. However, because the displacements are small, the particle positions are quite close and the solution at one set of coordinates serves as an excellent initial guess for the solution at the other set. The second solve converges very rapidly ( $O(1)$ iterations are required in our tests, compared to $O(10)$ for the first solve), and therefore contributes a negligible amount to the overall computational cost of the algorithm.

Using RFD, the mean drift of the Brownian displacement and stresslet can be computed from saddle point solves,

(4.12) $$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}{\mathcal{M}}(\boldsymbol{q}\pm \unicode[STIX]{x0394}\boldsymbol{q}) & \boldsymbol{B}\\ \boldsymbol{B}^{\text{T}} & -\unicode[STIX]{x1D64D}_{FU}^{nf}(\boldsymbol{q}\pm \unicode[STIX]{x0394}\boldsymbol{q})\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}c@{}}{\mathcal{F}}_{\pm }\\ \boldsymbol{U}_{\pm }\end{array}\right]=\left[\begin{array}{@{}c@{}}\mathbf{0}\\ -\unicode[STIX]{x0394}\boldsymbol{ q}\end{array}\right],\end{eqnarray}$$

where ${\mathcal{F}}_{\pm }$ and $\boldsymbol{U}_{\pm }$ are the solutions at the positive and negative displacements for the central difference. The far-field hydrodynamic stresslet is $\boldsymbol{S}_{\pm }^{ff}=[\mathbf{0}~~\unicode[STIX]{x1D644}]\boldsymbol{\cdot }{\mathcal{F}}_{\pm }$ and the total stresslet is

(4.13) $$\begin{eqnarray}\boldsymbol{S}_{\pm }=\boldsymbol{S}_{\pm }^{ff}+\unicode[STIX]{x1D64D}_{SU}^{nf}(\boldsymbol{q}\pm \unicode[STIX]{x0394}\boldsymbol{q})\boldsymbol{\cdot }\boldsymbol{U}_{\pm }.\end{eqnarray}$$

The mean drift of the Brownian displacement is

(4.14) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64D}_{FU}^{-1}=\frac{1}{\unicode[STIX]{x1D716}}\langle \boldsymbol{U}_{+}-\boldsymbol{U}_{-}\rangle +O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

and the mean drift of the Brownian stresslet is

(4.15) $$\begin{eqnarray}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D64D}_{SU}\boldsymbol{\cdot }\unicode[STIX]{x1D64D}_{FU}^{-1})=\frac{1}{\unicode[STIX]{x1D716}}\langle \boldsymbol{S}_{+}-\boldsymbol{S}_{-}\rangle +O(\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$

Figure 3.

Transport coefficients for simple cubic lattices as a function of volume fraction. (a) Spin viscosity, $\unicode[STIX]{x1D701}$ , computed using the FSD algorithm (circles) compared to ASD results of Sierou & Brady (Reference Sierou and Brady2001) (line). (b) Independent components, $\unicode[STIX]{x1D6FC}$ (circles) and $\unicode[STIX]{x1D6FD}$ (squares), of the high-frequency shear viscosity, $\unicode[STIX]{x1D702}^{S}$ . The solid and dashed lines are the values of $\unicode[STIX]{x1D6FC}$ and  $\unicode[STIX]{x1D6FD}$ , respectively, computed by Brady et al. (Reference Brady, Phillips, Lester and Bossis1988).

5.1 Algorithm validation

A specific implementation of the FSD algorithm for monodisperse spheres is given in appendix B. We validate our FSD algorithm against known static and dynamic results for monodisperse dispersions. Figure 3 shows the spin viscosity and high-frequency shear viscosity of simple cubic lattices of particles. The spin viscosity $\boldsymbol{\unicode[STIX]{x1D702}}^{L}$ for a cubic lattice defines the hydrodynamic torque $\boldsymbol{L}^{H}$ experienced by the particles as they rotate with angular velocity  $\unicode[STIX]{x1D734}$ ,

(5.1) $$\begin{eqnarray}\boldsymbol{L}^{H}=-\boldsymbol{\unicode[STIX]{x1D702}}^{L}\boldsymbol{\cdot }\unicode[STIX]{x1D734}.\end{eqnarray}$$

Because the cubic lattice is isotropic, the spin viscosity can be represented by a single scalar value,

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{ij}^{L}=\unicode[STIX]{x1D701}\unicode[STIX]{x1D6FF}_{ij}.\end{eqnarray}$$

Similarly, the shear viscosity defines the relationship between the hydrodynamic stresslet $\boldsymbol{S}^{H}$ acquired by particles in response to a rate of strain  $\boldsymbol{E}$ ,

(5.3) $$\begin{eqnarray}\boldsymbol{S}^{H}=-\boldsymbol{\unicode[STIX]{x1D702}}^{S}\boldsymbol{ : }\boldsymbol{E}.\end{eqnarray}$$

Due to symmetry, the fourth-order tensor $\boldsymbol{\unicode[STIX]{x1D702}}^{S}$ describing the shear viscosity can be decomposed in terms of two scalar coefficients (Brady et al. Reference Brady, Phillips, Lester and Bossis1988),

(5.4) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{ijkl}^{S}={\textstyle \frac{1}{2}}(1+\unicode[STIX]{x1D6FD})\left(\unicode[STIX]{x1D6FF}_{ik}\unicode[STIX]{x1D6FF}_{jl}+\unicode[STIX]{x1D6FF}_{il}\unicode[STIX]{x1D6FF}_{jk}-{\textstyle \frac{2}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D6FF}_{kl}\right)+(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD})\left(\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D6FF}_{ik}\unicode[STIX]{x1D6FF}_{il}\unicode[STIX]{x1D6FF}_{jk}\unicode[STIX]{x1D6FF}_{jl}\unicode[STIX]{x1D6FF}_{kl}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D6FF}_{kl}\right)\!.\end{eqnarray}$$

In general, these viscosity expressions require an ensemble average over all particles, but all particles are identical in a cubic lattice so the averages are omitted here. The scalar coefficients of the spin and shear viscosities were computed using a relative error tolerance of $10^{-3}$ in all approximations. The scalar spin viscosity $\unicode[STIX]{x1D701}$ is plotted as a function of volume fraction in figure 3(a) and is compared with the values computed by Sierou & Brady (Reference Sierou and Brady2001) using the accelerated Stokesian dynamics (ASD) algorithm. The two scalar components to the shear viscosity, $\unicode[STIX]{x1D6FC}$ and  $\unicode[STIX]{x1D6FD}$ , are plotted as a function of volume fraction in figure 3(b) and are compared with the values obtained by Brady et al. (Reference Brady, Phillips, Lester and Bossis1988) using the original SD algorithm. In all cases, excellent quantitative agreement is observed between our FSD result and the standard Stokesian dynamics approaches.

Figure 4.

Equilibrium and transport characteristics for random hard sphere suspensions. (a) Static structure factor for $N=1000$ particles computed from dynamic simulations using the FSD algorithm (lines) compared to the analytical Percus–Yevick result for hard spheres (Regnaut & Ravey Reference Regnaut and Ravey1989, Reference Regnaut and Ravey1990). (b) Short-time self-diffusion coefficient computed from the constrained stochastic displacements $D_{s}=(\unicode[STIX]{x0394}t/6k_{B}T)\langle \boldsymbol{U}^{B}\boldsymbol{\cdot }\boldsymbol{U}^{B}\rangle$ generated by the FSD algorithm (circles) compared to the analytical results of Ladd (Reference Ladd1990). Error bars represent standard deviations of the measurements.

To validate the proposed sampling and integration methods, we performed dynamic simulations of hard sphere particles with a relative error tolerance of $10^{-3}$ for all hydrodynamics calculations and a time step of $\unicode[STIX]{x0394}t=10^{-3}$ . To avoid unphysical overlap of the hard sphere particles, we used a short-ranged repulsive potential (Sierou & Brady Reference Sierou and Brady2001). We computed the structure factor of the resulting equilibrium microstructures from 100 configurations extracted from the particle trajectory every 10 bare particle diffusion times. Figure 4(a) shows the structure factor computed from our simulations compared to analytical results derived using the Percus–Yevick closure (Regnaut & Ravey Reference Regnaut and Ravey1989, Reference Regnaut and Ravey1990). The quantitative agreement between our simulations and the analytical theory indicates that the stochastic displacements are consistent with fluctuation–dissipation and produce the correct equilibrium particle configurations. Figure 4(b) shows the short-time self-diffusion coefficient computed from the ensemble mean of the constrained Brownian velocity,  $\boldsymbol{U}^{B}$ ,

(5.5) $$\begin{eqnarray}D_{s}=\frac{1}{3}\frac{\unicode[STIX]{x0394}t}{2k_{B}T}\langle \boldsymbol{U}^{B}\boldsymbol{\cdot }\boldsymbol{U}^{B}\rangle ,\end{eqnarray}$$

compared to the high-resolution calculations of Ladd (Reference Ladd1990). Again, quantitative agreement is observed between our calculations and established literature results, indicating that the simulation produces the correct dynamic displacements.

5.2 Algorithm performance

The SD algorithm presented here, which we call fast Stokesian dynamics, is designed for efficient scaling and computations in large system. As detailed in § 4 and appendix B, each component of the FSD algorithm has a computational cost that scales linearly with the number of particles. The most time-consuming calculation is the saddle point inversion because no more effective preconditioners for solution of the stresslet constraint are known. Further acceleration to the algorithm could be achieved by the development of such a preconditioner, should one exist, and would reduce the number of iterations required for the solve. However, the improvement would be modest, as the solution is already evaluated in $O(10)$ iterations, as outlined in § 3.3. In this section, we describe the performance characteristics of the FSD algorithm and compare them to our previous work with simpler approximations of the HI.

The expected performance of the FSD algorithm can be estimated from the previously observed performance of our RPY code. This value will serve as a performance baseline. At the RPY level of approximation for the hydrodynamic interactions, the PSE algorithm has a peak performance of around $3\times 10^{6}$ particle time steps per second (PTPS) (Fiore et al. Reference Fiore, Usabiaga, Donev and Swan2017), which is independent of the system size. PTPS is equal to the total number of particles in the simulation divided by the time (in seconds) required to performed a single simulation time step’s worth of calculations. The RPY performance can be used to estimate the optimal performance of the PSE algorithm including lubrication by comparing the operation counts between the two methods. At the RPY level of approximation, there are $O(3N)$ degrees of freedom associated with the three independent components of the hydrodynamic force on each particle, and a single matrix-vector multiply is required to evaluate the hydrodynamic interactions. Therefore, computing the RPY approximation for the hydrodynamics requires $O(3N)$ operations. In contrast, the FSD approximation includes $O(17N)$ degrees of freedom: $3N$ from the hydrodynamic force, $8N$ from the traceless hydrodynamic couplet, and $6N$ from the near-field hydrodynamic force and torque. Additionally, three solves of the saddle point problem need to be evaluated to compute the particle displacements: one solve for the deterministic contribution to particle motion, and two for the RFD approximation to the Brownian drift. However, as previously discussed, the solution of the first RFD solve can be used as the initial guess in the second solve to reduce the second solve to a negligible cost. Therefore, from an operation cost perspective, only two solves are effectively required. Each solve requires approximately 30 iterations regardless of system size, so the total operation count for the FSD algorithm is $O(17N\times 30\times 2)=O(1020N)$ , which is 340 times that of the RPY approximation. With this estimate, the expected peak performance of the FSD algorithm is $10^{4}$ PTPS.

Figure 5.

Performance of the PSE algorithm at three levels of hydrodynamic approximation: RPY (light grey circles), RPY plus stresslets (medium grey squares) and FSD (solid black triangles). Performance of the SEASD algorithm of Wang & Brady (Reference Wang and Brady2016) is shown in open black triangles. The calculations are performed with a relative error tolerance of $10^{-3}$ . The performance is measured in particle time steps per second, which is equal to the total number of particles divided by the time required to computed one time steps worth of calculations. The performance of the RPY and RPY+stresslet algorithms is independent of  $\unicode[STIX]{x1D719}$ , so only one representative curve is shown. Because the FSD performance is a function of  $\unicode[STIX]{x1D719}$ , the performance for $\unicode[STIX]{x1D719}=\{0.1,0.2,0.3,0.4\}$ is computed separately and denoted by the solid, dashed, dashed-dot and dotted lines, respectively.

Figure 5 shows the performance of the PSE algorithm applied to dynamic simulations of hard sphere dispersions performed with a global relative error tolerance of $10^{-3}$ , which is typical in many applications (Wang & Brady Reference Wang and Brady2016). Performance curves are shown for the FSD algorithm presented here, as well as for the RPY level of hydrodynamic approximation both with and without the stresslet constraints. In all cases, the computational cost includes the cost of performing the Brownian sampling. The RPY performance, with and without the stresslet constraint, has been previously published and was evaluated on an NVIDIA Tesla K40 graphics card (Fiore et al. Reference Fiore, Usabiaga, Donev and Swan2017; Fiore & Swan Reference Fiore and Swan2018). The FSD performance presented here was evaluated on an NVIDIA GTX 980Ti graphics card. The GTX has a maximum memory bandwidth of $337~\text{GB}~\text{s}^{-1}$ , compared to $288~\text{GB}~\text{s}^{-1}$ for the Tesla, and a theoretical maximum performance of 6 teraflops for single precision calculations, compared to 5 teraflops for the Tesla. Whether the application is memory-bound or not, the GTX is expected to be approximately 20 % faster than the Tesla. In our case, we believe the computations are not memory bound. As previously shown, the performance of the RPY stresslet-constrained RPY algorithms is independent of particle volume, so only one curve is plotted for each of these methods.

In contrast to our previous work, the FSD algorithm displays a modest dependence on volume fraction, which is expected based on the convergence rates shown in figure 1(b). As shown in figure 5, there is a factor of 1.3 difference in performance between the $\unicode[STIX]{x1D719}=0.1$ and $\unicode[STIX]{x1D719}=0.4$ results with $N=10^{5}$ . The largest difference between volume fractions is a factor of 3 for $N=10^{3}$ . For the largest system sizes studied, $N=10^{5}$ , the FSD algorithm is about two times slower than RPY with stresslets. This is substantially better than what would be expected purely from an operation count perspective, even accounting for the difference in maximum theoretical performance between the GPUs. This discrepancy is likely due to the fact that the full computational bandwidth of the GPU is not being used by the RPY or stresslet systems. We have previously noted that the PTPS measured with the RPY algorithm continues to increase up to $4\times 10^{6}$ particles (Fiore et al. Reference Fiore, Usabiaga, Donev and Swan2017), which is evidence of bandwidth underutilization. Furthermore, it is obvious in figure 5 that the throughput of the FSD simulation increases systematically from $N=10^{2}$ to $N=10^{5}$ , which again indicates that more of the available GPU hardware is being utilized by the larger systems. Additionally, some minor improvements have been made to the far-field algorithms since the method for RPY with stresslet constraints was published, which may contribute slightly to the relative performance differences observed.

Nevertheless, when considered with the performance data for the RPY method, figure 5 clearly demonstrates that the algorithm implementation is optimal and exhibits linear scaling in the total computational complexity with the number of particles modelled. Our algorithm reaches $6\times 10^{4}$ PTPS with $N=10^{5}$ particles, which for the first time enables such large scale simulations of colloidal materials with full hydrodynamic interactions, including lubrication. This is made plain by comparing the relative performance of our FSD algorithm with the spectral Ewald accelerated Stokesian dynamics (SEASD) algorithm (Wang & Brady Reference Wang and Brady2016), which also leverages GPUs, the non-uniform fast Fourier transform (NUFFT) acceleration of Lindbo & Tornberg (Reference Lindbo and Tornberg2010), and a modestly effective preconditioner for the far-field matrix inversion. The SEASD performance data were obtained for bidisperse suspensions of spheres, for which the performance is expected to be only slightly lower than the monodisperse case considered here. Wang & Brady (Reference Wang and Brady2016) report that their GPU implementation achieves an order of magnitude speedup over an equivalent CPU implementation, which means that the two orders of magnitude difference between FSD and SEASD shown in figure 5 are due primarily, if not exclusively, to the algorithmic improvements presented in this work (Wang & Brady used an NVIDIA GTX 680 GPU, which has a maximum memory bandwidth of $192~\text{GB}~\text{s}^{-1}$ and a theoretical maximum performance of 3.25 teraflops for single precision calculations. Therefore, hardware differences amount to at most a factor of 2 difference between their performance and the FSD performance presented here). Specifically, the saddle point hydrodynamic formulation, which avoids nested iterative matrix inversions and for which highly effective preconditioners are developed, and the Brownian force decomposition, which leverages the PSE sampling technique in the far field and a novel preconditioner applied to Krylov subspace techniques in the near field, drastically improve calculation speeds. Furthermore, FSD maintains linear scaling of the computation cost up to $10^{5}$ particles, while SEASD performance deteriorates for systems larger than $10^{3}$ particles.