Convection-diffusion equation#

One more linear operator, built the same way as the Poisson example: diffusion plus a constant convective term, in 1D, 2D and 3D. Every number and every plot below except the interactive panel comes from the code shown.

Problem#

-\nabla \cdot (\epsilon \nabla u) - b \, \mathbf{1} \cdot \nabla u = g \text{ in } \Omega, \qquad u = u_{\text{exact}} \text{ on the far sides}, \qquad \partial_n u = h \text{ on the outflow sides}, \qquad \Omega = (0,1)^D

with a constant diffusion coefficient $\epsilon$ and a constant convective speed $$b$$ along every coordinate direction at once: the vector field $-b(1,\dots,1)$ . That field points towards the sides $$x_i = 0$$ , so the fluid leaves through xmin, ymin and zmin, where the natural condition is a given flux $\partial_n u = h$ . The other sides, where the fluid enters, keep the Dirichlet condition $u = u_{\text{exact}}$ . The manufactured solution is the same as the other two examples, $u_{\text{exact}}(x) = e^{\sum_i x_i}$ , so $g = -D\,(b+\epsilon)\,u_{\text{exact}}$ , and on an outflow side the outward normal is $$-e_i$$ , so $h = -u_{\text{exact}}$ is exact.

Adding the convective term#

using Bramble
using Random

ϵ, b = 1.0, 0.1
sol(x) = exp(sum(x))
rhs(x) = -length(x) * sol(x) * (b + ϵ)

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

Seeded as on the Poisson page.

The outflow sides carry the Neumann relation through ghost points: gridspace(Ωₕ; ghosts = (:xmin, :ymin)) adds one layer of ghosts outside those two sides and constrained fixes each ghost from the relation ∂ₙ(u) = h, where h is a function of the boundary point. The far sides stay Dirichlet and are handed to assemble as side => data. The convective term is now the strong operator ∇̽ₕ, the second-order centred gradient on a non-uniform grid, summed over the directions, because a ghost value enters a difference stencil and not the staggered average of the weak form. In 1D ∇̽ₕ is a bare node rather than a one-element tuple, so the sum is only taken for D > 1:

outflow_flux(x) = -sol(x)

convection(u, ::Val{1}) = ∇̽ₕ(u)
convection(u, ::Val) = sum(∇̽ₕ(u))

function convdiff_system(Ωₕ, D::Int)
    inflow = (:xmax, :ymax, :zmax)[1:D]
    outflow = (:xmin, :ymin, :zmin)[1:D]
    Cₕ = constrained(gridspace(Ωₕ; ghosts = outflow),
        u -> map(side -> side => ∂ₙ(u) => outflow_flux, outflow))
    a = form(Cₕ, Cₕ,
        (u, v) -> innerₕ(-ϵ * Δₕ(u) - b * convection(u, Val(D)), v))
    l = form(Cₕ, v -> innerₕ(rhs, v))
    A, F = assemble(a, l; dirichlet = map(side -> side => sol, inflow))
    return Cₕ, A, F
end

Cₕ, A, F = convdiff_system(Ωₕ, 2)

uₕ = element(Wₕ)
uₕ .= 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 Neumann 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 an outflow side, the three-point forward difference of the lifted solution along the normal axis, with the non-uniform spacings, negated for the outward normal, must equal $$h$$ up to the truncation error, of order $$h^2$$ . flux_residual returns the largest residual over the nodes of one side, axis being 1 for xmin and 2 for ymin:

function flux_residual(Ωₕ, Cₕ, u, axis)
    U = parent(Bramble.lift(Cₕ, u))
    G = Bramble.mesh(Cₕ)
    n = Bramble.npoints(Ωₕ, Tuple)
    xs = ntuple(d -> collect(points(Ωₕ(d))), 2)
    h₁, h₂ = xs[axis][2] - xs[axis][1], xs[axis][3] - xs[axis][2]
    lin = LinearIndices(Bramble.npoints(G, Tuple))
    at(i, j) = U[lin[Bramble.padded_index(G, CartesianIndex(i, j))]]
    along(k, t) = axis == 1 ? at(k, t) : at(t, k)
    return maximum(1:n[3 - axis]) do t
        dₓ = -(2h₁ + h₂) / (h₁ * (h₁ + h₂)) * along(1, t) +
             (h₁ + h₂) / (h₁ * h₂) * along(2, t) - h₁ / (h₂ * (h₁ + h₂)) * along(3, t)
        point = axis == 1 ? (xs[1][1], xs[2][t]) : (xs[1][t], xs[2][1])
        abs(-dₓ - outflow_flux(point))
    end
end

residuals = [flux_residual(Ωₕ, Cₕ, A \ F, axis) for axis in (1, 2)]

2-element Vector{Float64}:
 0.0002478185736558558
 0.0010144662635300783

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

coarse_to_fine = let
    Random.seed!(20260903)
    Ωc = mesh(Ω, (12, 12), (false, false))
    out = Float64[]
    for level in 1:2
        Cc, Ac, Fc = convdiff_system(Ωc, 2)
        push!(out, flux_residual(Ωc, Cc, Ac \ Fc, 1))
        level < 2 && iterative_refinement!(Ωc)
    end
    out
end

2-element Vector{Float64}:
 0.030187781962067728
 0.0072637706835516624

coarse_to_fine[2] < coarse_to_fine[1]

true

Péclet number and flow direction#

The panel below solves a different problem: homogeneous convection-diffusion driven by the boundary values, with no manufactured solution, so the field changes with the parameters.

-\nabla \cdot (\epsilon \nabla u) + \mathbf{\beta} \cdot \nabla u = 0 \text{ in } \Omega = (0,1)^2, \qquad u = 1 \text{ on } \{x = 0\}, \qquad u = 0 \text{ on the rest of } \partial\Omega

with $\epsilon = 1$ and $\mathbf{\beta} = \text{Pe}\,(\cos\theta, \sin\theta)$ ; the sliders set the Péclet number Pe and the flow angle $\theta$ . The panel runs its own JavaScript solver on this problem; only the starting resolution comes from uₕ. With the centred stencil, raise Pe until the cell Péclet number Pe_h passes 1 and the grid-scale oscillation appears; the upwind stencil removes it at the cost of first-order accuracy.

Convection-diffusion, 2D

Convergence with a convective term#

The same pattern as the other two examples: one random coarse mesh per dimension, refined in place with iterative_refinement!. Every dimension has the same mixed problem, Neumann on the sides xmin, ymin, zmin that exist and Dirichlet on the opposite ones:

function convdiff_series(D::Int; n0::Int = 5, levels::Int)
    Ω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
        Wc = gridspace(Ωc)
        _, A_c, F_c = convdiff_system(Ωc, D)
        u_c = element(Wc)
        u_c .= A_c \ F_c

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

Second order in every dimension (fitted orders above 1.9, 1.9 and 1.8). Neither the Neumann relation nor the convective term changes the rate.

Where to go next#

Linear Poisson for the pure-diffusion case this one extends.