The art of coarse Stokes: Richardson extrapolation improves the accuracy and efficiency of the method of regularized stokeslets.

The method of regularized stokeslets is widely used in microscale biological fluid dynamics due to its ease of implementation, natural treatment of complex moving geometries, and removal of singular functions to integrate. The standard implementation of the method is subject to high computational cost due to the coupling of the linear system size to the numerical resolution required to resolve the rapidly varying regularized stokeslet kernel. Here, we show how Richardson extrapolation with coarse values of the regularization parameter is ideally suited to reduce the quadrature error, hence dramatically reducing the storage and solution costs without loss of accuracy. Numerical experiments on the resistance and mobility problems in Stokes flow support the analysis, confirming several orders of magnitude improvement in accuracy and/or efficiency.

Introduction: the method of regularized stokeslets

Flow problems associated with flagellar propulsion of cells, cilia-driven fluid transport, and synthetic microswimmers, are characterized by the inertialess regime of approximately zero Reynolds number flow, described mathematically—in Newtonian flow—by the Stokes flow equations,Typically, these conditions are associated with no-flux, no-penetration conditions on complex-shaped moving boundaries modelling cell surfaces and motile appendages. For a detailed introduction to the subject, see the recent text [1]. A range of mathematical and computational techniques are available to approach this problem; a computational method that has seen significant uptake and development over the last two decades is the method of regularized stokeslets, first described by Cortez [2] and subsequently elaborated for three-dimensional flow [3,4].

This technique can be viewed as a modification of the method of fundamental solutions and/or the boundary integral method for Stokes flow [5], the basis for which is the stokeslet [6] or Oseen tensor [7]:andThe pair of tensors S, P provide the solutions and p = (8π)−1P to the singularly forced Stokes flow equations,andwhere is the three-dimensional Dirac delta distribution and is a unit basis vector pointing in the k-direction. Equations (1.2) and (1.3) are singular when the source point and field point coincide. To facilitate numerical computation, the method of regularized stokeslets instead considers the Stokes flow equation with spatially smoothed point force,andwhere is a family of ‘blob’ functions approximating as ε → 0.

Several different choices for and associated regularized stokeslets have been studied; the most extensively used was presented in the original three-dimensional formulation of Cortez et al. [3],Developments focusing on the use of alternative blob functions to improve convergence include [8] (near-field) and, more recently, [9] (far-field).

The pressure and velocity as ε → 0; moreover the corresponding single layer boundary integral equation iswhere p = 1 for on or near B and p = 2 otherwise [3]. In equation (1.11) and below, summation over repeated indices in j = 1, 2, 3 or k = 1, 2, 3 is implied. The reduction to the single-layer potential is discussed by e.g. [3,5,10]; in brief, this equation can describe flow due to motion of a rigid body, or with suitable adjustment to f, the flow exterior to a body which does not change volume. A feature common to both standard and regularized stokeslet versions of the boundary integral equation is non-uniqueness of the solution f. This non-uniqueness occurs due to incompressibility of the stokeslet, i.e. provided the interior of B maintains its volume, then so that if f is a solution of equation (1.11) then so is f + an for any constant a. From the perspective of the original partial differential equation system, the non-uniqueness follows from the fact that the pressure part of the solution to equations (1.1) with velocity-only boundary conditions is determined only up to an additive constant. This issue is not dynamically important, and moreover the discretized approximations to the system described below result in invertible matrices.

Boundary integral methods have the major advantage of removing the need for a volumetric mesh, which both reduces computational cost, and moreover avoids the need for complex meshing and mesh movement. The key strength of the method of regularized stokeslets is in enabling the boundary integral method to be implemented in a particularly simple way: by replacing the integral by a numerical quadrature rule (abscissae, weight and surface metric), equation (1.11) may be approximated by,As is standard terminology in numerical methods for integral equations, we will refer to this as the Nyström discretization [11]. By allowing m = 1, …, N and j = 1, 2, 3, a dense system of 3N linear equations in 3N unknowns is formed. The diagonal entries when j = k and m = n are finite but numerically on the order of 1/ε, leading to (by the Gershgorin circle theorem) a well-conditioned matrix system.

The approach outlined above can be used to solve the resistance problem in Stokes flow, which involves prescribing a rigid body motion and calculating the force distribution, and hence total force and moment on the body. Once the force and moment associated with each of the six rigid body modes (unit velocity translation in the x direction, unit angular velocity rotation about x axis, for j = 1, 2, 3) are calculated, the grand resistance matrix A can be formed [5], which by linearity of the Stokes flow equations relates the force and moment to the velocity and angular velocity for any rigid body motion;For example, for a sphere of radius a centred at the origin, the matrix blocks are A = 6πμaI, and where I is the 3 × 3 identity matrix.

A closely related problem is the two-step calculation of the flow field due to a prescribed boundary motion; starting with prescribed surface velocities , first, the discrete force distribution F[n] is found by inversion of the Nyström matrix system; the velocity field at any point in the fluid can then be found through the summation,

The mobility problem is formulated by prescribing the total force and moment on the body (yielding six scalar equations) and augmenting the system with unknown velocity and angular velocity , which adds six scalar unknowns, so that a (3N + 6) × (3N + 6) system is formed. At a given time, these unknowns can be related to the evolution of the body trajectories (in terms of position and two basis vectors and ), through a system of nine ordinary differential equationswhich can be solved using available packages such as MATLAB’s ode45.

Finally, the swimming problem further prescribes the motion of cilia or flagella with respect to a body frame (typically, a frame in which the cell body is stationary), and often assumes zero total force and moment (neglecting gravity and other forces such as charge), again resulting in a (3N + 6) × (3N + 6) system. The key numerical features and challenges of the method of regularized stokeslets are exhibited by the resistance and mobility problems, which will therefore be our primary focus.

Convergence properties of the Nyström discretization

Equation (1.12) is subject to the O(ε) regularization error in the boundary integral equation, and the discretization error in the approximation of the integral. The integrand consists of a product: the slowly varying traction and the stokeslet kernel which is rapidly varying when . The error associated with discretization of the traction is at worst O(h), where h is the characteristic spacing between points. The dominant error for the stokeslet kernel can be shown to be(see [12], contained case, equation (2.7)).

Reducing the O(ε) regularization error by reducing ε therefore increases the O(ε−1h2) stokeslet quadrature error, necessitating refinement of the discretization length h. To reduce ε by a factor of R requires indicatively reducing h by a factor of , hence increasing the number of surface points and therefore degrees of freedom N by a factor of R. The cost of assembling the dense linear system then increases by a factor of R2, and the cost of a direct linear solver by a factor of R3. This calculation shows that, for example, improving from a 10% relative error to a 1% relative error may indicatively incur a cost increase of 1000 times. There are several approaches available already to address this issue, which involve a range of computational complexities: the fast multipole method [13], boundary element regularized stokeslet method [14] and the nearest-neighbour discretization [15], for example. In the next section, we will describe and analyse a very simple technique which alone, or potentially in combination with the above, improves the order of the regularization error, thereby enabling a coarser ε and hence alleviating the quadrature error. We will then briefly review an alternative ‘coarse’ approach, the nearest-neighbour method, a benchmark with similar implementational simplicity. Numerical experiments will be shown in the Results (§5), and we close with brief Discussion (§6).

Richardson extrapolation in regularization error

Consider the approximation of a physical quantity (e.g. moment on a rotating body) which has exact value M*. The value of this quantity calculated with discretization of size h and regularization parameter ε is denotedwhere E(ε) is the regularization error associated with the (undiscretized) integral equation, and E(h; ε) is the discretization error, which as indicated also has an indirect dependence on ε via the quadrature.

