Mercurial > repos > public > sbplib_julia
changeset 674:621460cf8279
Merge feature/stencil_convenience.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Sun, 07 Feb 2021 21:16:40 +0100 |
| parents | 7d7c1d636de3 (current diff) e5e935c95471 (diff) |
| children | 1ce3a104afc8 1accc3e051d0 d52902f36868 fdd48f6ace1c |
| files | |
| diffstat | 5 files changed, 50 insertions(+), 34 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/readoperator.jl Sun Jan 31 20:57:58 2021 +0100 +++ b/src/SbpOperators/readoperator.jl Sun Feb 07 21:16:40 2021 +0100 @@ -26,13 +26,13 @@ closureStencils = (closureStencils..., get_stencil(operators, "D2", "closure_stencils", i; center=i)) end # TODO: Get rid of the padding here. Any padding should be handled by the consturctor accepting the stencils. - eClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, e["closure"]), boundarySize), center=1) - dClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, d1["closure"]), boundarySize), center=1) + eClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, e["closure"]), boundarySize)..., center=1) + dClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, d1["closure"]), boundarySize)..., center=1) q_tuple = pad_tuple(toml_string_array_to_tuple(Float64, H["closure"]), boundarySize) quadratureClosure = Vector{typeof(innerStencil)}() for i ∈ 1:boundarySize - quadratureClosure = (quadratureClosure..., Stencil((q_tuple[i],), center=1)) + quadratureClosure = (quadratureClosure..., Stencil(q_tuple[i], center=1)) end d2 = SbpOperators.D2( @@ -99,7 +99,7 @@ center = div(width,2)+1 end - return Stencil(Tuple(stencil_weights), center=center) + return Stencil(stencil_weights..., center=center) end """
--- a/src/SbpOperators/stencil.jl Sun Jan 31 20:57:58 2021 +0100 +++ b/src/SbpOperators/stencil.jl Sun Feb 07 21:16:40 2021 +0100 @@ -1,3 +1,5 @@ +export CenteredStencil + struct Stencil{T<:Real,N} range::Tuple{Int,Int} weights::NTuple{N,T} @@ -13,13 +15,24 @@ Create a stencil with the given weights with element `center` as the center of the stencil. """ -function Stencil(weights::NTuple; center::Int) +function Stencil(weights::Vararg{Number}; center::Int) N = length(weights) range = (1, N) .- center return Stencil(range, weights) end +function CenteredStencil(weights::Vararg) + if iseven(length(weights)) + throw(ArgumentError("a centered stencil must have an odd number of weights.")) + end + + r = length(weights) ÷ 2 + + return Stencil((-r, r), weights) +end + + """ scale(s::Stencil, a)
--- a/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl Sun Jan 31 20:57:58 2021 +0100 +++ b/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl Sun Feb 07 21:16:40 2021 +0100 @@ -32,7 +32,7 @@ the closure stencils are those of the quadrature operator (and not the inverse). """ function InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {T,M} - inv_inner_stencil = Stencil(Tuple{T}(1),center=1) + inv_inner_stencil = Stencil(one(T), center=1) inv_closure_stencils = reciprocal_stencil.(closure_stencils) return InverseQuadrature(grid, inv_inner_stencil, inv_closure_stencils) end
--- a/src/SbpOperators/volumeops/quadratures/quadrature.jl Sun Jan 31 20:57:58 2021 +0100 +++ b/src/SbpOperators/volumeops/quadratures/quadrature.jl Sun Feb 07 21:16:40 2021 +0100 @@ -30,7 +30,7 @@ closure stencils (i.e the operator is diagonal) """ function DiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {M,T} - inner_stencil = Stencil(Tuple{T}(1),center=1) + inner_stencil = Stencil(one(T), center=1) return Quadrature(grid, inner_stencil, closure_stencils) end export DiagonalQuadrature
--- a/test/testSbpOperators.jl Sun Jan 31 20:57:58 2021 +0100 +++ b/test/testSbpOperators.jl Sun Feb 07 21:16:40 2021 +0100 @@ -24,9 +24,12 @@ @test eltype(s) == Float64 @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.)) - @test Stencil((1,2,3,4), center=1) == Stencil((0, 3),(1,2,3,4)) - @test Stencil((1,2,3,4), center=2) == Stencil((-1, 2),(1,2,3,4)) - @test Stencil((1,2,3,4), center=4) == Stencil((-3, 0),(1,2,3,4)) + @test Stencil(1,2,3,4; center=1) == Stencil((0, 3),(1,2,3,4)) + @test Stencil(1,2,3,4; center=2) == Stencil((-1, 2),(1,2,3,4)) + @test Stencil(1,2,3,4; center=4) == Stencil((-3, 0),(1,2,3,4)) + + @test CenteredStencil(1,2,3,4,5) == Stencil((-2, 2), (1,2,3,4,5)) + @test_throws ArgumentError CenteredStencil(1,2,3,4) end @testset "parse_rational" begin @@ -67,40 +70,40 @@ parsed_toml = TOML.parse(toml_str) @testset "get_stencil" begin - @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil((-1/2, 0., 1/2), center=2) - @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil((-1/2, 0., 1/2); center=1) - @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil((-1/2, 0., 1/2); center=3) + @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil(-1/2, 0., 1/2, center=2) + @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil(-1/2, 0., 1/2; center=1) + @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil(-1/2, 0., 1/2; center=3) - @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil((1.,), center=1) + @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil(1.; center=1) @test_throws AssertionError get_stencil(parsed_toml, "meta", "type") @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils") end @testset "get_stencils" begin - @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil((-1., 1.), center=1),) - @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil((-1., 1.), center=2),) - @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil((-1., 1.), center=2),) + @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil(-1., 1., center=1),) + @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil(-1., 1., center=2),) + @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil(-1., 1., center=2),) @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == ( - Stencil(( 2., -5., 4., -1., 0., 0.), center=1), - Stencil(( 1., -2., 1., 0., 0., 0.), center=1), - Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=1), - Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=1), + Stencil( 2., -5., 4., -1., 0., 0., center=1), + Stencil( 1., -2., 1., 0., 0., 0., center=1), + Stencil( -4/43, 59/43, -110/43, 59/43, -4/43, 0., center=1), + Stencil( -1/49, 0., 59/49, -118/49, 64/49, -4/49, center=1), ) @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == ( - Stencil(( 2., -5., 4., -1., 0., 0.), center=4), - Stencil(( 1., -2., 1., 0., 0., 0.), center=2), - Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=3), - Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=1), + Stencil( 2., -5., 4., -1., 0., 0., center=4), + Stencil( 1., -2., 1., 0., 0., 0., center=2), + Stencil( -4/43, 59/43, -110/43, 59/43, -4/43, 0., center=3), + Stencil( -1/49, 0., 59/49, -118/49, 64/49, -4/49, center=1), ) @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == ( - Stencil(( 2., -5., 4., -1., 0., 0.), center=1), - Stencil(( 1., -2., 1., 0., 0., 0.), center=2), - Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=3), - Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=4), + Stencil( 2., -5., 4., -1., 0., 0., center=1), + Stencil( 1., -2., 1., 0., 0., 0., center=2), + Stencil( -4/43, 59/43, -110/43, 59/43, -4/43, 0., center=3), + Stencil( -1/49, 0., 59/49, -118/49, 64/49, -4/49, center=4), ) @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3)) @@ -116,8 +119,8 @@ end @testset "VolumeOperator" begin - inner_stencil = Stencil(1/4 .* (1.,2.,1.),center=2) - closure_stencils = (Stencil(1/2 .* (1.,1.),center=1),Stencil((0.,1.),center=2)) + inner_stencil = CenteredStencil(1/4, 2/4, 1/4) + closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(0.,1.; center=2)) g_1D = EquidistantGrid(11,0.,1.) g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.)) g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.)) @@ -290,7 +293,7 @@ @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2 end end - + @testset "2D" begin l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2)); binomials = () @@ -400,7 +403,7 @@ op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) @testset "1D" begin H = DiagonalQuadrature(g_1D,op.quadratureClosure) - inner_stencil = Stencil((1.,),center=1) + inner_stencil = CenteredStencil(1.) @test H == Quadrature(g_1D,inner_stencil,op.quadratureClosure) @test H isa TensorMapping{T,1,1} where T end @@ -484,7 +487,7 @@ op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) @testset "1D" begin Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure); - inner_stencil = Stencil((1.,),center=1) + inner_stencil = CenteredStencil(1.) closures = () for i = 1:length(op.quadratureClosure) closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights))