diff --git a/src/glgl.jl b/src/glgl.jl
index 47d47cd5e93664ea5c664d66877502772089ecc7..54207a060c001df34c369f22da71452ddc3f9b9e 100644
--- a/src/glgl.jl
+++ b/src/glgl.jl
@@ -1,22 +1,23 @@
module GLGL
-using LinearAlgebra, Jacobi, ToggleableAsserts
-import dg1d: purge_zeros!, lagrange_poly
+using Jacobi, LinearAlgebra, ToggleableAsserts
+import dg1d: purge_zeros!, lagrange_poly, first_derivative_matrix, second_derivative_matrix
"""
-Computes quadrature nodes, weights and derivative weights
-for the Generalized-Gauss-Lobatto quadrature rule.
-n_pts ... number of quadrature nodes
-Return values
-x ... quadrature nodes, vector
-w ... quadrature weights, vector
-v ... derivative weight, number
+ rule(Npts::Integer)
+
+Computes quadrature nodes and weights for the Generalized-Legendre-Gauss-Lobatto quadrature rule.
"""
-function rule(n_pts::Integer)
+function rule(Npts::Integer)
+ x, w, v, D = inner_product_rule(Npts)
+ @. w += v * (D[1,:] - D[end,:])
+ return x, w
+end
+function inner_product_rule(Npts::Integer)
+ N = Npts - 2
- N = n_pts - 2
gb_x, gb_w = [], []
if N > 0
# Integration nodes and weights for Gauss-Gegenbauer quadrature
@@ -24,229 +25,127 @@ function rule(n_pts::Integer)
gb_x = zgj(N, α, β)
gb_w = wgj(gb_x, α, β)
end
-
- x = zeros(n_pts)
- w = zeros(n_pts)
+ x = zeros(Npts)
+ w = zeros(Npts)
x[2:end-1] .= gb_x
x[1] = -1.0
x[end] = 1.0
-
w[2:end-1] .= @. gb_w / (1.0 - gb_x^2)^2
-
w[1] = w[end] = 8.0 / 3.0 * (2.0 * N^2 + 10.0 * N + 9.0) / ((N+1.0) * (N+2.0) * (N+3.0) * (N+4.0))
v = 8.0 / ((N+1.0) * (N+2.0) * (N+3.0) * (N+4.0))
- return x, w, v
-end # function
-
-
-
-"""
- integrate(f, x, w, v, D)
-
-Integrate `f` sampled on GLGL quadrature nodes `x` with quadrature weights `w`,
-boundary weight `v` and derivative matrix `D` based on `x`.
-Use `rule` to obtain `x, w, v` and use `first_derivative_matrix` to obtain `D`.
-"""
-function integrate(f::AbstractVector, x::AbstractVector, w::AbstractVector, v::Real, D::AbstractVector)
-
- N = length(w)
- @toggled_assert (N,N) == size(D)
- @toggled_assert N == length(f)
-
- sum = dot(f, w)
- sum += v * dot(view(D,1,:), f)
- sum -= v * dot(view(D,N,:), f)
+ D = first_derivative_matrix(x)
- return sum
+ return x, w, v, D
end
"""
- integrate(f1, f2, w, v, D)
+ integrate(f, x, w)
-Similar to `integrate(f, w)`, but computes the inner product of `f1*f2` where both
-`f1, f2` must be sampled on the GLGL quadrature nodes `x` associated with the quadrature
-weights `w`, boundary weight `v` and derivative matrix `D` based on `x`.
-Use `rule` to obtain `x, w, v` and use `first_derivative_matrix` to obtain `D`.
+Integrate `f` sampled on GLGL quadrature nodes `x` with quadrature weights `w`.
+Use `rule` to obtain `x, w`.
"""
-function integrate(f1::AbstractVector, f2::AbstractVector, x::AbstractVector, w::AbstractVector,
- v::Real, D::AbstractMatrix)
-
- N = length(w)
- @toggled_assert (N,N) == size(D)
- @toggled_assert N == length(f1)
- @toggled_assert N == length(f2)
-
- sum = dot(f1.*f2, w)
- vD1, vDN = view(D,1,:), view(D,N,:)
- sum += v * (f1[1] * dot(vD1, f2) + f2[1] * dot(vD1, f1))
- sum -= v * (f1[end] * dot(vDN, f2) + f2[end] * dot(vDN, f1))
-
- return sum
+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 GLGL quadrature nodes `x` with quadrature weights `w`,
-boundary weight `v` and derivative matrix `D` based on `x`.
+Integrate `f` sampled on GLGL 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, v` and use `first_derivative_matrix` to obtain `D`.
+Use `rule` to obtain `x, w`.
"""
-function integrate(fn::Function, idxs::AbstractVector{<:Integer}, x::AbstractVector, w::AbstractVector,
- v::Real, D::AbstractMatrix)
-
- N = length(w)
- @toggled_assert (N,N) == size(D)
- @toggled_assert N == length(idxs)
-
+function integrate(fn::Function, idxs::AbstractVector{<:Integer}, x::AbstractVector)
+ @toggled_assert length(w) == length(idxs)
sum = 0.0
- for i in 1:N
- fi = fn(idxs[i])
- sum += (w[i] + v * (D[1,i] - D[N,i])) * fi
+ for i in 1:length(w)
+ sum += w[i] * f(idxs[i])
end
-
return sum
end
"""
- integrate(f::Function, x, w, v, D)
+ integrate(f::Function, x, w)
-Integrate `f` sampled on GLGL quadrature nodes `x` with quadrature weights `w`,
-boundary weight `v` and derivative matrix `D` based on `x`.
+Integrate `f` sampled on GLGL 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, v` and use `first_derivative_matrix` to obtain `D`.
+Use `rule` to obtain `x, w`.
-This is similar to `integrate(f, idxs, x, w, v, D)`, but is better suited if the integrand is
+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, v::Real, D::AbstractMatrix)
- @toggled_assert length(w) == length(x) == size(D, 1) == size(D, 2)
+function integrate(f::Function, x::AbstractVector, w::AbstractVector)
+ @toggled_assert length(w) == length(x)
sum = 0.0
N = length(w)
for i in 1:N
fi = f(x[i])
- sum += (w[i] + v * (D[1,i] - D[N,i])) * fi
+ sum += w[i] * fi
end
return sum
end
"""
-M_ij = (li, lj)
-x ... GLGL nodes
-w ... GLGL weights
-v ... GLGL derivative weight
-D ... Derivative matrix based on x.
+ inner_product(f1, f2, x, w)
+
+Computed the inner product between `f1` and `f2` sampled on GLGL quadrature nodes `x`
+with quadrature weights `w`, boundary weights `v` and first derivative matrix `D`
+Use `inner_product_rule` to obtain `x, w, v, D`.
"""
-function mass_matrix(x::AbstractVector, w::AbstractVector, v::Real, D::AbstractMatrix)
+function inner_product(f1::AbstractVector, f2::AbstractVector, x::AbstractVector, w::AbstractVector,
+ v::Real, D::AbstractMatrix)
+ N = length(w)
+ @toggled_assert (N,N) == size(D)
+ @toggled_assert N == length(f1)
+ @toggled_assert N == length(f2)
+ sum = dot(f1.*f2, w)
+ vD1, vDN = view(D,1,:), view(D,N,:)
+ sum += v * (f1[1] * dot(vD1, f2) + f2[1] * dot(vD1, f1))
+ sum -= v * (f1[end] * dot(vDN, f2) + f2[end] * dot(vDN, f1))
+ return sum
+end
- N = length(x)
- @assert N == length(w)
- @assert (N,N) == size(D)
+"""
+ mass_matrix(Npts::Integer)
+
+Mass matrix for Lagrange polynomials based on GLGL quadrature nodes.
+M_ij = (li, lj)
+"""
+function mass_matrix(Npts::Integer)
+ x, w, v, D = inner_product_rule(Npts)
M = diagm(w)
M[1,:] += v * D[1,:]
M[:,1] += v * D[1,:]
M[end,:] -= v * D[end,:]
M[:,end] -= v * D[end,:]
-
return M
-end # function
-
+end
"""
+ stiffness_matrix(Npts::Integer)
+
+Stiffness matrix for Lagrange polynomials based on GLGL quadrature nodes.
S_ij = (li, dlj/dx)
-Exactly representable, because product is of order 2N - 1 which can be integrated
-exactly using LGL quadrature.
-x ... GLGL nodes
-w ... GLGL weights
-v ... GLGL derivative weight
-D ... First derivative matrix based on x.
-D2 ... Second derivative matrix based on x.
"""
-function stiffness_matrix(x::AbstractVector, w::AbstractVector, v::Real,
- D::AbstractMatrix, D2::AbstractMatrix)
-
- N = length(x)
- @assert N == length(w)
- @assert (N,N) == size(D)
- @assert (N,N) == size(D2)
-
+function stiffness_matrix(Npts::Integer)
+ x, w, v, D = inner_product_rule(Npts)
+ D2 = second_derivative_matrix(x)
S = w .* D
S[1,:] += v * D2[1,:]
S[end,:] -= v * D2[end,:]
S .+= v * D[1,:] .* transpose(D[1,:])
S .-= v * D[end,:] .* transpose(D[end,:])
-
return S
end
-# TODO Add versions of mass_matrix, stiffness_matrix that take
-# a single argument n_pts.
-
-
-"""
-L_ij = (P_i, P_j)
-Analogon to mass matrix of Lagrange polynomials, expresses the linear dependency
-of the numerical representation of the Legendre basis. Ideally, it is a
-diagonal matrix (which should be ensured by the GLGL quadrature).
-x ... GLGL nodes
-w ... GLGL weights
-v ... GLGL derivative weight
-D ... First derivative matrix based on x.
-"""
-function linear_dependency_matrix(x::AbstractVector, w::AbstractVector, v::Real, D::AbstractMatrix)
- N = length(x)
- L = zeros(N,N)
- for i in 1:N
- Pi = map(x->legendre(x, i-1), x)
- for j in i:N
- Pj = map(x->legendre(x, j-1), x)
- lij = integrate(Pi, Pj, x, w, v, D)
- L[i,j] = lij
- if i != j
- L[j,i] = L[i,j]
- end
- end
- end
- purge_zeros!(L)
-
- return L
end
-
-
-"""
-E_ij = (l_i, P_j)
-Commonly used to project either Lagrange polynomials onto Legendres or
-reverse.
-x ... GLGL nodes
-w ... GLGL weights
-v ... GLGL derivative weight
-D ... First derivative matrix based on x.
-"""
-function expansion_matrix(x::AbstractVector, w::AbstractVector, v::Real, D::AbstractMatrix)
- N = length(x)
- T = zeros(N,N)
- for i in 1:N
- li = map(xi->lagrange_poly(i, xi, x), x)
- for j in 1:N
- Pj = map(x->legendre(x, j-1), x)
- tij = integrate(li, Pj, x, w, v, D)
- T[i,j] = tij
- end
- end
- purge_zeros!(T)
-
- return T
-end
-
-
-
-end # module
diff --git a/src/lgl.jl b/src/lgl.jl
index 31309a4173066a78e23a5f26f5b05577a94ab6ff..6a7c97e859347df251d63d181ce136bb110301eb 100644
--- a/src/lgl.jl
+++ b/src/lgl.jl
@@ -5,21 +5,19 @@ using Jacobi, LinearAlgebra, ToggleableAsserts
import dg1d: lagrange_poly, first_derivative_matrix, purge_zeros!
-
"""
-Quadrature nodes and weights for Gauss-Legendre rule.
-n_pts ... number of quadrature points
-Exact up to order 2*n_pts - 1.
+ rule(Npts::Integer)
+
+Computes quadrature nodes and weights for the Legendre-Gauss-Lobatto-Legendre quadrature rule.
"""
-function rule(n_pts::Integer)
+function rule(Npts::Integer)
α = 0.0; β = 0.0
- x = zglj(n_pts, α, β)
+ x = zglj(Npts, α, β)
w = wglj(x, α, β)
return x, w
end
-
"""
integrate(f, x, w)
@@ -27,30 +25,11 @@ Integrate `f` sampled on LGL 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)
+ @toggled_assert length(w) == length(f) == length(x)
return dot(w, f)
end
-"""
- integrate(f1, f2, x, w)
-
-Similar to `integrate(f, w)`, but computes the inner product of `f1*f2` where both
-`f1, f2` must be sampled on the LGL quadrature nodes `x` associated with the quadrature
-weights `w`.
-Use `rule` to obtain `x, w`.
-"""
-function integrate(f1::AbstractVector, f2::AbstractVector, x::AbstractVector, w::AbstractVector)
- @toggled_assert length(w) == length(f1)
- @toggled_assert length(w) == length(f2)
- I = 0.0
- for i = 1:length(w)
- I += w[i] * f1[i] * f2[i]
- end
- return I
-end
-
-
"""
integrate(f::Function, idxs::AbstractVector{<:Integer}, x, w)
@@ -88,6 +67,23 @@ function integrate(f::Function, x::AbstractVector, w::AbstractVector)
end
+"""
+ inner_product(f1, f2, x, w)
+
+Computed the inner product between `f1` and `f2` sampled on LGL 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
+
+
"""
Inverse of exact full mass matrix, computed using Sherman-Morrison formula,
cf. PhD thesis Marcus Bugner, Notes spectral elements Bernd Bruegmann,
@@ -105,176 +101,64 @@ function inverse_mass_matrix_exact(n_pts::Integer)
end
-
"""
n_pts ... number of qudrature nodes = number of lagrange polynomials
"""
-function inverse_mass_matrix_lumped(n_pts::Integer)
- _, w = rule(n_pts)
+function inverse_mass_matrix_lumped(Npts::Integer)
+ _, w = rule(Npts)
invM = diagm(1.0 ./ w)
return invM
end
-
"""
M_ij = (li, lj)
n_pts ... number of qudrature nodes = number of lagrange polynomials
"""
-function mass_matrix_exact(n_pts::Integer)
- z, w = rule(n_pts)
- M = zeros(n_pts, n_pts)
- aux_z, aux_w = rule(n_pts+2)
+function mass_matrix_exact(Npts::Integer)
+ z, w = rule(Npts)
+ M = zeros(Npts, Npts)
+ aux_z, aux_w = rule(Npts+2)
L(i, x) = lagrange_poly(i, x, z)
- M = zeros(n_pts, n_pts)
- for i=1:n_pts, j=1:n_pts
+ M = zeros(Npts, Npts)
+ for i=1:Npts, j=i:Npts
li = L.(i, aux_z)
- for j = i:n_pts
+ for j = i:Npts
lj = L.(j, aux_z)
- M[i,j] = integrate(li, lj, aux_w)
- if (i != j)
+ M[i,j] = inner_product(li, lj, aux_z, aux_w)
+ if j > i
M[j,i] = M[i,j]
end
end
end
-
return M
end
-
"""
M_ij = (li, lj)
Approximated mass matrix due to under integration of LGL quadrature.
n_pts ... number of qudrature nodes = number of lagrange polynomials
"""
-function mass_matrix_lumped(n_pts::Integer)
- z, w = rule(n_pts)
+function mass_matrix_lumped(Npts::Integer)
+ z, w = rule(Npts)
M = diagm(w)
return M
end
-
"""
S_ij = (li, dlj/dx)
Exactly representable, because product is of order 2N - 1 which can be integrated
exactly using LGL quadrature.
n_pts ... number of quadrature nodes = number of lagrange polynomials
"""
-function stiffness_matrix(n_pts::Integer)
- z, w = rule(n_pts)
+function stiffness_matrix(Npts::Integer)
+ z, w = rule(Npts)
D = first_derivative_matrix(z)
S = w .* D
return S
end
-
-"""
-L_ij = (P_i, P_j)
-Analogon to mass matrix of Lagrange polynomials, expresses the linear dependency
-of the numerical representation of the Legendre basis.
-x ... GLGL nodes
-w ... GLGL weights
-"""
-function linear_dependency_matrix_lumped(x::AbstractVector, w::AbstractVector)
- N = length(x)
- L = zeros(N,N)
- for i in 1:N
- Pi = map(x->legendre(x, i-1), x)
- for j in i:N
- Pj = map(x->legendre(x, j-1), x)
- lij = integrate(Pi, Pj, x, w)
- L[i,j] = lij
- if i != j
- L[j,i] = L[i,j]
- end
- end
- end
- purge_zeros!(L)
-
- return L
-end
-
-
-"""
-E_ij = (l_i, P_j)
-Commonly used to project either Lagrange polynomials onto Legendres or
-reverse.
-x ... GLGL nodes
-w ... GLGL weights
-"""
-function expansion_matrix_lumped(x::AbstractVector, w::AbstractVector)
- N = length(x)
- T = zeros(N,N)
- for i in 1:N
- li = map(xi->lagrange_poly(i, xi, x), x)
- for j in 1:N
- Pj = map(x->legendre(x, j-1), x)
- tij = integrate(li, Pj, x, w)
- T[i,j] = tij
- end
- end
- purge_zeros!(T)
-
- return T
-end
-
-
-"""
-L_ij = (P_i, P_j)
-Analogon to mass matrix of Lagrange polynomials, expresses the linear dependency
-of the numerical representation of the Legendre basis.
-x ... GLGL nodes
-w ... GLGL weights
-"""
-function linear_dependency_matrix_exact(x::AbstractVector, w::AbstractVector)
- N = length(x)
- N2 = N+2
- x2, w2 = rule(N2)
- L = zeros(N,N)
- for i in 1:N
- Pi = map(x->legendre(x, i-1), x2)
- for j in i:N
- Pj = map(x->legendre(x, j-1), x2)
- lij = integrate(Pi, Pj, x2, w2)
- L[i,j] = lij
- if i != j
- L[j,i] = L[i,j]
- end
- end
- end
- purge_zeros!(L)
-
- return L
-end
-
-
-"""
-E_ij = (l_i, P_j)
-Commonly used to project either Lagrange polynomials onto Legendres or
-reverse.
-x ... GLGL nodes
-w ... GLGL weights
-"""
-function expansion_matrix_exact(x::AbstractVector, w::AbstractVector)
- N = length(x)
- N2 = N + 2
- x2, w2 = rule(N2)
- T = zeros(N,N)
- for i in 1:N
- li = map(xi->lagrange_poly(i, xi, x), x2)
- for j in 1:N
- Pj = map(x->legendre(x, j-1), x2)
- tij = integrate(li, Pj, x2, w2)
- T[i,j] = tij
- end
- end
- purge_zeros!(T)
-
- return T
end
-
-
-end # module LGL
diff --git a/src/spectralelement.jl b/src/spectralelement.jl
index dcf0b4186767f86b92bce7675df695bf51f1a907..81b56b75597663f8e3b9e77288cba8baa86d59b3 100644
--- a/src/spectralelement.jl
+++ b/src/spectralelement.jl
@@ -3,10 +3,9 @@ Base.@kwdef struct SpectralElement
Npts::Int64 # number of quadrature points = N + 1
z::Vector{Float64} # quadrature points
w::Vector{Float64} # quadrature weights
- v::Float64 # derivative weight (only needed for GLGL quadrature rule)
- bw::Vector{Float64} # barycentric weights
D::Matrix{Float64} # derivative matrix
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
Ml_lhs::Vector{Float64} # = M^-1 l(-1)
@@ -16,7 +15,7 @@ Base.@kwdef struct SpectralElement
quadr::Symbol # :LGL, :LGL_mass_lumped, :GLGL
F::Matrix{Float64} # spectral filtering matrix
- function SpectralElement(N::Int, quadr::Symbol;
+ function SpectralElement(N::Integer, quadr::Symbol;
# filter parameters
ηc::Float64=0.5, sh::Float64=16.0, α::Float64=36.0)
Npts = N + 1
@@ -25,35 +24,26 @@ Base.@kwdef struct SpectralElement
z, w = LGL.rule(Npts)
v = -1.0 # unused here
purge_zeros!(z) # ensures that the centered node for even N basis sits exactly at x = 0.
- bw = barycentric_weights(z)
D = first_derivative_matrix(z)
S = LGL.stiffness_matrix(Npts)
+ M = LGL.mass_matrix_exact(Npts)
invM = LGL.inverse_mass_matrix_exact(Npts)
- # TODO Figure out how much the aliasing error influences the accuracy of the inner product.
- # Could be totally dealised using the Orszag 2/3rds rule.
- L = LGL.linear_dependency_matrix_exact(z, w)
- E = LGL.expansion_matrix_exact(z, w)
elseif quadr === :LGL_lumped
z, w = LGL.rule(Npts)
v = -1.0 # unused here
purge_zeros!(z) # ensures that the centered node for even N basis sits exactly at x = 0.
- bw = barycentric_weights(z)
D = first_derivative_matrix(z)
S = LGL.stiffness_matrix(Npts)
+ M = LGL.mass_matrix_lumped(Npts)
invM = LGL.inverse_mass_matrix_lumped(Npts)
- L = LGL.linear_dependency_matrix_lumped(z, w)
- E = LGL.expansion_matrix_lumped(z, w)
elseif quadr === :GLGL
- z, w, v = GLGL.rule(Npts)
+ z, w = GLGL.rule(Npts)
purge_zeros!(z) # ensures that the centered node for even N basis sits exactly at x = 0.
- bw = barycentric_weights(z)
D = first_derivative_matrix(z)
D2 = second_derivative_matrix(z)
- M = GLGL.mass_matrix(z, w, v, D)
- S = GLGL.stiffness_matrix(z, w, v, D, D2)
+ M = GLGL.mass_matrix(Npts)
+ S = GLGL.stiffness_matrix(Npts)
invM = inv(M)
- L = GLGL.linear_dependency_matrix(z, w, v, D)
- E = GLGL.expansion_matrix(z, w, v, D)
else
error("Invalid quadrature rule specified! Use one of :LGL_lumped, :LGL, :GLGL.")
end
@@ -72,14 +62,14 @@ Base.@kwdef struct SpectralElement
Ml_lhs = invM * l_lhs
Ml_rhs = invM * l_rhs
- return new(N, Npts, z, w, v, bw, D, S, invM, MDM, Ml_lhs, Ml_rhs,
+ return new(N, Npts, z, w, D, S, M, invM, MDM, Ml_lhs, Ml_rhs,
V, invV, quadr, F)
end
end
# (5.16) from Hesthaven, Warburton 2007
-function spectral_filter(η, prms)
+function spectral_filter(η::Real, prms)
ηc, sh, α = prms
return if 0 ≤ η ≤ ηc
1
@@ -89,24 +79,23 @@ function spectral_filter(η, prms)
end
-@inline function inner_product(u1, u2, el::SpectralElement)
- @unpack quadr = el
- @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
- if quadr === :LGL || quadr === :LGL_lumped
- return LGL.integrate(u1, u2, el.z, el.w)
- else # quadr === :GLGL
- return GLGL.integrate(u1, u2, el.z, el.w, el.v, el.D)
+@inline function inner_product(u1::AbstractVector, u2::AbstractVector, el::SpectralElement)
+ @unpack Npts, M = el
+ result = 0.0
+ @turbo for i in 1:Npts, j in 1:Npts
+ result += u1[i] * M[i,j] * u2[j]
end
+ return result
end
-@inline function integrate(u, el::SpectralElement)
+@inline function integrate(u::AbstractVector, el::SpectralElement)
@unpack quadr = el
@toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
if quadr === :LGL || quadr === :LGL_lumped
return LGL.integrate(u, el.z, el.w)
else # quadr === :GLGL
- return GLGL.integrate(u, el.z, el.w, el.v, el.D)
+ return GLGL.integrate(u, el.z, el.w)
end
end
@@ -117,15 +106,16 @@ end
if quadr === :LGL || quadr === :LGL_lumped
return LGL.integrate(fn, idxs, el.z, el.w)
else # quadr === :GLGL
- return GLGL.integrate(fn, idxs, el.z, el.w, el.v, el.D)
+ return GLGL.integrate(fn, idxs, el.z, el.w)
end
end
-@inline function differentiate(u, el::SpectralElement)
+@inline function differentiate(u::AbstractArray, el::SpectralElement)
return el.D * u
end
-@inline function differentiate!(du, u, el::SpectralElement)
+
+@inline function differentiate!(du::AbstractArray, u::AbstractArray, el::SpectralElement)
mul!(du, el.D, u)
end
diff --git a/test/IntegrationTests/refs/burgers_sine_avmda/output1d.h5 b/test/IntegrationTests/refs/burgers_sine_avmda/output1d.h5
index d2e27152e4d1ac944e2d38cf87060dd1e7100f3f..1a041f434509447ed292d10e1db86b39f72d8094 100644
Binary files a/test/IntegrationTests/refs/burgers_sine_avmda/output1d.h5 and b/test/IntegrationTests/refs/burgers_sine_avmda/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/euler_sod_shock_tube_bernstein/output1d.h5 b/test/IntegrationTests/refs/euler_sod_shock_tube_bernstein/output1d.h5
index 7d2cb4eb0cb6a458e63086ae2c59b38af0b3f9e4..387eebaa35970d5fb6599bdae781512889b74e08 100644
Binary files a/test/IntegrationTests/refs/euler_sod_shock_tube_bernstein/output1d.h5 and b/test/IntegrationTests/refs/euler_sod_shock_tube_bernstein/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5 b/test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5
index 2c258cd6e73e717289c01b4a4a953373f070538d..add2a996056d73e37a32e3de6aa665dd59fc14b9 100644
Binary files a/test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5 and b/test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5 differ
diff --git a/test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5 b/test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5
index 740e709f24c311ee88cb745632d6f3fb0f3db056..0bafcfb8acc690abf8d822611a6f20431289162b 100644
Binary files a/test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5 and b/test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/grhd_bondi_infall_avmda/output1d.h5 b/test/IntegrationTests/refs/grhd_bondi_infall_avmda/output1d.h5
index 6f9227a415288f4c49e200f88648973114f66f4f..220135c6a0abbf3b14330cca8d8708b86cf1097e 100644
Binary files a/test/IntegrationTests/refs/grhd_bondi_infall_avmda/output1d.h5 and b/test/IntegrationTests/refs/grhd_bondi_infall_avmda/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5 b/test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5
index 065bba8f1975121e930d94ae25151e5ccb47a7e3..8d7c701d82b877ec4bdcea885527db8d2f2716ea 100644
Binary files a/test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5 and b/test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5 differ
diff --git a/test/UnitTests/src/test_spectral.jl b/test/UnitTests/src/test_spectral.jl
index 3a2a38601b4d6e46d8767170b3bf5bf527f0c86d..cdb3f6f47d17015d4325fe53ed2179e7b4dcaa3a 100644
--- a/test/UnitTests/src/test_spectral.jl
+++ b/test/UnitTests/src/test_spectral.jl
@@ -37,18 +37,14 @@ end
f(x) = x^3 - 5.0 * x^2 + 1.0
f2_ana = 118.0/21.0 # according to wolfram alpha
- # use quadrature order N+1 = 3+1 = 4 here, because :LGL and :LGL_lumped are only exact up
- # to and including 2N-1
+ # LGL quadrature is only exact up to order 2*3-1=5, so use N=4 instead
N = 4
L = 2.0
- quadrs = [ :LGL, :GLGL, :LGL_lumped ]
-
- for quadr in quadrs
- el = SpectralElement(N, quadr)
- fx = f.(el.z)
- f2 = dg1d.inner_product(fx, fx, el)
- @test isapprox(f2, f2_ana)
- end
+ quadr = :LGL
+ el = SpectralElement(N, quadr)
+ fx = f.(el.z)
+ f2 = dg1d.inner_product(fx, fx, el)
+ @test isapprox(f2, f2_ana)
# check that :GLGL quadrature is exact up to and including order 2N
N = 3
@@ -64,8 +60,8 @@ end
@testset "Barycentric interpolation" begin
# setup
- n_pts = 6
- x, _ = dg1d.LGL.rule(n_pts)
+ Npts = 6
+ x, _ = dg1d.LGL.rule(Npts)
bw = dg1d.barycentric_weights(x)
# test function
@@ -83,7 +79,7 @@ end
fxintrp = [ dg1d.barycentric_interp(xi, x, bw, fx) for xi in xsmpl ]
@test all(isapprox.(fxintrp, fxsmpl))
- # test interpolation using matrix multiplication
+ # test interpolation using matrix multiplication
I = dg1d.barycentric_interp_matrix(xsmpl, x)
fxintrp = I * fx
@test all(isapprox.(fxintrp, fxsmpl))
@@ -101,62 +97,83 @@ end
end
-@testset "Operators GLGL" begin
+@testset "Test orthogonality relation GLGL" begin
- # set up quadrature rule
- N = 5
- N_pts = N + 1
- x, w, v = dg1d.GLGL.rule(N_pts)
- D = dg1d.first_derivative_matrix(x)
+ Npts = 6
+
+ x, w, v, D = dg1d.GLGL.inner_product_rule(Npts)
+ L = zeros(Npts,Npts)
+ for i in 1:Npts
+ Pi = map(x->dg1d.legendre(x, i-1), x)
+ for j in i:Npts
+ Pj = map(x->dg1d.legendre(x, j-1), x)
+ lij = dg1d.GLGL.inner_product(Pi, Pj, x, w, v, D)
+ L[i,j] = lij
+ if i != j
+ L[j,i] = L[i,j]
+ end
+ end
+ end
+ dg1d.purge_zeros!(L)
# analytical result
- L_ana = zeros(N_pts, N_pts)
- for i in 1:N_pts
+ L_ana = zeros(Npts, Npts)
+ for i in 1:Npts
L_ana[i,i] = 2.0 / (2.0 * i - 1.0)
end
- L = dg1d.GLGL.linear_dependency_matrix(x, w, v, D)
-
@test isapprox(L, L_ana)
end
-@testset "Operators LGL" begin
-
- # set up quadrature rule
- n_pts = 6
- x, w = dg1d.LGL.rule(n_pts)
- quad(f1, f2) = dg1d.LGL.inner_product(w, f1, f2)
-
- # analytical result
- L_ana = zeros(n_pts, n_pts)
- for i in 1:n_pts
- L_ana[i,i] = 2.0 / (2.0 * i - 1.0)
+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)
- L = dg1d.LGL.linear_dependency_matrix_exact(x, w)
-
- @test isapprox(L, L_ana)
-
+ return L
end
@testset "Test orthogonality relation" begin
- # set up quadrature rule
- n_pts = 6
- x, w = dg1d.LGL.rule(n_pts)
- quad(f1, f2) = dg1d.LGL.inner_product(w, f1, f2)
+ Npts = 6
+
+ Npts2 = Npts+2
+ x2, w2 = dg1d.LGL.rule(Npts2)
+ L = zeros(Npts,Npts)
+ for i in 1:Npts
+ Pi = map(x->dg1d.legendre(x, i-1), x2)
+ for j in i:Npts
+ Pj = map(x->dg1d.legendre(x, j-1), x2)
+ lij = dg1d.LGL.inner_product(Pi, Pj, x2, w2)
+ L[i,j] = lij
+ if i != j
+ L[j,i] = L[i,j]
+ end
+ end
+ end
+ dg1d.purge_zeros!(L)
# analytical result
- L_ana = zeros(n_pts, n_pts)
- for i in 1:n_pts
+ L_ana = zeros(Npts, Npts)
+ for i in 1:Npts
L_ana[i,i] = 2.0 / (2.0 * i - 1.0)
end
- L = dg1d.LGL.linear_dependency_matrix_exact(x, w)
-
@test isapprox(L, L_ana)
end
@@ -164,16 +181,15 @@ end
@testset "More integrator tests" begin
- n_pts = 6
+ Npts = 6
p = x -> (x-4)^2 * (x+1)^3
- x, w = dg1d.LGL.rule(n_pts)
+ x, w = dg1d.LGL.rule(Npts)
int = dg1d.LGL.integrate(p, x, w)
@test int ≈ 140/3
- x, w, v = dg1d.GLGL.rule(n_pts)
- D = dg1d.first_derivative_matrix(x)
- int = dg1d.GLGL.integrate(p, x, w, v, D)
+ x, w = dg1d.GLGL.rule(Npts)
+ int = dg1d.GLGL.integrate(p, x, w)
@test int ≈ 140/3
end