PaperlayerView exact revision on arXiv

Source-linked AI summary

Tolerance-driven close evaluation of the Stokes double layer potential on axisymmetric surfaces

David Krantz, Pritpal Matharu, Anna-Karin Tornberg

arXiv:2609.11778v1math.NAcs.NA

TL;DR

Close to a particle surface, standard quadrature loses accuracy and ever-finer upsampling becomes impractical. The paper presents a tolerance-driven classifier for choosing among standard, upsampled, and special quadrature; across tested interactions, the prescribed tolerance is met nearly always.

  • Problem

    Near-surface targets make the double-layer integrand sharply peaked, so standard quadrature loses accuracy and increasingly large upsampling factors become impractical.

    Source: §1 · p. 2

  • Method

    A fast classifier selects the least costly standard, upsampled, or special quadrature estimated to meet tolerance, using tabulated geometry indicators and an inexpensive density modifier.

    Source: §1.1 · p. 4, §6 · p. 29

  • Results

    The prescribed tolerance is met for nearly all tested target–particle interactions, with errors within a modest factor of tolerance in the few remaining cases.

    Source: §6 · p. 29

  • Takeaways & Limitations

    The framework is integrated into a Stokes mobility and resistance solver for off-surface interactions during both the boundary integral solve and subsequent flow evaluation.

    Source: §6 · p. 29

    The error indicators retain only the closest singularity pair, although additional pairs can contribute to quadrature error.

    Source: §3.1 · p. 9

Abstract

from arXiv · show

We consider boundary integral methods for Stokes mobility and resistance problems involving smooth axisymmetric particles. A primary numerical challenge is the accurate and efficient evaluation of layer potentials at off-surface points close to particle surfaces. We present a tolerance-driven workflow for evaluating the Stokes double layer potential at such evaluation points (targets) to prescribed accuracy while avoiding unnecessary computational cost. For each target--particle interaction, a fast classifier selects the least costly option estimated to meet the tolerance among standard, upsampled, and special quadrature. Geometry-dependent unit-density error indicators are precomputed and tabulated in reduced cylindrical coordinates, then combined on the fly with a local layer-density modifier, making its cost negligible relative to evaluating the potential. Targets requiring special quadrature are treated using a stabilized version of singularity swap surface quadrature: the periodic azimuthal integral is evaluated first using translated singularity swap quadrature to prevent severe cancellation near the surface, followed by adaptive Gauss--Legendre quadrature in the polar direction guided by error indicators. We integrate this workflow into a boundary integral solver with precomputed quadrature by expansion for on-surface self-interactions and demonstrate the workflow's performance for challenging configurations of spheroidal particles. Numerical results show that target classification is highly accurate. The prescribed tolerance is met for nearly all target--particle interactions and the error remains within a modest factor of the tolerance in the few remaining cases. Although the experiments focus on the Stokes double layer potential for spheroids, the off-surface framework applies to general smooth axisymmetric surfaces and can be easily adapted to other Stokes layer potentials.

1 Introduction

Close evaluation of the Stokes double layer potential is difficult because its off-surface integrand becomes sharply peaked near a particle, while quadrature choices trade accuracy against cost. The section introduces a tolerance-driven classifier and fast error indicators, plus stabilized special quadrature, for smooth axisymmetric particles.

1 Introduction

  • Boundary integral methods discretize only particle surfaces, but accurate evaluation of singular and nearly singular layer potentials remains a central numerical challenge.The paper focuses on the Stokes double layer potential at off-surface targets close to particle surfaces.Source: §1 · p. 2
  • Near a particle surface, the double-layer integrand becomes sharply peaked, so standard quadrature loses accuracy and upsampling eventually becomes computationally impractical.The appropriate method depends on target location, particle geometry, layer density, and the prescribed tolerance for each target–particle interaction.Source: §1 · p. 2
  • Earlier distance-based selection uses calibrated thresholds, while error indicators can estimate quadrature error directly but require costly target-dependent computations.Reliable classification matters because missed corrections compromise accuracy and unnecessary corrections add cost.Source: §1 · p. 3

1.1 Contributions and outline

  • Geometry-dependent unit-density error indicators are precomputed over reduced target coordinates and combined with local density information online at very low cost.Axisymmetry reduces the storage needed for the indicator tables.Source: §1.1 · p. 4
  • The special-quadrature component stabilizes S3Q with TSSQ in the azimuthal treatment to reduce cancellation as targets approach the surface.The workflow also precomputes SSQ weights and selects S3Q parameters adaptively from the tolerance.Source: §1.1 · p. 4
  • The framework meets the prescribed tolerance for nearly all tested target–particle interactions, with errors within a modest factor of tolerance in the few remaining cases.The complete workflow is demonstrated for spheroids and includes precomputed QBX for on-surface self-interactions.Source: §1.1 · p. 4
  • The off-surface framework applies to general smooth axisymmetric surfaces and can be adapted to other Stokes layer potentials.The complete solver workflow is demonstrated for spheroids, while the classifier is also tested on a capsule-shaped geometry.Source: §1.1 · p. 4

2 Problem setting

The section formulates exterior Stokes flow around rigid axisymmetric particles as a boundary integral problem and describes its surface discretization. It uses QBX for on-surface self-interactions, with precomputed weights for repeated evaluations.

2.1 The Stokes problem and boundary integral formulation

  • Mobility and resistance formulations yield a Fredholm second-kind equation for the density using a completed Stokes double-layer representation.The completion flow accounts for nonzero net particle forces and torques.Source: §2.1 · pp. 5–6
  • Nyström discretization enforces the equation at surface quadrature nodes, producing a dense system solved iteratively with GMRES.After solving for the density, particle quantities and the fluid velocity are recovered from its functionals and the representation formula.Source: §2.1 · p. 6

2.2 Surface discretization and quadrature

  • The base surface grid combines Gauss–Legendre nodes in the polar direction with equispaced trapezoidal nodes in the azimuthal direction.Upsampling interpolates the density and geometry onto a finer grid; it assumes they are resolved on the base grid and increases direct evaluation cost by κ^2.Source: §2.2 · pp. 6–7

2.3 QBX for on-surface evaluation

  • For each surface node, QBX places expansion centers on both sides of the boundary, and averaging the one-sided limits gives the principal-value integral.The two-sided construction is used for on-surface evaluations.Source: §2.3 · p. 7
  • Precomputed target-specific QBX weights reduce online evaluation to applying the weights to the current density at O(N) cost per target.Axisymmetry reduces storage to O(nθN), and weights can be reused for particles differing only by rigid-body motion.Source: §2.3 · p. 7
  • QBX handles each particle’s on-surface self-contribution, while off-surface interactions use the quadrature-selection procedure.For arbitrary off-surface targets, the complete density-to-potential map cannot be precomputed.Source: §2.3 · p. 7

3 Fast precomputed error indicators and target classification

This section develops inexpensive, axisymmetry-based error indicators for classifying off-surface target–particle interactions among standard, upsampled, and special quadrature, using quadrature-error estimates adapted to the Stokes double layer. The indicators closely track measured quadrature error at low evaluation cost, with local rather than global density scaling yielding sharper classifications.

3 Fast precomputed error indicators and target classification

  • The classifier selects standard quadrature, an upsampling level, or special quadrature according to an inexpensive estimate of error magnitude.Its purpose is to distinguish among these accuracy regimes for each target, particle, and prescribed tolerance.Source: §3 · p. 8

3.1 Error indicators for line potentials on curves

  • The leading nearby singularity yields asymptotic quadrature-error indicators that predict accuracy loss at moderate orders, while retaining only the closest singularity pair can capture error magnitude within a modest factor.Other singularity pairs decay exponentially faster with quadrature order under comparable numerator and geometry prefactors.Source: §3.1 · p. 9

3.2 Error indicators for layer potentials on surfaces

  • Surface error indicators combine estimated polar and azimuthal quadrature errors, separating precomputable unit-density terms from local density information.The unit-density indicators depend on kernel, geometry, discretization, and target, but not on layer density.Source: §3.2 · pp. 10–11

3.3 Error indicators for the tensor-valued stresslet kernel

  • For the Stokes double layer, the stresslet kernel uses singularity exponent p = 5/2, and componentwise density scaling is collapsed to one classification indicator per input density component.This avoids storing and interpolating the full 3 × 3 collection of unit-density indicators at runtime.Source: §3.3 · p. 12

3.4 Tabulation in reduced target coordinates

  • Axisymmetry enables precomputation of unit-density indicators in reduced target coordinates, with tables reusable across targets and densities for a fixed geometry and discretization.The rapidly varying indicator magnitudes are tabulated in logarithmic form and interpolated at arbitrary targets.Source: §3.4 · pp. 12–13
  • For tolerance ϵ, the classifier selects the smallest available upsampling factor κ whose estimated error is at most ϵ, using special quadrature if none qualifies.κ = 1 denotes standard quadrature; the fast indicator combines precomputed unit-density indicators with a local density modifier.Source: §3.4 · p. 14, §3.5 · p. 15

3.5 Performance of the fast error indicator

  • The evaluation uses median timings over three runs, and measured error is the Euclidean norm of the velocity difference from a large-upsampling reference.Experiments were run in MATLAB on a 3.4 GHz quad-core Intel i7-6700 CPU.Source: §3.5 · p. 15

3.5.1 Accuracy and efficiency

  • For the capsule test, the polar root required by the density modifier is recovered by interpolation because no closed-form root is available.The experiment uses R = 1, L = 6, C = 3, and a 50 × 50 surface discretization, with a smooth density varying by approximately five orders of magnitude.Source: §3.5.1 · p. 15
  • With a strongly localized density, the fast indicator’s contours closely follow measured quadrature-error levels on the capsule geometry.Figure 4 compares measured error for standard quadrature with indicator contours at levels from 10^-10 to 10^-2.Source: §3.5.1 · p. 16, Figure 4 · p. 16
  • For 10^5 targets, indicator evaluation takes less than 0.1 seconds, exceeding 10^6 targets per second with cost growing linearly after low-count overhead.The authors report that this cost is negligible compared with the upsampled or special quadrature the indicator activates.Source: §3.5.1 · p. 16, Figure 4 · p. 16

3.5.2 Effect of the density modifier

  • For the Stokes double layer potential (3) for two close-to-touching spheroids, the local-density-modified indicator closely follows measured quadrature error, while global L∞-based scaling overestimates it by approximately two orders of magnitude and unnecessarily enlarges regions assigned to upsampled or special quadrature.Source: §3.5.2 · pp. 16–17, Figure 5 · p. 17

4 The stabilized S3Q method

The method evaluates close-target interactions with stabilized S3Q, using TSSQ for azimuthal evaluation and adaptive Gauss–Legendre quadrature for the polar integral. The implementation uses a global polar discretization and does not exploit S3Q’s possible panel localization.

4 The stabilized S3Q method

  • For targets classified as requiring special quadrature, the Stokes double layer potential is evaluated with stabilized S3Q.The method overview includes density subtraction before evaluation.Source: §4 · p. 17
  • TSSQ stabilizes the azimuthal evaluation against severe cancellation, while composite Gauss–Legendre quadrature adaptively refines the polar integral.The procedure also constructs target-specific azimuthal weights through precomputation and one-dimensional interpolation.Source: §4 · p. 17
  • The implementation uses the global polar discretization required by its current QBX implementation, so it does not localize S3Q to panels near the target.S3Q can be applied to a multi-panel base discretization and restricted to nearby panels.Source: §4 · p. 17

4.1 Overview of the close-evaluation procedure

The close-evaluation procedure subtracts the density at the closest surface-grid node, then evaluates the resulting potential as an iterated integral. It applies azimuthal quadrature first and adaptively refines the polar quadrature.

