10: Upwind FD SBP schemes
General tensor product SBP methods are supported via the DGMulti solver in a reasonably complete way, see the previous tutorial. Nevertheless, there is also experimental support for SBP methods with other solver and mesh types.
The first step is to set up an SBP operator. A classical (central) SBP operator can be created as follows.
using Trixi
D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
derivative_order = 1, accuracy_order = 2,
xmin = 0.0, xmax = 1.0, N = 11)SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] using 11 nodes
and coefficients of Mattsson, Nordström (2004)
Summation by parts operators for finite difference approximations of second
derivatives.
Journal of Computational Physics 199, pp. 503-540.Instead of prefixing the source of coefficients MattssonNordström2004(), you can also load the package SummationByPartsOperators.jl. Either way, this yields an object representing the operator efficiently. If you want to compare it to coefficients presented in the literature, you can convert it to a matrix.
Matrix(D_SBP)11×11 Matrix{Float64}:
-10.0 10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-5.0 0.0 5.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 -5.0 0.0 5.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 -5.0 0.0 5.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 -5.0 0.0 5.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 -5.0 0.0 5.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 -5.0 0.0 5.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 -5.0 0.0 5.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 -5.0 0.0 5.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -5.0 0.0 5.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -10.0 10.0Upwind SBP operators are a concept introduced in 2017 by Ken Mattsson. You can create them as follows.
D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
derivative_order = 1, accuracy_order = 2,
xmin = 0.0, xmax = 1.0, N = 11)Upwind SBP first-derivative operators of order 2 on a grid in [0.0, 1.0] using 11 nodes
and coefficients of Mattsson2017Upwind operators are derivative operators biased towards one direction. The "minus" variants has a bias towards the left side, i.e., it uses values from more nodes to the left than from the right to compute the discrete derivative approximation at a given node (in the interior of the domain). In matrix form, this means more non-zero entries are left from the diagonal.
Matrix(D_upw.minus)11×11 Matrix{Float64}:
-10.0 10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-10.0 10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
5.0 -20.0 15.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 5.0 -20.0 15.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 5.0 -20.0 15.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 5.0 -20.0 15.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 5.0 -20.0 15.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 5.0 -20.0 15.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 5.0 -20.0 15.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.0 -16.0 10.0 2.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.0 -50.0 30.0Analogously, the "plus" variant has a bias towards the right side.
Matrix(D_upw.plus)11×11 Matrix{Float64}:
-30.0 50.0 -20.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-2.0 -10.0 16.0 -4.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 -15.0 20.0 -5.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 -15.0 20.0 -5.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 -15.0 20.0 -5.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 -15.0 20.0 -5.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 -15.0 20.0 -5.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 -15.0 20.0 -5.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -15.0 20.0 -5.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -10.0 10.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -10.0 10.0For more information on upwind SBP operators, please refer to the documentation of SummationByPartsOperators.jl and references cited there.
The basic idea of upwind SBP schemes is to apply a flux vector splitting and use appropriate upwind operators for both parts of the flux. In 1D, this means to split the flux
\[f(u) = f^-(u) + f^+(u)\]
such that $f^-(u)$ is associated with left-going waves and $f^+(u)$ with right-going waves. Then, we apply upwind SBP operators $D^-, D^+$ with an appropriate upwind bias, resulting in
\[\partial_x f(u) \approx D^+ f^-(u) + D^- f^+(u)\]
Note that the established notations of upwind operators $D^\pm$ and flux splittings $f^\pm$ clash. The right-going waves from $f^+$ need an operator biased towards their upwind side, i.e., the left side. This upwind bias is provided by the operator $D^-$.
Many classical flux vector splittings have been developed for finite volume methods and are described in the book "Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction" of Eleuterio F. Toro (2009), DOI: 10.1007/b79761. One such a well-known splitting provided by Trixi.jl is splitting_steger_warming.
Trixi.jl comes with several example setups using upwind SBP methods with flux vector splitting, e.g.,
Package versions
These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["Trixi", "SummationByPartsOperators"],
mode = PKGMODE_MANIFEST)Julia Version 1.10.6
Commit 67dffc4a8ae (2024-10-28 12:23 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 7763 64-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/Trixi.jl/Trixi.jl/docs/Manifest.toml`
[9f78cca6] SummationByPartsOperators v0.5.70
[a7f1ee26] Trixi v0.9.3 `~/work/Trixi.jl/Trixi.jl`This page was generated using Literate.jl.