Example: approximating derivatives using DG
MeshData
can be used to compute DG derivatives. Suppose $f$ is a differentiable function and the domain $\Omega$ can be decomposed into non-overlapping elements $D^k$. The approximation of $\frac{\partial f}{\partial x}$ can be approximated using the following formulation: find piecewise polynomial $u$ such that for all piecewise polynomials $v$
\[\int_{\Omega} u v = \sum_k \left( \int_{D^k} \frac{\partial u}{\partial x}v + \int_{\partial D^k} \frac{1}{2} \left[u\right]n_x v \right)\]
Here, $\left[u\right] = u^+ - u$ denotes the jump across an element interface, and $n_x$ is the $x$-component of the outward unit normal on $D^k$.
Discretizing the left-hand side of this formulation yields a mass matrix. Inverting this mass matrix to the right hand side yields the DG derivative. We show how to compute it for a uniform triangular mesh using MeshData
and StartUpDG.jl
.
We first construct the triangular mesh and initialize md::MeshData
.
using StartUpDG
using Plots
N = 3
K1D = 8
rd = RefElemData(Tri(), N)
VXY, EToV = uniform_mesh(Tri(), K1D)
md = MeshData(VXY, EToV, rd)
We can approximate a function $f(x, y)$ using interpolation
f(x, y) = exp(-5 * (x^2 + y^2)) * sin(1 + pi*x) * sin(2 + pi*y)
(; x, y) = md
u = @. f(x, y)
or using quadrature-based projection
(; Pq ) = rd
(; x, y, xq, yq ) = md
u = Pq * f.(xq, yq)
We can use scatter
in Plots.jl to quickly visualize the approximation. This is not intended to create a high quality image (see other libraries, e.g., Makie.jl
,VTK.jl
, or Triplot.jl
for publication-quality images).
(; Vp ) = rd
xp, yp, up = Vp * x, Vp * y, Vp * u # interp to plotting points
scatter(xp, yp, uxp, zcolor=uxp, msw=0, leg=false, ratio=1, cam=(0, 90))
Both interpolation and projection create a matrix u
of size $N_p \times K$ which contains coefficients (nodal values) of the DG polynomial approximation to $f(x, y)$. We can approximate the derivative of $f(x, y)$ using the DG derivative formulation
function dg_deriv_x(u, rd::RefElemData, md::MeshData)
(; Vf, Dr, Ds, LIFT ) = rd
(; rxJ, sxJ, J, nxJ, mapP ) = md
uf = Vf * u
ujump = uf[mapP] - uf
# derivatives using chain rule + lifted flux terms
ux = rxJ .* (Dr * u) + sxJ .* (Ds * u)
dudxJ = ux + LIFT * (.5 * ujump .* nxJ)
return dudxJ ./ J
end
We can visualize the result as follows:
dudx = dg_deriv_x(u, rd, md)
uxp = Vp * dudx
scatter(xp, yp, uxp, zcolor=uxp, msw=0, leg=false, ratio=1, cam=(0,90))
Plots of the polynomial approximation $u(x,y)$ and the DG approximation of $\frac{\partial u}{\partial x}$ are given below
⠀