4.1 Overview of the close-evaluation procedure

  • The resulting double layer potential is evaluated as an iterated integral, with the azimuthal integral taken first.The surface Jacobian is included in the integrand.Source: §4.1 · pp. 17–18
  • Subtracting the density at the closest surface-grid node reduces the sharp local kernel peak without changing the exact exterior potential.The subtracted density vanishes at that node.Source: §4.1 · pp. 17–18
  • At each polar node, the trapezoidal rule is used when accurate enough; otherwise SSQ or TSSQ evaluates the azimuthal integral.The polar values are combined by composite Gauss–Legendre quadrature on adaptively refined panels, with surface data and density interpolated from the original grid.Source: §4.1 · p. 18

4.2 Stabilized azimuthal evaluation

The azimuthal evaluation uses singularity swap quadrature, with a modified Fourier basis and stable constant-coefficient computation to address close-target cancellation. Indicators and a cancellation estimate select the azimuthal rule at each polar node.

4.2 Stabilized azimuthal evaluation

  • SSQ pairs the distance denominator with a simpler function sharing its complex-plane zeros, isolating the near singularity for quadrature.The smooth numerator is approximated in a suitable basis, with basis integrals evaluated for the chosen setting.Source: §4.2 · pp. 18–19

4.2.1 Standard Fourier formulation

  • The standard S3Q azimuthal formulation uses a Fourier basis and the selected azimuthal root to construct SSQ weights.This is the azimuthal SSQ formulation used in the original S3Q method [29, Section 4].Source: §4.2.1 · p. 19

4.2.2 Stabilization using the modified Fourier basis

  • For close targets, near-vanishing stresslet numerators and large basis integrals can cause severe cancellation in standard SSQ.TSSQ uses a modified Fourier basis whose nonconstant functions vanish at the root’s real part, and computes the constant coefficient using the known kernel and geometry with interpolated density.Source: §4.2.2 · pp. 19–20
  • At each polar node, a trapezoidal-rule indicator determines whether correction is needed; a cancellation estimate then selects SSQ or TSSQ.The component-wise error and cancellation estimates are combined by taking their respective maxima to choose a common rule for all output components.Source: §4.2.2 · p. 20

4.3 Efficient evaluation of azimuthal SSQ weights

Target-specific azimuthal weights let the method apply the integral directly to sampled data rather than explicitly forming basis coefficients. Precomputed Chebyshev interpolants reduce repeated basis-integral evaluation.

4.3 Efficient evaluation of azimuthal SSQ weights

  • Target-specific quadrature weights evaluate the azimuthal integral directly from sampled data, without explicitly forming the basis coefficients.The modified Fourier basis uses a stabilized adjoint formulation and FFT acceleration.Source: §4.3 · p. 20
  • Chebyshev interpolants of the basis-integral quantities are precomputed and evaluated with the second barycentric formula.The quantities depend smoothly on χ, and only nonnegative modes are tabulated.Source: §4.3 · p. 20
  • The interpolation tables keep absolute error below 10^-10 and are reused across requested quadrature tolerances.Reducing interpolation accuracy provides only negligible savings in total evaluation time.Source: §4.3 · p. 20

4.4 Adaptive polar Gauss–Legendre quadrature

The polar integral is evaluated with composite Gauss–Legendre quadrature on an adaptive panel grid, refined according to target-dependent error indicators. A minimum grid resolves the base geometry and density before refinement, while panel order balances accuracy control against computational cost.

4.4 Adaptive polar Gauss–Legendre quadrature

  • Azimuthal integrals use the trapezoidal rule, SSQ, or TSSQ, while smooth χ-dependent basis quantities support efficient one-dimensional interpolation.The fixed panel order trades fewer panels at higher order against more azimuthal evaluations per panel.Source: §4.4 · p. 21, Figure 6 · p. 21, Effect of panel order · p. 22

Minimum refinement depth.

  • The initial grid is uniformly bisected until it contains at least nθ nodes, resolving the base-grid geometry and density before target-specific refinement.Geometry and density values at panel nodes are interpolated from the original grid using precomputed barycentric Lagrange matrices.Source: Minimum refinement depth · p. 21

Target-dependent refinement.

  • Composite Gauss–Legendre panels use a target-dependent indicator combining square-root and reduced-order near-singular factors, with scaling for the remaining integrand factors.The logarithmic factor is not assigned a separate indicator because its resolution requirements resemble those of the square-root factor; the shared refinement criterion uses the maximum across output components.Source: §4.4 · p. 20, Target-dependent refinement · pp. 21–22

Effect of panel order.

  • At nGL = 8, 16, and 32, the indicators track quadrature error, but lower-order rules provide finer error control, motivating nGL = 8 for later experiments.For the closer target, higher-order rules can reduce error by several orders of magnitude in one bisection, potentially overshooting the requested tolerance.Source: Effect of panel order · p. 22, Figure 7 · p. 22

5 Numerical experiments

The experiments test the tolerance-driven quadrature workflow on spheroids, from single-particle evaluations to close particle gaps and a 100-particle periodic configuration. Across these settings, special quadrature controls errors near surfaces, and the classifier generally selects the least costly treatment that meets tolerance.

5 Numerical experiments

  • The experiments combine a fast classifier with stabilized S3Q and test post-processing evaluations against off-surface QBX references on spheroids.The sequence ranges from a single spheroid and two-particle tests to close-to-touching particles and a larger periodic system.Source: §5 · pp. 22–23

5.1 Setup

  • The common hierarchy uses QBX for on-surface self-interactions and selects standard quadrature, the smallest adequate upsampling factor, or S3Q for off-surface interactions.The experiments use two prolate spheroid types; mobility problems are considered in Sections 5.2–5.5 and a resistance problem in Section 5.6.Source: §5.1 · p. 23

5.2 Single spheroid and validation

  • For a single spheroid, S3Q stays below the prescribed tolerance across target sets, while standard and upsampled quadrature lose essentially all digits for d < 10^-3.Increasing κ extends accurate evaluation closer to the surface but does not prevent eventual error growth, even at κ = 6.Source: §5.2 · p. 23, Figure 8 · p. 24

5.3 Stabilization of S3Q near the surface

  • Across distances 10^-8 to 10^-3, stabilized S3Q removes standard S3Q’s O(1/d^3) error growth and achieves errors near or below tolerance.At tolerance 10^-6, the maximum error is slightly above tolerance in the closest regime, but about 10 digits more accurate than standard S3Q.Source: §5.3 · p. 24, Figure 9 · p. 25

5.4 Evaluation of target classification and error control

  • Shifting assigned quadrature down one level raises the fraction of affected targets exceeding tolerance from 0.02% to more than 95%, showing that classifications are nearly cost-minimal.For the tolerance 10^-6 example, reducing κ from 4 to 3 puts 4246 of 4447 points (95.5%) above tolerance.Source: §5.4 · pp. 25–26, Table 1 · p. 26, Figure 10 · p. 26

5.5 Close-to-touching spheroids

  • In the close-to-touching case, the classifier activates S3Q for targets as near as 1.9 × 10^-5 to a surface; mean and maximum errors are 8.46 × 10^-8 and 9.4 × 10^-6.The spheroids have a tip gap of 2.5 × 10^-3, and their opposing motion produces sharply varying density near the facing tips.Source: §5.5 · pp. 26–27, Figure 11 · p. 27

5.6 Triply periodic flow through a spheroid cluster

  • For 100 periodic spheroids, the maximum quadrature error remains just below the requested tolerance, while the area-normalized L2 error is smaller.On five check surfaces, the maximum no-slip error is 1.37 × 10^-5, corresponding to about five digits in the maximum norm.Source: Figure 13 · p. 28, §5.6 · p. 28, Figure 12 · p. 28

6 Conclusions

The paper presents a tolerance-driven framework that chooses among quadrature treatments for off-surface Stokes double-layer evaluation near smooth axisymmetric particles. Tests meet the prescribed tolerance for nearly all interactions, while the solver's on-surface QBX implementation is scoped to spheroids with a global tensor-product discretization.

