Software – Gonçalo Pena

Bramble.jl

Supraconvergent finite difference discretizations on nonuniform Cartesian grids.

v3.44.1MIT licenceJulia 1.13+Released 14 October 2026DOI 10.5281/zenodo.14230821

GitHub repository Documentation Zenodo DOI Related papers: Convergence of FDM/FEM

Poisson problem on a nonuniform mesh

Find $u$ such that

$$\begin{cases}-\Delta u = f & \text{in } \Omega=(0,1)^2,\\ \hphantom{-\Delta}\,u = 0 & \text{on } \partial\Omega,\end{cases}$$

with $u=\sin(\pi x)\sin(\pi y)$. On nodes $x_0<\dots<x_N$ and $y_0<\dots<y_M$, let $h_i=x_i-x_{i-1}$ and $h_{i+1/2}=(h_i+h_{i+1})/2$, with $k_j$, $k_{j+1/2}$ likewise in $y$. The inner products are

$$(u_H,v_H)_H=\sum_{i,j}h_{i+1/2}\,k_{j+1/2}\,u_H(x_i,y_j)v_H (x_i,y_j),$$ $$\begin{aligned}(p_H,q_H)_+={}&\sum_{i,j}h_i\,k_{j+1/2}\,p_H(x_i,y_j)\,q_H(x_i,y_j)\\&+\sum_{i,j}h_{i+1/2}\,k_j\,p_H(x_i,y_j)\,q_H(x_i,y_j),\end{aligned}$$

and the discrete gradient is ($D_{-y}$ likewise)

$$\nabla_H u=(D_{-x}u_H,\,D_{-y}u_H),\qquad D_{-x}u_H(x_i,y_j)=\frac{u_H(x_i,y_j)-u_H(x_{i-1},y_j)}{h_i}.$$

The problem now reads as: find $u_H$ vanishing on $\partial\Omega$ such that, for every such $v_H$,

$$(\nabla_H u_H,\nabla_H v_H)_+=(f_H,v_H)_H,$$

where $f_H$ is $f$ on the nodes. On a random mesh the truncation error is first order, yet the error decreases like $H^2_{\max}$. Run the convergence study to see it.

The same problem in Bramble.jl

Code and output from the Bramble.jl documentation.

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))
Ωₕ = mesh(Ω, (33, 33), (true, true))
Wₕ = gridspace(Ωₕ)

gₕ = element(Wₕ)
Rₕ!(gₕ, g)

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

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

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

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

0.00040178883968709

Convergence

Continuing from the code above:

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.00040178883968709, 0.00010041090485760582, 2.0005215336872038)

Halving the spacing cuts the error by four: second order.

Features

Nonuniform meshes, second order

1D, 2D and 3D Cartesian meshes with arbitrary spacing that keep second-order accuracy through supraconvergence, using staggered dual meshes.

Weak forms, sparse matrices

Write a bilinear form as an expression in u and v; assemble walks the mesh and returns a SparseMatrixCSC.

Discrete calculus

Metric-weighted D₊, D₋ and centred differences with matching discrete inner products and norms.

Allocation-free kernels

Runtime kernels allocate nothing, so assembly can be repeated inside time loops.

Coupled systems

Composite spaces Wₕ^Val(N) handle systems such as coupled reaction–diffusion.

Differentiable

Forward- and reverse-mode automatic differentiation through ForwardDiff, ReverseDiff and Enzyme.

Export

VTK for ParaView and PGFPlots/TikZ for papers.

Extensions

Optional Metal GPU, Makie/Plots, MUMPS, AlgebraicMultigrid and SciML integrations.