Streamline-Upwind/Petrov-Galerkin (SUPG) advection operator

This is an example of the 1D scalar SUPG advection problem.

Problem formulation

The advection problem is given by

\[\dfrac{\partial u}{\partial t} + c \dfrac{\partial u}{\partial x} = 0\]

where $u = u(x,t)$ and $c$ is its associated phase speed.

The SUPG advection weak form is given by

\[\int_\Omega v \left(\dfrac{\partial u}{\partial t} \right) dV - \int_\Omega \dfrac{\partial v}{\partial x} \left(c u - c τ \left( \dfrac{\partial u}{\partial t} + c \dfrac{\partial u}{\partial x} \right) \right) dV = 0, \forall v \in V\]

for an appropriate test space $V \subseteq H^1 \left( \Omega \right)$ on the domain. In this weak formulation, boundary terms have been omitted, as they are not present on the infinite grid for Local Fourier Analysis. The SUPG stabilization is controlled by the parameter $τ$, where $τ = 0$ gives the classical Galerkin formulation and $τ = \dfrac{h}{2}$ gives a nodally exact solution to the steady advection equation with source when using linear elements (this can be extended to advection-diffusion with a further scaling that depends on the cell Péclet number).

For discussion on SUPG, see Thomas JR Hughes (1979), Alexander N Brooks, Thomas JR Hughes (1982), and Christian H Whiting, Kenneth E Jansen, Saikat Dey (2003).

LFAToolkit code

The symbol of the continuous advection operator $u_t + c u_x = 0$ applied to the Fourier mode $e^{i\theta x}$ is $i\theta$. The finite element discretization yields $M u_t + A u = 0$, and thus we are interested in the symbol of $-M^{-1} A$. One may compare the continuous spectrum with the discrete symbol, which is necessarily periodic on $[-\pi, \pi)$, to understand the behavior for high wave numbers, including the high frequencies that will limit stable time steps. To understand dispersion within the resolved frequencies, we instead plot the phase speed $\lambda/\theta$, which should be very close to $c$ through the resolved frequencies. Here we show the SUPG advection operator on $H^1$ Lagrange basis.

# ---------------------------------------------------------------
# SUPG advection example following Melvin, Staniforth and Thuburn
# Q.J:R. Meteorol. Soc. 138: 1934-1947, Oct (2012)
# non polynomial advection example using transplanted (transformed)
# quadrature formulas from conformal maps following Hale and
# Trefethen (2008) SIAM J. NUMER. ANAL. Vol. 46, No. 2
# In this example, we can call other mapping options available
# i.e, sausage_transformation(9), kosloff_talezer_transformation(0.98)
# otherwise, run without mapping
# ---------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
Δx = 1.0
mesh = Mesh1D(Δx)
dxdξ = Δx / 2 # 2 comes from quadrature domain of [-1,1]
dξdx = 1 / dxdξ
det_dxdξ = dxdξ # Determinant of mapping Jacobian

p = 2;
q = p;
collocate = false
mapping = nothing
basis =
    TensorH1LagrangeBasis(p, q, 1, 1, collocatedquadrature = collocate, mapping = mapping)

# frequency set up
numbersteps = 100
θ_min = 0
θ_max = (p - 1) * π
θ = LinRange(θ_min, θ_max, numbersteps)

# associated phase speed
c = 1.0

# Tau scaling for SUPG
τ = 0.5 / (p - 1) # 0 returns Galerkin method; empirically, 0.25 is too big for p = 3 (quadratic)

# weak form for SUPG advection
function supgadvectionweakform(U::Matrix{Float64}, w::Array{Float64})
    u = U[1, :]
    du = U[2, :]
    dv = dξdx * (c * u - c * τ * (c * du * dξdx)) * w[1] * det_dxdξ
    return [dv]
end

# supg advection operator
inputs = [
    OperatorField(
        basis,
        [EvaluationMode.interpolation, EvaluationMode.gradient],
        "advected field",
    ),
    OperatorField(basis, [EvaluationMode.quadratureweights], "quadrature weights"),
]
outputs = [OperatorField(basis, [EvaluationMode.gradient])]

supgadvection = Operator(supgadvectionweakform, mesh, inputs, outputs)

# supg mass operator
function supgmassweakform(udot::Array{Float64}, w::Array{Float64})
    v = udot * w[1] * det_dxdξ
    dv = dξdx * c * τ * udot * w[1] * det_dxdξ
    return ([v; dv],)
end
supgmass = Operator(
    supgmassweakform,
    mesh,
    [
        OperatorField(basis, [EvaluationMode.interpolation], "uₜ"),
        OperatorField(basis, [EvaluationMode.quadratureweights], "quadrature weights"),
    ],
    [OperatorField(basis, [EvaluationMode.interpolation, EvaluationMode.gradient])],
)

# compute operator symbols
function advection_supg_symbol(θ)
    A = computesymbols(supgadvection, [θ]) # transform from reference to physical on dx=1 grid
    M = computesymbols(supgmass, [θ]) # mass matrix
    return sort(imag.(eigvals(-M \ A)))
end

eigenvalues = hcat(advection_supg_symbol.(θ)...)'

# ---------------------------------------------------------------