Difference, jump and average operators

Bramble provides the finite difference building blocks that discrete schemes are written in: differences, jumps, averages, and their algebraic structures. This tutorial covers:

1. The three operator families and how their names are built.
2. Applying an operator to a grid function, and what happens at the boundary.
3. The same operator as a sparse matrix.
4. Gradients and the other vectorial forms.
5. Summation by parts and skew-symmetry.
6. A convergence study, and the boundary effect that will otherwise spoil it.
7. Interpolating a grid function onto a different mesh, numerically and symbolically.

Every number below was produced by the code shown.

1. The operator families

There are three families:

The jump is the one family with no backward form. A jump belongs to the interface between two cells rather than to a direction of travel across it, so $\llbracket u \rrbracket = u_{i+1} - u_i$ at the interface between $x_i$ and $x_{i+1}$ is a single quantity; a backward jump would name that same interface from the other side and give the same numbers shifted by one index.

1.1 How the names are built

A name is a stem, a direction, and a coordinate:

So D₋ₓ is the backward finite difference along $x$, M₊ᵧ the forward average along $y$, and ∇₋ₕ the backward finite difference in every coordinate, which is the discrete gradient and has that extra name for it.

jump takes no direction, for the reason given above: it is jumpₓ, jumpᵧ, jump₂ and jumpₕ.

2. Applying an operator

An operator takes a VectorElement and returns a new one on the same space.

julia> using Bramble
julia> Ωₕ = mesh(domain(interval(0.0, 1.0)), 5, true);
julia> Wₕ = gridspace(Ωₕ);
julia> points(Ωₕ)5-element Vector{Float64}:
 0.0
 0.25
 0.5
 0.75
 1.0
julia> spacings(Ωₕ)5-element Vector{Float64}:
 0.25
 0.25
 0.25
 0.25
 0.25
julia> uₕ = Rₕ(Wₕ, x -> x^2);
julia> parent(uₕ)5-element Vector{Float64}:
 0.0
 0.0625
 0.25
 0.5625
 1.0

The two backward operators on that grid function:

julia> parent(D₋ₓ(uₕ))5-element Vector{Float64}:
 0.0
 0.25
 0.75
 1.25
 1.75
julia> parent(M₋ₓ(uₕ))5-element Vector{Float64}:
 0.0
 0.03125
 0.15625
 0.40625
 0.78125

Reading the second entry of each: the plain difference is $u_2 - u_1 = 0.0625$, so D₋ₓ divides that by $h_2 = 0.25$ to get $0.25$, and M₋ₓ averages $(u_1 + u_2)/2 = 0.03125$.

The jump has no backward form; forward, it is that plain difference, undivided:

julia> parent(jumpₓ(uₕ))5-element Vector{Float64}:
  0.0625
  0.1875
  0.3125
  0.4375
 -1.0

3. What happens at the boundary

Every operator has one slice where its stencil runs off the grid: the first point for a backward operator, the last for a forward one. There is no neighbour there, so the stencil is truncated, and the finite difference and the jump truncate differently.

The finite difference is zero on its truncated slice, because there is no one-sided stencil to divide by a spacing: the backward one is truncated at $x_1$ and the forward one at $x_5$. The jump instead behaves as if the missing neighbour were zero, which is what makes it agree with its matrix:

julia> parent(D₋ₓ(uₕ))[1]0.0
julia> parent(D₊ₓ(uₕ))[end]0.0
julia> parent(jumpₓ(uₕ))[end]   # -u₅, not 0-1.0

Section 9 shows why this matters in practice.

In two or more dimensions, directional operators apply along the coordinate lines of the tensor grid, and each directional family truncates along its corresponding boundary slice:

4. Operators as matrices

Passing a mesh or a grid space, rather than a grid function, returns the operator itself as a sparse matrix:

julia> A = D₋ₓ(Wₕ);
julia> typeof(A)SparseArrays.SparseMatrixCSC{Float64, Int64}
julia> A * parent(uₕ) ≈ parent(D₋ₓ(uₕ))true

Both routes give the same answer. Applying the operator directly to uₕ is the fast path and is what a time-stepping loop should use; the matrix is what to reach for when assembling a linear system, and it is also how the test suite checks the fast path.

5. Gradients and the other vectorial forms

The ₕ suffix applies the operator along every coordinate and returns a tuple with one entry per dimension. On a one-dimensional mesh it returns the single element itself rather than a one-tuple.

