From 801db563bf29667fef80618c04a66aec39106bd5 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 16 Jan 2024 10:20:05 -0600 Subject: [PATCH 01/13] add subcell limiting operators --- src/RefElemData_SBP.jl | 86 ++++++++++++++++++++++++++++++++++++++++++ src/StartUpDG.jl | 5 ++- 2 files changed, 89 insertions(+), 2 deletions(-) diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index dbc780c6..35055dad 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -376,3 +376,89 @@ function sparse_low_order_SBP_operators(rd::RefElemData{3, Hex, <:Union{<:SBP, < return sparse.((Qr, Qs, Qt)), sparse(E) end + +# builds subcell operators D * R such that for r such that sum(r) = 0, +# D * Diagonal(θ) * R * r is also diagonal for any choice of θ. This is +# useful in subcell limiting (see, for example +# https://doi.org/10.1016/j.compfluid.2022.105627 for a 1D example) +function subcell_limiting_operators(Qr::AbstractMatrix; tol = 100 * eps()) + Qr = sparse(Qr) + Sr = droptol!(Qr - Qr', tol) + Br = droptol!(Qr + Qr', tol) + + num_interior_fluxes = nnz(Sr) ÷ 2 + num_boundary_indices = nnz(Br) + num_fluxes = num_interior_fluxes + num_boundary_indices + + interior_indices = findall(Sr .> tol) + boundary_indices = findall(abs.(Br) .> tol) + + matrix_indices = zeros(Int, size(Sr)) + for (i, ij) in enumerate(boundary_indices) + matrix_indices[ij] = i + matrix_indices[ij.I[2], ij.I[1]] = i + end + + for (i, ij) in enumerate(interior_indices) + matrix_indices[ij] = i + length(boundary_indices) + matrix_indices[ij.I[2], ij.I[1]] = i + length(boundary_indices) + end + + D = zeros(size(Qr, 2), num_fluxes) + for i in axes(matrix_indices, 1), j in axes(matrix_indices, 2) + if matrix_indices[i, j] !== 0 + D[i, matrix_indices[i, j]] = Qr[i, j] + end + end + + d = (ones(size(D, 1))' * D)' + ids = findall(@. abs(d) > 1e2 * eps()) + + Z = nullspace(D)[ids, :] + A = pinv(D)[ids,:] + X = randn(size(Z, 2), size(D, 1)) + y = randn(length(ids)) + e = ones(size(D, 1)) + + # compute R via linear algebra + m, n = size(A) + z = size(Z, 2) + C = hcat(kron(I(n), Z), kron(ones(n), I(m))) + sol = pinv(C) * -vec(A) + X = reshape(sol[1:n*z], (z, n)) + R = pinv(D) + nullspace(D) * X + + # check if R satisfies our pseudoinverse condition + @assert norm(D * R - I(size(D,1))) < tol * length(D) + + return D, R +end + +""" + Δrst, Rrst = subcell_limiting_operators(rd::RefElemData) + +Returns tuples of subcell limiting operators Drst = (Δr, Δs, ...) and R = (Rr, Rs, ...) +such that for r where sum(r) = 0, sum(D * Diagonal(θ) * R * r) = 0 for any choice of θ. +These operators are useful for conservative subcell limiting techniques (see +https://doi.org/10.1016/j.compfluid.2022.105627 for an example of such an approach on +tensor product elements). + +Sparse SBP operators used in an intermediate step when buidling these subcell +limiting operators; by default, these operators are constructed using +`sparse_low_order_SBP_operators`. To construct subcell limiting operators for a +general SBP operator, one can use the following: + + Δ, R = subcell_limiting_operators(Q::AbstractMatrix; tol = 100 * eps()) +""" +function subcell_limiting_operators(rd::RefElemData) + Qrst, _ = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = subcell_limiting_operators.(Qrst) + return zip(Δrst, Rrst) +end + +# specialize for single dimension +function subcell_limiting_operators(rd::RefElemData{1}) + Qrst, _ = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = subcell_limiting_operators(Qrst[1]) + return Δrst, Rrst +end diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 41b098f0..13681022 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -7,7 +7,7 @@ using FillArrays: 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 + cond, diagm, eigvals, Diagonal, UniformScaling, I, mul!, norm, qr, ColumnNorm, Symmetric, nullspace, pinv using NodesAndModes: meshgrid, find_face_nodes, face_vertices @reexport using NodesAndModes # for basis functions using PathIntersections: PathIntersections @@ -17,7 +17,7 @@ using RecipesBase: RecipesBase using StaticArrays: SVector, SMatrix using Setfield: setproperties, @set # for "modifying" structs (setproperties) @reexport using SimpleUnPack: @unpack -using SparseArrays: sparse, droptol!, blockdiag +using SparseArrays: sparse, droptol!, blockdiag, nnz using Triangulate: Triangulate, TriangulateIO, triangulate @reexport using WriteVTK @@ -37,6 +37,7 @@ include("RefElemData_SBP.jl") export SBP, DefaultSBPType, TensorProductLobatto, Hicken, Kubatko # types for SBP node dispatch export LobattoFaceNodes, LegendreFaceNodes # type parameters for SBP{Kubatko{...}} export hybridized_SBP_operators, sparse_low_order_SBP_operators +export subcell_limiting_operators export inverse_trace_constant, face_type include("ref_elem_utils.jl") From ede4c696151d06e814c2092933474943e60e3557 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 16 Jan 2024 10:20:17 -0600 Subject: [PATCH 02/13] refactor tests --- test/reference_elem_tests.jl | 49 ---------------------- test/runtests.jl | 29 ++++++------- test/sparse_SBP_operator_tests.jl | 68 +++++++++++++++++++++++++++++++ 3 files changed, 83 insertions(+), 63 deletions(-) create mode 100644 test/sparse_SBP_operator_tests.jl diff --git a/test/reference_elem_tests.jl b/test/reference_elem_tests.jl index 7d487087..73a23000 100644 --- a/test/reference_elem_tests.jl +++ b/test/reference_elem_tests.jl @@ -336,55 +336,6 @@ end end end -@testset "Sparse SBP operators" begin - tol = 5e2*eps() - N = 3 - @testset "Line" begin - rd = RefElemData(Line(), Polynomial{Gauss}(), N) - (Q,), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Q, dims=2)) < tol - @test Q + Q' ≈ E' * Diagonal([-1,1]) * E - end - - @testset "Tri" begin - rd = RefElemData(Tri(), N) - (Qr, Qs), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Qr, dims=2)) < tol - @test norm(sum(Qs, dims=2)) < tol - @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E - @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E - end - - @testset "Quad (SBP)" begin - rd = RefElemData(Quad(), SBP(), N) - (Qr, Qs), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Qr, dims=2)) < tol - @test norm(sum(Qs, dims=2)) < tol - @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E - @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E - end - - @testset "Quad (Polynomial{Gauss})" begin - rd = RefElemData(Quad(), Gauss(), N) - (Qr, Qs), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Qr, dims=2)) < tol - @test norm(sum(Qs, dims=2)) < tol - @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E - @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E - end - - @testset "Hex (Polynomial{Gauss})" begin - rd = RefElemData(Hex(), Gauss(), N) - (Qr, Qs, Qt), E = sparse_low_order_SBP_operators(rd) - @test norm(sum(Qr, dims=2)) < tol - @test norm(sum(Qs, dims=2)) < tol - @test norm(sum(Qt, dims=2)) < tol - @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E - @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E - @test Qt + Qt' ≈ E' * Diagonal(rd.wf .* rd.ntJ) * E - end -end - @testset "RefElemData SBP keyword arguments" begin for elem in [Line(), Tri(), Quad(), Hex()] # make sure Nplot is accepted as a kwarg diff --git a/test/runtests.jl b/test/runtests.jl index af4b0d34..9be78b5f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,17 +8,18 @@ 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("multidim_sbp_tests.jl") -include("SummationByPartsOperatorsExt_tests.jl") -include("MeshData_tests.jl") -include("boundary_util_tests.jl") -include("hybrid_mesh_tests.jl") -include("noncon_mesh_tests.jl") -include("cut_mesh_tests.jl") -include("misc_tests.jl") -include("gmsh_parse_tests.jl") -include("hohqmesh_tests.jl") +# include("write_vtk_tests.jl") +# include("named_array_partition_tests.jl") +# include("triangulate_tests.jl") +# include("reference_elem_tests.jl") +# include("multidim_sbp_tests.jl") +include("sparse_SBP_operator_tests.jl") +# include("SummationByPartsOperatorsExt_tests.jl") +# include("MeshData_tests.jl") +# include("boundary_util_tests.jl") +# include("hybrid_mesh_tests.jl") +# include("noncon_mesh_tests.jl") +# include("cut_mesh_tests.jl") +# include("misc_tests.jl") +# include("gmsh_parse_tests.jl") +# include("hohqmesh_tests.jl") diff --git a/test/sparse_SBP_operator_tests.jl b/test/sparse_SBP_operator_tests.jl new file mode 100644 index 00000000..dff68673 --- /dev/null +++ b/test/sparse_SBP_operator_tests.jl @@ -0,0 +1,68 @@ +@testset "Sparse SBP operators for approx_type <: Polynomial" begin + tol = 5e2*eps() + N = 3 + @testset "Line" begin + rd = RefElemData(Line(), Polynomial{Gauss}(), N) + (Q,), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Q, dims=2)) < tol + @test Q + Q' ≈ E' * Diagonal([-1,1]) * E + end + + @testset "Tri" begin + rd = RefElemData(Tri(), N) + (Qr, Qs), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Qr, dims=2)) < tol + @test norm(sum(Qs, dims=2)) < tol + @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E + @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E + end + + @testset "Quad (SBP)" begin + rd = RefElemData(Quad(), SBP(), N) + (Qr, Qs), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Qr, dims=2)) < tol + @test norm(sum(Qs, dims=2)) < tol + @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E + @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E + end + + @testset "Quad (Polynomial{Gauss})" begin + rd = RefElemData(Quad(), Gauss(), N) + (Qr, Qs), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Qr, dims=2)) < tol + @test norm(sum(Qs, dims=2)) < tol + @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E + @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E + end + + @testset "Hex (Polynomial{Gauss})" begin + rd = RefElemData(Hex(), Gauss(), N) + (Qr, Qs, Qt), E = sparse_low_order_SBP_operators(rd) + @test norm(sum(Qr, dims=2)) < tol + @test norm(sum(Qs, dims=2)) < tol + @test norm(sum(Qt, dims=2)) < tol + @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E + @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E + @test Qt + Qt' ≈ E' * Diagonal(rd.wf .* rd.ntJ) * E + end +end + +@testset "Subcell limiting operators for approx_type <: SBP" begin + @testset "$elem_type" for elem_type in [Tri()] # [Tri(), Quad(), Hex()] + N = 3 + rd = RefElemData(elem_type, SBP(), N) + Qrst, E = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = subcell_limiting_operators(rd) + for dim in eachindex(Rrst) + @test Δrst[dim] * Rrst[dim] ≈ I + + # check conservation for a random set of limiting factors + r = randn(size(Rrst[dim], 2)) + r = r .- sum(r) / length(r) # so that sum(r) = 0 + θ = rand(size(Rrst[dim], 1)) + sum(Δrst[dim] * Diagonal(θ) * Rrst[dim] * r) + @show d + @test all(d .≈ d[1]) # test that d is a constant vector + end + end +end \ No newline at end of file From babc3a0862d2c6bf3761b48b55f15f4e445e2fbc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 25 Mar 2024 13:51:21 -0500 Subject: [PATCH 03/13] add factor --- src/RefElemData_SBP.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index 35055dad..139eda70 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -247,7 +247,7 @@ end """ - function sparse_low_order_SBP_operators(rd) + function sparse_low_order_SBP_operators(rd; factor=1.01) Constructs sparse low order SBP operators given a `RefElemData`. Returns operators `Qrst..., E ≈ Vf * Pq` that satisfy a generalized @@ -255,7 +255,7 @@ summation-by-parts (GSBP) property: `Q_i + Q_i^T = E' * B_i * E` """ -function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}) where {NDIMS} +function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}; factor=1.01) where {NDIMS} (; Pq, Vf, rstq, wf, nrstJ) = rd # if element is a simplex, convert nodes to an equilateral triangle for symmetry @@ -269,7 +269,7 @@ function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}) where {NDIMS} # heuristic cutoff - we take NDIMS + 1 neighbors, but the smallest distance = 0, # so we need to access the first NDIMS + 2 sorted distance entries. dist = sort(view(D, i, :))[NDIMS + 2] - for j in findall(view(D, i, :) .< 1.01 * dist) + for j in findall(view(D, i, :) .< factor * dist) A[i, j] = 1 end end From 04e139c8c4df3db8e7455cac6a1c15fa395a7af8 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 25 Mar 2024 14:01:56 -0500 Subject: [PATCH 04/13] improve docs and error message --- src/RefElemData_SBP.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index 139eda70..15f85fb4 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -254,6 +254,9 @@ Returns operators `Qrst..., E ≈ Vf * Pq` that satisfy a generalized summation-by-parts (GSBP) property: `Q_i + Q_i^T = E' * B_i * E` + +`factor` is a scaling which determines how close a node must be to +another node to be considered a neighbor. """ function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}; factor=1.01) where {NDIMS} (; Pq, Vf, rstq, wf, nrstJ) = rd @@ -278,7 +281,10 @@ function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}; factor=1.01) whe sorted_eigvals = sort(eigvals(L)) # check that there's only one zero null vector - @assert sorted_eigvals[2] > 100 * eps() + if sorted_eigvals[2] < 10 * size(A, 1) * eps() + error("Warning: the connectivity between nodes yields a graph with " * + "more than one connected component. Try increasing the value of `factor`.") + end E_dense = Vf * Pq E = zeros(size(E_dense)) From 7d794ff6b2f4a348760bbc6189e3fc58db14f78f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 25 Mar 2024 21:24:44 -0500 Subject: [PATCH 05/13] fix a test --- test/sparse_SBP_operator_tests.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/sparse_SBP_operator_tests.jl b/test/sparse_SBP_operator_tests.jl index dff68673..ccd1310b 100644 --- a/test/sparse_SBP_operator_tests.jl +++ b/test/sparse_SBP_operator_tests.jl @@ -60,9 +60,7 @@ end r = randn(size(Rrst[dim], 2)) r = r .- sum(r) / length(r) # so that sum(r) = 0 θ = rand(size(Rrst[dim], 1)) - sum(Δrst[dim] * Diagonal(θ) * Rrst[dim] * r) - @show d - @test all(d .≈ d[1]) # test that d is a constant vector + @test abs(sum(Δrst[dim] * Diagonal(θ) * Rrst[dim] * r)) < 10 * length(r) * eps() end end end \ No newline at end of file From afa5d521fd41ae9280fd414de048420ec00ffdc8 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 25 Mar 2024 21:33:05 -0500 Subject: [PATCH 06/13] comments --- src/RefElemData_SBP.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index 15f85fb4..f8ef4620 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -286,6 +286,8 @@ function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}; factor=1.01) whe "more than one connected component. Try increasing the value of `factor`.") end + # note: this part doesn't do anything for a nodal SBP operator. In that case, E_dense = E = Vf + # and E should just reduce to Vf, which is a face extraction operator. E_dense = Vf * Pq E = zeros(size(E_dense)) for i in axes(E, 1) From 0eca7246392982246c9e8d959d0d4a5c82159001 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 19 Apr 2024 08:06:48 -0500 Subject: [PATCH 07/13] clean up subcell ops --- src/RefElemData_SBP.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index f8ef4620..8fe47bc1 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -419,14 +419,11 @@ function subcell_limiting_operators(Qr::AbstractMatrix; tol = 100 * eps()) end end - d = (ones(size(D, 1))' * D)' + d = vec(ones(size(D, 1))' * D) ids = findall(@. abs(d) > 1e2 * eps()) Z = nullspace(D)[ids, :] A = pinv(D)[ids,:] - X = randn(size(Z, 2), size(D, 1)) - y = randn(length(ids)) - e = ones(size(D, 1)) # compute R via linear algebra m, n = size(A) From 8275a097b46b511a4b94096e358b17b4338fec9f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 21 May 2024 16:14:45 -0500 Subject: [PATCH 08/13] refactor and organization d --- src/RefElemData_SBP.jl | 244 ----------------------------------------- src/StartUpDG.jl | 6 +- src/low_order_sbp.jl | 224 +++++++++++++++++++++++++++++++++++++ src/ref_elem_utils.jl | 24 ++++ 4 files changed, 252 insertions(+), 246 deletions(-) create mode 100644 src/low_order_sbp.jl diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index 14d580f3..14fba196 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -194,249 +194,5 @@ function hybridized_SBP_operators(rd) return Qrsth, VhP, Ph, Vh end -# default to doing nothing -map_nodes_to_symmetric_element(element_type, rst...) = rst -# for triangles, map to an equilateral triangle -function map_nodes_to_symmetric_element(::Tri, r, s) - # biunit right triangular vertices - v1, v2, v3 = SVector{2}.(zip(nodes(Tri(), 1)...)) - denom = (v2[2] - v3[2]) * (v1[1] - v3[1]) + (v3[1] - v2[1]) * (v1[2] - v3[2]) - L1 = @. ((v2[2] - v3[2]) * (r - v3[1]) + (v3[1] - v2[1]) * (s - v3[2])) / denom - L2 = @. ((v3[2] - v1[2]) * (r - v3[1]) + (v1[1] - v3[1]) * (s - v3[2])) / denom - L3 = @. 1 - L1 - L2 - - # equilateral vertices - v1 = SVector{2}(2 * [-.5, -sqrt(3) / 6]) - v2 = SVector{2}(2 * [.5, -sqrt(3)/6]) - v3 = SVector{2}(2 * [0, sqrt(3)/3]) - - x = @. v1[1] * L1 + v2[1] * L2 + v3[1] * L3 - y = @. v1[2] * L1 + v2[2] * L2 + v3[2] * L3 - - return x, y -end - - -""" - function sparse_low_order_SBP_operators(rd; factor=1.01) - -Constructs sparse low order SBP operators given a `RefElemData`. -Returns operators `Qrst..., E ≈ Vf * Pq` that satisfy a generalized -summation-by-parts (GSBP) property: - - `Q_i + Q_i^T = E' * B_i * E` - -`factor` is a scaling which determines how close a node must be to -another node to be considered a neighbor. -""" -function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}; factor=1.01) where {NDIMS} - (; Pq, Vf, rstq, wf, nrstJ) = rd - - # if element is a simplex, convert nodes to an equilateral triangle for symmetry - rstq = map_nodes_to_symmetric_element(rd.element_type, rstq...) - - # build distance and adjacency matrices - D = [norm(SVector{NDIMS}(getindex.(rstq, i)) - SVector{NDIMS}(getindex.(rstq, j))) - for i in eachindex(first(rstq)), j in eachindex(first(rstq))] - A = zeros(Int, length(first(rstq)), length(first(rstq))) - for i in axes(D, 1) - # heuristic cutoff - we take NDIMS + 1 neighbors, but the smallest distance = 0, - # so we need to access the first NDIMS + 2 sorted distance entries. - dist = sort(view(D, i, :))[NDIMS + 2] - for j in findall(view(D, i, :) .< factor * dist) - A[i, j] = 1 - end - end - A = (A + A' .> 0) # symmetrize - L = Diagonal(vec(sum(A, dims=2))) - A - sorted_eigvals = sort(eigvals(L)) - - # check that there's only one zero null vector - if sorted_eigvals[2] < 10 * size(A, 1) * eps() - error("Warning: the connectivity between nodes yields a graph with " * - "more than one connected component. Try increasing the value of `factor`.") - end - - # note: this part doesn't do anything for a nodal SBP operator. In that case, E_dense = E = Vf - # and E should just reduce to Vf, which is a face extraction operator. - E_dense = Vf * Pq - E = zeros(size(E_dense)) - for i in axes(E, 1) - # find all j such that E[i,j] ≥ 0.5, e.g., points which positively contribute to at least half of the - # interpolation. These seem to be associated with volume points "j" that are close to face point "i". - ids = findall(E_dense[i, :] .>= 0.5) - E[i, ids] .= E_dense[i, ids] ./ sum(E_dense[i, ids]) # normalize so sum(E, dims=2) = [1, 1, ...] still. - end - Brst = (nJ -> diagm(wf .* nJ)).(nrstJ) - - # compute potential - e = ones(size(L, 2)) - right_hand_sides = map(B -> 0.5 * sum(E' * B, dims=2), Brst) - psi_augmented = map(b -> [L e; e' 0] \ [b; 0], right_hand_sides) - psi = map(x -> x[1:end-1], psi_augmented) - - # construct sparse skew part - function construct_skew_matrix_from_potential(ψ) - S = zeros(length(ψ), length(ψ)) - for i in axes(S, 1), j in axes(S, 2) - if A[i,j] > 0 - S[i,j] = ψ[j] - ψ[i] - end - end - return S - end - - Srst = construct_skew_matrix_from_potential.(psi) - Qrst = map((S, B) -> S + 0.5 * E' * B * E, Srst, Brst) - return sparse.(Qrst), sparse(E) -end - -function sparse_low_order_SBP_1D_operators(rd::RefElemData) - E = zeros(2, rd.N+1) - E[1, 1] = 1 - E[2, end] = 1 - - # create volume operators - Q = diagm(1 => ones(rd.N), -1 => -ones(rd.N)) - Q[1,1] = -1 - Q[end,end] = 1 - Q = 0.5 * Q - - return (sparse(Q),), sparse(E) -end - -sparse_low_order_SBP_operators(rd::RefElemData{1, Line, <:Union{<:SBP, <:Polynomial{Gauss}}}) = - sparse_low_order_SBP_1D_operators(rd) - -function diagonal_1D_mass_matrix(N, ::SBP) - _, w1D = gauss_lobatto_quad(0, 0, N) - return Diagonal(w1D) -end - -function diagonal_1D_mass_matrix(N, ::Polynomial{Gauss}) - _, w1D = gauss_quad(0, 0, N) - return Diagonal(w1D) -end - -function sparse_low_order_SBP_operators(rd::RefElemData{2, Quad, <:Union{<:SBP, <:Polynomial{Gauss}}}) - (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd) - - # permute face node ordering for the first 2 faces - ids = reshape(1:(rd.N+1) * 2, :, 2) - Er = zeros((rd.N+1) * 2, rd.Np) - Er[vec(ids'), :] .= kron(I(rd.N+1), E1D) - Es = kron(E1D, I(rd.N+1)) - E = vcat(Er, Es) - - M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type) - Qr = kron(M1D, Q1D) - Qs = kron(Q1D, M1D) - - return sparse.((Qr, Qs)), sparse(E) -end - -function sparse_low_order_SBP_operators(rd::RefElemData{3, Hex, <:Union{<:SBP, <:Polynomial{Gauss}}}) - (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd) - - # permute face nodes for the faces in the ±r and ±s directions - ids = reshape(1:(rd.N+1)^2 * 2, rd.N+1, :, 2) - Er, Es, Et = ntuple(_ -> zeros((rd.N+1)^2 * 2, rd.Np), 3) - Er[vec(permutedims(ids, [2, 3, 1])), :] .= kron(I(rd.N+1), E1D, I(rd.N+1)) - Es[vec(permutedims(ids, [3, 2, 1])), :] .= kron(I(rd.N+1), I(rd.N+1), E1D) - Et .= kron(E1D, I(rd.N+1), I(rd.N+1)) - - # create boundary extraction matrix - E = vcat(Er, Es, Et) - - M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type) - Qr = kron(M1D, Q1D, M1D) - Qs = kron(M1D, M1D, Q1D) - Qt = kron(Q1D, M1D, M1D) - - return sparse.((Qr, Qs, Qt)), sparse(E) -end - -# builds subcell operators D * R such that for r such that sum(r) = 0, -# D * Diagonal(θ) * R * r is also diagonal for any choice of θ. This is -# useful in subcell limiting (see, for example -# https://doi.org/10.1016/j.compfluid.2022.105627 for a 1D example) -function subcell_limiting_operators(Qr::AbstractMatrix; tol = 100 * eps()) - Qr = sparse(Qr) - Sr = droptol!(Qr - Qr', tol) - Br = droptol!(Qr + Qr', tol) - - num_interior_fluxes = nnz(Sr) ÷ 2 - num_boundary_indices = nnz(Br) - num_fluxes = num_interior_fluxes + num_boundary_indices - - interior_indices = findall(Sr .> tol) - boundary_indices = findall(abs.(Br) .> tol) - - matrix_indices = zeros(Int, size(Sr)) - for (i, ij) in enumerate(boundary_indices) - matrix_indices[ij] = i - matrix_indices[ij.I[2], ij.I[1]] = i - end - - for (i, ij) in enumerate(interior_indices) - matrix_indices[ij] = i + length(boundary_indices) - matrix_indices[ij.I[2], ij.I[1]] = i + length(boundary_indices) - end - - D = zeros(size(Qr, 2), num_fluxes) - for i in axes(matrix_indices, 1), j in axes(matrix_indices, 2) - if matrix_indices[i, j] !== 0 - D[i, matrix_indices[i, j]] = Qr[i, j] - end - end - - d = vec(ones(size(D, 1))' * D) - ids = findall(@. abs(d) > 1e2 * eps()) - - Z = nullspace(D)[ids, :] - A = pinv(D)[ids,:] - - # compute R via linear algebra - m, n = size(A) - z = size(Z, 2) - C = hcat(kron(I(n), Z), kron(ones(n), I(m))) - sol = pinv(C) * -vec(A) - X = reshape(sol[1:n*z], (z, n)) - R = pinv(D) + nullspace(D) * X - - # check if R satisfies our pseudoinverse condition - @assert norm(D * R - I(size(D,1))) < tol * length(D) - - return D, R -end - -""" - Δrst, Rrst = subcell_limiting_operators(rd::RefElemData) - -Returns tuples of subcell limiting operators Drst = (Δr, Δs, ...) and R = (Rr, Rs, ...) -such that for r where sum(r) = 0, sum(D * Diagonal(θ) * R * r) = 0 for any choice of θ. -These operators are useful for conservative subcell limiting techniques (see -https://doi.org/10.1016/j.compfluid.2022.105627 for an example of such an approach on -tensor product elements). - -Sparse SBP operators used in an intermediate step when buidling these subcell -limiting operators; by default, these operators are constructed using -`sparse_low_order_SBP_operators`. To construct subcell limiting operators for a -general SBP operator, one can use the following: - - Δ, R = subcell_limiting_operators(Q::AbstractMatrix; tol = 100 * eps()) -""" -function subcell_limiting_operators(rd::RefElemData) - Qrst, _ = sparse_low_order_SBP_operators(rd) - Δrst, Rrst = subcell_limiting_operators.(Qrst) - return zip(Δrst, Rrst) -end - -# specialize for single dimension -function subcell_limiting_operators(rd::RefElemData{1}) - Qrst, _ = sparse_low_order_SBP_operators(rd) - Δrst, Rrst = subcell_limiting_operators(Qrst[1]) - return Δrst, Rrst -end diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 7fddd7a6..544737a4 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -24,7 +24,6 @@ using Triangulate: Triangulate, TriangulateIO, triangulate @inline mean(x) = sum(x) / length(x) -# reference element utility functions include("RefElemData.jl") include("RefElemData_polynomial.jl") @@ -38,7 +37,10 @@ export TensorProductWedge include("RefElemData_SBP.jl") export SBP, DefaultSBPType, TensorProductLobatto, Hicken, Kubatko # types for SBP node dispatch export LobattoFaceNodes, LegendreFaceNodes # type parameters for SBP{Kubatko{...}} -export hybridized_SBP_operators, sparse_low_order_SBP_operators +export hybridized_SBP_operators + +include("low_order_sbp.jl") +export sparse_low_order_SBP_operators export subcell_limiting_operators export inverse_trace_constant, face_type diff --git a/src/low_order_sbp.jl b/src/low_order_sbp.jl new file mode 100644 index 00000000..142455c7 --- /dev/null +++ b/src/low_order_sbp.jl @@ -0,0 +1,224 @@ +""" + function sparse_low_order_SBP_operators(rd; factor=1.01) + +Constructs sparse low order SBP operators given a `RefElemData`. +Returns operators `Qrst..., E ≈ Vf * Pq` that satisfy a generalized +summation-by-parts (GSBP) property: + + `Q_i + Q_i^T = E' * B_i * E` + +`factor` is a scaling which determines how close a node must be to +another node to be considered a neighbor. +""" +function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}; factor=1.01) where {NDIMS} + (; Pq, Vf, rstq, wf, nrstJ) = rd + + # if element is a simplex, convert nodes to an equilateral triangle for symmetry + rstq = map_nodes_to_symmetric_element(rd.element_type, rstq...) + + # build distance and adjacency matrices + D = [norm(SVector{NDIMS}(getindex.(rstq, i)) - SVector{NDIMS}(getindex.(rstq, j))) + for i in eachindex(first(rstq)), j in eachindex(first(rstq))] + A = zeros(Int, length(first(rstq)), length(first(rstq))) + for i in axes(D, 1) + # heuristic cutoff - we take NDIMS + 1 neighbors, but the smallest distance = 0, + # so we need to access the first NDIMS + 2 sorted distance entries. + dist = sort(view(D, i, :))[NDIMS + 2] + for j in findall(view(D, i, :) .< factor * dist) + A[i, j] = 1 + end + end + A = (A + A' .> 0) # symmetrize + L = Diagonal(vec(sum(A, dims=2))) - A + sorted_eigvals = sort(eigvals(L)) + + # check that there's only one zero null vector + if sorted_eigvals[2] < 10 * size(A, 1) * eps() + error("Warning: the connectivity between nodes yields a graph with " * + "more than one connected component. Try increasing the value of `factor`.") + end + + # note: this part doesn't do anything for a nodal SBP operator. In that case, E_dense = E = Vf + # and E should just reduce to Vf, which is a face extraction operator. + E_dense = Vf * Pq + E = zeros(size(E_dense)) + for i in axes(E, 1) + # find all j such that E[i,j] ≥ 0.5, e.g., points which positively contribute to at least half of the + # interpolation. These seem to be associated with volume points "j" that are close to face point "i". + ids = findall(E_dense[i, :] .>= 0.5) + E[i, ids] .= E_dense[i, ids] ./ sum(E_dense[i, ids]) # normalize so sum(E, dims=2) = [1, 1, ...] still. + end + Brst = (nJ -> diagm(wf .* nJ)).(nrstJ) + + # compute potential + e = ones(size(L, 2)) + right_hand_sides = map(B -> 0.5 * sum(E' * B, dims=2), Brst) + psi_augmented = map(b -> [L e; e' 0] \ [b; 0], right_hand_sides) + psi = map(x -> x[1:end-1], psi_augmented) + + # construct sparse skew part + function construct_skew_matrix_from_potential(ψ) + S = zeros(length(ψ), length(ψ)) + for i in axes(S, 1), j in axes(S, 2) + if A[i,j] > 0 + S[i,j] = ψ[j] - ψ[i] + end + end + return S + end + + Srst = construct_skew_matrix_from_potential.(psi) + Qrst = map((S, B) -> S + 0.5 * E' * B * E, Srst, Brst) + return sparse.(Qrst), sparse(E) +end + +function sparse_low_order_SBP_1D_operators(rd::RefElemData) + E = zeros(2, rd.N+1) + E[1, 1] = 1 + E[2, end] = 1 + + # create volume operators + Q = diagm(1 => ones(rd.N), -1 => -ones(rd.N)) + Q[1,1] = -1 + Q[end,end] = 1 + Q = 0.5 * Q + + return (sparse(Q),), sparse(E) +end + +sparse_low_order_SBP_operators(rd::RefElemData{1, Line, <:Union{<:SBP, <:Polynomial{Gauss}}}) = + sparse_low_order_SBP_1D_operators(rd) + +function diagonal_1D_mass_matrix(N, ::SBP) + _, w1D = gauss_lobatto_quad(0, 0, N) + return Diagonal(w1D) +end + +function diagonal_1D_mass_matrix(N, ::Polynomial{Gauss}) + _, w1D = gauss_quad(0, 0, N) + return Diagonal(w1D) +end + +function sparse_low_order_SBP_operators(rd::RefElemData{2, Quad, <:Union{<:SBP, <:Polynomial{Gauss}}}) + (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd) + + # permute face node ordering for the first 2 faces + ids = reshape(1:(rd.N+1) * 2, :, 2) + Er = zeros((rd.N+1) * 2, rd.Np) + Er[vec(ids'), :] .= kron(I(rd.N+1), E1D) + Es = kron(E1D, I(rd.N+1)) + E = vcat(Er, Es) + + M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type) + Qr = kron(M1D, Q1D) + Qs = kron(Q1D, M1D) + + return sparse.((Qr, Qs)), sparse(E) +end + +function sparse_low_order_SBP_operators(rd::RefElemData{3, Hex, <:Union{<:SBP, <:Polynomial{Gauss}}}) + (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd) + + # permute face nodes for the faces in the ±r and ±s directions + ids = reshape(1:(rd.N+1)^2 * 2, rd.N+1, :, 2) + Er, Es, Et = ntuple(_ -> zeros((rd.N+1)^2 * 2, rd.Np), 3) + Er[vec(permutedims(ids, [2, 3, 1])), :] .= kron(I(rd.N+1), E1D, I(rd.N+1)) + Es[vec(permutedims(ids, [3, 2, 1])), :] .= kron(I(rd.N+1), I(rd.N+1), E1D) + Et .= kron(E1D, I(rd.N+1), I(rd.N+1)) + + # create boundary extraction matrix + E = vcat(Er, Es, Et) + + M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type) + Qr = kron(M1D, Q1D, M1D) + Qs = kron(M1D, M1D, Q1D) + Qt = kron(Q1D, M1D, M1D) + + return sparse.((Qr, Qs, Qt)), sparse(E) +end + +# builds subcell operators D * R such that for r such that sum(r) = 0, +# D * Diagonal(θ) * R * r is also diagonal for any choice of θ. This is +# useful in subcell limiting (see, for example +# https://doi.org/10.1016/j.compfluid.2022.105627 for a 1D example) +function subcell_limiting_operators(Qr::AbstractMatrix; tol = 100 * eps()) + Qr = sparse(Qr) + Sr = droptol!(Qr - Qr', tol) + Br = droptol!(Qr + Qr', tol) + + num_interior_fluxes = nnz(Sr) ÷ 2 + num_boundary_indices = nnz(Br) + num_fluxes = num_interior_fluxes + num_boundary_indices + + interior_indices = findall(Sr .> tol) + boundary_indices = findall(abs.(Br) .> tol) + + matrix_indices = zeros(Int, size(Sr)) + for (i, ij) in enumerate(boundary_indices) + matrix_indices[ij] = i + matrix_indices[ij.I[2], ij.I[1]] = i + end + + for (i, ij) in enumerate(interior_indices) + matrix_indices[ij] = i + length(boundary_indices) + matrix_indices[ij.I[2], ij.I[1]] = i + length(boundary_indices) + end + + D = zeros(size(Qr, 2), num_fluxes) + for i in axes(matrix_indices, 1), j in axes(matrix_indices, 2) + if matrix_indices[i, j] !== 0 + D[i, matrix_indices[i, j]] = Qr[i, j] + end + end + + d = vec(ones(size(D, 1))' * D) + ids = findall(@. abs(d) > 1e2 * eps()) + + Z = nullspace(D)[ids, :] + A = pinv(D)[ids,:] + + # compute R via linear algebra + m, n = size(A) + z = size(Z, 2) + C = hcat(kron(I(n), Z), kron(ones(n), I(m))) + sol = pinv(C) * -vec(A) + X = reshape(sol[1:n*z], (z, n)) + R = pinv(D) + nullspace(D) * X + + # check if R satisfies our pseudoinverse condition + @assert norm(D * R - I(size(D,1))) < tol * length(D) + + return D, R +end + +""" + Δrst, Rrst = subcell_limiting_operators(rd::RefElemData) + +Returns tuples of subcell limiting operators Drst = (Δr, Δs, ...) and R = (Rr, Rs, ...) +such that for r where sum(r) = 0, sum(D * Diagonal(θ) * R * r) = 0 for any choice of θ. +These operators are useful for conservative subcell limiting techniques (see +https://doi.org/10.1016/j.compfluid.2022.105627 for an example of such an approach on +tensor product elements). + +Sparse SBP operators used in an intermediate step when buidling these subcell +limiting operators; by default, these operators are constructed using +`sparse_low_order_SBP_operators`. To construct subcell limiting operators for a +general SBP operator, one can use the following: + + Δ, R = subcell_limiting_operators(Q::AbstractMatrix; tol = 100 * eps()) +""" +function subcell_limiting_operators(rd::RefElemData) + Qrst, _ = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = subcell_limiting_operators.(Qrst) + return zip(Δrst, Rrst) +end + +# specialize for single dimension +function subcell_limiting_operators(rd::RefElemData{1}) + Qrst, _ = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = first(subcell_limiting_operators.(Qrst)) + return (Δrst,), (Rrst,) +end + +# TODO: specialize for quadrilaterals and hex? + diff --git a/src/ref_elem_utils.jl b/src/ref_elem_utils.jl index df9cab16..1fa75801 100644 --- a/src/ref_elem_utils.jl +++ b/src/ref_elem_utils.jl @@ -86,6 +86,30 @@ function init_face_data(elem::Tet; quad_rule_face=quad_nodes(Tri(), N)) return rf, sf, tf, wf, nrJ, nsJ, ntJ end +# default to doing nothing +map_nodes_to_symmetric_element(element_type, rst...) = rst + +# for triangles, map to an equilateral triangle +function map_nodes_to_symmetric_element(::Tri, r, s) + # biunit right triangular vertices + v1, v2, v3 = SVector{2}.(zip(nodes(Tri(), 1)...)) + + denom = (v2[2] - v3[2]) * (v1[1] - v3[1]) + (v3[1] - v2[1]) * (v1[2] - v3[2]) + L1 = @. ((v2[2] - v3[2]) * (r - v3[1]) + (v3[1] - v2[1]) * (s - v3[2])) / denom + L2 = @. ((v3[2] - v1[2]) * (r - v3[1]) + (v1[1] - v3[1]) * (s - v3[2])) / denom + L3 = @. 1 - L1 - L2 + + # equilateral vertices + v1 = SVector{2}(2 * [-.5, -sqrt(3) / 6]) + v2 = SVector{2}(2 * [.5, -sqrt(3)/6]) + v3 = SVector{2}(2 * [0, sqrt(3)/3]) + + x = @. v1[1] * L1 + v2[1] * L2 + v3[1] * L3 + y = @. v1[2] * L1 + v2[2] * L2 + v3[2] * L3 + + return x, y +end + """ function inverse_trace_constant(rd::RefElemData) From e1b425eab10177457ac9fdd9986d78ec9ae964d9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 21 May 2024 21:37:10 -0500 Subject: [PATCH 09/13] improve comments --- src/low_order_sbp.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/low_order_sbp.jl b/src/low_order_sbp.jl index 142455c7..7e0d1b2f 100644 --- a/src/low_order_sbp.jl +++ b/src/low_order_sbp.jl @@ -213,12 +213,12 @@ function subcell_limiting_operators(rd::RefElemData) return zip(Δrst, Rrst) end -# specialize for single dimension +# specialize for single dimension to return tuples of operators function subcell_limiting_operators(rd::RefElemData{1}) Qrst, _ = sparse_low_order_SBP_operators(rd) Δrst, Rrst = first(subcell_limiting_operators.(Qrst)) return (Δrst,), (Rrst,) -end +end # TODO: specialize for quadrilaterals and hex? From cfd695719fb91a0e9e1c9e13016c2b15c29768fc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 21 May 2024 21:37:14 -0500 Subject: [PATCH 10/13] load nnz --- src/StartUpDG.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 0ed77702..59664a0b 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -17,7 +17,7 @@ using RecipesBase: RecipesBase @reexport using RecursiveArrayTools: NamedArrayPartition using StaticArrays: SVector, SMatrix using Setfield: setproperties, @set # for "modifying" structs (setproperties) -using SparseArrays: sparse, droptol!, blockdiag +using SparseArrays: sparse, droptol!, blockdiag, nnz using Triangulate: Triangulate, TriangulateIO, triangulate @reexport using WriteVTK From 678e61c225e71ffd5c81ed80c92c7f92198520e4 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 21 May 2024 21:37:24 -0500 Subject: [PATCH 11/13] add tensor product subcell operator tests --- test/sparse_SBP_operator_tests.jl | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/test/sparse_SBP_operator_tests.jl b/test/sparse_SBP_operator_tests.jl index ccd1310b..758cbb72 100644 --- a/test/sparse_SBP_operator_tests.jl +++ b/test/sparse_SBP_operator_tests.jl @@ -48,7 +48,7 @@ end @testset "Subcell limiting operators for approx_type <: SBP" begin - @testset "$elem_type" for elem_type in [Tri()] # [Tri(), Quad(), Hex()] + @testset "$elem_type" for elem_type in [Tri()] N = 3 rd = RefElemData(elem_type, SBP(), N) Qrst, E = sparse_low_order_SBP_operators(rd) @@ -63,4 +63,26 @@ end @test abs(sum(Δrst[dim] * Diagonal(θ) * Rrst[dim] * r)) < 10 * length(r) * eps() end end + + @testset "Tensor product elements: $elem_type" for elem_type in [Quad(), Hex()] + N = 2 + rd = RefElemData(elem_type, SBP(), N) + Qrst, E = sparse_low_order_SBP_operators(rd) + Δrst, Rrst = subcell_limiting_operators(rd) + + for dim in eachindex(Rrst) + conservation_vectors = nullspace(Matrix(Qrst[dim])) + Δ = Δrst[dim] + R = Rrst[dim] + + # make random conservative vector + r = randn(length(rd.r)) + r = r - conservation_vectors * (conservation_vectors' * r) + + # create limiting factors + θ = rand(size(R, 1)) + + @test norm(conservation_vectors' * (Δ * (θ .* (R * r)))) < 100 * eps() + end + end end \ No newline at end of file From f3879293740f4c26e778fca130e3266808a9e9ac Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 21 May 2024 21:39:09 -0500 Subject: [PATCH 12/13] reactivate all tests --- test/runtests.jl | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 9be78b5f..63caac7f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,18 +8,17 @@ 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("multidim_sbp_tests.jl") +include("write_vtk_tests.jl") +include("triangulate_tests.jl") +include("reference_elem_tests.jl") +include("multidim_sbp_tests.jl") include("sparse_SBP_operator_tests.jl") -# include("SummationByPartsOperatorsExt_tests.jl") -# include("MeshData_tests.jl") -# include("boundary_util_tests.jl") -# include("hybrid_mesh_tests.jl") -# include("noncon_mesh_tests.jl") -# include("cut_mesh_tests.jl") -# include("misc_tests.jl") -# include("gmsh_parse_tests.jl") -# include("hohqmesh_tests.jl") +include("SummationByPartsOperatorsExt_tests.jl") +include("MeshData_tests.jl") +include("boundary_util_tests.jl") +include("hybrid_mesh_tests.jl") +include("noncon_mesh_tests.jl") +include("cut_mesh_tests.jl") +include("misc_tests.jl") +include("gmsh_parse_tests.jl") +include("hohqmesh_tests.jl") From c2e3654572caa877246c9bd47bbcc4d0b7bf0584 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 21 May 2024 21:52:10 -0500 Subject: [PATCH 13/13] add dropped wedge/pyr tests back --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 63caac7f..cfdeabbd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,6 +15,7 @@ include("multidim_sbp_tests.jl") include("sparse_SBP_operator_tests.jl") include("SummationByPartsOperatorsExt_tests.jl") include("MeshData_tests.jl") +include("MeshData_wedge_pyr_tests.jl") include("boundary_util_tests.jl") include("hybrid_mesh_tests.jl") include("noncon_mesh_tests.jl")