Recall thatandwhere E(h) is the error associated with the force discretization and E(h; ε) is the quadrature error. The analysis below will focus on the situation in which the regularization parameter ε is not excessively small, so that the quadrature error (h2/ε) is subleading and hence the discretization error has minimal dependence on ε, thus E(h; ε) ≈ E(h; ε0) for some representative value ε0. Writingwe may then expand,Evaluation of M(εℓ, h) for three values of εℓ in this range results in a linear system,Applying the matrix inverse,Hence, the estimateprovides an approximation to M* that has errorThis improvement in order of accuracy comes at a small multiplicative cost associated with solving the problem three times; however, as these are three independent calculations they are ideally placed to exploit parallel computing architecture, thus reducing the additional computational cost.

Comparison with the nearest-neighbour regularized stokeslet method

Before carrying out numerical experiments, we will briefly recap a different strategy to address the ε-dependence of the linear system size which we have developed and described recently, in order to provide a benchmark with similar implementational simplicity. The nearest-neighbour version of the regularized stokeslet method [16] aims to remove the ε-dependence of the linear system size. This change is achieved by separating the degrees of freedom for traction from the quadrature by using two discretizations: a ‘coarse force’ set for the traction and a finer set for the quadrature. If these sets are identical, the method reduces to the familiar Nyström discretization. In general, choosing N < Q leverages the fact that the traction is more slowly varying than the near-field of the regularized stokeslet kernel. Discretizing the integral equation (1.11) on the fine set givesBased on the observation that the traction and associated weighting are slowly varying, the method employs degrees of freedom F[n] in the neighbourhood of each point of the coarse discretization, so thatwhere ν[q, n] is a sparse matrix defined so that ν[q, n] = 1 if the closest coarse point to is , and ν[q, n] = 0 otherwise.

A detail that was not addressed in our recent papers ([15,17], for example) is that the closest coarse point to a given quadrature point may not be uniquely defined. Moreover, it is occasionally possible that, for sufficiently distorted discretizations, a coarse point may have no quadrature points associated with it at all, resulting in a singular matrix. In the former case, the weighting may be split between two or more coarse points, so that the sum of each row of ν[q, n] is still equal to 1. In the latter case, the coarse point may be removed from the problem, or (better) the quadrature discretization refined.

The approximation (4.2) leads to the linear system,The computational complexity of the system is given by the 3N × 3Q function evaluations required to assemble the stokeslet matrix, followed by the O(N3) solution of the dense linear system (for direct methods).

The nearest-neighbour method is subject to similar O(ε) regularization error and O(h) discretization error (where h is characteristic of the force point spacing) as the Nyström method. Analysis of the quadrature error associated with collocation [12] identifies two dominant contributions: For contained force and quadrature discretizations (i), the cost of quadrature is still an important consideration. Reducing ε by a factor of R, necessitates reducing by a factor of R, and hence increasing the number of quadrature points—and associated matrix assembly cost—by a factor of R. Therefore, any improvement to the order of convergence of the regularization error will result in a corresponding improvement in the reduction of quadrature error.

Contained case: Quadrature centred about a force point which is also contained in the quadrature set is subject to a dominant error term , where h is the spacing of the quadrature points; the Nyström method described above is a special case of this, with h = h;

Disjoint case: Quadrature centred about a force point which is not contained in the quadrature set is subject to a dominant error term O(h/δ)2 h), where δ > 0 is the minimum distance between the force point and quadrature points. This term does not appear in the Nyström method error analysis. The term is written in this form because δ is typically similar in size to h for a given quadrature set, so with a little care, h/δ behaves as a multiplicative constant.

However, when disjoint force and quadrature discretizations (ii) are employed, the nearest-neighbour method is able to entirely decouple the strong dependence of the degrees of freedom (tied only to h) on the regularization parameter ε and quadrature discretization h. The nearest-neighbour method, therefore, provides a relatively efficient and accurate implementation of the regularized stokeslet method that, with minor care in the construction of the discretization sets, can be used as a benchmark. In the following section, we will assess the Richardson extrapolation approach against analytic solutions for two examples of the resistance problem, and against the nearest-neighbour method for an example mobility problem.

