Mercurial > repos > public > sbplib_julia

changeset 990:b6238afd3bb0 feature/stencil_set_type

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add methods for creating derivative operators in 1D from stencil sets without providing directions
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 18 Mar 2022 13:02:46 +0100
parents 7bf3121c6864
children 37fd8c1cadb2
files src/SbpOperators/readoperator.jl src/SbpOperators/volumeops/derivatives/first_derivative.jl src/SbpOperators/volumeops/derivatives/second_derivative.jl test/SbpOperators/volumeops/derivatives/first_derivative_test.jl test/SbpOperators/volumeops/derivatives/second_derivative_test.jl test/SbpOperators/volumeops/laplace/laplace_test.jl
diffstat 6 files changed, 49 insertions(+), 30 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/readoperator.jl	Thu Mar 17 21:31:20 2022 +0100
+++ b/src/SbpOperators/readoperator.jl	Fri Mar 18 13:02:46 2022 +0100
@@ -14,12 +14,7 @@
 
 A `StencilSet` contains a set of associated stencils. The stencils
 are are stored in a table, and can be accesed by indexing into the `StencilSet`.
-"""
-struct StencilSet
-    table
-end
 
-"""
     StencilSet(filename; filters)
 
 Creates a `StencilSet` from a TOML file based on some key-value
@@ -38,7 +33,13 @@
 
 See also [`sbp_operators_path`](@ref), [`get_stencil_set`](@ref), [`parse_stencil`](@ref), [`parse_scalar`](@ref), [`parse_tuple`](@ref),.
 """
-StencilSet(filename; filters...) = StencilSet(get_stencil_set(TOML.parsefile(filename); filters...))
+struct StencilSet
+    table
+    function StencilSet(filename; filters...)
+        return new(get_stencil_set(TOML.parsefile(filename); filters...))
+    end
+end
+
 Base.getindex(set::StencilSet,I...) = set.table[I...]
 
 """
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Thu Mar 17 21:31:20 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Fri Mar 18 13:02:46 2022 +0100
@@ -16,11 +16,18 @@
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd, direction)
 end
-first_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = first_derivative(grid,inner_stencil,closure_stencils,1)
 
 
 """
-    first_derivative(grid, stencil_set, direction)
+    first_derivative(grid, inner_stencil, closure_stencils)
+
+Creates a `first_derivative` operator on a 1D `grid` given `inner_stencil` and `closure_stencils`.
+"""
+first_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = first_derivative(grid, inner_stencil, closure_stencils, 1)
+
+
+"""
+    first_derivative(grid, stencil_set::StencilSet, direction)
 
 Creates a `first_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`.
 """
@@ -30,6 +37,10 @@
     first_derivative(grid,inner_stencil,closure_stencils,direction);
 end
 
-# TODO: Not possible to remove ::Stencil from first_derivative(grid::EquidistantGrid{1},...) due to type deduction failing. 
-# Is this due to the type ambuguity of StencilSet? Possible to address in some other way? 
-# If not, should we drop ::StencilSet from the method above (it is not required)?
\ No newline at end of file
+
+"""
+    first_derivative(grid, stencil_set)
+
+Creates a `first_derivative` operator on a 1D `grid` given a `stencil_set`.
+"""
+first_derivative(grid::EquidistantGrid{1}, stencil_set) = first_derivative(grid, stencil_set, 1)
--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Thu Mar 17 21:31:20 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Fri Mar 18 13:02:46 2022 +0100
@@ -16,7 +16,15 @@
     h_inv = inverse_spacing(grid)[direction]
     return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction)
 end
