We consider the flow of dense suspensions of rigid bodies in a Stokesian fluid. Such flows are difficult to compute numerically due to the presence of close-to-touching interactions, which may require a large number of unknowns to resolve sharply peaked surface forces, a large number of GMRES iterations to solve the discretized PDE, and an extremely small time step. A common way of dealing with these difficulties is to introduce a repulsion force between particles to prevent them from getting too close. However, this additional repulsion force is non-physical and may fundamentally alter the results of a simulation.

For suspensions of identical discs in 2D, we present a fast and accurate boundary integral method that mitigates these challenges without introducing artificial forces. Through precomputation, compression and interpolation of the close-to-touching part of the interaction operator, our method—termed interpolated compressed inverse preconditioning—efficiently handles close-to-touching interactions down to distances of \(10^{-10}\) with only a coarse discretization of the boundary. Additionally, we present a preconditioner that significantly reduces the number of GMRES iterations required to solve the Stokes mobility problem at each time step by effectively reusing the Krylov subspace from previous time steps. Coupled with high-order, adaptive time-stepping using spectral deferred correction, we are able to take larger time steps, mitigating the temporal stiffness resulting from close-to-touching interactions.

For a graphical description of this work, see: https://danfortunato.com/talks/ICIP-poster.pdf.

1 Stokes mobility problem

We consider \(\Nobj \) rigid discs \(\Omega = \{\Omega _1, \cdots , \Omega _{\Nobj }\}\) embedded in a Stokesian fluid. The fluid velocity in the exterior of \(\Omega \) is governed by the Stokes equations,

\begin{align} -\Delta \vec {u} + \nabla p &= 0 & \text {in } \Real ^2 \setminus \Omega , \label {e:stokes-eq0} \\ \nabla \cdot \vec {u} &= 0 & \text {in } \Real ^2 \setminus \Omega , \label {e:stokes-eq1} \end{align} where \(\vec {u}\) is the fluid velocity and \(p\) is the fluid pressure. Equations (1) and (2) denote the momentum balance and incompressibility constraints, respectively. In addition, we also assume that the fluid velocity at infinity decays to zero,

\begin{align*} \vec {u}(\vec {x}) &\rightarrow 0 & \text {as } |\vec {x}| \rightarrow \infty . \end{align*} Each disc has a net force \(\vec {F}_k\) and a net torque \(T_k\) acting about a point \(\vec {x}^c_k\). The discs undergo rigid body motion with the velocity \(\RigidBodyVel \) given by,

\begin{align*} \RigidBodyVel (x) &= \vec {v}_k + \omega _k (\vec {x} - \vec {x}^c_k)^{\perp } & \text {for all } \vec {x} \in \Omega _k, \end{align*} where \(\vec {v}_k\) is the translational velocity and \(\omega _k\) is the angular velocity of \(\Omega _k\) about the point \(\vec {x}^c_k\). A slip velocity boundary condition \(\vec {u}_s\) between the rigid bodies and the fluid is prescribed. Therefore, the fluid velocity on the boundary \(\partial \Omega \) is given by,

\begin{align*} \vec {u} &= \RigidBodyVel + \vec {u}_s & \text {on } \partial \Omega . \end{align*} In the mobility problem, we are given \(\vec {u}_s\), \(\vec {F}_k\), and \(T_k\) about \(\vec {x}^c_k\) for each \(\Omega _k\). The rigid body motion \(\RigidBodyVel \) (i.e., \(\vec {v}_k\) and \(\omega _k\) for each \(\Omega _k\)) is not known and must be determined.

Using the Stokes single- and double-layer potentials to represent the fluid velocity \(\vec {u}\) in terms of an unknown surface density \(\vec {\sigma }\), a boundary integral equation (BIE) for the Stokes mobility problem can be formulated as given in [1]:

\begin{align} \mathcal {K} \vec {\sigma } &= g & \text {on } \partial \Omega , \label {e:bie} \end{align} where \(\mathcal {K}\) is a second-kind boundary integral operator and \(g\) encodes the given slip velocity, net force, and net torque.

2 Close-to-touching interactions

Consider the model problem of two discs separated by a distance \(d\), with each disc discretized into a set of high-order Gauss–Legendre panels. The two disc problem serves as an effective pairwise preconditioner in a simulation with many close-touching discs. When the distance \(d\) between two discs gets small, the solution \(\sigma \) to the BIE in (3) becomes highly peaked. This requires an extremely fine discretization of the boundary in the close-to-touching region. We label the close-to-touching region as \(\Gamma _{2}\) and the remaining boundary as \(\Gamma _1 = \partial \Omega \setminus \Gamma _2\). Then, (3) can be discretized as a block linear system,

\begin{align} \begin{pmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end {pmatrix} \begin{pmatrix} \sigma _1 \\ \sigma _2 \end {pmatrix} = \begin{pmatrix} g_1 \\ g_2 \end {pmatrix}, \label {e:bie-block-lin} \end{align} where \(g_1\) and \(\sigma _1\) are the boundary conditions and unknowns on \(\Gamma _1\), \(g_2\) and \(\sigma _2\) are the boundary conditions and unknowns on \(\Gamma _2\), and \(K_{ij}\) represents a sub-block of the discretized operator \(\mathcal {K}\) that computes interactions from sources on \(\Gamma _j\) to targets on \(\Gamma _i\). Right preconditioning (4) with the block diagonal preconditioner \(\smash {\big (\begin {smallmatrix} I & 0 \\ 0 & K_{22}^{-1} \end {smallmatrix}\big )}\) yields the system

\begin{align} \begin{pmatrix} K_{11} & K_{12} K_{22}^{-1} \\ K_{21} & I \end {pmatrix} \begin{pmatrix} \sigma _1 \\ \overline {\sigma }_2 \end {pmatrix} = \begin{pmatrix} g_1 \\ g_2 \end {pmatrix} , \label {e:bie-right-precond} \end{align} where \(\overline {\sigma }_{2} = K_{22} \sigma _{2}\) is a new unknown on \(\Gamma _2\). While (5) may require fewer GMRES iterations to solve than (4), it still requires an excessively fine discretization in the close-to-touching region \(\Gamma _2\). Additionally, computing \(K_{22}^{-1}\) on the fly can be expensive, especially for problems with moving boundaries. However, one may show that \(\overline {\sigma }_{2}\) can be discretized on a coarse mesh and that the off-diagonal block \(K_{12} K_{22}^{-1}\) is low rank.

2.1 Compressing close-to-touching interactions

Since \(\Gamma _1\) and \(\Gamma _2\) are disjoint, the discretized boundary integral operators \(K_{12}\) and \(K_{21}\) are low rank, with the numerical rank independent of the distance \(d\). Hence, the column space of \(K_{21}\) is comprised of smooth functions that can be discretized using piecewise polynomials on a coarse mesh. From (5), we have \(\overline {\sigma }_{2} = g_{2} - K_{21} \sigma _{1}\); therefore, \(\overline {\sigma }_2\) is smooth whenever \(g_{2}\) is smooth and it can be discretized on a coarse mesh. Since \(K_{12}\) is low rank, so is \(K_{12} K_{22}^{-1}\) (with numerical rank independent of \(d\)). There are several ways of constructing a compressed representation for \(K_{12} K_{22}^{-1}\). To retain the boundary integral structure and allow for acceleration by the fast multipole method (FMM), we use a representation of the form \(K_{12} R\) where \(R\) is a low-rank operator such that \(K_{12} R \approx K_{12} K_{22}^{-1}\), up to a given numerical tolerance. We now describe the numerical construction of \(R\).

Consider two different panelizations of \(\Gamma _2\): a fine mesh where the panels on each disc are refined dyadically towards the closest point between the discs, and a coarse mesh with a small number of uniformly sized panels on each disc. We denote quantities on the coarse mesh with a superscript “\(c\)”; all other quantities are assumed to live on the fine mesh. Define the prolongation operator \(P\) that interpolates data from the coarse mesh to the fine mesh, and diagonal matrices \(W_f\) and \(W_c\) containing the weights for smooth integration on the fine and coarse meshes, respectively. Then, \(W_c^{-1} P^{T} W_f\) computes an \(L^2\) projection from the fine mesh to the coarse mesh.

Assuming that the boundary data \(g_2\) is smooth and therefore representable on the coarse mesh, we have

\begin{align} g_{2} &= P \, g^{c}_{2}, \label {e:coarse-g} \\ \overline {\sigma }_{2} &= P \, \overline {\sigma }^{c}_{2}. \label {e:coarse-sigma} \end{align} Since \(K_{12}\) and \(K_{21}\) are low rank, they can be approximated accurately by their coarse discretizations,

\begin{align} K_{12} &= K^{c}_{12} W_c^{-1} P^{T} W_f, \label {e:coarse-k12} \\ K_{21} &= P K^{c}_{21}, \label {e:coarse-k21} \end{align} Substituting (6)–(9) in (5), we obtain

\begin{align} \begin{pmatrix} K_{11} & K^{c}_{12} R \\ K^{c}_{21} & I \end {pmatrix} \begin{pmatrix} \sigma _1 \\ \overline {\sigma }^{c}_2 \end {pmatrix} = \begin{pmatrix} g_1 \\ g^{c}_2 \end {pmatrix} \end{align} where \(R = W_c^{-1} P^{T} W_f K_{22}^{-1} P\). This definition of \(R\) is used in the RCIP (recursively compressed inverse preconditioning) method [2]. For an order-\(p\) fine mesh with \(\BigO {\log d}\) levels of refinement, direct construction of \(R\) takes \(\BigO {p^3 (\log d)^3}\) operations; the RCIP method provides a faster algorithm to construct \(R\) without directly computing \(K_{22}^{-1}\), taking \(\BigO {p^3 \log d}\) operations. Our main result—termed ICIP (interpolated compressed inverse preconditioning)—instead constructs \(R\) through precomputation and interpolation, requiring only \(\BigO {p^2}\) work.

2.2 Interpolated compressed inverse preconditioning (ICIP)

Constructing \(R = R(d)\) each time for a different value of \(d\) can be expensive since it requires computing \(K_{22}^{-1}\). Instead, we construct a polynomial interpolant for \(R(d)\) over a range \(d \in [d_\text {min}, \, d_\text {max}]\) (where \(0<d_\text {min}<d_\text {max}\)). Then for any value of \(d\) in the interval, we construct \(R(d)\) through entrywise interpolation. We use Chebyshev polynomials in \(\log d\) as our interpolation basis, i.e.,

\[ [R(d)]_{ij} \approx \sum _{k=0}^q [R_k]_{ij} T_k(\log d) \]

where \([R_k]_{ij}\) is the \(k\)th Chebyshev coefficient for the \(ij\)th entry of \(R(d)\). For accurate interpolation of \(R(d)\) over a large dynamic range of \(10^{-10} < d < 10^{-1}\), only a moderate interpolation order of \(q=32\) is required. After an offline precomputation to generate \(\{R_k\}_{k=0}^q\), constructing \(R(d)\) at each time step costs \(\BigO {p^2 q}\) operations.

3 Accelerating timestepping with subspace recycling

While the two-disc preconditioner is effective at lowering the number of GMRES iterations induced by close-to-touching interactions, a significant number of GMRES iterations may still be required at each time step for problems with many discs. To ameliorate this effect, we propose a preconditioner which effectively reuses the Krylov subspace from previous time steps.

After the \(k\)th iteration of GMRES, the Krylov matrix is given by \(X = [b \; Ab \; \cdots \; A^{k-1}b]\). Let \(QR = AX\) be the QR decomposition of \(AX\). Then the matrix \(P\) given by

\[ P = I - QQ^T + XR^{-1}Q^T \]

has the following properties:

\[ \begin {aligned} PAx = x & \text { for all } x \in \operatorname {span}(X), \\ Py = y & \text { for all } y \perp \operatorname {span}(X). \end {aligned} \]

Hence, \(P\) effectively reuses the given Krylov subspace \(X\) when used as a preconditioner in a Krylov method. In a high-order time-stepping scheme based on spectral deferred corrections, this preconditioner can drastically reduce the number GMRES iterations by recycling the Krylov subspace between time steps.

References

• [1]  C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow, Cambridge University Press, 1992.

• [2]  J. Helsing, Solving integral equations on piecewise smooth boundaries using the RCIP method: a tutorial, 2022.