From 8d21c88ada19b2b0574bd9204f0ce97ca0fca6c0 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Tue, 21 May 2024 22:03:57 -0500 Subject: [PATCH] Add subcell finite volume differencing and reconstruction operators (#163) * add subcell limiting operators * refactor tests * add factor * improve docs and error message * fix a test * comments * clean up subcell ops * refactor and organization d * improve comments * load nnz * add tensor product subcell operator tests * reactivate all tests * add dropped wedge/pyr tests back --- src/RefElemData_SBP.jl | 153 -------------------- src/StartUpDG.jl | 11 +- src/low_order_sbp.jl | 224 ++++++++++++++++++++++++++++++ src/ref_elem_utils.jl | 24 ++++ test/reference_elem_tests.jl | 49 ------- test/runtests.jl | 1 + test/sparse_SBP_operator_tests.jl | 88 ++++++++++++ 7 files changed, 344 insertions(+), 206 deletions(-) create mode 100644 src/low_order_sbp.jl create mode 100644 test/sparse_SBP_operator_tests.jl diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl index b335ca75..14fba196 100644 --- a/src/RefElemData_SBP.jl +++ b/src/RefElemData_SBP.jl @@ -194,158 +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) - -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` -""" -function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}) 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, :) .< 1.01 * 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 - @assert sorted_eigvals[2] > 100 * eps() - - 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 diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index e7ad42df..59664a0b 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,13 +17,12 @@ 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 @inline mean(x) = sum(x) / length(x) -# reference element utility functions include("RefElemData.jl") include("RefElemData_polynomial.jl") @@ -37,7 +36,11 @@ 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 include("ref_elem_utils.jl") diff --git a/src/low_order_sbp.jl b/src/low_order_sbp.jl new file mode 100644 index 00000000..7e0d1b2f --- /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 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 + +# 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) diff --git a/test/reference_elem_tests.jl b/test/reference_elem_tests.jl index 66ad58f8..46f7a419 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 91d32de4..cfdeabbd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,6 +12,7 @@ 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("MeshData_wedge_pyr_tests.jl") diff --git a/test/sparse_SBP_operator_tests.jl b/test/sparse_SBP_operator_tests.jl new file mode 100644 index 00000000..758cbb72 --- /dev/null +++ b/test/sparse_SBP_operator_tests.jl @@ -0,0 +1,88 @@ +@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()] + 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)) + @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