Boundary integral operators

A central piece of integral equation methods is the efficient and accurate computation of integral operators. In the first part of this tutorial we will cover how to assemble and manipulate the four integral operators of Calderón calculus, namely the single-layer, double-layer, hypersingular, and adjoint operators [1, 2], for some predefined kernels in Inti.jl. In the second part we will show how to extend the package to handle custom kernels.

Predefined kernels and integral operators

To simplify the construction of integral operators for some commonly used PDEs, Inti.jl defines a few AbstractDifferentialOperators types. For each of these PDEs, the package provides a SingleLayerKernel, DoubleLayerKernel, HyperSingularKernel, and AdjointDoubleLayerKernel that can be used to construct the corresponding kernel functions, e.g.:

using Inti, StaticArrays, LinearAlgebra
op = Inti.Helmholtz(; dim = 2, k = 2π)
G   = Inti.SingleLayerKernel(op)

Inti.SingleLayerKernel{ComplexF64, Inti.Helmholtz{2, Float64}}(Helmholtz operator in 2 dimensions: -Δu - k²u)

Typically, we are not interested in the kernels themselves, but in the integral operators they define. Two functions, single_double_layer and adj_double_layer_hypersingular, are provided as a high-level syntax to construct the four integral operators of Calderón calculus:

Γ = Inti.parametric_curve(s -> SVector(cos(s), sin(s)), 0, 2π) |> Inti.Domain
Q = Inti.Quadrature(Γ; meshsize = 0.1, qorder = 5)
S, D = Inti.single_double_layer(;
    op,
    target = Q,
    source = Q,
    compression = (method = :none,),
    correction = (method = :dim,)
)
K, N = Inti.adj_double_layer_hypersingular(;
    op,
    target = Q,
    source = Q,
    compression = (method = :none,),
    correction = (method = :dim,)
)

Much goes on under the hood in the function above, and the sections on correction and compression methods will provide more details on the options available. The important thing to keep in mind is that S, D, K, and N are discrete approximations of the following (linear) operators:

\[\begin{aligned} S[\sigma](\boldsymbol{x}) &:= \int_{\Gamma} G(\boldsymbol{x}, \boldsymbol{y}) \sigma(\boldsymbol{y}) \mathrm{d} s_{\boldsymbol{y}}, \quad &&D[\sigma](\boldsymbol{x}) := \mathrm{p.v.} \int_{\Gamma} \frac{\partial G}{\partial \nu_{\boldsymbol{y}}}(\boldsymbol{x}, \boldsymbol{y}) \sigma(\boldsymbol{y}) \mathrm{d} s_{\boldsymbol{y}} \\ K[\sigma](\boldsymbol{x}) &:= \mathrm{p.v.} \int_{\Gamma} \frac{\partial G}{\partial \nu_{\boldsymbol{x}}}(\boldsymbol{x}, \boldsymbol{y}) \sigma(\boldsymbol{y}) \mathrm{d} s_{\boldsymbol{y}}, \quad &&N[\sigma](\boldsymbol{x}) := \mathrm{f.p.} \int_{\Gamma} \frac{\partial^2 G}{\partial \nu_{\boldsymbol{x}} \partial \nu_{\boldsymbol{y}}}(\boldsymbol{x}, \boldsymbol{y}) \sigma(\boldsymbol{y}) \mathrm{d} s_{\boldsymbol{y}} \end{aligned}\]

The actual type of S, D, K, and N depends on the compression and correction methods. In the simple case above, these are simply matrices:

map(typeof, (S, D, K, N))

(Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64})

If we turn on a compression method, such as :fmm, the types may change into something different:

using FMM2D # will load the extension
Sfmm, Dfmm = Inti.single_double_layer(;
    op,
    target = Q,
    source = Q,
    compression = (method = :fmm, tol = 1e-10),
    correction = (method = :dim, )
)
Kfmm, Nfmm = Inti.adj_double_layer_hypersingular(;
    op,
    target = Q,
    source = Q,
    compression = (method = :fmm, tol = 1e-10),
    correction = (method = :dim,)
)
typeof(Sfmm)

