Document the layer potential usage

docs/src/index.md

@@ -57,7 +57,7 @@ dependencies at once:
 
 ```@example weakdeps
 using Inti
-Inti.stack_weakdeps_env!(; verbose = false, update = false)
+Inti.stack_weakdeps_env!(; verbose = false, update = true)
 ```
 
 Note that the first time you run this command, it may take a while to download

docs/src/tutorials/compression_methods.md

@@ -4,6 +4,11 @@
 CurrentModule = Inti
 ```
 
+!!! note "Important points covered in this tutorial"
+    - Overview of the compression methods available in Inti.jl
+    - Details and limitations of the various compression methods
+    - Guideline on how to choose a compression method
+
 Inti.jl wraps several external libraries providing acceleration routines for
 integral operators. In general, acceleration routines have the signature
 `assemble_*(iop, args...; kwargs...)`, and take an [`IntegralOperator`](@ref) as

docs/src/tutorials/correction_methods.md

@@ -3,3 +3,40 @@
 ```@meta
 CurrentModule = Inti
 ```
+
+!!! warning "Work in progress"
+    This tutorial is still a work in progress. We will update it with more
+    details and examples in the future.
+
+!!! note "Important points covered in this tutorial"
+    - Overview of the correction methods available in Inti.jl
+    - Details and limitations of the various correction methods
+    - Guideline on how to choose a correction method
+
+When the underlying kernel is singular, a *correction* is usually necessary in
+order to obtain accurate results in the approximation of the underlying integral
+operator by a quadrature. At present, Inti.jl provides the following functions
+to correct for singularities:
+
+- [`adaptive_correction`](@ref)
+- [`bdim_correction`](@ref)
+- [`vdim_correction`](@ref)
+
+They have different strengths and weaknesses, and we will discuss them in the
+following sections.
+
+!!! note "High-level API"
+    Note that the [`single_double_layer`](@ref),
+    [`adj_double_layer_hypersingular`](@ref), and [`volume_potential`](@ref)
+    functions have high-level API with a `correction` keyword argument that
+    allows you to specify the correction method to use when constructing the
+    integral operators; see the documentation of these functions for more
+    details.
+
+## Adaptive correction
+
+## Boundary density interpolation method
+
+## Volume density interpolation method
+
+## Martensen-Kussmaul method

docs/src/tutorials/geo_and_meshes.md

@@ -173,6 +173,7 @@ for (n,shape) in enumerate(Inti.PREDEFINED_SHAPES)
       msh = Inti.meshgen(Γ; meshsize = 0.1)
       i,j = (n-1) ÷ ncols + 1, (n-1) % ncols + 1
       ax = Axis3(fig[i,j]; aspect = :data, title = shape)
+      hidedecorations!(ax)
       viz!(msh; showsegments = true)
 end
 fig # hide

docs/src/tutorials/integral_operators.md

@@ -154,6 +154,7 @@ 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)
+@assert e₊ < 1e-5 && e₋ < 1e-5 # hide
 println("projection error for C₊: $e₊")
 println("projection error for C₋: $e₋")
 ```

docs/src/tutorials/layer_potentials.md

@@ -6,18 +6,90 @@ CurrentModule = Inti
 
 !!! note "Important points covered in this tutorial"
     - Nearly singular evaluation of layer potentials
