From 909e3d509892fa7ab8d87bfdc98d176e441355af Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Mon, 19 Aug 2024 19:58:02 +0000
Subject: [PATCH] dg1d: Add Bsplines
(https://git.tpi.uni-jena.de/dg/dg1d.jl/-/merge_requests/210)
This adds multiple things:
- A `Bspline` module for working with 2nd order Bsplines.
- Extends `DGElement` with the `kind=:modal_bspline2` option.
- Updates docstrings and refactors some tests.
- Adds a `LG` module for Legendre-Gauss quadrature.
- Adds a plotting script `plot/bspline.jl` to tinker with the `Bspline2` approximation.
Also adds a new option to the `Mesh` parameters section:
- `kind = "modal_bspline2"`: Utilize a 2nd order modal Bspline Ansatz for the DG evolution.
---
plot/bspline2.jl | 145 +++++++++++++++++++
src/bernstein.jl | 2 +-
src/bspline.jl | 209 +++++++++++++++++++++++++++
src/dg1d.jl | 2 +
src/dgelement.jl | 98 ++++++++++---
src/legendre.jl | 2 +-
src/lg.jl | 86 +++++++++++
src/lgl.jl | 2 +-
src/parameters.jl | 62 +++++---
test/UnitTests/src/test_dgelement.jl | 76 ++++++----
test/UnitTests/src/test_utils.jl | 3 +-
11 files changed, 613 insertions(+), 74 deletions(-)
create mode 100644 plot/bspline2.jl
create mode 100644 src/bspline.jl
create mode 100644 src/lg.jl
diff --git a/plot/bspline2.jl b/plot/bspline2.jl
new file mode 100644
index 00000000..75baad63
--- /dev/null
+++ b/plot/bspline2.jl
@@ -0,0 +1,145 @@
+using dg1d
+using dg1d.Jacobi
+using GLMakie
+using LinearAlgebra
+
+
+# bs = dg1d.Bspline.Bspline2(5)
+x, _ = dg1d.LGL.rule(5)
+bs = dg1d.Bspline.Bspline2(x)
+xx = collect(range(-2.0,2.0,length=300))
+Bs = [ dg1d.Bspline.polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:bs.Npts+1 ]
+∂Bs = [ dg1d.Bspline.derivative_polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:bs.Npts+1 ]
+
+Vs = [ dg1d.Bspline.polynomial(bs.xs[2+i], j, bs) for i in 1:bs.Npts, j in 1:bs.Npts+1 ]
+# display(Vs)
+∂Vs = [ dg1d.Bspline.derivative_polynomial(bs.xs[2+i], j, bs) for i in 1:bs.Npts, j in 1:bs.Npts+1 ]
+# display(∂Vs)
+
+M = dg1d.Bspline.mass_matrix(bs)
+S = dg1d.Bspline.stiffness_matrix(bs)
+ST = transpose(S)
+BB = similar(M); fill!(BB, 0.0)
+o = 2
+for i in 1:bs.Npts+1, j in 1:bs.Npts+1
+ for k in 1:bs.Npts-1
+ xl, xr = bs.xs[o+k], bs.xs[o+k+1]
+ Nil, Nir = dg1d.Bspline.polynomial(xl, i, bs), dg1d.Bspline.polynomial(xr, i, bs)
+ Njl, Njr = dg1d.Bspline.polynomial(xl, j, bs), dg1d.Bspline.polynomial(xr, j, bs)
+ BB[i,j] += Nir*Njr - Nil*Njl
+ end
+end
+RBR = dg1d.purge_zeros(S+ST)
+display(M)
+display(S)
+display(RBR)
+display(BB)
+
+if any(isnan, Bs)
+ error("polynomial produces NaNs")
+end
+if any(isnan, ∂Bs)
+ error("derivative_polynomial produces NaNs")
+end
+if any(isnan, Vs)
+ error("polynomial produces NaNs")
+end
+if any(isnan, ∂Vs)
+ error("derivative_polynomial produces NaNs")
+end
+
+fig = Figure()
+ax = Axis(fig[1,1])
+vlines!(ax, bs.xs, color=(:black,0.5), linestyle=:dash)
+
+for (i,(B,V)) in enumerate(zip(eachcol(Bs),eachcol(Vs)))
+ lines!(ax, xx, B, label="$i", color=Makie.wong_colors()[i])
+ scatter!(ax, bs.xs[3:2+bs.Npts], V, color=Makie.wong_colors()[i])
+end
+PoU = vec(sum(Bs, dims=2))
+lines!(ax, xx, PoU, label=L"\Sigma_n B_n", linestyle=:dot)
+
+# for (i,(∂B,∂V)) in enumerate(zip(eachcol(∂Bs),eachcol(∂Vs)))
+# lines!(ax, xx, ∂B, label="$i", linestyle=:dash, color=Makie.wong_colors()[i])
+# scatter!(ax, bs.xs[3:2+bs.Npts], ∂V, color=Makie.wong_colors()[i])
+# end
+
+# FD_∂B = diff(Bs, dims=1) ./ diff(xx)
+# xx_ctr = xx[1:end-1] .+ (xx[2]-xx[1])/2
+# for (i,∂B) in enumerate(eachcol(FD_∂B))
+# lines!(ax, xx_ctr, ∂B, label="$i", color=Makie.wong_colors()[i])
+# end
+
+axislegend(ax)
+display(fig)
+
+
+# # fn = x -> x < 0.15 ? -x+1 : 0.0
+# # fn = x -> x < 0.0 ? -x : 0.0
+Npts = 4
+xs, _ = dg1d.LGL.rule(Npts)
+xx = collect(range(extrema(xs)...,length=100))
+fn = x -> x < -0.5 ? -x+1 : 0.0
+fx = fn.(xs)
+
+knots = vcat(xs[1], (xs[1:end-1].+xs[2:end])./2, xs[end])
+avg_fx = vcat(fx[1], (fx[1:end-1].+fx[2:end])./2, fx[end])
+bs = dg1d.Bspline.Bspline2(xs)
+Vs = dg1d.Bspline.vandermonde_matrix(knots, bs)
+intrp_Bs = [ dg1d.Bspline.polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:bs.Npts+1 ]
+A_lsqr = transpose(Vs)*Vs
+b_lsqr = transpose(Vs)*avg_fx
+
+# nodal Bspline
+nodal_Bs = Vs \ avg_fx
+fn_nodal_Bs = intrp_Bs * nodal_Bs
+
+# quasi-nodal Bspline
+quasi_nodal_Bs = avg_fx
+fn_quasi_nodal_Bs = intrp_Bs * quasi_nodal_Bs
+
+# corrected nodal Bspline
+c_nodal_Bs = max.(nodal_Bs, 0.0)
+fn_c_nodal_Bs = intrp_Bs * c_nodal_Bs
+
+# least-square Bspline
+lsqr_Bs = A_lsqr \ b_lsqr
+fn_lsqr_Bs = intrp_Bs * lsqr_Bs
+
+# corrected least-sqaure Bspline
+c_lsqr_Bs = max.(lsqr_Bs, 0.0)
+fn_c_lsqr_Bs = intrp_Bs * c_lsqr_Bs
+
+# nodal Lagrange
+intrp_La = dg1d.lagrange_interp_matrix(xx, xs)
+fn_La = intrp_La * fx
+
+# modal Legendre
+Ln = dg1d.Legendre.projection(fn, Npts)
+fn_Ln = dg1d.Legendre.interpolate.(xx, Ref(Ln))
+
+fig = Figure()
+ax = Axis(fig[1,1])
+
+vlines!(ax, xs, color=(:black,0.4), linestyle=:dash)
+# lines!(ax, xx, fn_Ln, color=:blue, label="Legendre")
+lines!(ax, xx, fn_La, color=:orange, label="Lagrange")
+scatterlines!(ax, xs, fx, color=:orange, markersize=22, linestyle=:dash)
+
+# lines!(ax, xx, fn_nodal_Bs, color=:blue, label="nodal Bspline")
+# scatterlines!(ax, knots, nodal_Bs, color=:blue, markersize=18, linestyle=:dash)
+
+lines!(ax, xx, fn_quasi_nodal_Bs, color=:green, label="quasi-nodal Bspline")
+scatterlines!(ax, knots, quasi_nodal_Bs, color=:green, markersize=14, linestyle=:dash)
+
+# lines!(ax, xx, fn_c_nodal_Bs, color=:red, label="corrected nodal Bspline")
+# scatterlines!(ax, knots, c_nodal_Bs, color=:red, markersize=10, linestyle=:dash)
+
+lines!(ax, xx, fn_lsqr_Bs, color=:violet, label="least-square Bspline")
+scatterlines!(ax, knots, lsqr_Bs, color=:violet, markersize=14, linestyle=:dash)
+
+lines!(ax, xx, fn_c_lsqr_Bs, color=:purple, label="corrected least-square Bspline")
+scatterlines!(ax, knots, c_lsqr_Bs, color=:purple, markersize=14, linestyle=:dash)
+
+axislegend(ax)
+display(fig)
diff --git a/src/bernstein.jl b/src/bernstein.jl
index 5ccaec4a..f7311e03 100644
--- a/src/bernstein.jl
+++ b/src/bernstein.jl
@@ -110,7 +110,7 @@ end
interpolate(x::Real, fn::AbstractVector)
Interpolate a Bernstein polynomial with quasi-nodal coefficients `fn` onto `x`.
-`x` is assumed to lie within `[-1,1`].
+`x` is assumed to lie within `[-1,1]`.
The `fn` are assumed to be sampled on the nodes `xs = [ -1.0 + 2*i/N for i in 0:N ]`
where `N = length(fn)-1`.
"""
diff --git a/src/bspline.jl b/src/bspline.jl
new file mode 100644
index 00000000..2ac90140
--- /dev/null
+++ b/src/bspline.jl
@@ -0,0 +1,209 @@
+module Bspline
+
+
+using dg1d
+using Jacobi
+
+
+struct Bspline2
+ Npts::Int64 # number of unique knots
+ xs::Vector{Float64} # the knots, start and end point have multiplicity three
+end
+
+
+function Bspline2(z::Vector{Float64})
+ Npts = length(z)
+ if Npts < 2
+ error("z must contain at least two points, but you provided $Npts")
+ end
+ if any(d -> d < 0.0, diff(z))
+ error("z must be a strictly increasing point set, but you provided $z")
+ end
+ xs = vcat([z[1],z[1]], z, [z[end],z[end]])
+ return Bspline2(Npts, xs)
+end
+
+
+function polynomial(x::Float64, i::Int64, bs::Bspline2)
+ @unpack Npts, xs = bs
+ 1 ≤ i ≤ Npts+1 || return 0.0
+ @inbounds begin
+ xi, xi1, xi2, xi3 = xs[i], xs[i+1], xs[i+2], xs[i+3]
+ end
+ if x < xi || x > xi3
+ return 0.0
+ end
+ if i == 1 && x == xi2
+ # compute derivative on lhs bdry from the right
+ @goto sub3
+ end
+ if i == 2 && x == xi1
+ # compute derivative on lhs bdry from the right
+ @goto sub2
+ end
+ if x ≤ xi1
+ return (x-xi)^2/((xi2-xi)*(xi1-xi))
+ end
+ if x ≤ xi2
+ @label sub2
+ return (x-xi)*(xi2-x)/((xi2-xi)*(xi2-xi1)) + (xi3-x)*(x-xi1)/((xi3-xi1)*(xi2-xi1))
+ end
+ if x ≤ xi3
+ @label sub3
+ return (xi3-x)^2/((xi3-xi1)*(xi3-xi2))
+ end
+ # unreachable
+ return NaN
+end
+
+
+function derivative_polynomial(x::Float64, i::Int64, bs::Bspline2)
+ @unpack Npts, xs = bs
+ 1 ≤ i ≤ Npts+1 || return 0.0
+ @inbounds begin
+ xi, xi1, xi2, xi3 = xs[i], xs[i+1], xs[i+2], xs[i+3]
+ end
+ if x < xi || x > xi3
+ return 0.0
+ end
+ if i == 1 && x == xi2
+ # compute derivative on lhs bdry from the right
+ @goto sub3
+ end
+ if i == 2 && x == xi1
+ # compute derivative on lhs bdry from the right
+ @goto sub2
+ end
+ if x ≤ xi1
+ return 2*(x-xi)/((xi2-xi)*(xi1-xi))
+ end
+ if x ≤ xi2
+ @label sub2
+ return 2*(xi2-x)/((xi2-xi)*(xi2-xi1)) - 2*(x-xi1)/((xi3-xi1)*(xi2-xi1))
+ end
+ if x ≤ xi3
+ @label sub3
+ return -2*(xi3-x)/((xi3-xi1)*(xi3-xi2))
+ end
+ # unreachable
+ return NaN
+end
+
+
+"""
+ interpolate(x::Real, fn::AbstractVector, bs::Bspline2)
+
+Interpolate a Bspline polynomial with quasi-nodal coefficients `fn` onto `x`.
+`x` is assumed to lie within `[-1,1]`.
+The `fn` are assumed to be sampled on the nodes `xs = [ -1.0 + 2*i/N for i in 0:N ]`
+where `N = length(fn)-1`.
+"""
+function interpolate(x::Real, coeffs::AbstractVector, bs::Bspline2)
+ @unpack Npts, xs = bs
+ if length(coeffs) != Npts+1
+ error("expected Npts+1 = $Npts coefficients, received $(length(coeffs))")
+ end
+ return unsafe_interpolate(x, coeffs, bs)
+end
+function unsafe_interpolate(x::Real, coeffs::AbstractVector, bs::Bspline2)
+ @unpack Npts, xs = bs
+ @inbounds xl, xr = xs[1], xs[end]
+ xl ≤ x ≤ xr || return 0.0
+ o = 2 # offset duplicated bdry knots
+ v_xs = view(xs, o .+ (1:Npts))
+ ir = searchsortedfirst(v_xs, x)
+ il = ir-1
+ il < 1 && return coeffs[1] # x == xl
+ is_supp = max(il-1,1):min(ir+1,Npts+1)
+ B = 0.0
+ @inbounds for i in is_supp
+ B += coeffs[i] * polynomial(x, i, bs)
+ end
+ return B
+end
+
+
+"""
+ mass_matrix(bs::Bspline2)
+
+Compute `M_ij = (Bi, Bj)` exactly where `Bi, Bj` are `i`th and `j`th Bspline of order 2.
+"""
+function mass_matrix(bs::Bspline2)
+ @unpack Npts, xs = bs
+ M = zeros(Float64, Npts+1, Npts+1)
+ zs, ws = dg1d.LG.rule(2)
+ o = 2 # offset bdry knots
+ for k in 1:Npts-1
+ il, ir = o+k, o+k+1
+ xl, xr = xs[il], xs[ir]
+ for i in 1:Npts+1, j in 1:Npts+1
+ mij = dg1d.LG.integrate(zs, ws) do z
+ x = ((xr-xl)*z + xr+xl)/2
+ dxdz = (xr-xl)/2
+ Bi = polynomial(x, i, bs)
+ Bj = polynomial(x, j, bs)
+ Bi*Bj*dxdz
+ end
+ M[i,j] += mij
+ end
+ end
+ return M
+end
+
+
+"""
+ stiffness_matrix(bs::Bspline2)
+
+Compute `S_ij = (Bi, dBj/dx)` exactly where `Bi, Bj` are `i`th and `j`th Bspline of order 2.
+"""
+function stiffness_matrix(bs::Bspline2)
+ @unpack Npts, xs = bs
+ S = zeros(Float64, Npts+1, Npts+1)
+ zs, ws = dg1d.LG.rule(3)
+ o = 2 # offset bdry knots
+ for k in 1:Npts-1
+ il, ir = o+k, o+k+1
+ xl, xr = xs[il], xs[ir]
+ for i in 1:Npts+1, j in 1:Npts+1
+ sij = dg1d.LG.integrate(zs, ws) do z
+ x = ((xr-xl)*z + xr+xl)/2
+ dxdz = (xr-xl)/2
+ Bi = polynomial(x, i, bs)
+ ∂Bj = derivative_polynomial(x, j, bs)
+ Bi*∂Bj*dxdz
+ end
+ S[i,j] += sij
+ end
+ end
+ return S
+end
+
+
+function vandermonde_matrix(z::AbstractVector, bs::Bspline2)
+ return [ Bspline.polynomial(z[i], j, bs) for i in 1:length(z), j in 1:bs.Npts+1 ]
+end
+
+
+# analytical integration formula from
+# Splines and PDEs: From Approximation Theory to Numerical Linear Algebra 2017, p. 13
+function integrate(f::AbstractVector, bs::Bspline2)
+ @toggled_assert length(f) == bs.Npts+1
+ sum = 0.0
+ @inbounds for i in bs.Npts+1
+ fi, xl, xr = f[i], bs.xs[i], bs.xs[i+3]
+ sum += fi * (xr-xl)/3
+ end
+ return sum
+end
+function integrate(f::Function, idxs::AbstractVector{<:Integer}, bs::Bspline2)
+ @toggled_assert length(idxs) == bs.Npts+1
+ sum = 0.0
+ @inbounds for i in bs.Npts+1
+ fi, xl, xr = f[idxs[i]], bs.xs[i], bs.xs[i+3]
+ sum += fi * (xr-xl)/3
+ end
+ return sum
+end
+
+
+end
diff --git a/src/dg1d.jl b/src/dg1d.jl
index 39ac3bfb..b9cbe9f4 100644
--- a/src/dg1d.jl
+++ b/src/dg1d.jl
@@ -51,7 +51,9 @@ include("output.jl")
include("lagrange.jl")
include("glgl.jl")
include("lgl.jl")
+include("lg.jl")
include("bernstein.jl")
+include("bspline.jl")
include("legendre.jl")
#######################################################################
diff --git a/src/dgelement.jl b/src/dgelement.jl
index dcb1914e..b5698cf4 100644
--- a/src/dgelement.jl
+++ b/src/dgelement.jl
@@ -1,3 +1,7 @@
+# TODO Make DGElement an abstract type and add the following subtypes
+# - NodalDGElement
+# - ModalBernsteinDGElement
+# - ModalBspline2DGElement
struct DGElement
N::Int64 # polynomial order
Npts::Int64 # number of quadrature points = N + 1
@@ -17,7 +21,7 @@ struct DGElement
S::Matrix{Float64} # stiffness matrix
M::Matrix{Float64} # mass matrix
invM::Matrix{Float64} # inverse mass matrix
- MDM::Matrix{Float64} # = (M^-1 D^T M) = M^-1 S^T
+ MDM::Matrix{Float64} # = (M^-1 D^T M) = M^-1 S^T (TODO: this identity only applies to nodal DG!!!)
Ml_lhs::Vector{Float64} # = M^-1 l(-1)
Ml_rhs::Vector{Float64} # = M^-1 l(+1)
@@ -25,13 +29,23 @@ struct DGElement
V::Matrix{Float64} # = P_(j-1)(z_i), Vandermonde matrix
invV::Matrix{Float64} # inv(V), inverse Vandermonde matrix
F::Matrix{Float64} # spectral filtering matrix
+ bs::Bspline.Bspline2 # only set for :modal_bspline
+ # TODO Replace N with Npts and remove N where ever possible.
function DGElement(N::Integer, quadr::Symbol, kind::Symbol;
# filter parameters
ηc::Float64=0.5, sh::Float64=16.0, α::Float64=36.0)
+ if quadr ∉ (:lgl, :lgl_lumped, :glgl)
+ error("Invalid quadrature rule specified for kind $kind: $quadr. Use one of lgl_lumped, lgl, glgl.")
+ end
+ if kind ∉ (:nodal_lagrange, :modal_bernstein, :modal_bspline2)
+ error("Invalid DG formulation specified: $kind. Use one of nodal_lagrange, modal_bernstein, modal_bspline2.")
+ end
+
Npts = N + 1
+ local bs
if kind === :nodal_lagrange
if quadr === :lgl
quadr_z, quadr_w = LGL.rule(Npts)
@@ -51,14 +65,12 @@ struct DGElement
M = GLGL.mass_matrix(Npts)
S = GLGL.stiffness_matrix(Npts)
invM = inv(M)
- else
- error("Invalid quadrature rule specified for kind $kind: $quadr. Use one of lgl_lumped, lgl, glgl.")
end
z = quadr_z
D = first_derivative_matrix(z)
V = vandermonde_matrix_legendre(z)
- l_lhs = zeros(Float64, N+1)
- l_rhs = zeros(Float64, N+1)
+ l_lhs = zeros(Float64, Npts)
+ l_rhs = zeros(Float64, Npts)
l_lhs[1] = 1.0
l_rhs[end] = 1.0
elseif kind === :modal_bernstein
@@ -66,25 +78,59 @@ struct DGElement
quadr_z, quadr_w = LGL.rule(Npts)
elseif quadr === :glgl
quadr_z, quadr_w = GLGL.rule(Npts)
- else
- error("Invalid quadrature rule specified for kind $kind: $quadr. Use one of lgl_lumped, lgl, glgl.")
end
purge_zeros!(quadr_z)
z = collect(range(-1.0,1.0,Npts))
purge_zeros!(z)
- D = [ Bernstein.derivative_bernstein_polynomial.(z[i], j-1, N) for i in 1:Npts, j in 1:Npts ]
+ D = [ Bernstein.derivative_bernstein_polynomial(z[i], j-1, N) for i in 1:Npts, j in 1:Npts ]
M = Bernstein.mass_matrix(Npts)
S = Bernstein.stiffness_matrix(Npts)
invM = inv(M)
V = Bernstein.vandermonde_matrix(z)
- l_lhs = zeros(Float64, N+1)
- l_rhs = zeros(Float64, N+1)
+ l_lhs = zeros(Float64, Npts)
+ l_rhs = zeros(Float64, Npts)
+ l_lhs[1] = 1.0
+ l_rhs[end] = 1.0
+ elseif kind === :modal_bspline2
+ if N < 2
+ error("Requested N+1=$(N+1) DOFs for DG formulation of kind $kind, but require at least 3.")
+ end
+ # User asked for Nth order polynomial sampled on Npts=N+1 points,
+ # but we can only provide N 2nd order Bsplines, which would give a total of Npts+1 DOFs.
+ # This is more than requested, so we reduce it by one.
+ # bs_z is now the grid across which the Bsplines are C^1 continuous.
+ if quadr === :lgl || quadr === :lgl_lumped
+ bs_z, _ = LGL.rule(N)
+ elseif quadr === :glgl
+ bs_z, _ = GLGL.rule(N)
+ end
+ purge_zeros!(bs_z)
+ bs = Bspline.Bspline2(bs_z)
+ # The dual grid of bs_z is the one on which we define the Npts quasi-nodal DOFs.
+ # These are the nodes on which we evaluate residuals etc.
+ z = vcat(bs_z[1], (bs_z[1:end-1].+bs_z[2:end])./2, bs_z[end])
+ purge_zeros!(z)
+ # Ignore given quadrature rule since we have to evaluate integrals on the subintervals,
+ # for which a 4th order accurate LG rule should suffice.
+ # We use a Gauss rule to avoid evaluating the Bspline on subcell interfaces and bdrys.
+ quadr_z, quadr_w = LG.rule(2)
+ D = [ Bspline.derivative_polynomial(z[i], j, bs) for i in 1:Npts, j in 1:Npts ]
+ M = Bspline.mass_matrix(bs)
+ S = Bspline.stiffness_matrix(bs)
+ invM = inv(M)
+ V = Bspline.vandermonde_matrix(z, bs)
+ l_lhs = zeros(Float64, Npts)
+ l_rhs = zeros(Float64, Npts)
l_lhs[1] = 1.0
l_rhs[end] = 1.0
- else
- error("Invalid DG formulation specified: $kind. Use on of nodal_lagrange, modal_bernstein.")
end
+ if kind !== :modal_bspline2
+ # just set up something, its unused for the other kinds
+ bs = Bspline.Bspline2(z)
+ end
+
+ # TODO Rename MDM to MST or so, because the identity MDM = M^{-1} S^{T} only applies for nodal DG
MDM = invM * transpose(S)
Ml_lhs = invM * l_lhs
Ml_rhs = invM * l_rhs
@@ -98,13 +144,15 @@ struct DGElement
purge_zeros!(invV)
purge_zeros!(Ml_lhs)
purge_zeros!(Ml_rhs)
+
+ # TODO Filter operators should likely go into separate struct
Λ = diagm([ spectral_filter((i-1)/N, (ηc,sh,α)) for i = 1:Npts ])
F = V * Λ * invV
return new(N, Npts, kind,
quadr, quadr_z, quadr_w,
z, D, S, M, invM, MDM, Ml_lhs, Ml_rhs,
- V, invV, F)
+ V, invV, F, bs)
end
end
@@ -131,19 +179,27 @@ end
@inline function integrate(u::AbstractVector, el::DGElement)
- if el.quadr === :lgl || el.quadr === :lgl_lumped
- return LGL.integrate(u, el.quadr_z, el.quadr_w)
- else # el.quadr === :glgl
- return GLGL.integrate(u, el.quadr_z, el.quadr_w)
+ if el.kind === :modal_bspline2
+ return Bspline.integrate(u, el.bs)
+ else
+ if el.quadr === :lgl || el.quadr === :lgl_lumped
+ return LGL.integrate(u, el.quadr_z, el.quadr_w)
+ else # el.quadr === :glgl
+ return GLGL.integrate(u, el.quadr_z, el.quadr_w)
+ end
end
end
@inline function integrate(fn::Function, idxs::AbstractVector{<:Integer}, el::DGElement)
- if el.quadr === :lgl || el.quadr === :lgl_lumped
- return LGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
- else # el.quadr === :glgl
- return GLGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
+ if el.kind === :modal_bspline2
+ return Bspline.integrate(fn, idxs, el.bs)
+ else
+ if el.quadr === :lgl || el.quadr === :lgl_lumped
+ return LGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
+ else # el.quadr === :glgl
+ return GLGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
+ end
end
end
diff --git a/src/legendre.jl b/src/legendre.jl
index fab3a9d7..5e6302ae 100644
--- a/src/legendre.jl
+++ b/src/legendre.jl
@@ -74,7 +74,7 @@ end
"""
- legendre_projection(f::Function, Npts::Integer)
+ projection(f::Function, Npts::Integer)
Compute the modal Legendre coefficients `f` using a LGL quadrature rule.
If a smooth `f` is exactly representable as an `N`th polynomial,
diff --git a/src/lg.jl b/src/lg.jl
new file mode 100644
index 00000000..e0ca1c9a
--- /dev/null
+++ b/src/lg.jl
@@ -0,0 +1,86 @@
+module LG
+
+
+using Jacobi, LinearAlgebra, ToggleableAsserts
+
+
+"""
+ rule(Npts::Integer)
+
+Computes quadrature nodes and weights for the Legendre-Gauss quadrature rule.
+"""
+function rule(Npts::Integer)
+ α = 0.0; β = 0.0
+ x = zgj(Npts, α, β)
+ w = wgj(x, α, β)
+ return x, w
+end
+
+
+"""
+ integrate(f, x, w)
+
+Integrate `f` sampled on LG quadrature nodes `x` with quadrature weights `w`.
+Use `rule` to obtain `x, w`.
+"""
+function integrate(f::AbstractVector, x::AbstractVector, w::AbstractVector)
+ @toggled_assert length(w) == length(f) == length(x)
+ return dot(w, f)
+end
+
+
+"""
+ integrate(f::Function, idxs::AbstractVector{<:Integer}, x, w)
+
+Integrate `f` sampled on LG quadrature nodes `x` with quadrature weights `w`.
+`f(idxs[i])` shall return the value of the integrand evaluated at `x[i]`.
+Use `rule` to obtain `x, w`.
+"""
+function integrate(f::Function, idxs::AbstractVector{<:Integer}, x::AbstractVector, w::AbstractVector)
+ @toggled_assert length(w) == length(idxs)
+ sum = 0.0
+ for i in 1:length(w)
+ sum += w[i] * f(idxs[i])
+ end
+ return sum
+end
+
+
+"""
+ integrate(f::Function, x, w)
+
+Integrate `f` sampled on LG quadrature nodes `x` with quadrature weights `w`.
+`f(x[i])` shall return the value of the integrand evaluated at `x[i]`.
+Use `rule` to obtain `x, w`.
+
+This is similar to `integrate(f, idxs, x, w)`, but is better suited if the integrand is
+given in closed form.
+"""
+function integrate(f::Function, x::AbstractVector, w::AbstractVector)
+ @toggled_assert length(w) == length(x)
+ sum = 0.0
+ for i in 1:length(w)
+ sum += w[i] * f(x[i])
+ end
+ return sum
+end
+
+
+"""
+ inner_product(f1, f2, x, w)
+
+Computed the inner product between `f1` and `f2` sampled on LG quadrature nodes `x`
+with quadrature weights `w`.
+Use `rule` to obtain `x, w`.
+"""
+function inner_product(f1::AbstractVector, f2::AbstractVector, x::AbstractVector, w::AbstractVector)
+ @toggled_assert length(w) == length(f1) == length(f2) == length(x)
+ sum = 0.0
+ for i in 1:length(w)
+ sum += w[i] * f1[i] * f2[i]
+ end
+ return sum
+end
+
+
+end
diff --git a/src/lgl.jl b/src/lgl.jl
index 68a1a556..44670809 100644
--- a/src/lgl.jl
+++ b/src/lgl.jl
@@ -8,7 +8,7 @@ import dg1d: lagrange_poly, first_derivative_matrix, purge_zeros!
"""
rule(Npts::Integer)
-Computes quadrature nodes and weights for the Legendre-Gauss-Lobatto-Legendre quadrature rule.
+Computes quadrature nodes and weights for the Legendre-Gauss-Lobatto quadrature rule.
"""
function rule(Npts::Integer)
α = 0.0; β = 0.0
diff --git a/src/parameters.jl b/src/parameters.jl
index 31b2078f..2160f97c 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -489,17 +489,17 @@ end
# TODO Rename basis to quadrature
"""
Nodal basis type used for the polynomial approximation. Available options
- - `lgl_lumped` ... mass-lumped Legendre-Gauss-Lobatto grid
- - `lgl` ... Legendre-Gauss-Lobatto grid
- - `glgl` ... Generalized Legendre-Gauss-Lobatto grid
+ - `lgl_lumped`: mass-lumped Legendre-Gauss-Lobatto grid
+ - `lgl`: Legendre-Gauss-Lobatto grid
+ - `glgl`: Generalized Legendre-Gauss-Lobatto grid
"""
basis = "lgl_lumped"
@check basis in [ "lgl", "glgl", "lgl_lumped" ]
"""
Cellwise approximation scheme of solution. Available options
- - `DG` ... Discontinuous Galerkin
- - `FV` ... Finite Volume (central scheme)
+ - `DG`: Discontinuous Galerkin
+ - `FV`: Finite Volume (central scheme)
"""
scheme = "DG"
@check scheme in [ "DG", "FV" ]
@@ -507,14 +507,30 @@ end
# TODO Rename this to basis
"""
DG scheme used for the evolution. Available options
- - `nodal_lagrange` ... Nodal approximation using Lagrange polynomials collocated with
+ - `nodal_lagrange`: Nodal approximation using Lagrange polynomials collocated with
the quadrature nodes provided by `basis`. This is the standard method discussed
in Hesthaven, Warburton, 2007.
- - `modal_bernstein` ... A modal approximation using Bernstein polynomials.
+ See [1].
+ - `modal_bernstein`: A modal approximation using Bernstein polynomials.
The coefficients correspond to the nodal values sampled on an equidistant grid,
but the Bernstein polynomials are not interpolating.
- Hence, this Ansatz is also refered to as quasi-nodal.
- See Beisiegel, Behrens, Environ Earth Sci (2015) 74:7275–7284.
+ We refer to these coefficients as being quasi-nodal.
+ See [2,3,4].
+ - `modal_bspline2`: A modal approximation using 2nd order Bspline polynomials.
+ The coefficients correspond to the nodal values sampled on a grid that is dual (including the bdrys)
+ to the grid of a `nodal_lagrange` Ansatz of order `n-1`.
+ I.e. a `modal_bspline2` element has `n+1` DOF.
+ We adapt a quasi-nodal interpretation for the coefficients.
+ The bdry values correspond to the same bdry values as in a nodal Ansatz,
+ the cell interior values correspond to the cell-averages (determined by linear interpolation)
+ of a nodal Ansatz.
+ See [3,4].
+ Note: Spectral filtering with the `filter_*` parameters is not implemented correctly.
+
+ [1] Hesthaven, Warburton, Nodal Discontinuous Galerkin Methods (2008).
+ [2] Beisiegel, Behrens, Environ Earth Sci (2015) 74:7275–7284.
+ [3] Piegl, Tiller, The NURBS book (1997).
+ [4] Kunoth et al., Splines and PDEs: From Approximation Theory to Numerical Linear Algebra (2018).
"""
dg_kind = "nodal_lagrange"
@check scheme in [ "nodal_lagrange", "modal_bernstein" ]
@@ -556,13 +572,13 @@ end
"""
Time stepping method. Available options
- - `midpt` ... 2nd order explicit Runge Kutta method -- midpoint formula
- - `rk4` ... 4th order explicit Runge Kutta method
- - `rk4_twothirds` ... 4th order explicit Runge Kutta method -- two thirds version
- - `rk4_ralston` ... 4th order explicit Runge Kutta method -- Ralston's method, designed for minimal truncation error
- - `ssprk3` ... 3rd order strong stability preserving Runge Kutta method
- - `lserk4` ... 4th order low storage explicit Runge Kutta method
- - `rkf10` ... 10th order Runge Kutta Feagin method
+ - `midpt`: 2nd order explicit Runge Kutta method -- midpoint formula
+ - `rk4`: 4th order explicit Runge Kutta method
+ - `rk4_twothirds`: 4th order explicit Runge Kutta method -- two thirds version
+ - `rk4_ralston`: 4th order explicit Runge Kutta method -- Ralston's method, designed for minimal truncation error
+ - `ssprk3`: 3rd order strong stability preserving Runge Kutta method
+ - `lserk4`: 4th order low storage explicit Runge Kutta method
+ - `rkf10`: 10th order Runge Kutta Feagin method
"""
method = "lserk4"
@check method in [ "midpt", "lserk4", "rk4", "rk4_twothirds", "rk_ralston",
@@ -655,25 +671,25 @@ end
"""
Data format for 0d output. Available options:
- - h5 ... raw hdf5 output, all fields in one file
+ - `h5`: raw hdf5 output, all fields in one file
"""
format0d = "h5"
@check format0d in [ "h5" ]
"""
Data format for 1d output. Available options:
- - h5 ... raw hdf5 output, all fields in one file
+ - `h5`: raw hdf5 output, all fields in one file
"""
format1d = "h5"
@check format1d in [ "h5" ]
"""
Data format for 2d output. Available options:
- - h5 ... raw hdf5 output, all fields in one file
- - xdmf+bin ... xdmf data description and binary data dump
- - xdmf+h5 ... xdmf data description and hdf5 data dump
- - xdmf ... same as xdmf+h5
- - vtk ... vtk xml format utilizing unstructured grids (but without subcell resolution)
+ - `h5`: raw hdf5 output, all fields in one file
+ - `xdmf+bin`: xdmf data description and binary data dump
+ - `xdmf+h5`: xdmf data description and hdf5 data dump
+ - `xdmf`: same as xdmf+h5
+ - `vtk`: vtk xml format utilizing unstructured grids (but without subcell resolution)
"""
format2d = "vtk"
@check format2d in [ "h5", "xdmf", "xdmf+bin", "xdmf+h5", "vtk" ]
diff --git a/test/UnitTests/src/test_dgelement.jl b/test/UnitTests/src/test_dgelement.jl
index 028a9d33..210ac199 100644
--- a/test/UnitTests/src/test_dgelement.jl
+++ b/test/UnitTests/src/test_dgelement.jl
@@ -129,27 +129,6 @@ end
end
-function lgl_linear_dependency_matrix_exact(Npts::Integer)
- Npts2 = Npts+2
- x2, w2 = dg1d.LGL.rule(Npts2)
- L = zeros(N,N)
- for i in 1:N
- Pi = map(x->dg1d.legendre(x, i-1), x2)
- for j in i:N
- Pj = map(x->dg1d.legendre(x, j-1), x2)
- lij = dg1d.LGL.integrate(Pi, Pj, x2, w2)
- L[i,j] = lij
- if i != j
- L[j,i] = L[i,j]
- end
- end
- end
- dg1d.purge_zeros!(L)
-
- return L
-end
-
-
@testset "Test orthogonality relation" begin
Npts = 6
@@ -194,10 +173,14 @@ end
int = dg1d.GLGL.integrate(p, x, w)
@test int ≈ 140/3
+ x, w = dg1d.LG.rule(Npts)
+ int = dg1d.LG.integrate(p, x, w)
+ @test int ≈ 140/3
+
end
-@testset "verify some DG operator identities" begin
+@testset "verify summation-by-parts property of DG operators" begin
N = 10
@@ -218,12 +201,32 @@ end
TST = T*S_legendre*transpose(T) |> dg1d.purge_zeros
@test all(el.S.≈TST)
+ el = dg1d.DGElement(N, :lgl, :modal_bspline2)
+ bs = el.bs
+ M = el.M
+ S = el.S
+ ST = transpose(S)
+ # verify the sum of surface terms
+ BB = dg1d.purge_zeros(S+ST)
+ exact_BB = similar(M); fill!(exact_BB, 0.0)
+ o = 2
+ for i in 1:bs.Npts+1, j in 1:bs.Npts+1
+ for k in 1:bs.Npts-1
+ xl, xr = bs.xs[o+k], bs.xs[o+k+1]
+ Nil, Nir = dg1d.Bspline.polynomial(xl, i, bs), dg1d.Bspline.polynomial(xr, i, bs)
+ Njl, Njr = dg1d.Bspline.polynomial(xl, j, bs), dg1d.Bspline.polynomial(xr, j, bs)
+ exact_BB[i,j] += Nir*Njr - Nil*Njl
+ end
+ end
+ dg1d.purge_zeros!(exact_BB)
+ @test all(exact_BB .≈ BB)
+
end
-@testset "Legendre expansions and Bernstein polynomials" begin
+@testset "Bernstein consistency tests" begin
- # approximate the Bernstein basis through Legendre polynomials
+ ## approximate the Bernstein basis through Legendre polynomials
N = 10
Npts = N+1
Cs = zeros(Float64, Npts, Npts)
@@ -234,13 +237,36 @@ end
end
end
- xx = collect(range(-1.0,1.0,25))
+ xx = collect(range(-1.0,1.0,length=25))
for xi in xx
Bs = [ dg1d.Bernstein.bernstein_polynomial(xi, i-1, N) for i in 1:Npts ] |> dg1d.purge_zeros
approx_Bs = [ dg1d.Legendre.interpolate(xi, Cs[i,:]) for i in 1:Npts ] |> dg1d.purge_zeros
@test all(isapprox.(Bs, approx_Bs; atol=1e-10, rtol=1e-8))
end
+ ## partition of unity
+ Bs = [ dg1d.Bernstein.bernstein_polynomial(xx[i], j-1, N) for i in 1:length(xx), j in 1:Npts ]
+ POU = sum(Bs, dims=2)
+ @test all(POU .≈ ones(Float64, length(POU)))
+
+end
+
+
+@testset "Bspline consistency tests" begin
+
+ z, _ = dg1d.LGL.rule(5)
+ bs = dg1d.Bspline.Bspline2(z)
+ xx = collect(range(-1.0,1.0,length=25))
+ Bs = [ dg1d.Bspline.polynomial(xx[i], j, bs) for i in 1:length(xx), j in 1:bs.Npts+1 ]
+
+ ## partition of unity
+ POU = sum(Bs, dims=2)
+ @test all(POU .≈ ones(Float64, length(POU)))
+ # alternative test: evaluate a Bspline polynomial with all coefficients equal to one
+ coeffs = ones(Float64, bs.Npts+1)
+ POU = dg1d.Bspline.interpolate.(xx, Ref(coeffs), Ref(bs))
+ @test all(POU .≈ ones(Float64, length(POU)))
+
end
diff --git a/test/UnitTests/src/test_utils.jl b/test/UnitTests/src/test_utils.jl
index f60b9490..c4512da7 100644
--- a/test/UnitTests/src/test_utils.jl
+++ b/test/UnitTests/src/test_utils.jl
@@ -1,6 +1,5 @@
import dg1d: make_outputdir, flatten, dirname_path, query,
- splitallext, trim_keys, match_keys, prepend_keys,
- concretetypes, match_properties, indent
+ splitallext, trim_keys, match_keys, prepend_keys, indent
@testset "Utilities" begin
--
GitLab