Results

We now turn our attention to the application of Richardson extrapolation to a series of model problems, comprising the calculation of: For simulations (a) and (b), comparisons can be made to known exact solutions. For each test case (a–c), we compare the results of simulations using both the Nyström [Ny] and Nyström + Richardson [NyR] methods. For the latter, we choose extrapolation points . The choice of extrapolation rule is discussed further in appendix A.

the grand resistance matrix for a unit sphere;

the grand resistance matrix for a prolate spheroid; and

the motion of a torus sedimenting under gravity.

For each problem, we use the minimum distance between any two force points in the discretization as our comparative lengthscale h. For the [NyR] method, results are shown against the smallest value of the regularization parameter (ε1) used in the calculation. Simulations are performed with GPU acceleration (see [18]) using a Lenovo Thinkstation with an NVIDIA Quadro RTX 5000 GPU. Each of the test problems that we consider, however, are easily within the capabilities of more modest hardware.

The grand resistance matrix of a rigid sphere

Application of Stokes’ Law gives the force exerted by the translation of the unit sphere with velocity as , and the moment exerted by the unit sphere with rotational velocity as . From this, the grand resistance matrix A can be constructed as in equation (1.13). We solve equation (1.12) [Ny] and equations (1.12) and (3.8) [NyR] for unit translations and rotations about each axis to obtain the numerical approximation to A, . The relative error in the calculation is then given by the relationwhere ‖A‖2 denotes the 2-norm (‖A‖2 = sup‖Ax‖2/‖x‖2).

The unit sphere is discretized by projecting onto the six faces of a cube (figure 1a), with the number of scalar degrees of freedom (sDOF) shown plotted against the minimum spacing between points (h) in figure 1b. The relative error in calculating the grand resistance matrix as h and ε are varied is shown in figure 1c,d ([Ny] and [NyR], respectively). We report results for an identical range of ε (and h) for both methods, although as described in §3, the [NyR] method is specifically designed to exploit larger values of ε for which the quadrature error is small, so the [NyR] results with ε = 0.1–0.4 are most pertinent.

Figure 1

Relative error in calculating the grand resistance matrix for the unit sphere. (a) Sketch of the sphere discretization (orange dots). (b) The number of scalar degrees of freedom used in calculations as h is varied. (c) and (d) The relative error of the Nyström and Nyström + Richardson methods as ε and h are varied. (e) and (f) The same data plotted for each ε as h is varied.

The [Ny] method is found to achieve relative error for a select number of parameter pairs (ε, h). This is strongly dependent, however, on the ‘dip’ in error which appears as h is decreased for a given ε (evident in figure 1e) and is a consequence of the balance between the opposite-signed regularization and quadrature errors; the small h plateau remains above error for each choice of ε. By contrast, the [NyR] method is able to significantly reduce the error in the plateau (figure 1f ), resulting in sub- errors for ε as large as 0.2. Indeed with ε = 0.2, the range of values of h capable of producing acceptably accurate results extends from h = 0.00077 to h = 0.0076. As a result of the reduction in regularization error, brought about by the [NyR] extrapolation, this method is able to achieve a minimum relative error of compared with for the [Ny] method, and moreover, accurate performance no longer depends on a precise interplay between h and ε. In the simulations we performed, the [NyR] method was able to attain very accurate results ( error) in 250 s of walltime.

The grand resistance matrix of a rigid prolate spheroid

To assess the performance on a system involving a modest disparity of length scales, the second model problem is the calculation of the grand resistance matrix for a prolate spheroid of major axis length 5 and minor axis length 1. Moreover, prolate spheroids are often used as models for both entire microscopic swimming cells, and for their propulsive cilia and flagella, and so provide an informative test geometry. The exact solution in the absence of other bodies is well known (e.g. [19]). Details of the discretization of the prolate spheroid are provided in appendix B.1. A sketch of the discretization and plot of sDOF as h is varied are shown in figure 2a,b.

