From 3763b4028d7c1667e556c58d598355b18b02f369 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 6 Oct 2023 13:11:42 -0400 Subject: [PATCH 01/14] adding initial extension for SummationByPartsOperators --- Project.toml | 7 + ext/StartUpDGSummationByPartsOperatorsExt.jl | 343 +++++++++++++++++++ 2 files changed, 350 insertions(+) create mode 100644 ext/StartUpDGSummationByPartsOperatorsExt.jl diff --git a/Project.toml b/Project.toml index 6886d8fd..ab5a709f 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,7 @@ Requires = "1" Setfield = "1" SimpleUnPack = "1" StaticArrays = "1" +SummationByPartsOperators = "0.5" Triangulate = "2" WriteVTK = "1" julia = "1.7" diff --git a/ext/StartUpDGSummationByPartsOperatorsExt.jl b/ext/StartUpDGSummationByPartsOperatorsExt.jl new file mode 100644 index 00000000..daaf1f22 --- /dev/null +++ b/ext/StartUpDGSummationByPartsOperatorsExt.jl @@ -0,0 +1,343 @@ +module StartUpDGSummationByPartsOperatorsExt + using StartUpDG + using SummationByPartsOperators + using SummationByPartsOperators: AbstractDerivativeOperator, + AbstractNonperiodicDerivativeOperator, DerivativeOperator, + AbstractPeriodicDerivativeOperator, + PeriodicDerivativeOperator, grid + + using LinearAlgebra: LinearAlgebra, Diagonal, diag, norm, UniformScaling + using SparseArrays: sparse, droptol! + + + 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 \ No newline at end of file From af1c1e38123ebd0154c5fc1ce994eeeaccab6c9f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 6 Oct 2023 13:16:58 -0400 Subject: [PATCH 02/14] adding missing package imports --- ext/StartUpDGSummationByPartsOperatorsExt.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/ext/StartUpDGSummationByPartsOperatorsExt.jl b/ext/StartUpDGSummationByPartsOperatorsExt.jl index daaf1f22..919e738d 100644 --- a/ext/StartUpDGSummationByPartsOperatorsExt.jl +++ b/ext/StartUpDGSummationByPartsOperatorsExt.jl @@ -1,13 +1,16 @@ module StartUpDGSummationByPartsOperatorsExt + + using LinearAlgebra: LinearAlgebra, Diagonal, diag, norm, UniformScaling + using SparseArrays: sparse, droptol! + using StartUpDG + using SummationByPartsOperators using SummationByPartsOperators: AbstractDerivativeOperator, AbstractNonperiodicDerivativeOperator, DerivativeOperator, AbstractPeriodicDerivativeOperator, PeriodicDerivativeOperator, grid - using LinearAlgebra: LinearAlgebra, Diagonal, diag, norm, UniformScaling - using SparseArrays: sparse, droptol! function construct_1d_operators(D::AbstractDerivativeOperator, tol) From 3586f7ca8ef3616576a3225a3628e9c7a472f330 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 11:12:02 -0400 Subject: [PATCH 03/14] add SBP extension tests --- test/Project.toml | 2 + test/SummationByPartsOperatorsExt_tests.jl | 60 ++++++++++++++++++++++ test/runtests.jl | 1 + 3 files changed, 63 insertions(+) create mode 100644 test/SummationByPartsOperatorsExt_tests.jl 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..b052c92b --- /dev/null +++ b/test/SummationByPartsOperatorsExt_tests.jl @@ -0,0 +1,60 @@ +using SummationByPartsOperators + +@testset "StartUpDGSummationByPartsOperatorsExt tests" 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(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 "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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 76e6f3a9..7ab74a14 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ include("named_array_partition_tests.jl") include("triangulate_tests.jl") include("reference_elem_tests.jl") include("sbp_tests.jl") +include("SummationByPartsOperatorsExt_tests.jl") include("MeshData_tests.jl") include("boundary_util_tests.jl") include("hybrid_mesh_tests.jl") From afd636c81c27ef10d8cb08f6324480a8d93f5848 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 11:12:57 -0400 Subject: [PATCH 04/14] rename multidimensional SBP tests update name of multidim SBP tests d --- test/{sbp_tests.jl => multidim_sbp_tests.jl} | 3 +-- test/runtests.jl | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) rename test/{sbp_tests.jl => multidim_sbp_tests.jl} (98%) 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 7ab74a14..c8f03cf0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,7 @@ 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") From 787e796aea8ae6ad3eadf9024e29d2682126f9e7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 11:13:11 -0400 Subject: [PATCH 05/14] update docs --- docs/src/RefElemData.md | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) 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} From 9cec8a828f0a94c9d314ae0cf0abf415e8b39690 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 11:16:10 -0400 Subject: [PATCH 06/14] bump Julia compat to 1.9 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ab5a709f..4ef4b4e7 100644 --- a/Project.toml +++ b/Project.toml @@ -46,4 +46,4 @@ StaticArrays = "1" SummationByPartsOperators = "0.5" Triangulate = "2" WriteVTK = "1" -julia = "1.7" +julia = "1.9" From 293240470519966916a413e77f02a0c9c0603e26 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 11:25:05 -0400 Subject: [PATCH 07/14] revert back to Julia compat 1.7, add Requires --- Project.toml | 2 +- ext/StartUpDGSummationByPartsOperatorsExt.jl | 17 ++++++++++++----- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index 4ef4b4e7..ab5a709f 100644 --- a/Project.toml +++ b/Project.toml @@ -46,4 +46,4 @@ StaticArrays = "1" SummationByPartsOperators = "0.5" Triangulate = "2" WriteVTK = "1" -julia = "1.9" +julia = "1.7" diff --git a/ext/StartUpDGSummationByPartsOperatorsExt.jl b/ext/StartUpDGSummationByPartsOperatorsExt.jl index 919e738d..6b583463 100644 --- a/ext/StartUpDGSummationByPartsOperatorsExt.jl +++ b/ext/StartUpDGSummationByPartsOperatorsExt.jl @@ -5,11 +5,18 @@ module StartUpDGSummationByPartsOperatorsExt using StartUpDG - using SummationByPartsOperators - using SummationByPartsOperators: AbstractDerivativeOperator, - AbstractNonperiodicDerivativeOperator, DerivativeOperator, - AbstractPeriodicDerivativeOperator, - PeriodicDerivativeOperator, grid + # 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 From b9f220c8f94356b3bb7f7d5cc62140936e407ca3 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 11:27:39 -0400 Subject: [PATCH 08/14] add forgotten part --- src/StartUpDG.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 3e1db386..7d8ed700 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -109,5 +109,11 @@ function __init__() end 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 # module From c4454b53acf61418b44ffe7c2fe90863da06f9d8 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 11:31:48 -0400 Subject: [PATCH 09/14] add [extras] to Project.toml --- Project.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Project.toml b/Project.toml index ab5a709f..8c73ec3e 100644 --- a/Project.toml +++ b/Project.toml @@ -47,3 +47,6 @@ SummationByPartsOperators = "0.5" Triangulate = "2" WriteVTK = "1" julia = "1.7" + +[extras] +SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" From 5c281f1ccdd85def8e0d2f05882a51b1bda1c694 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 11:41:02 -0400 Subject: [PATCH 10/14] try to fix tests to work with Requires for pre 1.9 julia --- test/SummationByPartsOperatorsExt_tests.jl | 2 -- test/runtests.jl | 5 +++++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/test/SummationByPartsOperatorsExt_tests.jl b/test/SummationByPartsOperatorsExt_tests.jl index b052c92b..50276cf9 100644 --- a/test/SummationByPartsOperatorsExt_tests.jl +++ b/test/SummationByPartsOperatorsExt_tests.jl @@ -1,5 +1,3 @@ -using SummationByPartsOperators - @testset "StartUpDGSummationByPartsOperatorsExt tests" begin tol = 200*eps() diff --git a/test/runtests.jl b/test/runtests.jl index c8f03cf0..af4b0d34 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,11 @@ 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") From aa76ccc06b6472b37f48115106264695fae9e3df Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 12:00:41 -0400 Subject: [PATCH 11/14] fix Requires by moving ext into `init` function --- src/StartUpDG.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 7d8ed700..0758b40e 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -107,13 +107,13 @@ function __init__() @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin include("TriangulatePlots.jl") end -end -# Until Julia v1.9 is the minimum required version for StartUpDG.jl, we still support Requires.jl -@static if !isdefined(Base, :get_extension) + # 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 end # module From 38b9c1e88bb344d1138adced0cf5e8924624bbaf Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 13:57:49 -0400 Subject: [PATCH 12/14] add dropped import --- ext/StartUpDGSummationByPartsOperatorsExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/StartUpDGSummationByPartsOperatorsExt.jl b/ext/StartUpDGSummationByPartsOperatorsExt.jl index 6b583463..2a685465 100644 --- a/ext/StartUpDGSummationByPartsOperatorsExt.jl +++ b/ext/StartUpDGSummationByPartsOperatorsExt.jl @@ -1,7 +1,7 @@ module StartUpDGSummationByPartsOperatorsExt using LinearAlgebra: LinearAlgebra, Diagonal, diag, norm, UniformScaling - using SparseArrays: sparse, droptol! +using SparseArrays: sparse, droptol!, spzeros using StartUpDG From 15ecd51fb9520cf5d9ba10d5ded2108d3c9ab3e0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 13:57:55 -0400 Subject: [PATCH 13/14] formatting --- ext/StartUpDGSummationByPartsOperatorsExt.jl | 734 ++++++++++--------- src/StartUpDG.jl | 29 +- 2 files changed, 412 insertions(+), 351 deletions(-) diff --git a/ext/StartUpDGSummationByPartsOperatorsExt.jl b/ext/StartUpDGSummationByPartsOperatorsExt.jl index 2a685465..e036cda8 100644 --- a/ext/StartUpDGSummationByPartsOperatorsExt.jl +++ b/ext/StartUpDGSummationByPartsOperatorsExt.jl @@ -1,353 +1,413 @@ module StartUpDGSummationByPartsOperatorsExt - using LinearAlgebra: LinearAlgebra, Diagonal, diag, norm, UniformScaling +using LinearAlgebra: LinearAlgebra, Diagonal, diag, norm, UniformScaling using SparseArrays: sparse, droptol!, spzeros - using StartUpDG +using StartUpDG - # Required for visualization code - if isdefined(Base, :get_extension) - using SummationByPartsOperators: SummationByPartsOperators, DerivativeOperator, grid, - AbstractDerivativeOperator, AbstractNonperiodicDerivativeOperator, - PeriodicDerivativeOperator, AbstractPeriodicDerivativeOperator +# 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 - # 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 + 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))) - 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) + 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 - - 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 + 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))) - 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 + LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf)) + + rstf = (rf, sf) + nrstJ = (nrJ, nsJ) - # 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) + # 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 \ No newline at end of file +end diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 0758b40e..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,16 +105,16 @@ export ck45 # LSERK 45 using Requires: @require function __init__() - @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") + @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 - + end # module From a755b6af2063dc84ae7cb259451ae15a449f0d9b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 7 Oct 2023 13:57:58 -0400 Subject: [PATCH 14/14] add mroe tests --- test/SummationByPartsOperatorsExt_tests.jl | 71 +++++++++++++++++++--- 1 file changed, 63 insertions(+), 8 deletions(-) diff --git a/test/SummationByPartsOperatorsExt_tests.jl b/test/SummationByPartsOperatorsExt_tests.jl index 50276cf9..90ad63b3 100644 --- a/test/SummationByPartsOperatorsExt_tests.jl +++ b/test/SummationByPartsOperatorsExt_tests.jl @@ -1,25 +1,25 @@ -@testset "StartUpDGSummationByPartsOperatorsExt tests" begin +@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(Quad(), D) + 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)) ≈ 4 - @test abs(sum(rd.wf)) ≈ 8 - @test abs(sum(rd.wf .* rd.nrJ)) + abs(sum(rd.wf .* rd.nsJ)) < tol + @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 @@ -55,4 +55,59 @@ @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