Mercurial > repos > public > sbplib_julia
changeset 1099:05a25a5063bb refactor/sbpoperators/inflation
Try to remove volume_operator and boundary_operator methods
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Mon, 21 Mar 2022 12:51:39 +0100 |
| parents | 5be17f647018 |
| children | 157a78959e5d |
| files | src/SbpOperators/boundaryops/boundary_operator.jl src/SbpOperators/boundaryops/boundary_restriction.jl src/SbpOperators/boundaryops/normal_derivative.jl src/SbpOperators/volumeops/derivatives/first_derivative.jl src/SbpOperators/volumeops/derivatives/second_derivative.jl src/SbpOperators/volumeops/volume_operator.jl test/SbpOperators/boundaryops/boundary_operator_test.jl test/SbpOperators/volumeops/derivatives/first_derivative_test.jl test/SbpOperators/volumeops/volume_operator_test.jl |
| diffstat | 9 files changed, 116 insertions(+), 192 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/boundaryops/boundary_operator.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Mon Mar 21 12:51:39 2022 +0100 @@ -1,27 +1,3 @@ -""" - boundary_operator(grid,closure_stencil,boundary) - -Creates a boundary operator on a `Dim`-dimensional grid for the -specified `boundary`. The action of the operator is determined by `closure_stencil`. - -When `Dim=1`, the corresponding `BoundaryOperator` tensor mapping is returned. -When `Dim>1`, the `BoundaryOperator` `op` is inflated by the outer product -of `IdentityTensors` in orthogonal coordinate directions, e.g for `Dim=3`, -the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`. -""" -function boundary_operator(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary) - #TODO:Check that dim(boundary) <= Dim? - - d = dim(boundary) - op = BoundaryOperator(restrict(grid, d), closure_stencil, region(boundary)) - - # Create 1D IdentityTensors for each coordinate direction - one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid))) - Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids)) - - return LazyTensors.inflate(op, size(grid), d) -end - """ BoundaryOperator{T,R,N} <: LazyTensor{T,0,1}
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Mon Mar 21 12:51:39 2022 +0100 @@ -15,7 +15,9 @@ """ function boundary_restriction(grid, closure_stencil::Stencil, boundary) converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil) - return SbpOperators.boundary_operator(grid, converted_stencil, boundary) + + op = BoundaryOperator(restrict(grid, dim(boundary)), converted_stencil, region(boundary)) + return LazyTensors.inflate(op, size(grid), dim(boundary)) end """
--- a/src/SbpOperators/boundaryops/normal_derivative.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/src/SbpOperators/boundaryops/normal_derivative.jl Mon Mar 21 12:51:39 2022 +0100 @@ -13,7 +13,9 @@ function normal_derivative(grid, closure_stencil::Stencil, boundary) direction = dim(boundary) h_inv = inverse_spacing(grid)[direction] - return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary) + + op = BoundaryOperator(restrict(grid, dim(boundary)), scale(closure_stencil,h_inv), region(boundary)) + return LazyTensors.inflate(op, size(grid), dim(boundary)) end """
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl Mon Mar 21 12:51:39 2022 +0100 @@ -14,7 +14,9 @@ """ function first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction) h_inv = inverse_spacing(grid)[direction] - return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd, direction) + + D₁ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd) + return LazyTensors.inflate(D₁, size(grid), direction) end first_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = first_derivative(grid,inner_stencil,closure_stencils,1)
--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl Mon Mar 21 12:51:39 2022 +0100 @@ -14,7 +14,9 @@ """ function second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction) 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) + + D₂ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even) + return LazyTensors.inflate(D₂, size(grid), direction) end second_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = second_derivative(grid,inner_stencil,closure_stencils,1)
--- a/src/SbpOperators/volumeops/volume_operator.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/src/SbpOperators/volumeops/volume_operator.jl Mon Mar 21 12:51:39 2022 +0100 @@ -1,22 +1,3 @@ -""" - volume_operator(grid, inner_stencil, closure_stencils, parity, direction) - -Creates a volume operator on a `Dim`-dimensional grid acting along the -specified coordinate `direction`. The action of the operator is determined by -the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the -corresponding `VolumeOperator` tensor mapping is returned. When `Dim>1`, the -returned operator is the appropriate outer product of a one-dimensional -operators and `IdentityTensor`s, e.g for `Dim=3` the volume operator in the -y-direction is `I⊗op⊗I`. -""" -function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction) - #TODO: Check that direction <= Dim? - - op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity) - return LazyTensors.inflate(op, size(grid), direction) -end -# TBD: Should the inflation happen here or should we remove this method and do it at the caller instead? - """ VolumeOperator{T,N,M,K} <: TensorOperator{T,1} Implements a one-dimensional constant coefficients volume operator
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Mon Mar 21 12:51:39 2022 +0100 @@ -6,7 +6,8 @@ using Sbplib.RegionIndices import Sbplib.SbpOperators.Stencil import Sbplib.SbpOperators.BoundaryOperator -import Sbplib.SbpOperators.boundary_operator + +# TODO: What should happen to all the commented tests? Deleted? Replicated for user code? @testset "BoundaryOperator" begin closure_stencil = Stencil((0,2), (2.,1.,3.)) @@ -14,26 +15,26 @@ g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) @testset "Constructors" begin - @testset "1D" begin + @testset "1D" begin # TODO: Remove these testsets op_l = BoundaryOperator{Lower}(closure_stencil,size(g_1D)[1]) @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower()) - @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}()) @test op_l isa LazyTensor{T,0,1} where T - op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1]) + op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1]) # TBD: Is this constructor really needed? looks weird! @test op_r == BoundaryOperator(g_1D,closure_stencil,Upper()) - @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}()) @test op_r isa LazyTensor{T,0,1} where T end - @testset "2D" begin - e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}()) - @test e_w isa InflatedLazyTensor - @test e_w isa LazyTensor{T,1,2} where T - end + # @testset "2D" begin + # e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}()) + # @test e_w isa InflatedLazyTensor + # @test e_w isa LazyTensor{T,1,2} where T + # end end - op_l, op_r = boundary_operator.(Ref(g_1D), Ref(closure_stencil), boundary_identifiers(g_1D)) - op_w, op_e, op_s, op_n = boundary_operator.(Ref(g_2D), Ref(closure_stencil), boundary_identifiers(g_2D)) + + op_l = BoundaryOperator(g_1D, closure_stencil, Lower()) + op_r = BoundaryOperator(g_1D, closure_stencil, Upper()) + # op_w, op_e, op_s, op_n = boundary_operator.(Ref(g_2D), Ref(closure_stencil), boundary_identifiers(g_2D)) @testset "Sizes" begin @testset "1D" begin @@ -44,17 +45,17 @@ @test range_size(op_r) == () end - @testset "2D" begin - @test domain_size(op_w) == (11,15) - @test domain_size(op_e) == (11,15) - @test domain_size(op_s) == (11,15) - @test domain_size(op_n) == (11,15) + # @testset "2D" begin + # @test domain_size(op_w) == (11,15) + # @test domain_size(op_e) == (11,15) + # @test domain_size(op_s) == (11,15) + # @test domain_size(op_n) == (11,15) - @test range_size(op_w) == (15,) - @test range_size(op_e) == (15,) - @test range_size(op_s) == (11,) - @test range_size(op_n) == (11,) - end + # @test range_size(op_w) == (15,) + # @test range_size(op_e) == (15,) + # @test range_size(op_s) == (11,) + # @test range_size(op_n) == (11,) + # end end @testset "Application" begin @@ -76,43 +77,43 @@ @test (op_l'*u)[11] isa ComplexF64 end - @testset "2D" begin - v = rand(size(g_2D)...) - u = fill(3.124) - @test op_w*v ≈ 2*v[1,:] + v[2,:] + 3*v[3,:] rtol = 1e-14 - @test op_e*v ≈ 2*v[end,:] + v[end-1,:] + 3*v[end-2,:] rtol = 1e-14 - @test op_s*v ≈ 2*v[:,1] + v[:,2] + 3*v[:,3] rtol = 1e-14 - @test op_n*v ≈ 2*v[:,end] + v[:,end-1] + 3*v[:,end-2] rtol = 1e-14 + # @testset "2D" begin + # v = rand(size(g_2D)...) + # u = fill(3.124) + # @test op_w*v ≈ 2*v[1,:] + v[2,:] + 3*v[3,:] rtol = 1e-14 + # @test op_e*v ≈ 2*v[end,:] + v[end-1,:] + 3*v[end-2,:] rtol = 1e-14 + # @test op_s*v ≈ 2*v[:,1] + v[:,2] + 3*v[:,3] rtol = 1e-14 + # @test op_n*v ≈ 2*v[:,end] + v[:,end-1] + 3*v[:,end-2] rtol = 1e-14 - g_x = rand(size(g_2D)[1]) - g_y = rand(size(g_2D)[2]) + # g_x = rand(size(g_2D)[1]) + # g_y = rand(size(g_2D)[2]) - G_w = zeros(Float64, size(g_2D)...) - G_w[1,:] = 2*g_y - G_w[2,:] = g_y - G_w[3,:] = 3*g_y + # G_w = zeros(Float64, size(g_2D)...) + # G_w[1,:] = 2*g_y + # G_w[2,:] = g_y + # G_w[3,:] = 3*g_y - G_e = zeros(Float64, size(g_2D)...) - G_e[end,:] = 2*g_y - G_e[end-1,:] = g_y - G_e[end-2,:] = 3*g_y + # G_e = zeros(Float64, size(g_2D)...) + # G_e[end,:] = 2*g_y + # G_e[end-1,:] = g_y + # G_e[end-2,:] = 3*g_y - G_s = zeros(Float64, size(g_2D)...) - G_s[:,1] = 2*g_x - G_s[:,2] = g_x - G_s[:,3] = 3*g_x + # G_s = zeros(Float64, size(g_2D)...) + # G_s[:,1] = 2*g_x + # G_s[:,2] = g_x + # G_s[:,3] = 3*g_x - G_n = zeros(Float64, size(g_2D)...) - G_n[:,end] = 2*g_x - G_n[:,end-1] = g_x - G_n[:,end-2] = 3*g_x + # G_n = zeros(Float64, size(g_2D)...) + # G_n[:,end] = 2*g_x + # G_n[:,end-1] = g_x + # G_n[:,end-2] = 3*g_x - @test op_w'*g_y == G_w - @test op_e'*g_y == G_e - @test op_s'*g_x == G_s - @test op_n'*g_x == G_n - end + # @test op_w'*g_y == G_w + # @test op_e'*g_y == G_e + # @test op_s'*g_x == G_s + # @test op_n'*g_x == G_n + # end @testset "Regions" begin u = fill(3.124)
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Mon Mar 21 12:51:39 2022 +0100 @@ -77,6 +77,9 @@ @test D₁*v ≈ v rtol=tol end + + + # TODO: Different directions end end
--- a/test/SbpOperators/volumeops/volume_operator_test.jl Mon Mar 21 10:04:15 2022 +0100 +++ b/test/SbpOperators/volumeops/volume_operator_test.jl Mon Mar 21 12:51:39 2022 +0100 @@ -7,7 +7,6 @@ import Sbplib.SbpOperators.Stencil import Sbplib.SbpOperators.VolumeOperator -import Sbplib.SbpOperators.volume_operator import Sbplib.SbpOperators.odd import Sbplib.SbpOperators.even @@ -15,112 +14,68 @@ 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.)) @testset "Constructors" begin - @testset "1D" begin - op = VolumeOperator(inner_stencil,closure_stencils,(11,),even) - @test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,even) - @test op == volume_operator(g_1D,inner_stencil,closure_stencils,even,1) - @test op isa LazyTensor{T,1,1} where T - end - @testset "2D" begin - op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2) - Ix = IdentityTensor{Float64}((11,)) - Iy = IdentityTensor{Float64}((12,)) - @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy - @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even) - @test op_x isa LazyTensor{T,2,2} where T - @test op_y isa LazyTensor{T,2,2} where T - end - @testset "3D" begin - op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2) - op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3) - Ix = IdentityTensor{Float64}((11,)) - Iy = IdentityTensor{Float64}((12,)) - Iz = IdentityTensor{Float64}((10,)) - @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy⊗Iz - @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)⊗Iz - @test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),even) - @test op_x isa LazyTensor{T,3,3} where T - @test op_y isa LazyTensor{T,3,3} where T - @test op_z isa LazyTensor{T,3,3} where T - end + op = VolumeOperator(inner_stencil,closure_stencils,(11,),even) + @test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,even) + @test op isa LazyTensor{T,1,1} where T end @testset "Sizes" begin @testset "1D" begin - op = volume_operator(g_1D,inner_stencil,closure_stencils,even,1) + op = VolumeOperator(g_1D,inner_stencil,closure_stencils,even) @test range_size(op) == domain_size(op) == size(g_1D) end - - @testset "2D" begin - op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2) - @test range_size(op_y) == domain_size(op_y) == - range_size(op_x) == domain_size(op_x) == size(g_2D) - end - @testset "3D" begin - op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2) - op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3) - @test range_size(op_z) == domain_size(op_z) == - range_size(op_y) == domain_size(op_y) == - range_size(op_x) == domain_size(op_x) == size(g_3D) - end end - op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_2D,inner_stencil,closure_stencils,odd,2) - v = zeros(size(g_2D)) - Nx = size(g_2D)[1] - Ny = size(g_2D)[2] - for i = 1:Nx - v[i,:] .= i - end - rx = copy(v) - rx[1,:] .= 1.5 - rx[Nx,:] .= (2*Nx-1)/2 - ry = copy(v) - ry[:,Ny-1:Ny] = -v[:,Ny-1:Ny] + # op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1) + # op_y = volume_operator(g_2D,inner_stencil,closure_stencils,odd,2) + # v = zeros(size(g_2D)) + # Nx = size(g_2D)[1] + # Ny = size(g_2D)[2] + # for i = 1:Nx + # v[i,:] .= i + # end + # rx = copy(v) + # rx[1,:] .= 1.5 + # rx[Nx,:] .= (2*Nx-1)/2 + # ry = copy(v) + # ry[:,Ny-1:Ny] = -v[:,Ny-1:Ny] - @testset "Application" begin - @test op_x*v ≈ rx rtol = 1e-14 - @test op_y*v ≈ ry rtol = 1e-14 + # @testset "Application" begin + # @test op_x*v ≈ rx rtol = 1e-14 + # @test op_y*v ≈ ry rtol = 1e-14 - @test (op_x*rand(ComplexF64,size(g_2D)))[2,2] isa ComplexF64 - end + # @test (op_x*rand(ComplexF64,size(g_2D)))[2,2] isa ComplexF64 + # end - @testset "Regions" begin - @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14 - @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14 - @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14 - @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14 - @test (op_x*v)[Index(11,Upper),Index(3,Interior)] ≈ rx[11,3] rtol = 1e-14 + # @testset "Regions" begin + # @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14 + # @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14 + # @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14 + # @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14 + # @test (op_x*v)[Index(11,Upper),Index(3,Interior)] ≈ rx[11,3] rtol = 1e-14 - @test_throws BoundsError (op_x*v)[Index(3,Lower),Index(3,Interior)] - @test_throws BoundsError (op_x*v)[Index(9,Upper),Index(3,Interior)] + # @test_throws BoundsError (op_x*v)[Index(3,Lower),Index(3,Interior)] + # @test_throws BoundsError (op_x*v)[Index(9,Upper),Index(3,Interior)] - @test (op_y*v)[Index(3,Interior),Index(1,Lower)] ≈ ry[3,1] rtol = 1e-14 - @test (op_y*v)[Index(3,Interior),Index(2,Lower)] ≈ ry[3,2] rtol = 1e-14 - @test (op_y*v)[Index(3,Interior),Index(6,Interior)] ≈ ry[3,6] rtol = 1e-14 - @test (op_y*v)[Index(3,Interior),Index(11,Upper)] ≈ ry[3,11] rtol = 1e-14 - @test (op_y*v)[Index(3,Interior),Index(12,Upper)] ≈ ry[3,12] rtol = 1e-14 + # @test (op_y*v)[Index(3,Interior),Index(1,Lower)] ≈ ry[3,1] rtol = 1e-14 + # @test (op_y*v)[Index(3,Interior),Index(2,Lower)] ≈ ry[3,2] rtol = 1e-14 + # @test (op_y*v)[Index(3,Interior),Index(6,Interior)] ≈ ry[3,6] rtol = 1e-14 + # @test (op_y*v)[Index(3,Interior),Index(11,Upper)] ≈ ry[3,11] rtol = 1e-14 + # @test (op_y*v)[Index(3,Interior),Index(12,Upper)] ≈ ry[3,12] rtol = 1e-14 - @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)] - @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)] - end + # @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)] + # @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)] + # end - @testset "Inferred" begin - @test_skip @inferred apply(op_x, v,1,1) - @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower)) - @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower)) - @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower)) - @test_skip @inferred apply(op_y, v,1,1) - @inferred apply(op_y, v, Index(1,Lower),Index(1,Lower)) - @inferred apply(op_y, v, Index(1,Lower),Index(6,Interior)) - @inferred apply(op_y, v, Index(1,Lower),Index(11,Upper)) - end + # @testset "Inferred" begin + # @test_skip @inferred apply(op_x, v,1,1) + # @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower)) + # @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower)) + # @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower)) + # @test_skip @inferred apply(op_y, v,1,1) + # @inferred apply(op_y, v, Index(1,Lower),Index(1,Lower)) + # @inferred apply(op_y, v, Index(1,Lower),Index(6,Interior)) + # @inferred apply(op_y, v, Index(1,Lower),Index(11,Upper)) + # end end