Figure 2

Relative error in calculating the grand resistance matrix for a prolate spheroid with major axis length a = 5 and minor axis length c = 1. (a) Sketch of the discretization (orange dots). (b) The number of scalar degrees of freedom used in calculations as h is varied. (c) and (d) The relative error of the Nyström and Nyström + Richardson methods as ε and h are varied. (e) and (f) The same data plotted for each ε as h is varied.

Similarly to the case of the unit sphere, the [Ny] method is able to achieve a minimum error of for the smallest ε in this study and a specific choice of h within the error dip (figure 2c,e). For each choice of ε, the error plateau for small h is at least relative error. Relatively large ε = 0.2, 0.4 yield error plateaus of and , respectively.

The [NyR] method also exhibits this dip phenomenon; however, the reduction in regularization error provided by the Richardson extrapolation (figure 2d,f) results in significantly reduced error plateaus of and (ε = 0.2 and 0.4, respectively), again being more robustly maintained over a larger range of values of h. For this test problem, the [NyR] method achieved error in 390 s of walltime.

The motion of a torus sedimenting under gravity

As a final test case, we simulate the mobility problem of a torus sedimenting under the action of gravity (for detailed set-up and discretization, see appendix B.2). In the absence of an exact solution to this problem, we compare the distance travelled in the vertical direction after the system (equations (B 6)–(B 8)) are solved for t ∈ [0, 98.7]. We compare the results obtained with the [Ny] and [NyR] methods with those from a simulation using the nearest-neighbour method ([NEAREST]) with a refined force discretization, disjoint force and quadrature discretizations and ε = 10−6. Figure 3a–c shows, respectively, example discretizations for the [Ny]/[NyR], and [NEAREST] methods, and the number of sDOF used in the [Ny] and [NyR] methods as h is varied. For the [NEAREST] simulation, a highly resolved system is constructed with 14 667 sDOF and 231 744 quadrature points.

Figure 3

A torus, with central radius R = 2.5 and tube radius r = 1, sedimenting under gravity. (a) Sketch of the Nyström discretization (orange dots). (b) Sketch of the [NEAREST] force (large orange dots) and quadrature (small green dots) discretizations. (c) The number of scalar degrees of freedom used in Nyström and Nyström + Richardson calculations as h is varied. (d) and (e) The z-position of the torus at t = 98.7 calculated with the Nyström and Nyström + Richardson methods as ε and h are varied. (f) and (g) The same data plotted for each ε as h is varied, with a dotted line showing the result using the nearest-neighbour method for comparison. (h) and (i) The error in z-position at t = 98.7 relative to the nearest-neighbour calculation. (j) and (k) The same data plotted for each ε as h is varied. The cross in (e,i) denotes a parameter combination for which results could not be obtained due to near-singularity of the linear system.

Figure 3d–g shows the convergence in z-position of the torus at t = 98.7 as both ε and h are varied. The relative difference between these results and the [NEAREST] simulation are shown in figure 3h–k. The error behaves similarly to the previous cases: while [Ny] achieves accurate results with specific combinations of ε and h; by contrast [NyR] at relatively large values of ε = 0.1–0.4 attains sub- error over an extended range of h values.

As anticipated by the error analysis, the advantage of [NyR] appears in the range of relatively coarse values, i.e. ε = 0.1–0.4. A solution could not be obtained when ε = 0.4 and h < 0.039, due to the matrix system with ε3 = 2ε = 0.8 becoming close to singular. For the choice of ε = 0.4, the [NyR] method attained an error of (compared with the result using [NEAREST]) in 144 s of walltime.

The results for the smallest choice of regularization parameter, ε = 0.01, are not converged with h, consistent with our analysis in §3 focusing on moderate values of ε for which the quadrature error is subleading.

Discussion

This article considered the implementation of the regularized stokeslet method, a widely used approach in biological fluid dynamics for computational solution of the Stokes flow equations. An inherent challenge is the strong dependence of the degrees of freedom on the regularization parameter ε, which necessitates an inverse-cubic relationship between the linear solver cost and the regularization parameter.

