Bramble.jl GitHub

Neumann, Robin and pure Neumann boundaries#

Poisson's equation on a non-uniform grid with a flux or a mixed condition on the whole boundary, imposed through ghost points and constrained, and checked for second order. Every number and every plot below comes from the code shown.

Problem#

\[-\Delta u = g \text{ in } \Omega = (0,1)^2, \qquad \alpha\, u + \partial_n u = h \text{ on } \partial\Omega,\]

with \(\partial_n\) the outward normal derivative and \(\alpha \ge 0\) free to change from side to side. Three choices of \(\alpha\) on xmin, xmax, ymin, ymax: Robin and Neumann, \((2, 0, 2, 0)\); Robin, \((2, 1, 2, 1)\); pure Neumann, \((0, 0, 0, 0)\). The manufactured solution is \(u(x, y) = e^{x} \cos 2y\), so \(g = 3u\) and the data \(h = \alpha u + \partial_n u\) come from its gradient.

One boundary relation per side#

gridspace(Ωₕ; ghosts = :boundary) adds one layer of ghost points outside every side. constrained takes a function of the trial function u and returns one side => relation => data entry per side: ∂ₙ(u) on a Neumann side, α * u + ∂ₙ(u) on a Robin side. Each relation fixes a ghost value from the interior ones, so the physical nodes are the only unknowns. assemble(a, l) takes both forms on the constrained space, so the relations' data reach the right-hand side.

using Bramble
using LinearAlgebra: dot
using Random

sol(p) = exp(p[1]) * cos(2p[2])
∇sol(p) = (sol(p), -2 * exp(p[1]) * sin(2p[2]))
rhs(p) = 3 * sol(p)
rhs (generic function with 1 method)

Each side with its outward normal and the coordinate the normal points along:

const SIDES = ((:xmin, -1.0, 1), (:xmax, 1.0, 1), (:ymin, -1.0, 2), (:ymax, 1.0, 2))

function side_relation(u, α, j)
    side, normal, axis = SIDES[j]
    data = p -> α[j] * sol(p) + normal * ∇sol(p)[axis]
    relation = iszero(α[j]) ? ∂ₙ(u) : α[j] * u + ∂ₙ(u)
    return side => relation => data
end

function solve_case(Ωₕ, α)
    Cₕ = constrained(gridspace(Ωₕ; ghosts = :boundary),
        u -> ntuple(j -> side_relation(u, α, j), 4))
    A, F = assemble(form(Cₕ, Cₕ, (u, v) -> innerₕ(-1.0 * Δₕ(u), v)),
        form(Cₕ, v -> innerₕ(rhs, v)))
    Wₕ = gridspace(Ωₕ)
    uₕ = element(Wₕ)
    uₑ = Rₕ(Wₕ, sol)
    if all(iszero, α)
        w = assemble(form(Wₕ, v -> innerₕ(x -> 1.0, v)))
        uₕ .= pde_solve(A, F; mean_zero = w)
        uₑ .-= dot(w, uₑ) / sum(w)
    else
        uₕ .= pde_solve(A, F)
    end
    return uₕ, uₑ
end
solve_case (generic function with 1 method)

The Robin and Neumann case on one random grid, then the pointwise error:

Random.seed!(20260903)
Ω = domain(interval(0.0, 1.0) × interval(0.0, 1.0))
uₕ, uₑ = solve_case(mesh(Ω, (17, 17), (false, false)), (2.0, 0.0, 2.0, 0.0))

A pure Neumann problem fixes \(u\) only up to a constant, and the data must satisfy a compatibility condition, which the manufactured data satisfy in the continuum. The matrix is then singular, its null space the constants. solve_case passes pde_solve the keyword mean_zero, the vector of cell measures that is the discrete integral, which solves the system bordered by that constraint. The error is measured against the exact solution shifted to the same mean. The pure Neumann solution on the grid above:

The error is the discrete \(H^1\) norm norm₁ₕ, and the grid is a random non-uniform one refined in place by iterative_refinement!, so each level is the same random mesh split dyadically:

function convergence(α; levels = 5)
    Random.seed!(20260903)
    Ω = domain(interval(0.0, 1.0) × interval(0.0, 1.0))
    Ωₕ = mesh(Ω, (9, 9), (false, false))
    hs, errs = Float64[], Float64[]
    for level in 1:levels
        uₕ, uₑ = solve_case(Ωₕ, α)
        push!(hs, hₘₐₓ(Ωₕ))
        push!(errs, norm₁ₕ(uₕ .- uₑ))
        level < levels && iterative_refinement!(Ωₕ)
    end
    return hs, errs
end

order(hs, errs) = log(errs[end - 1] / errs[end]) / log(hs[end - 1] / hs[end])

hs_mixed, errs_mixed = convergence((2.0, 0.0, 2.0, 0.0))
hs_robin, errs_robin = convergence((2.0, 1.0, 2.0, 1.0))
hs_neumann, errs_neumann = convergence((0.0, 0.0, 0.0, 0.0))

order_mixed = order(hs_mixed, errs_mixed)
order_robin = order(hs_robin, errs_robin)
order_neumann = order(hs_neumann, errs_neumann)
(order_mixed, order_robin, order_neumann)
(2.000092631987504, 1.9999237630517444, 2.000173713957636)
order_mixed > 1.95 && order_robin > 1.95 && order_neumann > 1.95
true

Second order#

The error against \(h_{\max}\) on log-log axes, one curve per choice of \(\alpha\). The dashed line has slope 2 and is anchored at the finest point of the first curve.

convergence_plot(
    [(hs_mixed, errs_mixed, "Robin and Neumann", "#5B5FC7"),
        (hs_robin, errs_robin, "Robin", "#0E7C86"),
        (hs_neumann, errs_neumann, "pure Neumann", "#B26A00")];

All three curves run parallel to the slope-2 line. The points are random, so no symmetry cancels the one-sided error of the ghost layer. The Neumann and Robin relations keep second order.

Where to go next#