Boundary integral operators
- Define layer potentials and the four integral operators of Calderón calculus
- Construct block operators
- Set up a custom kernel
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 AbstractDifferentialOperator
s 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.4349129187055663e-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.1314846997420126e-8
projection error for C₋: 2.1314837674158142e-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.
If the integral operator of interest arises from a PDE, it is recommended to define a new AbstractDifferentialOperator
type, and implement the required methods for SingleLayerKernel
, DoubleLayerKernel
, AdjointDoubleLayerKernel
, and HyperSingularKernel
. This will enable the use of the high-level syntax for constructing boundary integral operators, as well as the use of the compression and correction methods specific to integral operators arising from PDEs.
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
Kernel functions passed to IntegralOperator
should always take two arguments, target
and source
, which are both of QuadratureNode
. This allows for extracting not only the coords
of the nodes, but also the normal
vector if needed (e.g. for double-layer or hypersingular kernels).
The approximation of Sop
now involves two steps:
- build a dense operator
S₀
that efficiently computes the matrix-vector productSop * x
for any vectorx
- 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.
Integral operators defined from custom kernel functions do not support all the features of the predefined ones. In particular, some singular integration methods (e.g. the Density Interpolation Method) and acceleration routines (e.g. Fast Multipole Method) used to correct for singular and nearly singular integral operators, and to accelerate the matrix vector products, are only available for specific kernels. Check the corrections and compression for more details concerning which methods are compatible with custom kernels.