Difference, jump and average operators#
You will apply a difference, a jump or an average to a discrete function, see what each does at the edge of the grid, and get the operator as a matrix. The mesh tutorial and the space tutorial come first, for meshes, grid spaces and discrete functions. The form tutorial assembles these operators into linear systems.
Apply an operator#
Take a grid of five points and the function \(u(x) = x^2\) restricted to it:
Ωₕ = mesh(domain(interval(0.0, 1.0)), 5)
Wₕ = gridspace(Ωₕ)
uₕ = Rₕ(Wₕ, x -> x^2)
points(Ωₕ), parent(uₕ)([0.0, 0.25, 0.5, 0.75, 1.0], [0.0, 0.0625, 0.25, 0.5625, 1.0])An operator takes a VectorElement such as uₕ and returns a new one on the same space. The backward difference ∇ₕ divides \(u_i - u_{i-1}\) by the spacing, and the backward average Mₕ takes \((u_{i-1} + u_i)/2\):
parent(∇ₕ(uₕ)), parent(Mₕ(uₕ))([0.0, 0.25, 0.75, 1.25, 1.75], [0.0, 0.03125, 0.15625, 0.40625, 0.78125])At the second entry the plain difference is \(u_2 - u_1 = 0.0625\) and the spacing is \(h_2 = 0.25\), so ∇ₕ gives \(0.25\). The average is \((u_1 + u_2)/2 = 0.03125\). The first entry of the difference is \(0\), which the next section explains.
The jump is the plain difference across an interface, with no division:
parent(jumpₕ(uₕ))5-element Vector{Float64}:
0.0625
0.1875
0.3125
0.4375
-1.0Here jumpₕ is the forward difference \(u_{i+1} - u_i\), so its first entry is \(0.0625\). The last entry, \(-1\), is explained below.
Try this
Replace x^2 by x. The derivative is \(1\), so ∇ₕ returns ones everywhere except at the first point, and jumpₕ returns the constant \(0.25\), the spacing, except at the last point.
Names are built from a stem and a direction (naming scheme).
| Family | Meaning | Backward form |
|---|---|---|
| finite difference | a difference divided by the spacing, so it approximates \(\partial u / \partial x\) | \(\dfrac{u_i - u_{i-1}}{h_i}\) |
| jump | the plain difference across an interface, undivided | \(u_{i+1} - u_i\) |
| average | the mean of a point and its neighbour | \(\dfrac{u_{i-1} + u_i}{2}\) |
| index shift | the neighbour's value, moved onto the point | \(u_{i-1}\) |
The jump has no backward form: it belongs to the interface between \(x_i\) and \(x_{i+1}\), so \(\llbracket u \rrbracket = u_{i+1} - u_i\) is one quantity, and a backward jump would name the same interface from the other side.
Four differences compared#
Four families combine the same one-sided differences in different ways:
| Operator | Stencil | Order, non-uniform | Order, uniform | Boundary |
|---|---|---|---|---|
D₋ |
\(\{i-1, i\}\) | 1 | 1 | truncated to 0 at the first point |
Dc |
\(\{i-1, i+1\}\), divided by the full span | 1 | 2 | truncated to 0 at both ends |
D̽ |
\(\{i-1, i, i+1\}\), a weighted mean of two backward differences | 2 | 2 | no convention of its own: falls back to D₊/D₋ |
D̃ |
\(\{i, i+1\}\), divided by the averaged spacing | 1 on smoothly graded grids, not consistent where neighbouring spacings differ by an O(1) ratio (random grids); summation by parts holds exactly | 1 | truncated to 0 at the last point, like D₊ |
On a uniform grid \(h_i = h_{i+1}\): D̽ and Dc collapse to the mean of D₋ and D₊, and D̃ collapses to D₊. They separate only where the spacing varies, which is why a uniform-grid benchmark cannot tell them apart. Use D̃ with innerₕ, for the summation by parts.
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.
The finite difference is zero on its truncated slice, because there is no one-sided stencil to divide by a spacing. The jump instead behaves as if the missing neighbour were zero.
parent(∇ₕ(uₕ))[1], parent(D₊ₓ(uₕ))[end], parent(jumpₕ(uₕ))[end](0.0, 0.0, -1.0)The last entry is \(-u_5 = -1\), which makes the jump agree with its matrix.
The index shifts follow the jump: an off-grid neighbour reads as zero, so S₊ₕ is zero at the last point and S₋ₕ at the first.
parent(S₊ₕ(uₕ)), parent(S₋ₕ(uₕ)), parent(S₊ₕ(uₕ)) - parent(uₕ) == parent(jumpₕ(uₕ))([0.0625, 0.25, 0.5625, 1.0, 0.0], [0.0, 0.0, 0.0625, 0.25, 0.5625], true)That convention gives three identities that hold at every point, the boundary included. The matrix of S₊ is the transpose of the matrix of S₋, S₊(u) - u is jump(u), and u - S₋(u) is the undivided backward difference. The exception is D₊: it is zero at the last point, so S₊(u) - u is not \(h\) times D₊(u) there. Inside a form the rule is the same, and a shift of a composed operand, such as S₊ₓ(D₋ₓ(u)), reads zero wherever the shifted stencil leaves the grid.
In two or more dimensions the directional operators apply along the coordinate lines of the tensor grid, and each directional family truncates along its own boundary slice:
Operators as matrices#
Pass a mesh or a grid space instead of a grid function and the operator comes back as a sparse matrix:
A = ∇ₕ(Wₕ)
typeof(A), A * parent(uₕ) ≈ parent(∇ₕ(uₕ))(SparseArrays.SparseMatrixCSC{Float64, Int64}, true)Both routes give the same answer. Applying the operator to uₕ is the fast path for time loops; the matrix is for assembling a linear system.
Gradients in two dimensions#
The ₕ suffix applies an operator along every coordinate and returns a tuple with one entry per dimension. On a one-dimensional mesh it returns the single element, which is why ∇ₕ(uₕ) above was not wrapped in a tuple.
Ω₂ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (4, 4))
W₂ = gridspace(Ω₂)
vₕ = Rₕ(W₂, x -> x[1] + 2x[2])
g = ∇ₕ(vₕ)
length(g), round.(reshape(parent(g[1]), 4, 4); digits = 8), round.(reshape(parent(g[2]), 4, 4); digits = 8)(2, [0.0 0.0 0.0 0.0; 1.0 1.0 1.0 1.0; 1.0 1.0 1.0 1.0; 1.0 1.0 1.0 1.0], [0.0 2.0 2.0 2.0; 0.0 2.0 2.0 2.0; 0.0 2.0 2.0 2.0; 0.0 2.0 2.0 2.0])The rows run along \(x\) and the columns along \(y\). Away from the truncated slices g[1] is \(1\) and g[2] is \(2\), the two partial derivatives of \(x + 2y\); g[1] is zero on the first row and g[2] on the first column. The suffix works for the other families as jumpₕ and Mₕ, and all of them accept a mesh, a grid space or a grid function.
Try this
Change the function to x -> x[1]^2 and print both matrices again. g[1] now grows down the rows, while g[2] is zero everywhere.
The vectorial operators of a composite space contract with LinearAlgebra's ⋅ and × to give the divergence and the curl:
using LinearAlgebra: ⋅, ×
uv = Rₕ(vector_gridspace(Ω₂), x -> (x[1] + 2x[2], 3x[1] - x[2]))
parent(∇ₕ ⋅ uv) == parent(divₕ(uv)), parent(∇ₕ × uv) == parent(curlₕ(uv))(true, true)Ghost points#
At a ghosted boundary the stencils read the ghost as a neighbour. The facet operators ∂ₙ, ∂ₙₙ and jumpₙ, the corner rules and the Laplacian inside a form are in the ghost point tutorial.