Getting started

A Poisson problem, end to end, in twenty lines. Six steps take a continuous problem to a solved grid function, and every tutorial afterwards is one of them in detail.

The problem

\[-\Delta u = g \text{ in } \Omega = (0,1)^2, \qquad u = u_{\text{exact}} \text{ on } \partial\Omega\]

with the manufactured solution $u_{\text{exact}}(x, y) = \sin(\pi x)\sin(\pi y)$, so $g = 2\pi^2 u_{\text{exact}}$. Knowing the answer in advance is what makes the error below checkable.

The code

using Bramble

uexact(x) = sinpi(x[1]) * sinpi(x[2])
g(x) = 2π^2 * uexact(x)

Ω = domain(interval(0.0, 1.0) × interval(0.0, 1.0))   # geometry plus boundary labels
Ωₕ = mesh(Ω, (33, 33), (true, true))                  # 33 × 33 points, uniformly spaced
Wₕ = gridspace(Ωₕ)                                    # one unknown per point

gₕ = element(Wₕ)
Rₕ!(gₕ, g)                                            # the source, sampled at the points

a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))    # the discrete Laplacian
l = form(Wₕ, v -> innerₕ(gₕ, v))                      # the load
bcs = dirichlet_constraints(Ω, :boundary => uexact)

A, F = assemble(a, l; dirichlet = bcs)

uₕ = element(Wₕ)
uₕ .= A \ F

normₕ(uₕ .- Rₕ(Wₕ, uexact))

0.0004017888396870855

Four parts in ten thousand. Halving the spacing should cut that by four, since the scheme is second order:

function poisson_error(n)
    Ωₕ = mesh(Ω, (n, n), (true, true))
    Wₕ = gridspace(Ωₕ)
    gₕ = element(Wₕ)
    Rₕ!(gₕ, g)
    A, F = assemble(
        form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v))),
        form(Wₕ, v -> innerₕ(gₕ, v));
        dirichlet = dirichlet_constraints(Ω, :boundary => uexact))
    uₕ = element(Wₕ)
    uₕ .= A \ F
    return normₕ(uₕ .- Rₕ(Wₕ, uexact))
end

e₁, e₂ = poisson_error(33), poisson_error(65)

e₁, e₂, log2(e₁ / e₂)

(0.0004017888396870855, 0.00010041090485689152, 2.0005215336974502)

The form is the step worth dwelling on. a and l are expressions in the trial and test functions, not matrices: form stores the expression, and assemble is what walks the mesh. That separation is what lets the same a be refilled every step of a time loop, handed to a nonlinear solver, or differentiated.

Where to go from here:

Discrete foundations covers the first three steps: domains, meshes and their metric, grid spaces and the operators that act on them.
Linear and bilinear forms covers the last three: writing a form, assembling it, imposing conditions, and coupled systems.
The worked examples run the whole chain on real problems, from nonlinear Poisson to 3D elasticity, a heat equation and an inverse problem.

Getting started

The problem

The code

What each step gives you