julia> Ω₂ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (4, 4), (true, true));
julia> W₂ = gridspace(Ω₂);
julia> vₕ = Rₕ(W₂, x -> x[1] + 2x[2]);
julia> g = ∇₋ₕ(vₕ);
julia> length(g)2

Away from the truncated slices, g[1] is 1.0 and g[2] is 2.0, the two partial derivatives of $x + 2y$. The same suffix works for the other families as jumpₕ and M₋ₕ, and all of them accept a mesh, a grid space or a grid function.

6. Summation by parts, and Dstar₊ₓ

Continuous integration by parts, $\int u' v = -\int u v'$ for $v$ vanishing on the boundary, has a discrete counterpart, and which forward difference it holds for is not the obvious one. The operator that satisfies it is Dstar₊ₓ: the forward difference divided by the averaged spacing rather than by the forward spacing,

\[\textrm{Dstar}_{+x}(u_h)(i) = \frac{u_{i+1} - u_i}{(h_i + h_{i+1})/2}\]

with the last point truncated to zero, as D₊ₓ is. On a uniform grid $h_i = h_{i+1}$ and it coincides with D₊ₓ; the two differ only where the spacing varies:

julia> Ωₙ = mesh(domain(interval(0.0, 1.0)), 5, true);
julia> set_points!(Ωₙ, [0.0, 0.1, 0.3, 0.7, 1.0])
julia> uₙ = Rₕ(gridspace(Ωₙ), x -> x^2);
julia> parent(D₊ₓ(uₙ))5-element Vector{Float64}:
 0.10000000000000002
 0.39999999999999997
 0.9999999999999999
 1.6999999999999997
 0.0
julia> parent(Dstar₊ₓ(uₙ))5-element Vector{Float64}:
 0.10000000000000002
 0.5333333333333333
 1.333333333333333
 1.4571428571428573
 0.0

The identity is

\[(\textrm{Dstar}_{+x} u_h,\, v_h)_h = -(u_h,\, D_{-x} v_h)_{+x}\]

for any vₕ that vanishes on the boundary. Note which product sits on each side: the left is innerₕ, weighted by the cell measures, and the right is inner₊ₓ, weighted by the staggered ones. Only vₕ has to vanish; uₕ is unconstrained, since the boundary term the identity discards is a product of the two.

julia> Ωᵣ = mesh(domain(interval(0.0, 1.0)), 21, false);   # a random, non-uniform grid
julia> Wᵣ = gridspace(Ωᵣ);
julia> aₕ = Rₕ(Wᵣ, x -> cos(x) + 0.7);                     # not zero at the boundary
julia> bₕ = Rₕ(Wᵣ, x -> sin(pi * x));                      # zero at both ends
julia> innerₕ(Dstar₊ₓ(aₕ), bₕ)-0.3224302783078163
julia> -inner₊ₓ(aₕ, D₋ₓ(bₕ))                              # equal to machine precision-0.32243027830781656
julia> innerₕ(D₊ₓ(aₕ), bₕ)                                # D₊ₓ does not agree-0.31334878829556

It holds per coordinate in two and three dimensions as well, with Dstar₊ᵧ, Dstar₊₂ and their inner products. Dstar₊ₕ returns all coordinates at once, as ∇₊ₕ does.

This is why the operator exists. Energy estimates for these schemes are derived by moving a difference from one factor to the other, and that step is exact only with this pairing: with D₊ₓ it leaves a residual that does not vanish under refinement, since it is a difference of quadrature weights and not a truncation error. Like the other difference families, Dstar₊ can also be had as a sparse matrix: Dstar₊ₓ(Wₕ) is diag(2/(hᵢ + hᵢ₊₁)) times the undivided forward difference $u_{i+1} - u_i$, with an empty last row.

7. The centered difference, Dcₓ

Both one-sided differences reach one point; the centered one reaches both ways, and divides by the whole span its stencil covers:

\[\textrm{Dc}_x(u_h)(i) = \frac{u_{i+1} - u_{i-1}}{h_i + h_{i+1}} = \frac{u_{i+1} - u_{i-1}}{x_{i+1} - x_{i-1}}\]

It is the only operator here that truncates on two slices, since neither the first nor the last point has a neighbour on both sides.

julia> parent(D₋ₓ(uₙ))5-element Vector{Float64}:
 0.0
 0.10000000000000002
 0.39999999999999997
 0.9999999999999999
 1.6999999999999997
julia> parent(D₊ₓ(uₙ))5-element Vector{Float64}:
 0.10000000000000002
 0.39999999999999997
 0.9999999999999999
 1.6999999999999997
 0.0
