sbp.jl

# Summation-by-parts (SBP) coefficients for the 2-1 stencil, taken
# from "coeffs/Derivatives_2_1". The notation "2-1" means that this
# operator is 2nd order accurate in the interior and 1st order
# accurate at the boundary.
#
# Generally, it is impossible to find SBP stencils that have the same
# order of accuracy everywhere. Luckily, since there is only a finite
# number of boundary points which does not increase when the
# resolution is increase, the overall observed order (in the continuum
# limit) is one more than the order at the boundary. That is, this
# stencil is overall 2nd order accurate.

"Integration weight at the boundary"
const a = [1.0 / 2.0]

"Derivative stencils near the boundary"
const q = [-1.0 -1.0/2.0
           1.0 0
           0 1.0/2.0]

# The above matrix `q` defines 2 stencils, each using 3 points. All
# stencils are to be evaluated with the leftmost points of the domain
# as input.
#
# In this case, the last (second) stencil happens to be the regular
# second-order finite difference. This is not true in general.

"Calculate a derivative"
function deriv(f::AbstractVector, h::Real)
    # Number of boundary points that need special treatment
    nb = size(q, 2)
    # We need at least twice as many points as there are boundary points
    @assert length(f) ≥ 2nb

    df = similar(f)
    # Left boundary
    df[1:nb] = q * f[1:size(q, 1)] / h
    # Right boundary (reverse grid points and change sign)
    df[end:-1:(end - nb + 1)] = -q * f[end:-1:(end - size(q, 1) + 1)] / h
    # Interior (standard second order derivative)
    for i in (nb + 1):(length(f) - nb)
        df[i] = (f[i - 1] + f[i + 1]) / 2h
    end

    return df
end

"Integrate (quadrature) over all points"
function integrate(f::AbstractVector, h::Real)
    # Number of boundary points that need special treatment
    nb = length(a)
    @assert length(f) ≥ 2nb

    s = zero(eltype(f))
    # Left boundary
    for i in 1:nb
        s += a[i] * f[i]
    end
    # Right boundary (reverse grid points)
    for i in 1:nb
        s += a[i] * f[length(f)-i+1]
    end
    # Interior
    s += sum(f[nb+1:end-nb])
    s *= h

    return s
end