Poisson Problem

Problem definition

In this example we will solve the Poisson equation in a domain $\Omega$ with Dirichlet boundary conditions on $\Gamma := \partial \Omega$:

\[ \begin{align*} -\Delta u &= f \quad \text{in } \quad \Omega\\ u &= g \quad \text{on } \quad \Gamma \end{align*}\]

where $f : \Omega \to \mathbb{R}$ and $g : \Gamma \to \mathbb{R}$ are given functions. To solve this problem using integral equations, we split the solution $u$ into a particular solution $u_p$ and a homogeneous solution $u_h$:

\[ u = u_p + u_h.\]

The function $u_p$ is given by

\[u_p(\boldsymbol{r}) = \int_{\Omega} G(\boldsymbol{r}, \boldsymbol{r'}) f(\boldsymbol{r'}) d\boldsymbol{r'}.\]

with $G$ the fundamental solution of $-\Delta$.

The function $u_h$ satisfies the homogeneous problem

\[ \begin{align*} \Delta u_h &= 0, \quad &&\text{in } \quad \Omega, \\ u_h &= g - u_p, \quad &&\text{on } \quad \Gamma, \end{align*}\]

which can be solved using the integral equation method. In particular, for this example, we employ a double-layer formulation:

\[u_h(\boldsymbol{r}) = \int_{\Gamma} G(\boldsymbol{r}, \boldsymbol{r'}) \sigma(\boldsymbol{r}') d\boldsymbol{r'},\]

where the density function $\sigma$ solves the integral equation

\[ -\frac{\sigma(\boldsymbol{x})}{2} + \int_{\Gamma} \partial_{\nu_{\boldsymbol{y}}}G(\boldsymbol{x}, \boldsymbol{y}) \sigma(\boldsymbol{y}) \ \mathrm{d} s_{\boldsymbol{y}} = g(\boldsymbol{x}) - u_p(\boldsymbol{x}).\]

In what follows we illustrate how to solve the problem in this manner.

Geometry and mesh

We use the Gmsh API to create a jellyfish-shaped domain and to generate a second order mesh of its interior and boundary:

using Inti, Gmsh
meshsize = 0.1
gmsh.initialize()
jellyfish = Inti.gmsh_curve(0, 2π; meshsize) do s
    r = 1 + 0.3 * cos(4 * s + 2 * sin(s))
    return r * Inti.Point2D(cos(s), sin(s))
end
cl = gmsh.model.occ.addCurveLoop([jellyfish])
surf = gmsh.model.occ.addPlaneSurface([cl])
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(2)
msh = Inti.import_mesh(; dim = 2)
gmsh.finalize()

Info    : Meshing 1D...
Info    : Meshing curve 1 (BSpline)
Info    : Done meshing 1D (Wall 0.0265784s, CPU 0.026579s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0752421s, CPU 0.07524s)
Info    : 566 nodes 1033 elements
Info    : Meshing order 2 (curvilinear on)...
Info    : [  0%] Meshing curve 1 order 2
Info    : [ 60%] Meshing surface 1 order 2
Info    : Surface mesh: worst distortion = 0.451559 (0 elements in ]0, 0.2]); worst gamma = 0.758919
Info    : Done meshing order 2 (Wall 0.0080392s, CPU 0.008039s)

We can now extract components of the mesh corresponding to the $\Omega$ and $\Gamma$ domains:

Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, msh)
Γ = Inti.boundary(Ω)
Ω_msh = view(msh, Ω)
Γ_msh = view(msh, Γ)

and visualize them:

using Meshes, GLMakie
viz(Ω_msh; showsegments = true)
viz!(Γ_msh; color = :red)

To conclude the geometric setup, we need a quadrature for the volume and boundary:

Ω_quad = Inti.Quadrature(Ω_msh; qorder = 4)
Γ_quad = Inti.Quadrature(Γ_msh; qorder = 6)

Integral operators

We can now assemble the required volume potential. To obtain the value of the particular solution $u_p$ on the boundary for the modified integral equation above we will need the volume integral operator mapping to points on the boundary, i.e. operator:

using FMM2D #to accelerate the maps
op = Inti.Laplace(; dim = 2)
# Newtonian potential mapping domain to boundary
V_d2b = Inti.volume_potential(;
    op,
    target = Γ_quad,
    source = Ω_quad,
    compression = (method = :fmm, tol = 1e-12),
    correction = (method = :dim, maxdist = 5 * meshsize, target_location = :on),
)

588×5100 LinearMaps.LinearCombination{Float64} with 2 maps:
  588×5100 LinearMaps.FunctionMap{Float64,true}(#3; issymmetric=false, ishermitian=false, isposdef=false)
  588×5100 LinearMaps.WrappedMap{Float64} of
    588×5100 SparseArrays.SparseMatrixCSC{Float64, Int64} with 3528 stored entries

We require also the boundary integral operators for the ensuing integral equation:

# Single and double layer operators on Γ
S_b2b, D_b2b = Inti.single_double_layer(;
    op,
    target = Γ_quad,
    source = Γ_quad,
    compression = (method = :fmm, tol = 1e-12),
    correction = (method = :dim,),
)

(588×588 LinearMaps.LinearCombination{Float64} with 2 maps:
  588×588 LinearMaps.FunctionMap{Float64,true}(#3; issymmetric=false, ishermitian=false, isposdef=false)
  588×588 LinearMaps.WrappedMap{Float64} of
    588×588 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4116 stored entries, 588×588 LinearMaps.LinearCombination{Float64} with 2 maps:
  588×588 LinearMaps.FunctionMap{Float64,true}(#4; issymmetric=false, ishermitian=false, isposdef=false)
  588×588 LinearMaps.WrappedMap{Float64} of
    588×588 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4116 stored entries)

Solving the linear system

We are now in a position to solve the original Poisson problem, but for that we need to specify the functions $f$ and $g$. In order to verify that our numerical approximation is correct, however, we will play a different game and specify instead a manufactured solution $u_e$ from which we will derive the functions $f$ and $g$:

# Create a manufactured solution
uₑ = (x) -> cos(2 * x[1]) * sin(2 * x[2])
fₑ = (x) -> 8 * cos(2 * x[1]) * sin(2 * x[2]) # -Δuₑ
g = map(q -> uₑ(q.coords), Γ_quad)
f = map(q -> fₑ(q.coords), Ω_quad)

With these, we can compute the right-hand-side of the integral equation for the homogeneous part of the solution:

rhs = g - V_d2b * f

and solve the integral equation for the integral density function $σ$:

using IterativeSolvers, LinearAlgebra
σ = gmres(-I / 2 + D_b2b, rhs; abstol = 1e-8, verbose = true, restart = 1000)

=== gmres ===
rest	iter	resnorm
  1	  1	2.08e+00
  1	  2	1.26e+00
  1	  3	1.76e-01
  1	  4	2.10e-02
  1	  5	1.34e-03
  1	  6	5.55e-04
  1	  7	2.34e-04
  1	  8	8.19e-05
  1	  9	1.35e-05
  1	 10	2.94e-06
  1	 11	9.76e-07
  1	 12	1.31e-07
  1	 13	2.50e-08

With the density function at hand, we can now reconstruct our approximate solution:

G = Inti.SingleLayerKernel(op)
dG = Inti.DoubleLayerKernel(op)
𝒱 = Inti.IntegralPotential(G, Ω_quad)
𝒟 = Inti.IntegralPotential(dG, Γ_quad)
u = (x) -> 𝒱[f](x) + 𝒟[σ](x)

#13 (generic function with 1 method)

and evaluate it at any point in the domain:

x = Inti.Point2D(0.1, 0.4)
println("error at $x: ", u(x) - uₑ(x))

error at [0.1, 0.4]: 0.00010371018731314408

Solution evaluation and visualization

Although we have "solved" the problem in the previous section, using the anonymous function u to evaluate the field is neither efficient nor accurate when there are either many points to evaluate, or when they lie close to the domain $\Omega$. The fundamental reason for this is the usual: the integral operators in the function u are dense matrices, and their evaluation inside or near to $\Omega$ suffers from inaccurate singular and near-singular quadrature.

To address this issue, we need to assemble accelerated and corrected versions of the integral operators. Let us suppose we wish to evaluate the solution $u$ at all the quadrature nodes of $\Omega$:

V_d2d = Inti.volume_potential(;
    op,
    target = Ω_quad,
    source = Ω_quad,
    compression = (method = :fmm, tol = 1e-8),
    correction = (method = :dim,),
)

5100×5100 LinearMaps.LinearCombination{Float64} with 2 maps:
  5100×5100 LinearMaps.FunctionMap{Float64,true}(#3; issymmetric=false, ishermitian=false, isposdef=false)
  5100×5100 LinearMaps.WrappedMap{Float64} of
    5100×5100 SparseArrays.SparseMatrixCSC{Float64, Int64} with 30600 stored entries

Likewise, we need operators mapping densities from our boundary quadrature to our mesh nodes:

We now evaluate the solution at all quadrature nodes and compare it to the manufactured:

u_quad = V_d2d * f + D_b2d * σ
er_quad = u_quad - map(q -> uₑ(q.coords), Ω_quad)
println("maximum error at all quadrature nodes: ", norm(er_quad, Inf))

maximum error at all quadrature nodes: 0.0012392645770104675

Lastly, let us visualize the solution and the error on the mesh nodes using quadrature_to_node_vals:

nodes = Inti.nodes(Ω_msh)
u_nodes = Inti.quadrature_to_node_vals(Ω_quad, u_quad)
er = u_nodes - map(uₑ, nodes)
colorrange = extrema(u_nodes)
fig = Figure(; size = (800, 300))
ax = Axis(fig[1, 1]; aspect = DataAspect())
viz!(Ω_msh; colorrange, color = u_nodes, interpolate = true)
cb = Colorbar(fig[1, 2]; label = "u", colorrange)
# plot error
log_er = log10.(abs.(er))
colorrange = extrema(log_er)
colormap = :inferno
ax = Axis(fig[1, 3]; aspect = DataAspect())
viz!(Ω_msh; colorrange, colormap, color = log_er, interpolate = true)
cb = Colorbar(fig[1, 4]; label = "log₁₀|u - uₑ|", colormap, colorrange)