julia> parent(Dcₓ(uₙ))5-element Vector{Float64}:
 0.0
 0.3
 0.7999999999999999
 1.3
 0.0

Writing the denominator as $x_{i+1} - x_{i-1}$ rather than as a pair of spacings buys two properties that hold on any grid, not only a uniform one.

First, it reproduces an affine function's derivative exactly, since numerator and denominator are then the same quantity:

julia> parent(Dcₓ(Rₕ(gridspace(Ωₙ), x -> 3x + 1)))5-element Vector{Float64}:
 0.0
 3.0
 2.9999999999999996
 3.0000000000000004
 0.0

Second, it is skew-symmetric in innerₕ for grid functions vanishing on the boundary:

julia> Ωₛ = mesh(domain(interval(0.0, 1.0)), 41, false);   # a random, non-uniform grid
julia> Wₛ = gridspace(Ωₛ);
julia> pₕ = Rₕ(Wₛ, x -> sin(pi * x));
julia> qₕ = Rₕ(Wₛ, x -> sin(2pi * x) * x * (1 - x));       # both zero at both ends
julia> innerₕ(Dcₓ(pₕ), qₕ)0.21082197004504277
julia> -innerₕ(pₕ, Dcₓ(qₕ))                               # equal to machine precision0.21082197004504272

The reason is the same cancellation that gives Dstar₊ₓ its identity in section 6: innerₕ weights point $i$ by the cell measure $(h_i + h_{i+1})/2$, which is exactly half the centered denominator. The weights cancel, and the left side collapses to

\[\tfrac{1}{2} \sum_i (u_{i+1} - u_{i-1})\, v_i\]

which shifting the index by one turns into minus the right side. Unlike the Dstar₊ₓ identity, which needs only vₕ to vanish, this one needs both: the discarded boundary term is symmetric in the two.

Accuracy follows the usual rule: the centered difference approximates the derivative at the midpoint of its stencil, which is $x_i$ only when the two spacings match. So it is second order on a uniform grid and first order otherwise, where the one-sided differences are first order on both. Like every other family, Dcₓ accepts a mesh or a grid space for the matrix and a grid function to apply it; Dcₕ gives every coordinate at once. Both end rows of the matrix are empty, which is the truncation.

8. Second order on a non-uniform grid, Dₕₓ

Dcₓ is second order only on a uniform grid. The fix is to take the same two one-sided differences and weight them by the opposite spacings:

\[\textrm{D}_{hx}(u_h)(i) = \frac{h_i}{h_i + h_{i+1}}\, D_{-x} u_h(x_{i+1}) + \frac{h_{i+1}}{h_i + h_{i+1}}\, D_{-x} u_h(x_i)\]

Compare that with Dcₓ, which is the same combination with the weights the other way round. When $h_i = h_{i+1}$ the two agree, and both reduce to the mean of $D_{-x}$ and $D_{+x}$; they part company only where the spacing varies.

The swap buys exactness on quadratics rather than only on affine functions, on any grid. With $u = x^2$ the weighted sum telescopes to $2 x_i (h_i + h_{i+1})$, and the denominator cancels:

julia> parent(Dcₓ(uₙ))5-element Vector{Float64}:
 0.0
 0.3
 0.7999999999999999
 1.3
 0.0
julia> parent(Dₕₓ(uₙ))5-element Vector{Float64}:
 0.10000000000000002
 0.2
 0.5999999999999999
 1.3999999999999997
 1.6999999999999997
julia> 2 .* points(Ωₙ)   # Dₕₓ hits this exactly in the interior5-element Vector{Float64}:
 0.0
 0.2
 0.6
 1.4
 2.0

That one order of extra exactness is one order of extra accuracy. Differencing $\sin$ against $\cos$ on a random grid, refined by halving every interval so the grids stay nested:

Dₕₓ is not skew-symmetric, so Dcₓ remains the one to reach for when the scheme needs that structure and Dₕₓ the one to reach for when it needs the order. Both accept a mesh or a grid space for the matrix and a grid function to apply it, and ∇ₕ gives every coordinate at once, the centered counterpart of ∇₋ₕ and ∇₊ₕ. They differ at the boundary: Dcₓ truncates both end rows to zero, while Dₕₓ has no truncated-boundary convention of its own and falls back to D₊ₓ/D₋ₓ there instead.

9. A convergence study, and the boundary

D₋ₓ is first order, so the error against a known derivative should fall by a factor of ten each time the grid is refined by ten. Measuring it naively does not show that:

