Mercurial > repos > public > sbplib_julia
changeset 765:fdd48f6ace1c operator_storage_array_of_table
Merge in default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Sun, 07 Feb 2021 21:32:21 +0100 |
| parents | d34b515b0ae7 (current diff) 621460cf8279 (diff) |
| children | 7624a1350ece |
| files | src/SbpOperators/readoperator.jl test/testSbpOperators.jl |
| diffstat | 8 files changed, 74 insertions(+), 37 deletions(-) [+] |
line wrap: on
line diff
--- a/TODO.md Sat Feb 06 20:45:34 2021 +0100 +++ b/TODO.md Sun Feb 07 21:32:21 2021 +0100 @@ -21,6 +21,8 @@ 3. [ ] Fix the rest of the library Should getregion also work for getregion(::Colon,...) - [ ] Add possibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`? + - [ ] Add custom pretty printing to LazyTensors/SbpOperators to enhance readability of e.g error messages. + See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing) ## Repo - [ ] Rename repo to Sbplib.jl
--- a/src/Grids/EquidistantGrid.jl Sat Feb 06 20:45:34 2021 +0100 +++ b/src/Grids/EquidistantGrid.jl Sun Feb 07 21:32:21 2021 +0100 @@ -46,10 +46,7 @@ The dimension of the grid. """ -function dimension(grid::EquidistantGrid) - return length(grid.size) -end - +dimension(grid::EquidistantGrid{Dim}) where Dim = Dim """ spacing(grid::EquidistantGrid) @@ -95,3 +92,15 @@ return EquidistantGrid(size, limit_lower, limit_upper) end export restrict + +""" + boundary_identifiers(::EquidistantGrid) + +Returns a tuple containing the boundary identifiers for the grid, stored as + (CartesianBoundary(1,Lower), + CartesianBoundary(1,Upper), + CartesianBoundary(2,Lower), + ...) +""" +boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),dimension(g)))...)...,) +export boundary_identifiers
--- a/src/SbpOperators/readoperator.jl Sat Feb 06 20:45:34 2021 +0100 +++ b/src/SbpOperators/readoperator.jl Sun Feb 07 21:32:21 2021 +0100 @@ -29,13 +29,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( @@ -136,7 +136,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 20:45:34 2021 +0100 +++ b/src/SbpOperators/stencil.jl Sun Feb 07 21:32:21 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 20:45:34 2021 +0100 +++ b/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl Sun Feb 07 21:32:21 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 20:45:34 2021 +0100 +++ b/src/SbpOperators/volumeops/quadratures/quadrature.jl Sun Feb 07 21:32:21 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/testGrids.jl Sat Feb 06 20:45:34 2021 +0100 +++ b/test/testGrids.jl Sun Feb 07 21:32:21 2021 +0100 @@ -1,5 +1,6 @@ using Sbplib.Grids using Test +using Sbplib.RegionIndices @testset "Grids" begin @@ -53,6 +54,15 @@ @test restrict(g, 2:3) == EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0)) @test restrict(g, [1,3]) == EquidistantGrid((2,3),(0.0,0.0),(2.0,3.0)) @test restrict(g, [2,1]) == EquidistantGrid((5,2),(0.0,0.0),(1.0,2.0)) + + @testset "boundary_identifiers" begin + g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) + bids = (CartesianBoundary{1,Lower}(),CartesianBoundary{1,Upper}(), + CartesianBoundary{2,Lower}(),CartesianBoundary{2,Upper}(), + CartesianBoundary{3,Lower}(),CartesianBoundary{3,Upper}()) + @test boundary_identifiers(g) == bids + @inferred boundary_identifiers(g) + end end end
--- a/test/testSbpOperators.jl Sat Feb 06 20:45:34 2021 +0100 +++ b/test/testSbpOperators.jl Sun Feb 07 21:32:21 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 @@ -105,40 +108,40 @@ end @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)) @@ -154,8 +157,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.)) @@ -438,7 +441,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 @@ -522,7 +525,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))