Here, we have investigated a simple modification of the widely used Nyström method, by employing Richardson extrapolation; performing calculations with three, coarse values of ε and extrapolating to significantly reduce the order of the regularization error. The method was compared with the original Nyström approach on three test problems: calculating the grand resistance matrices of the unit sphere and prolate spheroid, and simulating the motion of a torus sedimenting under gravity.

Investigation of these model problems has highlighted two significant phenomena, the first of which is well known but is worth repeating: (i) obtaining an acceptable level of error using the Nyström method is strongly dependent on being within the region where the (opposite-signed) regularization and quadrature errors exhibit significant cancellation, a phenomenon which has sensitive dependence on the discretization h as ε is varied. (ii) The improvement in the order of regularization error provided by Richardson extrapolation is able to significantly and robustly reduce errors for simulations with (relatively) large choices of ε, enabling highly accurate results with relatively modest computational resources. This advantage is (by design) only maintained for these coarse values of ε, so that the regularization error is subleading. Another approach which improves the order of convergence of the (important) local regularization error is given by Nguyen & Cortez [8], although the resulting regularized stokeslets may not be exactly divergence-free.

As discussed above, there are several existing approaches to improving the efficiency and accuracy of regularized stokeslet methods. The best approach in terms of strict computational complexity is the use of fast methods such as the kernel independent fast multipole method, which enables the approximation of the matrix-vector operation required for iterative solution of the linear problem [13,20], resulting in a O(N log N) method—although with somewhat greater implementational complexity. Another formulation is to borrow from the boundary element method developed for the standard singular stokeslet formulation [14], which has been applied to systems such as embryonic left-right symmetry breaking [21] and bacterial morphology [22]. The boundary element approach decouples the quadrature from the traction discretization and hence degrees of freedom of the system, enabling larger problems to be solved, although again at the expense of greater complexity through the need to construct a true surface mesh, with a mapping between elements and nodes. The nearest-neighbour discretization [15] retains much of the simplicity of the Nyström method, while separating the quadrature discretization from the degrees of freedom. Provided that the discretizations do not overlap, we still find this method to be an optimal combination of simplicity and efficiency. The Richardson approach does not avoid the need for the regularization parameter to not exceed the length scales characterizing the physical problem, for example the distance between objects. In this respect, the nearest-neighbour approach is advantageous because of its ability to accommodate smaller values of the regularization parameter.

In this work, we have focused on demonstrating how a numerically simple modification to the, already easy-to-implement, Nyström method can provide excellent improvements by employing coarse values of the regularization parameter ε. This approach can be considered complementary to the nearest-neighbour method in its coarse philosophy and style: both methods are figuratively coarse in their simplicity, and literally coarse in their approach of increasing numerical parameters. The Richardson approach allows increases in the regularization parameter; the nearest-neighbour approach allows increase the force discretization spacing h. Either method enables more accurate results to be achieved with greater robustness and for lower computational cost. Moreover, both have the advantage of being formulated in terms of basic linear algebra operations, and therefore can be further improved through the use of GPU parallelization with minimal modifications [18]. The choice of which method to use is a matter of preference; the Richardson approach has the advantage of being immediately adoptable by any group with a working Nyström code, alongside the repeated calculations being embarrassingly parallel; the nearest-neighbour approach has the advantage of completely removing the dependence of the system size on ε.

Accessible algorithmic improvements such as these provide the improved ability to solve a plethora of problems in very low Reynolds number hydrodynamics. Various potential application areas include microswimmers such as sperm [23,24], algae and bioconvection [25-29], mechanisms of flagellar mechanics [30,31], squirmers [32,33] and bio-inspired swimmers [34-36]. Stokeslet-based methods have been employed since the work of Gray and Hancock [6] in the 1950s; they continue to provide ease of implementation, efficiency and most importantly physical insight into biological systems.

Click here for additional data file.