for n in (11, 101, 1001, 10001)
    Ω = mesh(domain(interval(0.0, 1.0)), n, true)
    W = gridspace(Ω)
    u = Rₕ(W, sin)
    e = D₋ₓ(u) - Rₕ(W, cos)
    println(n, "  ", normₕ(e))
end

11  0.22500131738444287
101  0.07075844941899423
1001  0.022362202707526486
10001  0.007071116010228156

Over every point the observed order is one half, not one. The cause is section 3: D₋ₓ(uₙ)[1] is 0.0 while $\cos(0) = 1$, so that one point contributes an error of $1$ no matter how fine the grid is. It carries a weight of about $h/2$ in the discrete norm, so it alone contributes about $\sqrt{h/2}$, which is exactly the half order observed.

Excluding that single truncated point recovers the expected first order. So when measuring convergence, or assembling a scheme, treat the truncated slice explicitly: that is where the boundary condition belongs, and leaving the operator's truncated value in place silently halves the observed order.

10. Interpolation between different meshes

Every operator so far maps a grid space to itself. πₕ is the one that does not: it moves a grid function from one mesh to a genuinely different one, which is what makes a heterogeneous composite space — one whose leaves are built over different meshes — useful for more than indexing. Named after Rₕ/Rₕ!'s own convention: πₕ/πₕ! are the numeric pair here, and the same name πₕ — one argument fewer — is also the symbolic wrapper the next tutorial uses, told apart by argument count.

The idea is the standard piecewise (multi)linear interpolant: to read a value at a physical point $x$, find the source mesh's cell containing it (locate_cell) and blend the values at that cell's corners, weighted by how close $x$ is to each one.

The four weights always sum to $1$ — a partition of unity — so the interpolant never overshoots the range of the four corner values. interpolate_at computes this directly at one point; πₕ/πₕ! apply it at every point of a destination space, and are exactly Rₕ/Rₕ! applied to the interpolant as an ordinary function of position — restricting a continuous function and interpolating a discrete one are the same mechanism, πₕ is just the case where that function happens to be another grid function's own interpolant:

julia> Ωbig = mesh(domain(box((0.0, 0.0), (1.0, 1.0))), (10, 10), (true, true));
julia> Ωsmall = mesh(domain(box((0.0, 0.0), (1.0, 1.0))), (4, 4), (true, true));
julia> Wbig, Wsmall = gridspace(Ωbig), gridspace(Ωsmall);
julia> src = Rₕ(Wsmall, x -> x[1] + x[2]);   # affine, so the interpolant is exact
julia> dest = πₕ(Wbig, src);
julia> exact = Rₕ(Wbig, x -> x[1] + x[2]);
julia> maximum(abs, parent(dest) .- parent(exact))2.220446049250313e-16

Once πₕ returns an ordinary VectorElement, every operator above just applies to it as normal — D₋ₓ(dest), M₋ₓ(dest), a bilinear form, anything.

As a matrix

Like D₋ₓ(Wₕ) above, the interpolant is also available as a matrix — but between the two spaces rather than one, and rectangular rather than square, since Wdest and Wsrc generally carry a different number of degrees of freedom:

julia> P = interpolation_matrix(Wbig, Wsmall);
julia> size(P)(100, 16)
julia> P * parent(src) ≈ parent(dest)true

Each row of P has at most $2^D$ nonzero entries — one destination point's corner weights — so it is genuinely sparse, and it is always a SparseMatrixCSC regardless of either space's own backend matrix_type: unlike the shift-based matrices above, a destination point's source cell has no regular diagonal structure to exploit, so P is assembled directly from locate_cell rather than composed from shift.

πₕ!(dest, src) re-locates every destination point's cell on every call, which is wasted work when the same two meshes are interpolated between repeatedly (a time loop moving a coefficient across two composite leaves, say). Build P once and pass it to πₕ! instead: zero allocations, no locate_cell search, just mul! under the hood — the same "build the pattern once" split allocate_system_matrix/assemble! already use.

julia> dest2 = similar(dest);
julia> πₕ!(dest2, P, src);
julia> parent(dest2) ≈ parent(dest)true

Composing symbolically, inside a form

The same name, one argument fewer, is also the symbolic counterpart: πₕ(uₕ) wraps a grid function's interpolant as an AST source, so it composes with the same operators any other source does — D₋ₓ(πₕ(uₕ)), M₋ₓ(πₕ(uₕ)) — and can appear on the left of innerₕ inside a form. Dispatch tells the two πₕ apart by argument count, Wₕ/src against uₕ alone, not by a different name. See the forms tutorial for a worked example against a heterogeneous composite space.