changeset 1021:ee5a641a8277 feature/dissipation_operators

Add some tests and start implementing dissipation()
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 18 Mar 2022 15:44:03 +0100
parents 3f5137ce3aa1
children e74c41c4b60e
files src/SbpOperators/SbpOperators.jl src/SbpOperators/volumeops/derivatives/dissipation.jl test/SbpOperators/volumeops/derivatives/dissipation_test.jl
diffstat 3 files changed, 50 insertions(+), 2 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/SbpOperators.jl	Fri Mar 18 12:24:10 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Fri Mar 18 15:44:03 2022 +0100
@@ -9,6 +9,7 @@
 export normal_derivative
 export first_derivative
 export second_derivative
+export dissipation
 
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
--- a/src/SbpOperators/volumeops/derivatives/dissipation.jl	Fri Mar 18 12:24:10 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Fri Mar 18 15:44:03 2022 +0100
@@ -1,3 +1,12 @@
+function dissipation(g::EquidistantGrid, p, direction)
+    h_inv = inverse_spacing(g)[direction]
+
+    # D = volume_operator(g,CenteredStencil(1),(CenteredStencil(1)), )
+    return nothing, nothing
+end
+
+dissipation(g::EquidistantGrid{1}, p) = dissipation(g, p, 1)
+
 function dissipation_interior_weights(p)
    if p == 0
        return (1,)
--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Fri Mar 18 12:24:10 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Fri Mar 18 15:44:03 2022 +0100
@@ -1,8 +1,8 @@
 using Test
 
 using Sbplib.SbpOperators
-# using Sbplib.Grids
-# using Sbplib.LazyTensors
+using Sbplib.Grids
+using Sbplib.LazyTensors
 
 using Sbplib.SbpOperators: Stencil
 
@@ -13,6 +13,44 @@
 using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
 using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
 
+"""
+    monomial(x,k)
+
+Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
+Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
+"""
+function monomial(x,k)
+    if k < 0
+        return zero(x)
+    end
+    x^k/factorial(k)
+end
+
+@testset "dissipation" begin
+    g = EquidistantGrid(20, 0., 11.)
+    D,Dᵀ = dissipation(g, 1)
+
+    @test D isa TensorMapping{Float64,1,1} where T
+    @test Dᵀ isa TensorMapping{Float64,1,1} where T
+
+     @testset "Accuracy conditions" begin
+        N = 20
+        g = EquidistantGrid(N, 0//1,2//1)
+        @testset "D_$p" for p ∈ [1,2,3,4]
+            D,Dᵀ = dissipation(g, p)
+
+            @testset "x^$k" for k ∈ 0:1
+                v = evalOn(g, x->monomial(x,k))
+
+                x, = points(g)[10]
+                @test (D*v)[10] == monomial(x,k-1)
+            end
+
+            # Test Dᵀ works backwards and interior forwards
+        end
+    end
+end
+
 @testset "dissipation_interior_weights" begin
     @test dissipation_interior_weights(1) == (-1, 1)
     @test dissipation_interior_weights(2) == (1,-2, 1)