add tests for SBP property (#272)

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
ranocha · Jun 16, 2024 · 5402b15 · 5402b15
1 parent f09cbac
commit 5402b15
Show file tree

Hide file tree

Showing 5 changed files with 94 additions and 4 deletions.
diff --git a/test/SBP_operators_test.jl b/test/SBP_operators_test.jl
@@ -16,6 +16,9 @@ for source in D_test_list, T in (Float32,Float64)
     xmin = -one(T)
     xmax = 2*one(T)
     N = 101
+    B = zeros(T, N, N)
+    B[1, 1] = convert(T, -1)
+    B[end, end] = convert(T, 1)
     der_order = 1
 
     acc_order = 2
@@ -42,6 +45,9 @@ for source in D_test_list, T in (Float32,Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        # SBP property
+        M = mass_matrix(D)
+        @test M * Matrix(D) + Matrix(D)' * M ≈ B
         # interior and boundary
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 1000*eps(T), eachindex(res))
@@ -88,6 +94,9 @@ for source in D_test_list, T in (Float32,Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        # SBP property
+        M = mass_matrix(D)
+        @test M * Matrix(D) + Matrix(D)' * M ≈ B
         # interior and boundary
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 1000*eps(T), eachindex(res))
@@ -142,6 +151,9 @@ for source in D_test_list, T in (Float32,Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        # SBP property
+        M = mass_matrix(D)
+        @test M * Matrix(D) + Matrix(D)' * M ≈ B
         # interior and boundary
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 1000*eps(T), eachindex(res))
@@ -202,6 +214,9 @@ for source in D_test_list, T in (Float32,Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        # SBP property
+        M = mass_matrix(D)
+        @test M * Matrix(D) + Matrix(D)' * M ≈ B
         # interior and boundary
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 16000*eps(T), eachindex(res))
@@ -241,6 +256,8 @@ for source in D_test_list, T in (Float32,Float64)
     xmin = -one(T)
     xmax = 2*one(T)
     N = 101
+    eL = SummationByPartsOperators.OneHot(1, N)
+    eR = SummationByPartsOperators.OneHot(N, N)
     der_order = 2
 
     acc_order = 2
@@ -267,6 +284,11 @@ for source in D_test_list, T in (Float32,Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        # SBP property
+        M = mass_matrix(D)
+        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'
         # interior and boundary
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < eps(T), eachindex(res))
@@ -315,6 +337,11 @@ for source in D_test_list, T in (Float32,Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        # SBP property
+        M = mass_matrix(D)
+        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'
         # interior and boundary
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 2000*eps(T), eachindex(res))
@@ -369,6 +396,11 @@ for source in D_test_list, T in (Float32,Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        # SBP property
+        M = mass_matrix(D)
+        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'
         # interior and boundary
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 10000*eps(T), eachindex(res))
@@ -428,6 +460,11 @@ for source in D_test_list, T in (Float32,Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        # SBP property
+        M = mass_matrix(D)
+        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'
         # interior and boundary
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 10000*eps(T), eachindex(res))

diff --git a/test/fourier_operators_test.jl b/test/fourier_operators_test.jl
@@ -21,6 +21,8 @@ for T in (Float32, Float64)
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        M = mass_matrix(D)
+        @test isapprox(M * Matrix(D) + Matrix(D)' * M, zeros(T, N, N), atol = N*eps(T))
         u = compute_coefficients(zero, D)
         res = D*u
         for k in 0:(N÷2)-1

diff --git a/test/legendre_operators_test.jl b/test/legendre_operators_test.jl
@@ -16,13 +16,18 @@ for T in (Float32, Float64)
     xmax = one(T)
 
     for N in 2 .^ (1:4)
+        B = zeros(T, N, N)
+        B[1, 1] = convert(T, -1)
+        B[end, end] = convert(T, 1)
         D = legendre_derivative_operator(xmin, xmax, N)
         @test D == legendre_derivative_operator(xmin=xmin, xmax=xmax, N=N)
         println(devnull, D)
         @test SummationByPartsOperators.derivative_order(D) == 1
         @test issymmetric(D) == false
         @test SummationByPartsOperators.xmin(D) ≈ xmin
         @test SummationByPartsOperators.xmax(D) ≈ xmax
+        M = mass_matrix(D)
+        @test M * Matrix(D) + Matrix(D)' * M ≈ B
         u = compute_coefficients(zero, D)
         res = D*u
         for k in 1:N-1

diff --git a/test/periodic_operators_test.jl b/test/periodic_operators_test.jl
@@ -31,6 +31,8 @@ let T=Float32
     @test issymmetric(D) == false
     @test SummationByPartsOperators.xmin(D) ≈ xmin
     @test SummationByPartsOperators.xmax(D) ≈ xmax
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0)
     @test all(i->abs(res[i]) < eps(T), accuracy_order:length(res)-accuracy_order)
     mul!(res, D, x1)
@@ -48,6 +50,8 @@ let T=Float32
     @test issymmetric(D) == false
     @test SummationByPartsOperators.xmin(D) ≈ xmin
     @test SummationByPartsOperators.xmax(D) ≈ xmax
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0)
     @test all(i->abs(res[i]) < 50*eps(T), accuracy_order:length(res)-accuracy_order)
     mul!(res, D, x1)
@@ -69,6 +73,8 @@ let T=Float32
     @test issymmetric(D) == false
     @test SummationByPartsOperators.xmin(D) ≈ xmin
     @test SummationByPartsOperators.xmax(D) ≈ xmax
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0)
     @test all(i->abs(res[i]) < 50*eps(T), accuracy_order:length(res)-accuracy_order)
     mul!(res, D, x1)
@@ -96,6 +102,8 @@ let T=Float32
     @test issymmetric(D) == true
     @test SummationByPartsOperators.xmin(D) ≈ xmin
     @test SummationByPartsOperators.xmax(D) ≈ xmax
+    M = mass_matrix(D)
+    @test M * Matrix(D) - Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0)
     @test all(i->abs(res[i]) < eps(T), accuracy_order:length(res)-accuracy_order)
     mul!(res, D, x1)
@@ -115,6 +123,8 @@ let T=Float32
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == true
+    M = mass_matrix(D)
+    @test M * Matrix(D) - Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0)
     @test all(i->abs(res[i]) < 50N*eps(T), accuracy_order:length(res)-accuracy_order)
     mul!(res, D, x1)
@@ -136,6 +146,8 @@ let T=Float32
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == true
+    M = mass_matrix(D)
+    @test M * Matrix(D) - Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 50N*eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1000N*eps(T)
     mul!(res, D, x2); @test norm((res - 2 .* x0)[accuracy_order:end-accuracy_order], Inf) < 5000N*eps(T)
@@ -159,6 +171,8 @@ let T=Float32
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == true # because this operator is zero!
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1400*eps(T)
     mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 7000*eps(T)
@@ -172,6 +186,8 @@ let T=Float32
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == false
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 400*eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 10000N*eps(T)
     mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 50000N*eps(T)
@@ -187,6 +203,8 @@ let T=Float32
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == false
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 300*eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 80000N*eps(T)
     mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 300000N*eps(T)
@@ -237,6 +255,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == false
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0)
     @test all(i->abs(res[i]) < eps(T), accuracy_order:length(res)-accuracy_order)
     mul!(res, D, x1)
@@ -252,9 +272,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == false
-    accuracy_order = 4
-    D = periodic_central_derivative_operator(1, accuracy_order, xmin, xmax, N)
-    @test issymmetric(D) == false
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0)
     @test all(i->abs(res[i]) < 10*eps(T), accuracy_order:length(res)-accuracy_order)
     mul!(res, D, x1)
@@ -274,6 +293,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == false
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0)
     @test all(i->abs(res[i]) < 10*eps(T), accuracy_order:length(res)-accuracy_order)
     mul!(res, D, x1)
@@ -303,6 +324,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == true
+    M = mass_matrix(D)
+    @test M * Matrix(D) - Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1400*eps(T)
     mul!(res, D, x2); @test norm((res - 2 .* x0)[accuracy_order:end-accuracy_order], Inf) < 7000*eps(T)
@@ -314,6 +337,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == true
+    M = mass_matrix(D)
+    @test M * Matrix(D) - Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 400*eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 5000*eps(T)
     mul!(res, D, x2); @test norm((res - 2 .* x0)[accuracy_order:end-accuracy_order], Inf) < 11000*eps(T)
@@ -327,6 +352,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == true
+    M = mass_matrix(D)
+    @test M * Matrix(D) - Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1000*eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 5000*eps(T)
     mul!(res, D, x2); @test norm((res - 2 .* x0)[accuracy_order:end-accuracy_order], Inf) < 10020*eps(T)
@@ -348,6 +375,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == true # because this operator is zero!
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 1400*eps(T)
     mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 7000*eps(T)
@@ -359,6 +388,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == false
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 400*eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 20000*eps(T)
     mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 200000*eps(T)
@@ -372,6 +403,8 @@ let T = Float64
     @test SummationByPartsOperators.derivative_order(D) == derivative_order
     @test SummationByPartsOperators.accuracy_order(D) == accuracy_order
     @test issymmetric(D) == false
+    M = mass_matrix(D)
+    @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
     mul!(res, D, x0); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 300*eps(T)
     mul!(res, D, x1); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 70000*eps(T)
     mul!(res, D, x2); @test norm(res[accuracy_order:end-accuracy_order], Inf) < 200000*eps(T)
@@ -718,6 +751,8 @@ let T = Float32
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == false
+        M = mass_matrix(D)
+        @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < eps(T), stencil_width:length(res)-stencil_width)
         mul!(res, D, x1)
@@ -738,6 +773,8 @@ let T = Float32
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == false
+        M = mass_matrix(D)
+        @test M * Matrix(D) + Matrix(D)' * M ≈ zeros(T, N, N)
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 50*eps(T),  stencil_width:length(res)-stencil_width)
         mul!(res, D, x1)
@@ -770,6 +807,8 @@ let T = Float32
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == true
+        M = mass_matrix(D)
+        @test M * Matrix(D) - Matrix(D)' * M ≈ zeros(T, N, N)
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < eps(T),  stencil_width:length(res)-stencil_width)
         mul!(res, D, x1)
@@ -792,6 +831,8 @@ let T = Float32
         @test derivative_order(D) == der_order
         @test accuracy_order(D) == acc_order
         @test issymmetric(D) == true
+        M = mass_matrix(D)
+        @test M * Matrix(D) - Matrix(D)' * M ≈ zeros(T, N, N)
         mul!(res, D, x0)
         @test all(i->abs(res[i]) < 50N*eps(T),  stencil_width:length(res)-stencil_width)
         mul!(res, D, x1)

diff --git a/test/upwind_operators_test.jl b/test/upwind_operators_test.jl
@@ -63,7 +63,6 @@ end
   N = 21
   xmin = 0.
   xmax = Float64(N + 1)
-  interior = 10:N-10
   acc_order = 2
 
   D = upwind_operators(Mattsson2017, derivative_order=1, accuracy_order=acc_order,
@@ -97,6 +96,12 @@ end
   x = grid(D)
   @test integrate(x, D) == integrate(x, Dp)
 
+  B = zeros(N, N)
+  B[1, 1] = -1.0
+  B[end, end] = 1.0
+  M = mass_matrix(D)
+  @test M * Matrix(Dp) + Matrix(Dm)' * M ≈ B
+
   @test_throws ArgumentError UpwindOperators(
     derivative_operator(Mattsson2017(:minus  ), 1, acc_order, xmin, xmax, N),
     derivative_operator(Mattsson2017(:central), 1, acc_order, xmin, xmax, N+1),