Bramble.jl GitHub

Linear Poisson equation#

A worked problem, built from pieces the forms tutorial introduces one at a time. Every number and every plot below except the interactive panel comes from the code shown.

Problem#

\[-\Delta u = g \text{ in } \Omega = (0,1)^2, \qquad u = u_{\text{exact}} \text{ on } \texttt{xmin} \text{ and } \texttt{ymin}, \qquad u + \partial_n u = h \text{ on } \texttt{xmax} \text{ and } \texttt{ymax},\]

with \(\partial_n\) the outward normal derivative: a Dirichlet condition on the two sides at the origin and a Robin condition on the two others. The manufactured solution is \(u_{\text{exact}}(x, y) = e^{x+y}\), so \(g = -2 u_{\text{exact}}\), and the Robin data \(h = u_{\text{exact}} + \partial_n u_{\text{exact}} = 2 u_{\text{exact}}\) there, since the outward normal derivative of \(e^{x+y}\) on xmax and ymax is \(u_{\text{exact}}\) itself. A manufactured solution gives the exact answer the error is measured against.

Ghost points on the Robin sides#

The Robin sides go through ghost points: gridspace(Ωₕ; ghosts = (:xmax, :ymax)) adds one layer outside those two sides only. constrained takes a function of the trial function u and returns one side => relation => data entry per side, here u + ∂ₙ(u) with the data \(h\) as a function of the boundary point. Each relation fixes a ghost value from the interior ones, which leaves the physical nodes as the only unknowns.

The Dirichlet sides need no ghosts. There are two ways to impose them. The dirichlet keyword of assemble replaces the rows of those nodes, and symmetrize = true moves the known values to the right-hand side, which restores the symmetry of the matrix to rounding. The alternative is the relation u => g in constrained, which removes the nodes from the unknowns. The page solves with the keyword; the test suite checks that the relation gives the same solution.

Assembling and solving on a random grid#

using Bramble
using Random

sol(x) = exp(x[1] + x[2])
rhs(x) = -2 * sol(x)

Random.seed!(20260903)
Ω = domain(interval(0.0, 1.0) × interval(0.0, 1.0))
Ωₕ = mesh(Ω, (24, 24), (false, false))
Wₕ = gridspace(Ωₕ)
ScalarGridSpace {2D, Float64}:
  Mesh: MeshnD{2D, 576 pts}
  Dofs: 576

(false, false) builds a randomly perturbed non-uniform grid in both directions (the default is uniform). A uniform grid also converges at second order here, but its symmetry can hide an error that random points expose.

The seed fixes the draw: (false, false) takes the interior points from the global RNG, so the mesh is the same on every build and the suite can assert the printed numbers.

The bilinear form is the discrete Laplacian, written \(\langle -\Delta_h u, v \rangle_h\) because it reads the ghost values the relations fix; the forms tutorial builds the same operator in one dimension as inner₊(∇ₕ(u), ∇ₕ(v), :x). The source is the right-hand side sampled at the grid points and tested against v, innerₕ(rhs, v):

robin_data(x) = 2 * sol(x)

function poisson_system(Ωₕ)
    Cₕ = constrained(gridspace(Ωₕ; ghosts = (:xmax, :ymax)),
        u -> (:xmax => u + ∂ₙ(u) => robin_data, :ymax => u + ∂ₙ(u) => robin_data))
    a = form(Cₕ, Cₕ, (u, v) -> innerₕ(-1.0 * Δₕ(u), v))
    l = form(Cₕ, v -> innerₕ(rhs, v))
    A, F = assemble(a, l; dirichlet = (:xmin => sol, :ymin => sol), symmetrize = true)
    return Cₕ, A, F
end

Cₕ, A, F = poisson_system(Ωₕ)

uₕ = element(Wₕ)
uₕ .= pde_solve(A, F)

The solution with the grid lines drawn, then the pointwise error uₕ - sol: the colour scale of the second plot is the error's own.

The Robin relation on the lifted solution#

Bramble.lift(Cₕ, u) returns the solution on the padded grid, with the ghost values the relations fix filled in. Ghost values satisfy the relation by construction, so the check reads the physical nodes only: on xmax, the three-point backward difference of the lifted solution for the normal derivative, with the non-uniform spacings, must satisfy \(u + \partial_n u = h\) up to the truncation error, of order \(h^2\). robin_residual returns the largest residual over the nodes of the side:

function robin_residual(Ωₕ, Cₕ, u)
    U = parent(Bramble.lift(Cₕ, u))
    G = Bramble.mesh(Cₕ)
    nx, ny = Bramble.npoints(Ωₕ, Tuple)
    xs, ys = collect(points(Ωₕ(1))), collect(points(Ωₕ(2)))
    h₁, h₂ = xs[nx] - xs[nx - 1], xs[nx - 1] - xs[nx - 2]
    lin = LinearIndices(Bramble.npoints(G, Tuple))
    at(i, j) = U[lin[Bramble.padded_index(G, CartesianIndex(i, j))]]
    return maximum(1:ny) do j
        dₙ = (2h₁ + h₂) / (h₁ * (h₁ + h₂)) * at(nx, j) -
             (h₁ + h₂) / (h₁ * h₂) * at(nx - 1, j) + h₁ / (h₂ * (h₁ + h₂)) * at(nx - 2, j)
        abs(at(nx, j) + dₙ - robin_data((xs[nx], ys[j])))
    end
end

residual = robin_residual(Ωₕ, Cₕ, pde_solve(A, F))
0.00913610567344314

This is the truncation error of the one-sided difference, plus the error of the solution itself, so it is small but not round-off. It falls as \(h^2\): the same random mesh at 12 points per side, then split once, gives

residuals = let
    Random.seed!(20260903)
    Ωc = mesh(Ω, (12, 12), (false, false))
    out = Float64[]
    for level in 1:2
        Cc, Ac, Fc = poisson_system(Ωc)
        push!(out, robin_residual(Ωc, Cc, pde_solve(Ac, Fc)))
        level < 2 && iterative_refinement!(Ωc)
    end
    out
end
2-element Vector{Float64}:
 0.009464289630539469
 0.0023085554352881843
residuals[2] < residuals[1]
true

Browsing the solution#

An interactive panel, seeded from the resolution of uₕ above but solving its own problem in the browser: drag the resolution slider, switch between a uniform and a fully random mesh, or pick a different manufactured solution, and every plot and diagnostic below updates immediately.

Linear Poisson, 2D

Convergence in one, two and three dimensions#

One random coarse mesh per dimension, refined in place with iterative_refinement!, so every finer level is the same random mesh split dyadically; a fresh draw per level would add noise to the rate. Dirichlet on the sides at the origin, Robin on the others, in 1D, 2D and 3D.

function poisson_series(D::Int; n0::Int = 5, levels::Int)
    sol_d(x) = exp(sum(x))
    rhs_d(x) = -D * sol_d(x)
    robin_d(x) = 2 * sol_d(x)
    lo = (:xmin, :ymin, :zmin)[1:D]    # Dirichlet sides
    hi = (:xmax, :ymax, :zmax)[1:D]    # Robin sides, through ghost points

    Ωd = domain(reduce(×, ntuple(_ -> interval(0.0, 1.0), D)))
    Ωc = mesh(Ωd, ntuple(_ -> n0, D), ntuple(_ -> false, D))

    hs, errs = Float64[], Float64[]
    for level in 1:levels
        Cc = constrained(gridspace(Ωc; ghosts = hi),
            u -> ntuple(j -> hi[j] => u + ∂ₙ(u) => robin_d, D))
        a_c = form(Cc, Cc, (u, v) -> innerₕ(-1.0 * Δₕ(u), v))
        l_c = form(Cc, v -> innerₕ(rhs_d, v))
        A_c, F_c = assemble(a_c, l_c;
            dirichlet = ntuple(j -> lo[j] => sol_d, D), symmetrize = true)

        Wc = gridspace(Ωc)
        u_c = element(Wc)
        u_c .= pde_solve(A_c, F_c)
        uexact_c = Rₕ(Wc, sol_d)

        push!(hs, hₘₐₓ(Ωc))
        push!(errs, norm₁ₕ(u_c .- uexact_c))
        level < levels && iterative_refinement!(Ωc)
    end
    return hs, errs
end

Random.seed!(20260903)
hs1, errs1 = poisson_series(1; n0 = 6, levels = 7)  # 6 up to 321 points
Random.seed!(20260903)
hs2, errs2 = poisson_series(2; levels = 6)   # 5² up to 129² points
Random.seed!(20260903)
hs3, errs3 = poisson_series(3; levels = 4)   # 5³ up to 33³ points

order1 = log(errs1[end - 1] / errs1[end]) / log(hs1[end - 1] / hs1[end])
order2 = log(errs2[end - 1] / errs2[end]) / log(hs2[end - 1] / hs2[end])
order3 = log(errs3[end - 1] / errs3[end]) / log(hs3[end - 1] / hs3[end])
(order1, order2, order3)
(2.000007249133581, 1.9998135835124402, 1.9931291034460799)
order1 > 1.9 && order2 > 1.9 && order3 > 1.8
true
convergence_plot([(hs1, errs1, "1D", "#5B5FC7"), (hs2, errs2, "2D", "#0E7C86"), (hs3, errs3, "3D", "#B26A00")];

Second order in every dimension (fitted orders above 1.9, 1.9 and 1.8); each curve's markers follow the reference slope. The dyadic split multiplies the points by 2, 4 and 8 in 1D, 2D and 3D, so 3D runs 4 levels, 2D 6 and 1D 7.

Where to go next#