-    - Creating a smooth domain with splines using Gmsh.jl
+    - Creating a smooth domain with splines using [Gmsh](https://gmsh.info/)'s API
+    - Plotting values on a mesh
+
+In this tutorial we focus on **evaluating** the layer potentials given a source
+density. This is a common post-processing task in boundary integral equation
+methods, and while most of it is straightforward, some subtleties arise when the
+target points are close to the boundary (nearly-singular integrals). 
+
+## Integral potentials
+
+[`IntegralPotential`](@ref) represent the following mathematical objects:
+
+```math
+\mathcal{P}[\sigma](\boldsymbol{r}) = \int_{\Gamma} K(\boldsymbol{r}, \boldsymbol{r'}) \sigma(\boldsymbol{r'}) \, d\boldsymbol{r'}
+```
+
+where ``K`` is the kernel of the operator, ``\Gamma`` is the source's boundary,
+``\boldsymbol{r} \not \in \Gamma`` is a target point, and ``\sigma`` is the
+source density.
+
+Here is a simple example of how to create a kernel representing Laplace's
+double-layer potential:
+
+```@example layer_potentials
+using Inti, StaticArrays, LinearAlgebra
+# define a kernel function
+function K(target,source)
+    r = Inti.coords(target) - Inti.coords(source)
+    ny = Inti.normal(source)
+    return 1 / (2π * norm(r)^2) * dot(r, ny)
+end
+# define a domain
+Γ = Inti.parametric_curve(s -> SVector(cos(2π * s), sin(2π * s)), 0, 1) |> Inti.Domain
+# and a quadrature of Γ
+Q = Inti.Quadrature(Γ; meshsize = 0.1, qorder = 5)
+𝒮 = Inti.IntegralPotential(K, Q)
+```
+
+If you have a source density ``\sigma``, defined on the quadrature nodes of
+``\Gamma``, you can create a function that evaluates the layer potential at an
+arbitrary point:
+
+```@example layer_potentials
+σ = map(q -> 1.0, Q)
+u = 𝒮[σ]
+```
+
+`u` is now an anonymous function that evaluates the layer potential at any point:
+
+```@example layer_potentials
+r = SVector(0.1, 0.2)
+@assert u(r) ≈ -1 # hide
+u(r)
+```
+
+Although we created the single-layer potential for the Laplace kernel manually,
+it is often more convenient to use the `single_layer_potential` when working
+with a supported PDE, e.g.:
+
+```@example layer_potentials
+pde = Inti.Laplace(; dim = 2)
+𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q)
+```
+
+creates the single and double layer potentials for the Laplace equation in 2D.
 
 ## Direct evaluation of layer potentials
 
+We now show how to evaluate the layer potentials of an exact solution on a mesh
+created through the Gmsh API. Do to so, let us first define the PDE:g
+
 ```@example layer_potentials
 using Inti, StaticArrays, LinearAlgebra, Meshes, GLMakie, Gmsh
 # define the PDE
 k = 4π
 pde = Inti.Helmholtz(; dim = 2, k)
-meshsize = 2π / k / 10
-# create the domain and mesh using the Gmsh API
+```
+
+We will now use the [`gmsh_curve`](@ref) function to create a smooth domain of a
+kite using splines:
+
+```@example layer_potentials
 gmsh.initialize()
+meshsize = 2π / k / 4 
 kite = Inti.gmsh_curve(0, 1; meshsize) do s
     SVector(0.25, 0.0) + SVector(cos(2π * s) + 0.65 * cos(4π * s[1]) - 0.65, 1.5 * sin(2π * s))
 end
@@ -27,27 +99,58 @@ gmsh.model.occ.synchronize()
 gmsh.model.mesh.generate(2)
 msh = Inti.import_mesh(; dim = 2)
 gmsh.finalize()
+```
+
+!!! tip
+    The GMSH API is a powerful tool to create complex geometries and meshes
+    directly from Julia (the `gmsh_curve` function above is just a simple
+    wrapper around some spline functionality). For more information, see the
+    [official
+    documentation](https://gmsh.info/doc/texinfo/gmsh.html#Gmsh-application-programming-interface).
+
+We can visualize the triangular mesh using:
+
+```@example layer_potentials
+using Meshes, GLMakie
 # extract the domain Ω from the mesh entities
 ents = Inti.entities(msh)
 Ω = Inti.Domain(e->Inti.geometric_dimension(e) == 2, ents)
-# create a quadrature on the boundary
-Γ = Inti.boundary(Ω)
-Q = Inti.Quadrature(view(msh,Γ); qorder = 5)
+viz(msh[Ω]; showsegments = true, axis = (aspect = DataAspect(), ))
+```
+
+For the purpose of testing the accuracy of the layer potential evaluation, we
+will construct an exact solution of the Helmholtz equation on the interior
+domain and plot it:
+
+```@example layer_potentials
 # construct an exact interior solution as a sum of random plane waves
 dirs  = [SVector(cos(θ), sin(θ)) for θ in 2π*rand(10)]
 coefs = rand(ComplexF64, 10)
 u  =  (x)   -> sum(c*exp(im*k*dot(x, d)) for (c,d) in zip(coefs, dirs))
 du =  (x,ν) -> sum(c*im*k*dot(d, ν)*exp(im*k*dot(x, d)) for (c,d) in zip(coefs, dirs))
-# plot it 
+# plot the exact solution
 Ω_msh = view(msh, Ω)
 target = Inti.nodes(Ω_msh)
 viz(Ω_msh; showsegments = false, axis = (aspect = DataAspect(), ), color = real(u.(target)))
 ```
 
-Let us now compute the layer potentials of the exact solution on the boundary,
-and evaluate the error on the target nodes:
+Since `u` satisfies the Helmholtz equation, we know that the following
+representation holds:
+
+```math
+u(\boldsymbol{r}) = \mathcal{S}[\gamma_1 u](\boldsymbol{r}) - \mathcal{D}[\gamma_0 u](\boldsymbol{r}), \quad \boldsymbol{r} \in \Omega
+```
+
+where ``\gamma_0 u`` and ``\gamma_1 u`` are the Dirichlet and Neumann traces of
+``u``, and ``\mathcal{S}`` and ``\mathcal{D}`` are the single and double layer
+potentials over ``\Gamma := \partial \Omega``.
+
+Let's compare next the exact solution with the layer potential evaluation, based
+on a quadrature of ``\Gamma``:
 
 ```@example layer_potentials
+Γ = Inti.boundary(Ω)
+Q = Inti.Quadrature(view(msh,Γ); qorder = 5)
 # evaluate the layer potentials
 𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q)
 γ₀u = map(q -> u(q.coords), Q)
@@ -67,22 +170,43 @@ Colorbar(fig[1, 2]; label = "log₁₀(error)", colorrange)
 fig
 ```
 
+We see a common pattern of potential evaluation: the error is small away from
+the boundary, but grows near it. This is due to the nearly-singular nature of
+the layer potential integrals, which can be mitigated by using a correction
+method that accounts for the singularity of the kernel as ``\boldsymbol{r} \to
+\Gamma``.
+
 ## Near-field correction of layer potentials
 
+There are two cases where the direct evaluation of layer potentials is not
+recommended:
+
+1. When the target point is close to the boundary (nearly-singular integrals).
+2. When you wish to evaluate the layer potential at many target points and take
+   advantage of an acceleration routine.
+
+In such cases, it is recommended to use the `single_double_layer` function, (or
+to directly assemble an `IntegralOperator`) with a correction method. Here is an
+example of how to use the FMM acceleration with a near-field correction to
+evaluate the layer potentials::
+
 ```@example layer_potentials
+using FMM2D
 S, D = Inti.single_double_layer(; pde, target, source = Q,
-    compression = (method = :none, ),
+    compression = (method = :fmm, tol = 1e-12),
     correction = (method = :dim, target_location = :inside, maxdist = 0.2)
 )
 er_log10_cor = log10.(abs.(S*γ₁u - D*γ₀u - u.(target)))
-colorrange = extrema(er_log10_cor)
-fig, ax, pl = viz(Ω_msh;
-    color = er_log10_cor,
-    colormap = :viridis,
-    colorrange,
-    axis = (aspect = DataAspect(),), 
-    interpolate=true
-)
-Colorbar(fig[1, 2]; label = "log₁₀(error)", colorrange)
+colorrange = extrema(er_log10) # use scale without correction
+fig = Figure(resolution = (800, 400))
+ax1 = Axis(fig[1, 1], aspect = DataAspect(), title = "Naive evaluation")
+viz!(Ω_msh; color = er_log10, colormap = :viridis, colorrange,interpolate=true)
+ax2 = Axis(fig[1, 2], aspect = DataAspect(), title = "Nearfield correction")
+viz!(Ω_msh; color = er_log10_cor, colormap = :viridis, colorrange, interpolate=true)
+Colorbar(fig[1, 3]; label = "log₁₀(error)", colorrange)
 fig
 ```
+
+As can be seen, the near-field correction significantly reduces the error near
+the boundary, making if feasible to evaluate the layer potential near ``\Gamma``
+if necessary.

docs/src/tutorials/solvers.md

@@ -1,5 +1,16 @@
-# Solvers
+# Linear solvers
 
 ```@meta
 CurrentModule = Inti
 ```
+
+!!! warning "Work in progress"
+    This tutorial is still a work in progress. We will update it with more
+    details and examples in the future.
+
+Inti.jl does not provide its own linear solvers, but relies on external
+libraries such as
+[IterativeSolvers.jl](https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl)
+or the [LinearAlgebra](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/)
+standard library for the solving the linear systems that arise in the
+discretization of integral equations.

src/api.jl

@@ -225,6 +225,60 @@ function single_double_layer_potential(; pde, source)
     return 𝒮, 𝒟
 end
 
+"""
+    volume_potential(; pde, target, source::Quadrature, compression, correction)
+
+Compute the volume potential operator for a given PDE.
+
+## Arguments
+- `pde`: The PDE (Partial Differential Equation) to solve.
+- `target`: The target domain where the potential is computed.
+- `source`: The source domain where the potential is generated.
+- `compression`: The compression method to use for the potential operator.
+- `correction`: The correction method to use for the potential operator.
+
+## Returns
+
+The volume potential operator `V` that represents the interaction between the
+target and source domains.
+
+## Compression
+
+The `compression` argument is a named tuple with a `method` field followed by
+method-specific fields. It specifies how the dense linear operators should be
+compressed. The available options are:
+
+  - `(method = :none, )`: no compression is performed, the resulting matrices
+    are dense.
+  - `(method =:hmatrix, tol)`: the resulting operators are compressed using
+    hierarchical matrices with an absolute tolerance `tol` (defaults to `1e-8`).
+  - `(method = :fmm, tol)`: the resulting operators are compressed using the
+    fast multipole method with an absolute tolerance `tol` (defaults to `1e-8`).
+
+## Correction
+
+The `correction` argument is a named tuple with a `method` field followed by
+method-specific fields. It specifies how the singular and nearly-singular
+integrals should be computed. The available options are:
+
+  - `(method = :none, )`: no correction is performed. This is not recommented,
+    as the resulting approximation will be inaccurate if the source and target
+    are not sufficiently far apart.
+  - `(method = :dim, maxdist, target_location)`: use the density interpolation
+    method to compute the correction. `maxdist` specifies the distance between
+    source and target points above which no correction is performed (defaults to
+    `Inf`). `target_location` should be either `:inside`, `:outside`, or `:on`,
+    and specifies where the `target`` points lie relative to the to the
+    `source`'s boundary. When `target === source`, `target_location` is not
+    needed.
+
+## Details
+The volume potential operator is computed by assembling the integral operator
+`V` using the single-layer kernel `G`. The operator `V` is then compressed using
+the specified compression method. If no compression is specified, the operator
+is returned as is. If a correction method is specified, the correction is
+computed and added to the compressed operator.
+"""
 function volume_potential(; pde, target, source::Quadrature, compression, correction)
     correction = _normalize_correction(correction, target, source)
     compression = _normalize_compression(compression, target, source)