Extrapolated Regularization Method

• Extrapolated regularization method is a numerical technique that regularizes nearly singular kernels and employs extrapolation to achieve high-order accurate evaluations of boundary integrals.
• It uses multiple smoothing parameters with standard quadrature and local linear systems to eliminate leading errors and maintain O(h^5) accuracy even on challenging surfaces.
• Enhanced by parallel processing and treecode algorithms, the method significantly reduces computational cost and is well-suited for large-scale simulations in fluid dynamics and potential theory.

The Extrapolated Regularization Method encompasses a class of techniques in numerical analysis and scientific computing designed to accurately evaluate nearly singular boundary integrals, particularly for surface integrals arising in Stokes flow and related potential theory problems. These methods employ kernel regularization with a smoothing parameter and apply numerical extrapolation to obtain high-order accurate values for integrals evaluated at points very close to, or on, the integration surface. Recent developments have also focused on computational optimization, enabling these methods to efficiently handle large-scale problems and challenging geometries.

1. Principles of the Extrapolated Regularization Method

When evaluating single- and double-layer surface integrals for Stokes flow,

$$u_i(\mathbf{y}) = \frac{1}{8\pi}\int_{\Gamma} S_{ij}(\mathbf{y}, \mathbf{x}) f_j(\mathbf{x})\,dS(\mathbf{x}) \quad\text{and}\quad v_i(\mathbf{y}) = \frac{1}{8\pi}\int_{\Gamma} T_{ijk}(\mathbf{y}, \mathbf{x}) q_j(\mathbf{x}) n_k(\mathbf{x})\,dS(\mathbf{x}),$$

standard quadrature fails when the target $\mathbf{y}$ is on or near the surface $\Gamma$, due to the nearly singular nature of the kernels.

The extrapolated regularization method circumvents this by:

• Regularizing the kernel: Replace $1/r^p$ singularities with $s_i(r/\delta)/r^p$, where $r=|\mathbf{y} - \mathbf{x}|$, $\delta$ is the smoothing parameter, and $s_1, s_2, s_3$ are specific regularizing functions (e.g., involving the error function, $\operatorname{erf}$), ensuring smooth near-field behavior and rapid decay for $r/\delta > 6$.
• Applying numerical quadrature: Use a standard quadrature rule on the regularized kernel.
• Performing extrapolation: Carry out the integral for three different values of $\delta$ (typically $\delta_1=3h, \delta_2=4h, \delta_3=5h$, where $h$ is the grid size), and solve a local linear system (derived from a small parameter expansion of the regularized error) to eliminate the leading-order regularization errors and recover an $O(\delta^5)$ accurate estimate.

This approach maintains $O(h^5)$ accuracy uniformly near and on the surface, and works without the need for special near-singular quadrature rules or surface parameterizations.

2. Computational Accelerations: Parallelism and Fast Summation

For practical, large-scale applications, computational cost is a major concern. The method achieves scalability through several optimization strategies:

• OpenMP Parallelization: The evaluation at each target point is independent, allowing for straightforward parallelization over multiple threads. Empirically, this approach achieves near-ideal linear speedup (e.g., a 4x speedup with 4 threads).
• Local vs. Nonlocal Interaction Splitting: Since regularization only affects the near field ($r/\delta \le 6$), extrapolated sums (for all $\delta$) are computed only for neighboring source points. Far-field (nonlocal) contributions, being essentially unaffected by regularization, are computed only once and shared across the needed $\delta$ values.
• Kernel-Independent Treecode (KITC): The far-field summation, otherwise $O(N^2)$ for $N$ sources, is accelerated by a tree-based algorithm. KITC constructs a spatial hierarchy (octree), applies barycentric Lagrange interpolation at Chebyshev points for target-cluster interactions, and uses a multipole acceptance criterion (MAC) to balance accuracy and cost. This reduces the far-field complexity to $O(M\log N)$ for $M$ targets, allowing efficient evaluation even for tens or hundreds of thousands of surface points.

3. Selection and Effects of the Smoothing Parameter

The smoothing parameter $\delta$ plays a central role in controlling regularization error:

• Regularization error before extrapolation is $O(\delta)$; after extrapolation using three $\delta$ values, the error drops to $O(\delta^5)$.
• Practical recommendations: $\delta_i = \rho_i h$ with $\rho_i=3,4,5$. This ensures both accurate regularization and negligible quadrature error.
• The method is robust to the choice of $\delta/h$ so long as $\delta/h\geq 2$.
• For on-surface evaluation, specially designed smoothing functions allow $O(\delta^5)$ accuracy without extrapolation.

4. Tuning and Experimentally Optimal Computational Parameters

Systematic experiments have determined optimal choices for the KITC and overall algorithm:

• Treecode MAC parameter ($\theta$): $\theta=0.6$ balances accuracy and speed.
• Minimal leaf size ($N_0$): Start with $N_0=2000$ for $h=1/64$, double for each $h\rightarrow h/2$ refinement.
• Interpolation degree ($p$): Use $p=6$ for $h=1/64$, increment by $2$ for each grid refinement.
• This prescription preserves the $O(h^5)$ error scaling of the overall method.

Summarized parameter progression:

$h$	$N_0$	$p$
$1/64$	2000	6
$1/128$	4000	8
$1/256$	8000	10

5. Application: Stokes Flow Around Nearly Touching Spheres

The method is validated in the challenging example of Stokes flow around two spheres separated by a small gap ($\epsilon = 1/16^3$). Integrals involving points near the adjacent sphere's surface exhibit extreme near-singularity.

• The coupling integral equations between the spheres are solved using GMRES (typically requiring about 12 iterations for $10^{-10}$ tolerance).
• Regularization is localized: self-interactions are handled with a single, small $\delta=3h$; near interactions use the extrapolated regularization with three $\delta$ values.
• The method maintains $O(h^5)$ accuracy, with maximal error in the gap region but rapid error reduction as $h\to0$.
• Efficient parallel evaluation and treecode acceleration allow tractable computation without sacrificing accuracy.

6. Summary Table: Workflow and Optimal Practice

Component	Approach	Notes
Regularization	Evaluate using 3 $\delta$, extrapolate to eliminate leading errors	$\delta=3h,4h,5h$ recommended
Quadrature	Partition-of-unity, high-order for regularized kernels	$\delta/h \geq 2$ for stability
Near-field handling	Local region (within $r/\delta \leq 6$) computed anew per $\delta$	Efficient due to small neighborhood
Far-field summation	KITC computes clusters only once per target	Single computation reused for all $\delta$
Parallelization	OpenMP over targets	Linear speedup observed
Accuracy	Uniform $O(h^5)$ near/on surface; robust on challenging geometries	Confirmed in gap test

7. Significance and Broader Impact

The extrapolated regularization method, particularly in its optimized form with parallelization and fast summation, provides a reliable, high-order, and scalable technique for evaluating nearly singular integrals in computational fluid dynamics, electrostatics, and other potential problems. The method is simple to implement, requires no special parametrization or near-singular quadratures, and addresses both accuracy and computational efficiency even in challenging geometries such as nearly touching surfaces. Practical guidance on parameter selection and computational optimization makes this method suitable for deployment in large-scale scientific and engineering applications.

---

References to Main Formulas and Algorithms: See Beale and Tlupova, Adv. Comput. Math., 2024, for the original development; for the fast summation scheme, see Wang, Krasny, and Tlupova, 2020.

Key Regularized Kernel (example):

$$S_{ij}^\delta(\mathbf{y}, \mathbf{x}) = \frac{\delta_{ij}}{r} s_1\left(\frac{r}{\delta}\right) + \frac{(y_i-x_i)(y_j-x_j)}{r^3} s_2\left(\frac{r}{\delta}\right)$$

Treecode MAC condition:

$\frac{r_c}{R} \leq \theta$

Extrapolation system for fifth-order error cancellation:

$$u_i(\mathbf{y}) + c_1 \rho I_0(\lambda) + c_2 \rho^3 I_2(\lambda) = u^\delta_i(\mathbf{y}) + O(\delta^5)$$