Improve multidimensional matrix operators (#330)

* improve multidimensional matrix operators

* format

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* remove kron specialization again

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

Author: JoshuaLampert

src/SummationByPartsOperators.jl

@@ -152,7 +152,8 @@ export FilterCallback, ConstantFilter, ExponentialFilter
 export SafeMode, FastMode, ThreadedMode
 export derivative_order, accuracy_order, source_of_coefficients, grid, semidiscretize
 export mass_matrix, mass_matrix_boundary
-export integrate, integrate_boundary, restrict_boundary,
+export integrate, integrate_boundary,
+       restrict_interior, restrict_boundary,
        left_boundary_weight, right_boundary_weight,
        scale_by_mass_matrix!, scale_by_inverse_mass_matrix!,
        derivative_left, derivative_right,
@@ -165,6 +166,7 @@ export periodic_central_derivative_operator, periodic_derivative_operator,
        fourier_derivative_operator,
        legendre_derivative_operator, legendre_second_derivative_operator,
        upwind_operators, function_space_operator, tensor_product_operator_2D
+export normals, boundary_indices
 export UniformMesh1D, UniformPeriodicMesh1D
 export couple_continuously, couple_discontinuously
 export mul!

src/general_operators.jl

@@ -304,6 +304,15 @@ restrict_boundary(u, D::AbstractNonperiodicDerivativeOperator) = u[[begin, end]]
 
 restrict_boundary(u, D::AbstractPeriodicDerivativeOperator) = eltype(u)[]
 
+"""
+    restrict_interior(u, D::AbstractDerivativeOperator)
+
+Restrict the coefficients `u` to the interior nodes of the derivative operator `D`.
+"""
+restrict_interior(u, D::AbstractNonperiodicDerivativeOperator) = u[(begin + 1):(end - 1)]
+
+restrict_interior(u, D::AbstractPeriodicDerivativeOperator) = u
+
 """
     mass_matrix_boundary(D::AbstractDerivativeOperator)
 

src/multidimensional_matrix_operators.jl

@@ -74,10 +74,28 @@ Base.ndims(::AbstractMultidimensionalMatrixDerivativeOperator{Dim}) where {Dim}
 Base.getindex(D::AbstractMultidimensionalMatrixDerivativeOperator, i::Int) = D.Ds[i]
 Base.eltype(::AbstractMultidimensionalMatrixDerivativeOperator{Dim, T}) where {Dim, T} = T
 
+"""
+    normals(D::MultidimensionalMatrixDerivativeOperator)
+
+Return the normal vectors of the boundary nodes of a [`MultidimensionalMatrixDerivativeOperator`](@ref) `D`.
+"""
+normals(D::AbstractMultidimensionalMatrixDerivativeOperator) = D.normals
+
+"""
+    boundary_indices(D::MultidimensionalMatrixDerivativeOperator)
+
+Return the indices of the boundary nodes of a [`MultidimensionalMatrixDerivativeOperator`](@ref) `D`.
+"""
+boundary_indices(D::AbstractMultidimensionalMatrixDerivativeOperator) = D.boundary_indices
+
 function restrict_boundary(u, D::AbstractMultidimensionalMatrixDerivativeOperator)
     u[D.boundary_indices]
 end
 
+function restrict_interior(u, D::AbstractMultidimensionalMatrixDerivativeOperator)
+    u[setdiff(eachindex(u), D.boundary_indices)]
+end
+
 """
     integrate_boundary([func = identity,] u, D::MultidimensionalMatrixDerivativeOperator, dim)
 

test/SBP_operators_test.jl

@@ -45,6 +45,9 @@ for source in D_test_list, T in (Float32, Float64)
         # SBP property
         M = mass_matrix(D)
         @test M * Matrix(D) + Matrix(D)' * M ≈ mass_matrix_boundary(D)
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 1000 * eps(T), eachindex(res))
@@ -112,6 +115,9 @@ for source in D_test_list, T in (Float32, Float64)
         # SBP property
         M = mass_matrix(D)
         @test M * Matrix(D) + Matrix(D)' * M ≈ mass_matrix_boundary(D)
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 1000 * eps(T), eachindex(res))
@@ -189,6 +195,9 @@ for source in D_test_list, T in (Float32, Float64)
         # SBP property
         M = mass_matrix(D)
         @test M * Matrix(D) + Matrix(D)' * M ≈ mass_matrix_boundary(D)
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 1000 * eps(T), eachindex(res))
@@ -276,6 +285,9 @@ for source in D_test_list, T in (Float32, Float64)
         # SBP property
         M = mass_matrix(D)
         @test M * Matrix(D) + Matrix(D)' * M ≈ mass_matrix_boundary(D)
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 16000 * eps(T), eachindex(res))
@@ -438,6 +450,9 @@ for source in D_test_list, T in (Float32, Float64)
         dL = derivative_left(D, Val{1}())
         dR = derivative_right(D, Val{1}())
         @test M * Matrix(D) - Matrix(D)' * M ≈ eR * dR' - eL * dL' - dR * eR' + dL * eL'
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 2000 * eps(T), eachindex(res))
@@ -509,6 +524,9 @@ for source in D_test_list, T in (Float32, Float64)
         dL = derivative_left(D, Val{1}())
         dR = derivative_right(D, Val{1}())
         @test M * Matrix(D) - Matrix(D)' * M ≈ eR * dR' - eL * dL' - dR * eR' + dL * eL'
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 10000 * eps(T), eachindex(res))
@@ -587,6 +605,9 @@ for source in D_test_list, T in (Float32, Float64)
         dL = derivative_left(D, Val{1}())
         dR = derivative_right(D, Val{1}())
         @test M * Matrix(D) - Matrix(D)' * M ≈ eR * dR' - eL * dL' - dR * eR' + dL * eL'
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 10000 * eps(T), eachindex(res))
@@ -669,6 +690,9 @@ end
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == false
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < eps(T), eachindex(res))
@@ -730,6 +754,9 @@ end
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == false
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 100_000 * eps(T), eachindex(res))
@@ -800,6 +827,9 @@ end
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == false
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 1_000_000 * eps(T), eachindex(res))
@@ -884,6 +914,9 @@ end
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == false
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 10 * eps(T) / D.Δx^4, eachindex(res))
@@ -953,6 +986,9 @@ end
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == false
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 50_000_000 * eps(T), eachindex(res))
@@ -1033,6 +1069,9 @@ end
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == false
+
+        @test restrict_interior(x2, D) ≈ x2[(begin + 1):(end - 1)]
+        @test restrict_boundary(x2, D) ≈ [xmin^2, xmax^2]
         # interior and boundary
         mul!(res, D, x0)
         @test all(i -> abs(res[i]) < 50_000_000 * eps(T), eachindex(res))

test/multidimensional_matrix_operators_test.jl

@@ -121,8 +121,21 @@ end
         end
         @test eltype(D_t) == eltype(D_1) == eltype(D_2) == T
         @test real(D_t) == real(D_1) == real(D_2) == T
+        @test normals(D_t) == D_t.normals
+        @test boundary_indices(D_t) == D_t.boundary_indices
 
-        @test length(grid(D_t)) == N_x * N_y
+        nodes = grid(D_t)
+        @test length(nodes) == N_x * N_y
+        u = sin.(first.(nodes) .^ 2) .* cospi.(last.(nodes))
+        u_reference = copy(u)
+
+        SummationByPartsOperators.scale_by_mass_matrix!(u, D_t)
+        SummationByPartsOperators.scale_by_inverse_mass_matrix!(u, D_t)
+        @test u ≈ u_reference
+        @test length(restrict_interior(u, D_t)) == (N_x - 2) * (N_y - 2)
+        @test length(restrict_boundary(u, D_t)) == 2 * (N_x + N_y)
+
+        # derivative operators
         M = mass_matrix(D_t)
         D_x = D_t[1]
         @test D_x isa SparseMatrixCSC
@@ -148,6 +161,11 @@ end
         @test B_x ≈ Diagonal(kron(B_1D_1, M_1D_2))
         @test B_y ≈ Diagonal(kron(M_1D_1, B_1D_2))
 
+        # integrals
+        @test integrate(abs2, u, D_t) ≈ sum(M * abs2.(u))
+        @test integrate_boundary(abs2, u, D_t, 1) ≈ sum(B_x * abs2.(u))
+        @test integrate_boundary(abs2, u, D_t, 2) ≈ sum(B_y * abs2.(u))
+
         # accuracy test
         x = grid(D_t)
         x0 = ones(N_x * N_y)