Mercurial > repos > public > sbplib_julia

changeset 1035:ceda69b8f27a feature/dissipation_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add test for transpose equality and fix bugs found
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 22 Mar 2022 13:28:06 +0100
parents ed19c549c506
children 0e31b9901160
files src/SbpOperators/volumeops/derivatives/dissipation.jl test/SbpOperators/volumeops/derivatives/dissipation_test.jl
diffstat 2 files changed, 45 insertions(+), 18 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/derivatives/dissipation.jl	Tue Mar 22 12:33:58 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Tue Mar 22 13:28:06 2022 +0100
@@ -33,8 +33,16 @@
 midpoint(weights) = length(weights)÷2 + 1
 midpoint_transpose(weights) = length(weights)+1 - midpoint(weights)
 
-dissipation_interior_stencil(weights) =           Stencil(weights..., center=midpoint(weights))
-dissipation_transpose_interior_stencil(weights) = Stencil(weights..., center=midpoint_transpose(weights))
+function dissipation_interior_stencil(weights)
+    return Stencil(weights..., center=midpoint(weights))
+end
+function dissipation_transpose_interior_stencil(weights)
+    if iseven(length(weights))
+        weights = map(-, weights)
+    end
+
+    return Stencil(weights..., center=midpoint_transpose(weights))
+end
 
 dissipation_lower_closure_size(weights) = midpoint(weights) - 1
 dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights)
--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Tue Mar 22 12:33:58 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Tue Mar 22 13:28:06 2022 +0100
@@ -36,7 +36,7 @@
      @testset "Accuracy conditions" begin
         N = 20
         g = EquidistantGrid(N, 0//1,2//1)
-        h = spacing(g)[1]
+        h = only(spacing(g))
         @testset "D_$p" for p ∈ [1,2,3,4]
             D,Dᵀ = undevided_dissipation(g, p)
 
@@ -46,16 +46,35 @@
 
                 @test D*v == h^p * vₚₓ
             end
+        end
+    end
 
+    @testset "tanspose equality" begin
+        function get_matrix(D)
+            N = only(range_size(D))
+            M = only(domain_size(D))
 
-            @testset "x^$k" for k ∈ 0:p
-                v  = evalOn(g, x->monomial(x,k))
-                vₚₓ = evalOn(g, x->monomial(x,k-p))
+            Dmat = zeros(N,M)
+            e = zeros(M)
+            for i ∈ 1:M
+                if i > 1
+                    e[i-1] = 0.
+                end
+                e[i] = 1.
+                Dmat[:,i] = D*e
+            end
 
-                for i ∈ SbpOperators.lower_closure_size(Dᵀ)+1:N-SbpOperators.upper_closure_size(Dᵀ)
-                    @test (Dᵀ*v)[i] == h^p * vₚₓ[i]
-                end
-            end
+            return Dmat
+        end
+
+        g = EquidistantGrid(11, 0., 1.)
+        @testset "D_$p" for p ∈ [1,2,3,4]
+            D,Dᵀ = undevided_dissipation(g, p)
+
+            D̄  = get_matrix(D)
+            D̄ᵀ = get_matrix(Dᵀ)
+
+            @test D̄ == D̄ᵀ'
         end
     end
 end
@@ -68,17 +87,17 @@
 end
 
 @testset "dissipation_interior_stencil" begin
-    @test dissipation_interior_stencil(dissipation_interior_weights(1)) == Stencil(-1,1, center=2)
-    @test dissipation_interior_stencil(dissipation_interior_weights(2)) == Stencil(1,-2,1, center=2)
-    @test dissipation_interior_stencil(dissipation_interior_weights(3)) == Stencil(-1,3,-3,1, center=3)
-    @test dissipation_interior_stencil(dissipation_interior_weights(4)) == Stencil(1, -4, 6, -4, 1, center=3)
+    @test dissipation_interior_stencil(dissipation_interior_weights(1)) == Stencil(-1, 1, center=2)
+    @test dissipation_interior_stencil(dissipation_interior_weights(2)) == Stencil( 1,-2, 1, center=2)
+    @test dissipation_interior_stencil(dissipation_interior_weights(3)) == Stencil(-1, 3,-3, 1, center=3)
+    @test dissipation_interior_stencil(dissipation_interior_weights(4)) == Stencil( 1,-4, 6,-4, 1, center=3)
 end
 
 @testset "dissipation_transpose_interior_stencil" begin
-    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(1)) == Stencil(-1,1, center=1)
-    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(2)) == Stencil(1,-2,1, center=2)
-    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(3)) == Stencil(-1,3,-3,1, center=2)
-    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(4)) == Stencil(1, -4, 6, -4, 1, center=3)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(1)) == Stencil(1,-1, center=1)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(2)) == Stencil(1,-2, 1, center=2)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(3)) == Stencil(1,-3, 3,-1, center=2)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(4)) == Stencil(1,-4, 6,-4, 1, center=3)
 end
 
 @testset "midpoint" begin