LinearMaps.LinearCombination{ComplexF64, Tuple{LinearMaps.FunctionMap{ComplexF64, IntiFMM2DExt.var"#7#16"{Float64, Bool, Vector{Float64}, Matrix{Float64}, Matrix{Float64}}, Nothing, true}, LinearMaps.WrappedMap{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}}

This is because the FMM method is used to approximate the matrix-vector in a matrix-free way: the only thing guaranteed is that S and D can be applied to a vector:

x = map(q -> cos(q.coords[1] + q.coords[2]), Q)
norm(Sfmm*x - S*x, Inf)

3.4395341270496036e-15

The Sfmm object above in fact combines two linear maps:

Sfmm

378×378 LinearMaps.LinearCombination{ComplexF64} with 2 maps:
  378×378 LinearMaps.FunctionMap{ComplexF64,true}(#7; issymmetric=false, ishermitian=false, isposdef=false)
  378×378 LinearMaps.WrappedMap{ComplexF64} of
    378×378 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2268 stored entries

The FunctionMap computes a matrix-vector by performing a function call to the FMM2D library. The WrappedMap accounts for a sparse matrix used to correct for singular and nearly singular interactions. These two objects are added lazily using LinearMaps.

Operator composition

Effortlessly and efficiently composing operators is a powerful abstraction for integral equations, as it allows for the construction of complex systems from simple building blocks. To show this, let us show how to construct the Calderón projectors:

\[\begin{aligned} H = \begin{bmatrix} -D & S \\ -N & K \end{bmatrix} \end{aligned}\]

As is well-known [1, Theorem 3.1.3], the operators $C_\pm = I/2 \pm H$ are the projectors (i.e. $C_{\pm}^2 = C_{\pm}$):

using LinearMaps
# create the block operator
H = [-Dfmm Sfmm; -Nfmm Kfmm]
C₊ = I / 2 + H
C₋ = I / 2 - H
# define two density functions on Γ
u = map(q -> cos(q.coords[1] + q.coords[2]), Q)
v = map(q-> q.coords[1], Q)
x = [u; v]
# compute the error in the projector identity
e₊ = norm(C₊*(C₊*x) - C₊*x, Inf)
e₋ = norm(C₋*(C₋*x) - C₋*x, Inf)
println("projection error for C₊: $e₊")
println("projection error for C₋: $e₋")

projection error for C₊: 2.1314832112638394e-8
projection error for C₋: 2.131482645480564e-8

We see that the error in the projector identity is small, as expected. Note that such compositions are not limited to the Calderón projectors, and can be used e.g. to construct the combined field integral equation (CFIE), or to compose a formulation with an operator preconditioner.

Custom kernels

So far we have focused on problems for which Inti.jl provides predefined kernels, and used the high-level syntax of e.g. single_double_layer to construct the integral operators. We will now dig into the details of how to set up a custom kernel function, and how to build an integral operator from it.

For the sake of simplicity, let us consider the following kernel representing the half-space Dirichlet Green function for Helmholtz's equation in 2D:

\[ G_D(\boldsymbol{x}, \boldsymbol{y}) = \frac{i}{4} H^{(1)}_0(k |\boldsymbol{x} - \boldsymbol{y}|) - \frac{i}{4} H^{(1)}_0(k |\boldsymbol{x} - \boldsymbol{y}^*|),\]

where $\boldsymbol{y}^* = (y_1, -y_2)$. We can define this kernel as a

using SpecialFunctions # for hankelh1
function helmholtz_kernel(target, source, k)
    x, y  = Inti.coords(target), Inti.coords(source)
    yc = SVector(y[1], -y[2])
    d, dc  = norm(x-y), norm(x-yc)
    # the singularity at x = y needs to be handled separately, so just put a zero
    d == 0 ? zero(ComplexF64) : im / 4 * ( hankelh1(0, k * d) - hankelh1(0, k * dc))
end

helmholtz_kernel (generic function with 1 method)

Let us now consider the integral operator $S$ defined by:

\[ S[\sigma](\boldsymbol{x}) = \int_{\Gamma} G_D(\boldsymbol{x}, \boldsymbol{y}) \sigma(\boldsymbol{y}) \mathrm{d} s_{\boldsymbol{y}}, \quad \boldsymbol{x} \in \Gamma.\]

We can represent S by an IntegralOperator type:

k = 50π
λ = 2π/k
meshsize = λ / 10
geo = Inti.parametric_curve(s -> SVector(cos(s), 2 + sin(s)), 0, 2π)
Γ = Inti.Domain(geo)
msh = Inti.meshgen(Γ; meshsize)
Q = Inti.Quadrature(msh; qorder = 5)
# create a local scope to capture `k`
K = let k = k
    (t,q) -> helmholtz_kernel(t,q,k)
end
Sop = Inti.IntegralOperator(K, Q, Q)

9426×9426 Inti.IntegralOperator{ComplexF64, Main.var"#13#14"{Float64}, Inti.Quadrature{2, Float64}, Inti.Quadrature{2, Float64}}:
         0.0+0.0im         …  0.000148196+5.79811e-5im
  9.77388e-5+5.77601e-5im      8.87692e-5+5.77276e-5im
  5.90094e-5+5.69761e-5im      5.53068e-5+5.68505e-5im
   3.6894e-5+5.53946e-5im      3.45719e-5+5.51638e-5im
  2.37421e-5+5.3463e-5im       2.19418e-5+5.31446e-5im
  1.74977e-5+5.21219e-5im  …   1.58961e-5+5.17549e-5im
  1.59759e-5+5.17423e-5im      1.44192e-5+5.13629e-5im
  1.05873e-5+5.02037e-5im      9.18068e-6+4.9778e-5im
  2.47309e-6+4.71938e-5im      1.27443e-6+4.66931e-5im
   -5.367e-6+4.32253e-5im     -6.37251e-6+4.26473e-5im
            ⋮              ⋱             ⋮
 -7.33246e-7+4.86992e-5im      4.28174e-7+4.92372e-5im
  7.44013e-6+5.1517e-5im       8.80374e-6+5.19856e-5im
  1.28489e-5+5.29316e-5im      1.43597e-5+5.33569e-5im
  1.43724e-5+5.3277e-5im   …   1.59274e-5+5.36906e-5im
  2.06014e-5+5.44836e-5im      2.23526e-5+5.48511e-5im
  3.35692e-5+5.61654e-5im      3.58384e-5+5.6449e-5im
  5.47182e-5+5.74384e-5im      5.83651e-5+5.76196e-5im
  8.85522e-5+5.79443e-5im      9.74649e-5+5.80337e-5im
 0.000148196+5.79811e-5im  …          0.0+0.0im

The approximation of Sop now involves two steps:

• build a dense operator S₀ that efficiently computes the matrix-vector product Sop * x for any vector x
• correct for the inaccuracies of S₀ due to singular/nearly-singular interactions by adding to it a correction matrix δS

For the first step, we will use a hierarchical matrix:

using HMatrices
S₀ = Inti.assemble_hmatrix(Sop; rtol = 1e-4)

HMatrix of ComplexF64 with range 1:9426 × 1:9426
	 number of nodes in tree: 2893
	 number of leaves: 2170 (690 admissible + 1480 full)
	 min rank of sparse blocks : 4
	 max rank of sparse blocks : 30
	 min length of dense blocks : 625
	 max length of dense blocks : 3936
	 min number of elements per leaf: 625
	 max number of elements per leaf: 2468041
	 depth of tree: 0
	 compression ratio: 15.270194

The correction matrix δS will be constructed using adaptive_correction:

δS = Inti.adaptive_correction(Sop; tol = 1e-4, maxdist = 5*meshsize)

9426×9426 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 565560 stored entries:
⎡⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⎤
⎢⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥
⎣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦

How exactly one adds S₀ and δS to get the final operator depends on the intended usage. For instance, one can use the LinearMap type to simply add them lazily:

using LinearMaps
S = LinearMap(S₀) + LinearMap(δS)

9426×9426 LinearMaps.LinearCombination{ComplexF64} with 2 maps:
  9426×9426 LinearMaps.WrappedMap{ComplexF64} of
    9426×9426 HMatrix{ClusterTree{2, Float64}, ComplexF64}
  9426×9426 LinearMaps.WrappedMap{ComplexF64} of
    9426×9426 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 565560 stored entries

Or, one can add δS to S₀ to create a new object:

S = S₀ + δS

HMatrix of ComplexF64 with range 1:9426 × 1:9426
	 number of nodes in tree: 2893
	 number of leaves: 2170 (690 admissible + 1480 full)
	 min rank of sparse blocks : 4
	 max rank of sparse blocks : 30
	 min length of dense blocks : 625
	 max length of dense blocks : 3936
	 min number of elements per leaf: 625
	 max number of elements per leaf: 2468041
	 depth of tree: 0
	 compression ratio: 15.270194

or if performance/memory is a concern, one may want to directly add δS to S₀ in-place:

axpy!(1.0, δS, S₀)

HMatrix of ComplexF64 with range 1:9426 × 1:9426
	 number of nodes in tree: 2893
	 number of leaves: 2170 (690 admissible + 1480 full)
	 min rank of sparse blocks : 4
	 max rank of sparse blocks : 30
	 min length of dense blocks : 625
	 max length of dense blocks : 3936
	 min number of elements per leaf: 625
	 max number of elements per leaf: 2468041
	 depth of tree: 0
	 compression ratio: 15.270194

All of these should give an identical matrix-vector product, but the latter two allow e.g. for the use of direct solvers though an LU factorization.