From f73e661fb704a9c6958be126cba5889cbb429435 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Sat, 7 Oct 2023 13:59:16 -0500 Subject: [PATCH] Add SummationByPartsOperators.jl package extension (#152) * adding initial extension for SummationByPartsOperators * adding missing package imports * add SBP extension tests * rename multidimensional SBP tests update name of multidim SBP tests d * update docs * bump Julia compat to 1.9 * revert back to Julia compat 1.7, add Requires * add forgotten part * add [extras] to Project.toml * try to fix tests to work with Requires for pre 1.9 julia * fix Requires by moving ext into `init` function * add dropped import * formatting * add mroe tests --- Project.toml | 10 + docs/src/RefElemData.md | 17 +- ext/StartUpDGSummationByPartsOperatorsExt.jl | 413 +++++++++++++++++++ src/StartUpDG.jl | 23 +- test/Project.toml | 2 + test/SummationByPartsOperatorsExt_tests.jl | 113 +++++ test/{sbp_tests.jl => multidim_sbp_tests.jl} | 3 +- test/runtests.jl | 8 +- 8 files changed, 576 insertions(+), 13 deletions(-) create mode 100644 ext/StartUpDGSummationByPartsOperatorsExt.jl create mode 100644 test/SummationByPartsOperatorsExt_tests.jl rename test/{sbp_tests.jl => multidim_sbp_tests.jl} (98%) diff --git a/Project.toml b/Project.toml index 6886d8fd..8c73ec3e 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,12 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" +[weakdeps] +SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" + +[extensions] +StartUpDGSummationByPartsOperatorsExt = "SummationByPartsOperators" + [compat] ConstructionBase = "1" FillArrays = "0.13, 1" @@ -37,6 +43,10 @@ Requires = "1" Setfield = "1" SimpleUnPack = "1" StaticArrays = "1" +SummationByPartsOperators = "0.5" Triangulate = "2" WriteVTK = "1" julia = "1.7" + +[extras] +SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" diff --git a/docs/src/RefElemData.md b/docs/src/RefElemData.md index dfda3b96..ab997da2 100644 --- a/docs/src/RefElemData.md +++ b/docs/src/RefElemData.md @@ -72,9 +72,22 @@ By default, `RefElemData` is constructed for a nodal basis (in order to facilita ## `RefElemData` based on SBP finite differences -It is also possible to construct a [`RefElemData`](@ref) based on [multi-dimensional SBP finite difference operators](https://doi.org/10.1137/15M1038360). These utilize nodes constructed by [Tianheng Chen and Chi-Wang Shu](https://doi.org/10.1016/j.jcp.2017.05.025), [Ethan Kubatko](https://sites.google.com/site/chilatosu/ethan-bio), and [Jason Hicken](https://doi.org/10.1007/s10915-020-01154-8). +It is also possible to construct a [`RefElemData`](@ref) based on both traditional finite difference SBP operators and [multi-dimensional SBP finite difference operators](https://doi.org/10.1137/15M1038360). The traditional finite difference SBP operators are built using [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl). The multi-dimensional SBP operators utilize nodes constructed by [Tianheng Chen and Chi-Wang Shu](https://doi.org/10.1016/j.jcp.2017.05.025), [Ethan Kubatko](https://sites.google.com/site/chilatosu/ethan-bio), and [Jason Hicken](https://doi.org/10.1007/s10915-020-01154-8). -Some examples: +Some examples of traditional finite difference SBP operators on tensor product domains: +```julia +using StartUpDG +using SummationByPartsOperators # this package must be loaded to enable extension + +D = derivative_operator(MattssonNordström2004(), + derivative_order=1, accuracy_order=N+1, + xmin=-1.0, xmax=1.0, N=20) + +rd = RefElemData(Line(), D) +rd = RefElemData(Quad(), D) +rd = RefElemData(Hex(), D) +``` +Some examples of multi-dimensional SBP operators: ```julia N = 3 rd = RefElemData(Quad(), SBP(), N) # defaults to SBP{TensorProductLobatto} diff --git a/ext/StartUpDGSummationByPartsOperatorsExt.jl b/ext/StartUpDGSummationByPartsOperatorsExt.jl new file mode 100644 index 00000000..e036cda8 --- /dev/null +++ b/ext/StartUpDGSummationByPartsOperatorsExt.jl @@ -0,0 +1,413 @@ +module StartUpDGSummationByPartsOperatorsExt + +using LinearAlgebra: LinearAlgebra, Diagonal, diag, norm, UniformScaling +using SparseArrays: sparse, droptol!, spzeros + +using StartUpDG + +# Required for visualization code +if isdefined(Base, :get_extension) + using SummationByPartsOperators: + SummationByPartsOperators, + DerivativeOperator, + grid, + AbstractDerivativeOperator, + AbstractNonperiodicDerivativeOperator, + PeriodicDerivativeOperator, + AbstractPeriodicDerivativeOperator +else + # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl + using ..SummationByPartsOperators + using ..SummationByPartsOperators: + AbstractDerivativeOperator, + AbstractPeriodicDerivativeOperator, + AbstractNonperiodicDerivativeOperator, + DerivativeOperator, + PeriodicDerivativeOperator, + grid +end + +function construct_1d_operators(D::AbstractDerivativeOperator, tol) + nodes_1d = collect(grid(D)) + M = SummationByPartsOperators.mass_matrix(D) + if M isa UniformScaling + weights_1d = M * ones(Bool, length(nodes_1d)) + else + weights_1d = diag(M) + end + + # StartUpDG assumes nodes from -1 to +1. Thus, we need to re-scale everything. + # We can adjust the grid spacing as follows. + xmin = SummationByPartsOperators.xmin(D) + xmax = SummationByPartsOperators.xmax(D) + factor = 2 / (xmax - xmin) + @. nodes_1d = factor * (nodes_1d - xmin) - 1 + @. weights_1d = factor * weights_1d + + D_1d = droptol!(inv(factor) * sparse(D), tol) + I_1d = Diagonal(ones(Bool, length(nodes_1d))) + + return nodes_1d, weights_1d, D_1d, I_1d +end + +function StartUpDG.RefElemData(element_type::Line, + D::AbstractDerivativeOperator; + tol = 100 * eps(),) + approximation_type = D + N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree + + # 1D operators + nodes_1d, weights_1d, D_1d = construct_1d_operators(D, tol) + + # volume + rq = r = nodes_1d + wq = weights_1d + Dr = D_1d + M = Diagonal(wq) + Pq = LinearAlgebra.I + Vq = LinearAlgebra.I + + VDM = nothing # unused generalized Vandermonde matrix + + rst = (r,) + rstq = (rq,) + Drst = (Dr,) + + # face + face_vertices = StartUpDG.face_vertices(element_type) + face_mask = [1, length(nodes_1d)] + + rf = [-1.0; 1.0] + nrJ = [-1.0; 1.0] + wf = [1.0; 1.0] + if D isa AbstractPeriodicDerivativeOperator + # we do not need any face stuff for periodic operators + Vf = spzeros(length(wf), length(wq)) + else + Vf = sparse([1, 2], [1, length(nodes_1d)], [1.0, 1.0]) + end + LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf)) + + rstf = (rf,) + nrstJ = (nrJ,) + + # low order interpolation nodes + r1 = StartUpDG.nodes(element_type, 1) + V1 = StartUpDG.vandermonde(element_type, 1, r) / + StartUpDG.vandermonde(element_type, 1, r1) + + return RefElemData(element_type, + approximation_type, + N, + face_vertices, + V1, + rst, + VDM, + face_mask, + rst, + LinearAlgebra.I, # plotting + rstq, + wq, + Vq, # quadrature + rstf, + wf, + Vf, + nrstJ, # faces + M, + Pq, + Drst, + LIFT) +end + +function StartUpDG.RefElemData(element_type::Quad, + D::AbstractDerivativeOperator; + tol = 100 * eps(),) + approximation_type = D + N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree + + # 1D operators + nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol) + + # volume + s, r = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d)) # this is to match + # ordering of nrstJ + rq = r + sq = s + wr, ws = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d)) + wq = wr .* ws + Dr = kron(I_1d, D_1d) + Ds = kron(D_1d, I_1d) + M = Diagonal(wq) + Pq = LinearAlgebra.I + Vq = LinearAlgebra.I + + VDM = nothing # unused generalized Vandermonde matrix + + rst = (r, s) + rstq = (rq, sq) + Drst = (Dr, Ds) + + # face + face_vertices = StartUpDG.face_vertices(element_type) + face_mask = vcat(StartUpDG.find_face_nodes(element_type, r, s)...) + + rf, sf, wf, nrJ, nsJ = StartUpDG.init_face_data(element_type, + quad_rule_face = (nodes_1d, weights_1d)) + if D isa AbstractPeriodicDerivativeOperator + # we do not need any face stuff for periodic operators + Vf = spzeros(length(wf), length(wq)) + else + Vf = sparse(eachindex(face_mask), face_mask, ones(Bool, length(face_mask))) + end + LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf)) + + rstf = (rf, sf) + nrstJ = (nrJ, nsJ) + + # low order interpolation nodes + r1, s1 = StartUpDG.nodes(element_type, 1) + V1 = StartUpDG.vandermonde(element_type, 1, r, s) / + StartUpDG.vandermonde(element_type, 1, r1, s1) + + return RefElemData(element_type, + approximation_type, + N, + face_vertices, + V1, + rst, + VDM, + face_mask, + rst, + LinearAlgebra.I, # plotting + rstq, + wq, + Vq, # quadrature + rstf, + wf, + Vf, + nrstJ, # faces + M, + Pq, + Drst, + LIFT) +end + +function StartUpDG.RefElemData(element_type::Hex, + D::AbstractDerivativeOperator; + tol = 100 * eps(),) + approximation_type = D + N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree + + # 1D operators + nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol) + + # volume + # to match ordering of nrstJ + s, r, t = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d, nodes_1d)) + rq = r + sq = s + tq = t + wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d, weights_1d)) + wq = wr .* ws .* wt + Dr = kron(I_1d, I_1d, D_1d) + Ds = kron(I_1d, D_1d, I_1d) + Dt = kron(D_1d, I_1d, I_1d) + M = Diagonal(wq) + Pq = LinearAlgebra.I + Vq = LinearAlgebra.I + + VDM = nothing # unused generalized Vandermonde matrix + + rst = (r, s, t) + rstq = (rq, sq, tq) + Drst = (Dr, Ds, Dt) + + # face + face_vertices = StartUpDG.face_vertices(element_type) + face_mask = vcat(StartUpDG.find_face_nodes(element_type, r, s, t)...) + + rf, sf, tf, wf, nrJ, nsJ, ntJ = let + rf, sf = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d)) + wr, ws = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d)) + wf = wr .* ws + StartUpDG.init_face_data(element_type, quad_rule_face = (rf, sf, wf)) + end + Vf = sparse(eachindex(face_mask), face_mask, ones(Bool, length(face_mask))) + LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf)) + + rstf = (rf, sf, tf) + nrstJ = (nrJ, nsJ, ntJ) + + # low order interpolation nodes + r1, s1, t1 = StartUpDG.nodes(element_type, 1) + V1 = StartUpDG.vandermonde(element_type, 1, r, s, t) / + StartUpDG.vandermonde(element_type, 1, r1, s1, t1) + + return RefElemData(element_type, + approximation_type, + N, + face_vertices, + V1, + rst, + VDM, + face_mask, + rst, + LinearAlgebra.I, # plotting + rstq, + wq, + Vq, # quadrature + rstf, + wf, + Vf, + nrstJ, # faces + M, + Pq, + Drst, + LIFT) +end + +# specialized Hex constructor in 3D to reduce memory usage. +function StartUpDG.RefElemData(element_type::Hex, + D::AbstractPeriodicDerivativeOperator; + tol = 100 * eps(),) + approximation_type = D + N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree + + # 1D operators + nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol) + + # volume + # to match ordering of nrstJ + s, r, t = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d, nodes_1d)) + rq = r + sq = s + tq = t + wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d, weights_1d)) + wq = wr .* ws .* wt + Dr = kron(I_1d, I_1d, D_1d) + Ds = kron(I_1d, D_1d, I_1d) + Dt = kron(D_1d, I_1d, I_1d) + M = Diagonal(wq) + Pq = LinearAlgebra.I + Vq = LinearAlgebra.I + + VDM = nothing # unused generalized Vandermonde matrix + + rst = (r, s, t) + rstq = (rq, sq, tq) + Drst = (Dr, Ds, Dt) + + # face + # We do not need any face data for periodic operators. Thus, we just + # pass `nothing` to save memory. + face_vertices = ntuple(_ -> nothing, 3) + face_mask = nothing + wf = nothing + rstf = ntuple(_ -> nothing, 3) + nrstJ = ntuple(_ -> nothing, 3) + Vf = nothing + LIFT = nothing + + # low order interpolation nodes + V1 = nothing # do not need to store V1, since we specialize StartUpDG.MeshData to avoid using it. + + return RefElemData(element_type, + approximation_type, + N, + face_vertices, + V1, + rst, + VDM, + face_mask, + rst, + LinearAlgebra.I, # plotting + rstq, + wq, + Vq, # quadrature + rstf, + wf, + Vf, + nrstJ, # faces + M, + Pq, + Drst, + LIFT) +end + +function Base.show(io::IO, + mime::MIME"text/plain", + rd::RefElemData{NDIMS, ElementType, ApproximationType}) where { + NDIMS, + ElementType <: StartUpDG.AbstractElemShape, + ApproximationType <: AbstractDerivativeOperator, +} + @nospecialize rd + print(io, "RefElemData for an approximation using an ") + show(IOContext(io, :compact => true), rd.approximation_type) + print(io, " on $(rd.element_type) element") +end + +function Base.show(io::IO, + rd::RefElemData{NDIMS, ElementType, ApproximationType}) where { + NDIMS, + ElementType <: StartUpDG.AbstractElemShape, + ApproximationType <: AbstractDerivativeOperator, +} + @nospecialize rd + print(io, "RefElemData{", summary(rd.approximation_type), ", ", rd.element_type, "}") +end + +function StartUpDG.inverse_trace_constant(rd::RefElemData{ + NDIMS, + ElementType, + ApproximationType, +}) where { + NDIMS, + ElementType <: Union{Line, Quad, Hex}, + ApproximationType <: AbstractDerivativeOperator, +} + D = rd.approximation_type + + # the inverse trace constant is the maximum eigenvalue corresponding to + # M_f * v = λ * M * v + # where M_f is the face mass matrix and M is the volume mass matrix. + # Since M is diagonal and since M_f is just the boundary "mask" matrix + # (which extracts the first and last entries of a vector), the maximum + # eigenvalue is the inverse of the first or last mass matrix diagonal. + left_weight = SummationByPartsOperators.left_boundary_weight(D) + right_weight = SummationByPartsOperators.right_boundary_weight(D) + max_eigenvalue = max(inv(left_weight), inv(right_weight)) + + # For tensor product elements, the trace constant for higher dimensional + # elements is the one-dimensional trace constant multiplied by `NDIMS`. See + # "GPU-accelerated discontinuous Galerkin methods on hybrid meshes." + # Chan, Jesse, et al (2016), https://doi.org/10.1016/j.jcp.2016.04.003 + # for more details (specifically, Appendix A.1, Theorem A.4). + return NDIMS * max_eigenvalue +end + +# This is used in `estimate_dt`. `estimate_h` uses that `Jf / J = O(h^{NDIMS-1}) / O(h^{NDIMS}) = O(1/h)`. +# However, since we do not initialize `Jf` for periodic FDSBP operators, we specialize `estimate_h` +# based on the reference grid provided by SummationByPartsOperators.jl and information about the domain size +# provided by `md::MeshData``. +function StartUpDG.estimate_h(e, + rd::RefElemData{NDIMS, ElementType, ApproximationType}, + md::MeshData) where { + NDIMS, + ElementType <: StartUpDG.AbstractElemShape, + ApproximationType <: SummationByPartsOperators.AbstractPeriodicDerivativeOperator, +} + D = rd.approximation_type + x = grid(D) + + # we assume all SummationByPartsOperators.jl reference grids are rescaled to [-1, 1] + xmin = SummationByPartsOperators.xmin(D) + xmax = SummationByPartsOperators.xmax(D) + factor = 2 / (xmax - xmin) + + # If the domain has size L^NDIMS, then `minimum(md.J)^(1 / NDIMS) = L`. + # WARNING: this is not a good estimate on anisotropic grids. + return minimum(diff(x)) * factor * minimum(md.J)^(1 / NDIMS) +end + +end diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 3e1db386..fb98506e 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -1,4 +1,4 @@ -module StartUpDG +module StartUpDG using Reexport: @reexport @@ -6,7 +6,8 @@ using ConstructionBase: ConstructionBase using FillArrays: Ones, Zeros, Fill using HDF5: h5open # used to read in SBP triangular node data using Kronecker: kronecker # for Hex element matrix manipulations -using LinearAlgebra: cond, diagm, eigvals, Diagonal, UniformScaling, I, mul!, norm, qr, ColumnNorm, Symmetric +using LinearAlgebra: + cond, diagm, eigvals, Diagonal, UniformScaling, I, mul!, norm, qr, ColumnNorm, Symmetric using NodesAndModes: meshgrid, find_face_nodes, face_vertices @reexport using NodesAndModes # for basis functions using PathIntersections: PathIntersections @@ -77,8 +78,8 @@ export uniform_mesh include("mesh/gmsh_utilities.jl") export read_Gmsh_2D # unifies v2.2.8 and v4.1 mesh reading export readGmsh2D, readGmsh2D_v4 # TODO: deprecate -export read_Gmsh_2D_v2, read_Gmsh_2D_v4 -export MeshImportOptions +export read_Gmsh_2D_v2, read_Gmsh_2D_v4 +export MeshImportOptions include("mesh/hohqmesh_utilities.jl") export read_HOHQMesh @@ -104,10 +105,16 @@ export ck45 # LSERK 45 using Requires: @require function __init__() - @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin - include("TriangulatePlots.jl") - end + @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin + include("TriangulatePlots.jl") + end + + # Until Julia v1.9 is the minimum required version for StartUpDG.jl, we still support Requires.jl + @static if !isdefined(Base, :get_extension) + @require SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" begin + include("../ext/StartUpDGSummationByPartsOperatorsExt.jl") + end + end end - end # module diff --git a/test/Project.toml b/test/Project.toml index 013caf76..b65cc55d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" +SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -12,4 +13,5 @@ HOHQMesh = "0.2" Plots = "1" StaticArrays = "1" StructArrays = "0.6" +SummationByPartsOperators = "0.5" Suppressor = "0.2" diff --git a/test/SummationByPartsOperatorsExt_tests.jl b/test/SummationByPartsOperatorsExt_tests.jl new file mode 100644 index 00000000..90ad63b3 --- /dev/null +++ b/test/SummationByPartsOperatorsExt_tests.jl @@ -0,0 +1,113 @@ + +@testset "SummationByPartsOperators ext tests for nonperiodic operators" begin + tol = 200*eps() + N = 3 + f(N,r,s) = r^N + s^N # function for testing differentiation + + D = derivative_operator(MattssonNordström2004(), + derivative_order=1, accuracy_order=N+1, + xmin=-1.0, xmax=1.0, N=20) + @testset "Line" begin + rd = RefElemData(Line(), D) + @test rd.r == rd.rst[1] + @test rd.Np == length(rd.r) + @test rd.Nq == length(rd.rq) + @test abs(sum(rd.wq)) ≈ 2 + @test abs(sum(rd.wf .* rd.nrJ)) < tol + @test rd.Pq*rd.Vq ≈ I + @test all(vec.(rd.rstf) .≈ (x->getindex(x,rd.Fmask)).(rd.rst)) + @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd) + + md = MeshData(uniform_mesh(Line(), 2), rd) + @test estimate_h(1, rd, md) ≈ 0.5 + end + + @testset "Quad" begin + rd = RefElemData(Quad(), D) + @test rd.r == rd.rst[1] + @test rd.Np == length(rd.r) + @test rd.Nq == length(rd.rq) + @test abs(sum(rd.wq)) ≈ 4 + @test abs(sum(rd.wf)) ≈ 8 + @test abs(sum(rd.wf .* rd.nrJ)) + abs(sum(rd.wf .* rd.nsJ)) < tol + @test rd.Pq*rd.Vq ≈ I + @test all(vec.(rd.rstf) .≈ (x->getindex(x,rd.Fmask)).(rd.rst)) + @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd) + end + + @testset "Hex" begin + rd = RefElemData(Hex(), D) + @test propertynames(rd)[1] == :element_type + @test rd.t == rd.rst[3] + @test rd.tf == rd.rstf[3] + @test rd.tq == rd.rstq[3] + @test rd.rp == rd.rstp[1] + @test rd.sp == rd.rstp[2] + @test rd.tp == rd.rstp[3] + @test rd.Np == length(rd.r) + @test rd.Nq == length(rd.rq) + @test abs(sum(rd.wq)) ≈ 8 + @test abs(sum(rd.wf)) ≈ 6*4 + @test abs(sum(rd.wf .* rd.nrJ)) < tol + @test abs(sum(rd.wf .* rd.nsJ)) < tol + @test abs(sum(rd.wf .* rd.ntJ)) < tol + @test rd.Pq*rd.Vq ≈ I + @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd) + end + +end + +@testset "SummationByPartsOperators ext tests for periodic operators" begin + tol = 200*eps() + N = 3 + D = periodic_derivative_operator(derivative_order=1, accuracy_order=N+1, + xmin=0.0, xmax=2.0, N=20) + + @testset "Line" begin + rd = RefElemData(Line(), D) + @test rd.r == rd.rst[1] + @test rd.Np == length(rd.r) + @test rd.Nq == length(rd.rq) + @test abs(sum(rd.wq)) ≈ 2 + + # diff matrices should be skew symmetric + Qrst = (A -> rd.M * A).(rd.Drst) + @test all((A -> norm(A+A')).(Qrst) .< tol) + + @test rd.Pq*rd.Vq ≈ I + + end + + @testset "Quad" begin + rd = RefElemData(Quad(), D) + @test rd.r == rd.rst[1] + @test rd.Np == length(rd.r) + @test rd.Nq == length(rd.rq) + @test abs(sum(rd.wq)) ≈ 4 + + # diff matrices should be skew symmetric + Qrst = (A -> rd.M * A).(rd.Drst) + @test all((A -> norm(A+A')).(Qrst) .< tol) + + @test rd.Pq*rd.Vq ≈ I + end + + @testset "Hex" begin + rd = RefElemData(Hex(), D) + @test propertynames(rd)[1] == :element_type + @test rd.t == rd.rst[3] + @test rd.tq == rd.rstq[3] + @test rd.rp == rd.rstp[1] + @test rd.sp == rd.rstp[2] + @test rd.tp == rd.rstp[3] + @test rd.Np == length(rd.r) + @test rd.Nq == length(rd.rq) + @test abs(sum(rd.wq)) ≈ 8 + + # diff matrices should be skew symmetric + Qrst = (A -> rd.M * A).(rd.Drst) + @test all((A -> norm(A+A')).(Qrst) .< tol) + + @test rd.Pq*rd.Vq ≈ I + end +end \ No newline at end of file diff --git a/test/sbp_tests.jl b/test/multidim_sbp_tests.jl similarity index 98% rename from test/sbp_tests.jl rename to test/multidim_sbp_tests.jl index 164cf167..2ff50614 100644 --- a/test/sbp_tests.jl +++ b/test/multidim_sbp_tests.jl @@ -1,4 +1,4 @@ -@testset "SBP elements" begin +@testset "Multidimensional SBP tests" begin tol = 200*eps() @@ -64,5 +64,4 @@ @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd) end - end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 76e6f3a9..af4b0d34 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,13 +1,19 @@ using Test using Suppressor using LinearAlgebra + +# we need to load this first before StartUpDG for Julia pre-v1.9 +# since Requires.jl needs us to load this first. +using SummationByPartsOperators + using StartUpDG include("write_vtk_tests.jl") include("named_array_partition_tests.jl") include("triangulate_tests.jl") include("reference_elem_tests.jl") -include("sbp_tests.jl") +include("multidim_sbp_tests.jl") +include("SummationByPartsOperatorsExt_tests.jl") include("MeshData_tests.jl") include("boundary_util_tests.jl") include("hybrid_mesh_tests.jl")