WIP: first-order upwind SBP coefficients

src/SBP_coefficients/Mattsson2017.jl

@@ -35,7 +35,65 @@ end
 
 function first_derivative_coefficients(source::Mattsson2017, order::Int,
                                        T=Float64, mode=FastMode())
-    if order == 2
+    if order == 1
+        left_boundary_plus = (
+            DerivativeCoefficientRow{T,1,2}(SVector(T(-2),    # -2 for continuously coupled = same
+                                                    T(2), )), #  2 for continuously coupled = same
+        )
+        right_boundary_plus = (
+            DerivativeCoefficientRow{T,1,2}(SVector(T(0),
+                                                    T(0), )),
+        )
+        upper_coef_plus = SVector( T(1), )
+        central_coef_plus = T(-1)
+        lower_coef_plus = SVector{0,T}()
+        left_weights = SVector( T(1//2), ) # 1//2 for continuously coupled = same
+        right_weights = left_weights
+        left_boundary_derivatives = Tuple{}()
+        right_boundary_derivatives = left_boundary_derivatives
+
+        left_boundary_minus = (
+            DerivativeCoefficientRow{T,1,2}(SVector(T(0),
+                                                    T(0), )),
+        )
+        right_boundary_minus = (
+            DerivativeCoefficientRow{T,1,2}(SVector(T(2),      #  2 for continuously coupled = same
+                                                    T(-2), )), # -2 for continuously coupled = same
+        )
+        upper_coef_minus     = .- lower_coef_plus
+        central_coef_minus   = .- central_coef_plus
+        lower_coef_minus     = .- upper_coef_plus
+
+        left_boundary_central  = (left_boundary_plus  .+ left_boundary_minus)  ./ 2
+        right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2
+        upper_coef_central     = widening_plus(upper_coef_plus, upper_coef_minus) / 2
+        central_coef_central   = (central_coef_plus   + central_coef_minus) / 2
+        lower_coef_central     = widening_plus(lower_coef_plus, lower_coef_minus) / 2
+
+        if source.kind === :plus
+          left_boundary  = left_boundary_plus
+          right_boundary = right_boundary_plus
+          upper_coef     = upper_coef_plus
+          central_coef   = central_coef_plus
+          lower_coef     = lower_coef_plus
+        elseif source.kind === :minus
+          left_boundary  = left_boundary_minus
+          right_boundary = right_boundary_minus
+          upper_coef     = upper_coef_minus
+          central_coef   = central_coef_minus
+          lower_coef     = lower_coef_minus
+        elseif source.kind === :central
+          left_boundary  = left_boundary_central
+          right_boundary = right_boundary_central
+          upper_coef     = upper_coef_central
+          central_coef   = central_coef_central
+          lower_coef     = lower_coef_central
+        end
+        DerivativeCoefficients(left_boundary, right_boundary,
+                               left_boundary_derivatives, right_boundary_derivatives,
+                               lower_coef, central_coef, upper_coef,
+                               left_weights, right_weights, mode, 1, order, source)
+    elseif order == 2
         left_boundary_plus = (
             DerivativeCoefficientRow{T,1,3}(SVector(T(-3),
                                                     T(5),

test/upwind_operators_test.jl

@@ -2,52 +2,67 @@ using Test
 using LinearAlgebra
 using SummationByPartsOperators
 
+source_list = (Mattsson2017,)
+
 # check construction of interior part of upwind operators
 @testset "Check interior parts" begin
   N = 21
-  xmin = 0.
+  xmin = 0.0
   xmax = Float64(N + 1)
   interior = 10:N-10
 
+  for Source in source_list
+    for acc_order in 1:9
+      Dp_bounded, Dm_bounded, Dc_bounded = try
+        Dp_bounded = derivative_operator(Source(:plus   ), 1, acc_order, xmin, xmax, N)
+        Dm_bounded = derivative_operator(Source(:minus  ), 1, acc_order, xmin, xmax, N)
+        Dc_bounded = derivative_operator(Source(:central), 1, acc_order, xmin, xmax, N)
+        Dp_bounded, Dm_bounded, Dc_bounded
+      catch
+        Dp_bounded = nothing
+        Dm_bounded = nothing
+        Dc_bounded = nothing
+        Dp_bounded, Dm_bounded, Dc_bounded
+      end
+
+      Dp_bounded === nothing && continue
 
-  for acc_order in 2:9
-    Dp_bounded = derivative_operator(Mattsson2017(:plus   ), 1, acc_order, xmin, xmax, N)
-    Dm_bounded = derivative_operator(Mattsson2017(:minus  ), 1, acc_order, xmin, xmax, N)
-    Dc_bounded = derivative_operator(Mattsson2017(:central), 1, acc_order, xmin, xmax, N)
-    for compact in (true, false)
-      show(IOContext(devnull, :compact=>compact), Dp_bounded)
-      show(IOContext(devnull, :compact=>compact), Dm_bounded)
-      show(IOContext(devnull, :compact=>compact), Dc_bounded)
-      summary(IOContext(devnull, :compact=>compact), Dp_bounded)
-      summary(IOContext(devnull, :compact=>compact), Dm_bounded)
-      summary(IOContext(devnull, :compact=>compact), Dc_bounded)
-    end
-    M = mass_matrix(Dp_bounded)
-    @test M == mass_matrix(Dm_bounded)
-    @test M == mass_matrix(Dc_bounded)
-    Dp_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -(acc_order - 1) ÷ 2)
-    Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -acc_order + (acc_order - 1) ÷ 2)
-    Dp = Matrix(Dp_bounded)
-    Dm = Matrix(Dm_bounded)
-    @test Dp[interior,interior] ≈ Matrix(Dp_periodic)[interior,interior]
-    @test Dm[interior,interior] ≈ Matrix(Dm_periodic)[interior,interior]
-    res = M * Dp + Dm' * M
-    res[1,1] += 1
-    res[end,end] -= 1
-    @test norm(res) < N * eps()
-    x = grid(Dp_bounded)
-    for D in (Dp_bounded, Dm_bounded, Dc_bounded)
-      @test norm(D * x.^0) < N * eps()
-      for k in 1:acc_order÷2
-        @test D * x.^k ≈ k .* x.^(k-1)
+      for compact in (true, false)
+        show(IOContext(devnull, :compact => compact), Dp_bounded)
+        show(IOContext(devnull, :compact => compact), Dm_bounded)
+        show(IOContext(devnull, :compact => compact), Dc_bounded)
+        summary(IOContext(devnull, :compact => compact), Dp_bounded)
+        summary(IOContext(devnull, :compact => compact), Dm_bounded)
+        summary(IOContext(devnull, :compact => compact), Dc_bounded)
       end
-      for k in acc_order÷2+1:acc_order
-        @test (D * x.^k)[interior] ≈ (k .* x.^(k-1))[interior]
+
+      M = mass_matrix(Dp_bounded)
+      @test M == mass_matrix(Dm_bounded)
+      @test M == mass_matrix(Dc_bounded)
+      Dp_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -(acc_order - 1) ÷ 2)
+      Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -acc_order + (acc_order - 1) ÷ 2)
+      Dp = Matrix(Dp_bounded)
+      Dm = Matrix(Dm_bounded)
+      @test Dp[interior,interior] ≈ Matrix(Dp_periodic)[interior,interior]
+      @test Dm[interior,interior] ≈ Matrix(Dm_periodic)[interior,interior]
+      res = M * Dp + Dm' * M
+      res[1,1] += 1
+      res[end,end] -= 1
+      @test norm(res) < N * eps()
+      x = grid(Dp_bounded)
+      for D in (Dp_bounded, Dm_bounded, Dc_bounded)
+        @test norm(D * x.^0) < N * eps()
+        for k in 1:acc_order÷2
+          @test D * x.^k ≈ k .* x.^(k-1)
+        end
+        for k in acc_order÷2+1:acc_order
+          @test (D * x.^k)[interior] ≈ (k .* x.^(k-1))[interior]
+        end
       end
+      diss = M * (Dp - Dm)
+      @test diss ≈ diss'
+      @test maximum(eigvals(Symmetric(diss))) < N * eps()
     end
-    diss = M * (Dp - Dm)
-    @test diss ≈ diss'
-    @test maximum(eigvals(Symmetric(diss))) < N * eps()
   end
 end