6 Conclusions

  • The classifier selects the least costly estimated-adequate treatment—standard quadrature, selected upsampling, or special quadrature—using precomputed geometry indicators and an inexpensive local density modifier.Classification cost is negligible compared with subsequent potential evaluation.Source: §6 · p. 29
  • For special quadrature, TSSQ evaluates the azimuthal integral to prevent severe cancellation, while error indicators guide adaptive polar Gauss–Legendre refinement.S3Q's automatic parameter selection allows error control to the prescribed tolerance.Source: §6 · p. 29
  • At tolerance 10^-6, only 0.02% of 59,696 interactions assigned upsampling or special quadrature exceed tolerance, versus more than 95% under the next less expensive rule.Across the experiments, the tolerance is met for nearly all tested interactions, with remaining errors within a modest factor; the classifier almost always selects the least costly adequate treatment.Source: §6 · p. 29
  • The framework is integrated into a Stokes mobility and resistance solver, but its on-surface QBX implementation supports spheroids with a global tensor-product discretization.The classifier and S3Q extend to other smooth axisymmetric geometries, and the off-surface framework can be adapted to other Stokes layer potentials.Source: §6 · p. 29

Appendix A Numerical evaluation of one-dimensional error indicator integrals

Appendix A gives implementation formulas for the one-dimensional integrals used in surface error indicators. It approximates these integrals with locally corrected root approximations and low-order Gauss–Laguerre quadrature.

Appendix A Numerical evaluation of one-dimensional error indicator integrals

  • The implementation formulas follow [43, Section 6.1], building on the earlier treatment in [42, Section 6.1].Source: Appendix A · p. 29
  • The indicator integrals use one-dimensional quadrature factors and selected roots of the squared-distance function in the azimuthal and polar directions.The polar root is mapped to the Gauss–Legendre reference interval [−1, 1].Source: Appendix A · pp. 29–30
  • Because the integrals are generally unavailable in closed form, the formulas combine locally corrected roots with low-order Gauss–Laguerre quadrature.Each of the four half-line integrals uses an 8-point rule, sufficient in practice because the indicators serve only target classification.Source: Appendix A · p. 30
  • A local quadratic approximation estimates the polar root, then correction by the root at the closest surface node better captures the indicator's peak magnitude.The construction starts from the target's closest surface node and its local geometry.Source: Appendix A · p. 30
  • The azimuthal direction uses the same corrected-root construction, with K = ∥γ̄θ∥/∥γ̄ϕ∥ capturing exponential indicator decay away from the closest node.The azimuthal construction yields a corrected local approximation eϕ0.Source: Appendix A · p. 30

Appendix B Analytical formulas for roots of the squared-distance function

Appendix B collects analytical roots of the squared-distance function used by the error indicators. The azimuthal root is available for general axisymmetric surfaces, while analytical polar roots are given for spheres and spheroids.

Appendix B Analytical formulas for roots of the squared-distance function

  • For a general axisymmetric surface, the polar root is analytical only in special cases such as spheres and spheroids.The appendix formulates the surface parametrization and target in cylindrical coordinates.Source: Appendix B · pp. 30–31
  • For general axisymmetric surfaces, the azimuthal squared-distance equation has complex roots expressed using λ, with a specified square-root branch convention.The appendix defines atan2 as the argument of a complex number and specifies branch cuts for the square roots.Source: Appendix B · p. 31
  • For spheres, the polar roots are given in terms of atan2 and a logarithm, subject to the stated geometric condition.The sphere case assumes equal surface parameters a(ϑ) = c(ϑ) = a.Source: Appendix B · p. 31
  • For spheroids, the polar squared-distance equation is quartic, and the relevant root is the one with the smallest absolute imaginary part.The lemma assumes a ≠ c and provides expressions involving Δ, τ, and d².Source: Appendix B · p. 31
Loading 2609.11778v1…