Skip to content
Snippets Groups Projects
Commit a940bcc6 authored by Florian Atteneder's avatar Florian Atteneder
Browse files

dg1d: refactor lgl and glgl modules (!204) 

Unify interface for LGL and GLGL rules and update `SpectralElement` accordingly:
- most routine should just work with one argument `Npts`
- add more doc strings
- remove unused methods `linear_dependency_matrix` and `expansion_matrix`
- remove unused fields of `SpectralElement`
parent 70e66d2b
No related branches found
No related tags found
1 merge request!204dg1d: refactor lgl and glgl modules
Pipeline #7166 passed
Hide whitespace changes
Inline Side-by-side
Showing
with 193 additions and 404 deletions
src/glgl.jl
+ 67
− 168
View file @ a940bcc6
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
src/lgl.jl
+ 38
− 154
View file @ a940bcc6
......@@ -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
src/spectralelement.jl
+ 21
− 31
View file @ a940bcc6
......@@ -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
test/IntegrationTests/refs/burgers_sine_avmda/output1d.h5
+ 0
− 0
View file @ a940bcc6
No preview for this file type
test/IntegrationTests/refs/euler_sod_shock_tube_bernstein/output1d.h5
+ 0
− 0
View file @ a940bcc6
No preview for this file type
test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5
+ 0
− 0
View file @ a940bcc6
No preview for this file type
test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5
+ 0
− 0
View file @ a940bcc6
No preview for this file type
test/IntegrationTests/refs/grhd_bondi_infall_avmda/output1d.h5
+ 0
− 0
View file @ a940bcc6
No preview for this file type
test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5
+ 0
− 0
View file @ a940bcc6
No preview for this file type
test/UnitTests/src/test_spectral.jl
+ 67
− 51
View file @ a940bcc6
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment