Mercurial > repos > public > sbplib_julia
changeset 675:1ce3a104afc8 feature/boundary_quads
Merge in default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Sun, 07 Feb 2021 21:28:53 +0100 |
| parents | 538ccbaeb1f8 (current diff) 621460cf8279 (diff) |
| children | 1d3e29ffc6c6 |
| files | src/SbpOperators/volumeops/quadratures/quadrature.jl test/testSbpOperators.jl |
| diffstat | 5 files changed, 50 insertions(+), 34 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/readoperator.jl Sat Feb 06 12:04:06 2021 +0100 +++ b/src/SbpOperators/readoperator.jl Sun Feb 07 21:28:53 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 Sat Feb 06 12:04:06 2021 +0100 +++ b/src/SbpOperators/stencil.jl Sun Feb 07 21:28:53 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 Sat Feb 06 12:04:06 2021 +0100 +++ b/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl Sun Feb 07 21:28:53 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 Sat Feb 06 12:04:06 2021 +0100 +++ b/src/SbpOperators/volumeops/quadratures/quadrature.jl Sun Feb 07 21:28:53 2021 +0100 @@ -26,7 +26,7 @@ export quadrature function quadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T}}) where {M,T} - inner_stencil = Stencil(Tuple{T}(1),center=1) + inner_stencil = CenteredStencil(one(T)) return quadrature(grid, inner_stencil, closure_stencils) end @@ -52,7 +52,7 @@ end function boundary_quadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T}}, id::CartesianBoundary) where {M,T} - inner_stencil = Stencil(Tuple{T}(1),center=1) + inner_stencil = CenteredStencil(one(T)) return boundary_quadrature(grid,inner_stencil,closure_stencils,id) end
--- a/test/testSbpOperators.jl Sat Feb 06 12:04:06 2021 +0100 +++ b/test/testSbpOperators.jl Sun Feb 07 21:28:53 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.)) @@ -402,7 +405,7 @@ op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) @testset "1D" begin H = quadrature(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 @@ -519,7 +522,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))