-second_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
+
+
+"""
+    second_derivative(grid, inner_stencil, closure_stencils)
+
+Creates a `second_derivative` operator on a 1D `grid` given `inner_stencil` and `closure_stencils`.
+"""
+second_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = second_derivative(grid, inner_stencil, closure_stencils,1)
+
 
 """
     second_derivative(grid, stencil_set, direction)
@@ -27,8 +35,12 @@
     inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
     closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
     second_derivative(grid,inner_stencil,closure_stencils,direction);
-end 
+end
+
 
-# TODO: Not possible to remove ::Stencil from second_derivative(grid::EquidistantGrid{1},...) due to type deduction failing. 
-# Is this due to the type ambuguity of StencilSet? Possible to address in some other way? 
-# If not, should we drop ::StencilSet from the method above (it is not required)?
\ No newline at end of file
+"""
+    second_derivative(grid, stencil_set)
+
+Creates a `second_derivative` operator on a 1D `grid` given a `stencil_set`.
+"""
+second_derivative(grid::EquidistantGrid{1}, stencil_set) = second_derivative(grid, stencil_set, 1)
\ No newline at end of file
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Thu Mar 17 21:31:20 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Fri Mar 18 13:02:46 2022 +0100
@@ -29,12 +29,14 @@
         
         @test first_derivative(g₁, stencil_set, 1) isa TensorMapping{Float64,1,1}
         @test first_derivative(g₂, stencil_set, 2) isa TensorMapping{Float64,2,2}
+        @test first_derivative(g₁, stencil_set, 1) == first_derivative(g₁, stencil_set)
 
         interior_stencil = CenteredStencil(-1,0,1)
         closure_stencils = [Stencil(-1,1, center=1)]
 
         @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa TensorMapping{Float64,1,1}
         @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa VolumeOperator
+        @test first_derivative(g₁, interior_stencil, closure_stencils, 1) == first_derivative(g₁, interior_stencil, closure_stencils)
         @test first_derivative(g₂, interior_stencil, closure_stencils, 2) isa TensorMapping{Float64,2,2}
     end
 
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Thu Mar 17 21:31:20 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Fri Mar 18 13:02:46 2022 +0100
@@ -21,11 +21,12 @@
             Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils,1)
             @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils)
             @test Dₓₓ == second_derivative(g_1D,stencil_set,1)
+            @test Dₓₓ == second_derivative(g_1D,stencil_set)
             @test Dₓₓ isa VolumeOperator
         end
         @testset "2D" begin
             Dₓₓ = second_derivative(g_2D,inner_stencil,closure_stencils,1)
-            D2 = second_derivative(g_1D,inner_stencil,closure_stencils)
+            D2 = second_derivative(g_1D,inner_stencil,closure_stencils,1)
             I = IdentityMapping{Float64}(size(g_2D)[2])
             @test Dₓₓ == D2⊗I
             @test Dₓₓ == second_derivative(g_2D,stencil_set,1)
@@ -51,9 +52,7 @@
             # implies that L*v should be exact for monomials up to order 2.
             @testset "2nd order" begin
                 stencil_set = StencilSet(operator_path; order=2)
-                inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-			    closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-                Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
+                Dₓₓ = second_derivative(g_1D,stencil_set)
                 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
                 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
                 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
@@ -64,9 +63,7 @@
             # implies that L*v should be exact for monomials up to order 3.
             @testset "4th order" begin
                 stencil_set = StencilSet(operator_path; order=4)
-                inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-			    closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-                Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
+                Dₓₓ = second_derivative(g_1D,stencil_set)
                 # NOTE: high tolerances for checking the "exact" differentiation
                 # due to accumulation of round-off errors/cancellation errors?
                 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
@@ -92,9 +89,7 @@
             # implies that L*v should be exact for binomials up to order 2.
             @testset "2nd order" begin
                 stencil_set = StencilSet(operator_path; order=2)
-                inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-                closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-                Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,2)
+                Dyy = second_derivative(g_2D,stencil_set,2)
                 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
                 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
                 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
@@ -105,9 +100,7 @@
             # implies that L*v should be exact for binomials up to order 3.
             @testset "4th order" begin
                 stencil_set = StencilSet(operator_path; order=4)
-                inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
-                closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
-                Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,2)
+                Dyy = second_derivative(g_2D,stencil_set,2)
                 # NOTE: high tolerances for checking the "exact" differentiation
                 # due to accumulation of round-off errors/cancellation errors?
                 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Thu Mar 17 21:31:20 2022 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Fri Mar 18 13:02:46 2022 +0100
@@ -69,7 +69,7 @@
 @testset "laplace" begin
     @testset "1D" begin
         Δ = laplace(g_1D, inner_stencil, closure_stencils)
-        @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils)
+        @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils, 1)
         @test Δ isa TensorMapping{T,1,1}  where T
     end
     @testset "3D" begin