Correction methods

Important Points Covered in This Tutorial

Overview of the correction methods available in Inti.jl.
Details and limitations of the various correction methods.
Guidelines on how to choose a correction method.

When the underlying kernel is singular, a correction is usually necessary to obtain accurate results in the approximation of the underlying integral operator by a quadrature. Currently, Inti.jl provides the following functions to correct for singularities:

Each method has its own strengths and weaknesses, which will be discussed in the following sections.

High-Level API

The single_double_layer, adj_double_layer_hypersingular, and volume_potential functions provide a high-level API with a correction keyword argument. This allows users to specify the correction method to use when constructing the integral operators. See the documentation of these functions for more details.

Adaptive Correction

The adaptive_correction method combines adaptive quadrature for nearly singular integrals with a direct evaluation method for singular integrals, based on [7]. It is a robust method suitable for a wide range of kernels, as long as the singularities are no worse than a Hadamard finite-part (e.g., $1/r^3$ in 3D and $1/r^2$ in 2D). This makes it a good default choice for most problems.

Strengths

Robust method that works for a wide range of kernels.
Conceptually straightforward and easy to use.
Handles open surfaces.

Weaknesses

Can be slow for large problems and high accuracy requirements.
Sometimes difficult to tune parameters for optimal performance.
Round-off errors in certain cases can make achieving high accuracy challenging.

Docstrings

Inti.adaptive_correction — Function

adaptive_correction(iop::IntegralOperator; [maxdist, rtol, threads = true, kwargs...])
adaptive_correction(iop::IntegralOperator, maxdist, quads_dict::Dict, threads = true)

This function computes a sparse correction for the integral operator iop, addressing its singular or nearly singular entries.

The parameter maxdist specifies the maximum distance between target points and source elements to be considered for correction (only interactions within this distance are corrected).

The parameter rtol defines the tolerance for the adaptive quadrature used to compute the corrections for singular or nearly singular entries.

Additional kwargs arguments are passed to adaptive_quadrature; see its documentation for more information.

Selecting maxdist and rtol involves balancing accuracy and computational cost. A smaller maxdist reduces the number of corrections but may impact accuracy. Conversely, a smaller rtol improves correction accuracy but increases computational expense. The ideal values for maxdist and rtol depend on the kernel and the mesh/quadrature rule applied.

By default, maxdist and rtol are estimated using the local_correction_dist_and_tol, but it is often possible to improve performance by manually tunning these parameters.

Advanced usage

For finer control, you can provide a dictionary quads_dict that contains quadrature rules for each reference element type present in the mesh of source(iop). This allows you to fine-tune the quadrature rules for specific element types (e.g. use a fixed quadrature rule instead of an adaptive one).

The dictionary quads_dict must adhere to the following structure:

quads_dict[E].nearfield_quad: A function that integrates over the nearfield of the reference element type E. Used in the nearly-singular correction.
quads_dict[E].radial_quad: A function that integrates over the radial direction of the reference element type E. Used in the singular correction.
quads_dict[E].angular_quad: A function that integrates over the angular direction of the reference element type E. Used in the singular correction.

Here is an example of how to implement a custom quads_dict given an iop:

quads_dict = Dict()
msh = Inti.mesh(source(iop))
for E in Inti.element_types(msh)
    ref_domain = Inti.reference_domain(E)
    quads = (
        nearfield_quad = Inti.adaptive_quadrature(ref_domain; atol),
        radial_quad    = Inti.GaussLegendre(;order=5),
        angular_quad   = Inti.GuassLegendre(;order=20),
    )
    quads_dict[E] = quads
end

This will use an adaptive quadrature rule for the nearfield and fixed Gauss-Legendre quadrature rules for the radial and angular directions when computing the singular correction in polar coordinates on the reference domain. You can then call adaptive_correction(iop, maxdist, quads_dict) to use the custom quadrature.

source

Boundary Density Interpolation Method

The bdim_correction method implements the general-purpose version of the density interpolation method proposed in [8]. Is a global correction method that uses solutions of the underlying PDE, together with Green's identities, to interpolate the density on the boundary. It works best for low to moderate-order quadratures and is particularly useful for smooth boundaries when the PDE.

Strengths

Can be faster and more accurate for standard problems, such as scattering by closed surfaces.
Easier parameter tuning, as it only requires knowing whether the target surface is inside, outside, or on the boundary of the source.

Weaknesses

Only suitable for closed surfaces.
The underlying kernel must be related to the fundamental solution of a PDE.

Docstrings

Inti.bdim_correction — Function

bdim_correction(op,X,Y,S,D; green_multiplier, kwargs...)

Given a op and a (possibly inaccurate) discretizations of its single and double-layer operators S and D (taking a vector of values on Y and returning a vector on of values on X), compute corrections δS and δD such that S + δS and D + δD are more accurate approximations of the underlying single- and double-layer integral operators.

See [8] for more details on the method.

Arguments

Required:

op must be an AbstractDifferentialOperator
Y must be a Quadrature object of a closed surface
X is either inside, outside, or on Y
S and D are approximations to the single- and double-layer operators for op taking densities in Y and returning densities in X.
green_multiplier (keyword argument) is a vector with the same length as X storing the value of μ(x) for x ∈ X in the Green identity S\[γ₁u\](x) - D\[γ₀u\](x) + μ*u(x) = 0. See _green_multiplier.

Optional kwargs:

parameters::DimParameters: parameters associated with the density interpolation method
derivative: if true, compute the correction to the adjoint double-layer and hypersingular operators instead. In this case, S and D should be replaced by a (possibly innacurate) discretization of adjoint double-layer and hypersingular operators, respectively.
maxdist: distance beyond which interactions are considered sufficiently far so that no correction is needed. This is used to determine a threshold for nearly-singular corrections when X and Y are different surfaces. When X === Y, this is not needed.

source

Volume density interpolation method

TODO