Mercurial > repos > public > sbplib_julia
changeset 610:e40e7439d1b4 feature/volume_and_boundary_operators
Add a general boundary operator and make BoundaryRestriction a specialization of it.
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Sat, 05 Dec 2020 18:12:31 +0100 |
parents | 8f9b3eac128a |
children | e71f2f81b5f8 |
files | src/SbpOperators/SbpOperators.jl src/SbpOperators/boundaryops/boundary_operator.jl src/SbpOperators/boundaryops/boundary_restriction.jl test/testSbpOperators.jl |
diffstat | 4 files changed, 255 insertions(+), 144 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/SbpOperators.jl Wed Dec 02 09:35:14 2020 +0100 +++ b/src/SbpOperators/SbpOperators.jl Sat Dec 05 18:12:31 2020 +0100 @@ -14,6 +14,7 @@ include("quadrature/quadrature.jl") include("quadrature/inverse_diagonal_inner_product.jl") include("quadrature/inverse_quadrature.jl") +include("boundaryops/boundary_operator.jl") include("boundaryops/boundary_restriction.jl") end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Sat Dec 05 18:12:31 2020 +0100 @@ -0,0 +1,89 @@ +""" + 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 `IdentityMappings` 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{Dim,T}, closure_stencil::Stencil{T}, boundary::CartesianBoundary) where {Dim,T} + # Create 1D boundary operator + r = region(boundary) + d = dim(boundary) + op = BoundaryOperator(restrict(grid, d), closure_stencil, r) + + # Create 1D IdentityMappings for each coordinate direction + one_d_grids = restrict.(Ref(grid), Tuple(1:Dim)) + Is = IdentityMapping{T}.(size.(one_d_grids)) + + # Formulate the correct outer product sequence of the identity mappings and + # the boundary operator + parts = Base.setindex(Is, op, d) + return foldl(⊗, parts) +end +export boundary_operator + +""" + BoundaryOperator{T,R,N} <: TensorMapping{T,0,1} + +Implements the boundary operator `op` for 1D as a `TensorMapping` + +`op` is the restriction of a grid function to the boundary using some closure `Stencil{T,N}`. +The boundary to restrict to is determined by `R`. +`op'` is the prolongation of a zero dimensional array to the whole grid using the same closure stencil. +""" +struct BoundaryOperator{T,R<:Region,N} <: TensorMapping{T,0,1} + stencil::Stencil{T,N} + size::Int +end +export BoundaryOperator + +BoundaryOperator{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryOperator{T,R,N}(stencil, size) + +""" + BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil, region) + +Constructs the BoundaryOperator with stencil `closure_stencil` for a one-dimensional `grid`, restricting to +to the boundary specified by `region`. +""" +function BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil::Stencil{T,N}, region::Region) where {T,N} + return BoundaryOperator{T,typeof(region),N}(closure_stencil,size(grid)[1]) +end + +""" + closure_size(::BoundaryOperator) +The size of the closure stencil. +""" +closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N + +LazyTensors.range_size(op::BoundaryOperator) = () +LazyTensors.domain_size(op::BoundaryOperator) = (op.size,) + +function LazyTensors.apply(op::BoundaryOperator{T,Lower}, v::AbstractVector{T}) where T + apply_stencil(op.stencil,v,1) +end + +function LazyTensors.apply(op::BoundaryOperator{T,Upper}, v::AbstractVector{T}) where T + apply_stencil_backwards(op.stencil,v,op.size) +end + +function LazyTensors.apply_transpose(op::BoundaryOperator{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T + return op.stencil[Int(i)-1]*v[] +end + +function LazyTensors.apply_transpose(op::BoundaryOperator{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T + return op.stencil[op.size[1] - Int(i)]*v[] +end + +# Catch all combinations of Lower, Upper and Interior not caught by the two previous methods. +function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i::Index) where T + return zero(T) +end + +function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i) where T + r = getregion(i, closure_size(op), op.size) + apply_transpose(op, v, Index(i,r)) +end
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl Wed Dec 02 09:35:14 2020 +0100 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Sat Dec 05 18:12:31 2020 +0100 @@ -1,81 +1,15 @@ -""" - boundary_restriction(grid,closureStencil,boundary) - -Creates a boundary restriction operator on a `Dim`-dimensional grid for the -specified `boundary`. - -When `Dim=1`, the corresponding `BoundaryRestriction` tensor mapping is returned. -When `Dim>1`, the `BoundaryRestriction` `e` is inflated by the outer product -of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`, -the boundary restriction operator in the y-direction direction is `Ix⊗e⊗Iz`. -""" -function boundary_restriction(grid::EquidistantGrid{Dim,T}, closureStencil::Stencil{T,M}, boundary::CartesianBoundary) where {Dim,T,M} - # Create 1D boundary restriction operator - r = region(boundary) - d = dim(boundary) - e = BoundaryRestriction(restrict(grid, d), closureStencil, r) - - # Create 1D IdentityMappings for each coordinate direction - one_d_grids = restrict.(Ref(grid), Tuple(1:Dim)) - Is = IdentityMapping{T}.(size.(one_d_grids)) - - # Formulate the correct outer product sequence of the identity mappings and - # the boundary restriction operator - parts = Base.setindex(Is, e, d) - return foldl(⊗, parts) -end - -export boundary_restriction - -""" - BoundaryRestriction{T,R,N} <: TensorMapping{T,0,1} - -Implements the boundary operator `e` for 1D as a `TensorMapping` - -`e` is the restriction of a grid function to the boundary using some `closureStencil`. -The boundary to restrict to is determined by `R`. - -`e'` is the prolongation of a zero dimensional array to the whole grid using the same `closureStencil`. """ -struct BoundaryRestriction{T,R<:Region,N} <: TensorMapping{T,0,1} - stencil::Stencil{T,N} - size::Int -end -export BoundaryRestriction - -BoundaryRestriction{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryRestriction{T,R,N}(stencil, size) + BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) + BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) -function BoundaryRestriction(grid::EquidistantGrid{1}, closureStencil::Stencil{T,N}, region::Region) where {T,N} - return BoundaryRestriction{T,typeof(region),N}(closureStencil,size(grid)[1]) -end - -closure_size(::BoundaryRestriction{T,R,N}) where {T,R,N} = N - -LazyTensors.range_size(e::BoundaryRestriction) = () -LazyTensors.domain_size(e::BoundaryRestriction) = (e.size,) - -function LazyTensors.apply(e::BoundaryRestriction{T,Lower}, v::AbstractVector{T}) where T - apply_stencil(e.stencil,v,1) -end +Creates the boundary restriction operator `e` as a `TensorMapping` -function LazyTensors.apply(e::BoundaryRestriction{T,Upper}, v::AbstractVector{T}) where T - apply_stencil_backwards(e.stencil,v,e.size) -end - -function LazyTensors.apply_transpose(e::BoundaryRestriction{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T - return e.stencil[Int(i)-1]*v[] -end +`e` is the restriction of a grid function to the boundary specified by `boundary` or `region` using some `closure_stencil`. +`e'` is the prolongation of a grid function on the boundary to the whole grid using the same `closure_stencil`. +On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of +a `BoundaryOperator`. Also see the documentation of `boundary_operator(...)` for more details. +""" +BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = boundary_operator(grid, closure_stencil, boundary) +BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = BoundaryRestriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}()) -function LazyTensors.apply_transpose(e::BoundaryRestriction{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T - return e.stencil[e.size[1] - Int(i)]*v[] -end - -# Catch all combinations of Lower, Upper and Interior not caught by the two previous methods. -function LazyTensors.apply_transpose(e::BoundaryRestriction{T}, v::AbstractArray{T,0}, i::Index) where T - return zero(T) -end - -function LazyTensors.apply_transpose(e::BoundaryRestriction{T}, v::AbstractArray{T,0}, i) where T - r = getregion(i, closure_size(e), e.size) - apply_transpose(e, v, Index(i,r)) -end +export BoundaryRestriction
--- a/test/testSbpOperators.jl Wed Dec 02 09:35:14 2020 +0100 +++ b/test/testSbpOperators.jl Sat Dec 05 18:12:31 2020 +0100 @@ -181,58 +181,58 @@ @test Qinv*v == Qinv'*v end -@testset "BoundaryRestrictrion" begin - op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") +@testset "BoundaryOperator" begin + closure_stencil = SbpOperators.Stencil((0,2), (2.,1.,3.)) g_1D = EquidistantGrid(11, 0.0, 1.0) g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) @testset "Constructors" begin @testset "1D" begin - e_l = BoundaryRestriction{Lower}(op.eClosure,size(g_1D)[1]) - @test e_l == BoundaryRestriction(g_1D,op.eClosure,Lower()) - @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}()) - @test e_l isa TensorMapping{T,0,1} where T + 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 TensorMapping{T,0,1} where T - e_r = BoundaryRestriction{Upper}(op.eClosure,size(g_1D)[1]) - @test e_r == BoundaryRestriction(g_1D,op.eClosure,Upper()) - @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}()) - @test e_r isa TensorMapping{T,0,1} where T + op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1]) + @test op_r == BoundaryRestriction(g_1D,closure_stencil,Upper()) + @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}()) + @test op_r isa TensorMapping{T,0,1} where T end @testset "2D" begin - e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}()) + e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}()) @test e_w isa InflatedTensorMapping @test e_w isa TensorMapping{T,1,2} where T end end - e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}()) - e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}()) + op_l = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Lower}()) + op_r = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Upper}()) - e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}()) - e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}()) - e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}()) - e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}()) + op_w = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Lower}()) + op_e = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Upper}()) + op_s = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Lower}()) + op_n = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Upper}()) @testset "Sizes" begin @testset "1D" begin - @test domain_size(e_l) == (11,) - @test domain_size(e_r) == (11,) + @test domain_size(op_l) == (11,) + @test domain_size(op_r) == (11,) - @test range_size(e_l) == () - @test range_size(e_r) == () + @test range_size(op_l) == () + @test range_size(op_r) == () end @testset "2D" begin - @test domain_size(e_w) == (11,15) - @test domain_size(e_e) == (11,15) - @test domain_size(e_s) == (11,15) - @test domain_size(e_n) == (11,15) + @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(e_w) == (15,) - @test range_size(e_e) == (15,) - @test range_size(e_s) == (11,) - @test range_size(e_n) == (11,) + @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 @@ -241,6 +241,126 @@ @testset "1D" begin v = evalOn(g_1D,x->1+x^2) u = fill(3.124) + @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3] + @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2] + @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2] + @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)] + @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]] + 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 + + + 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_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_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 + + @testset "Regions" begin + u = fill(3.124) + @test (op_l'*u)[Index(1,Lower)] == 2*u[] + @test (op_l'*u)[Index(2,Lower)] == u[] + @test (op_l'*u)[Index(6,Interior)] == 0 + @test (op_l'*u)[Index(10,Upper)] == 0 + @test (op_l'*u)[Index(11,Upper)] == 0 + + @test (op_r'*u)[Index(1,Lower)] == 0 + @test (op_r'*u)[Index(2,Lower)] == 0 + @test (op_r'*u)[Index(6,Interior)] == 0 + @test (op_r'*u)[Index(10,Upper)] == u[] + @test (op_r'*u)[Index(11,Upper)] == 2*u[] + end + end + + @testset "Inferred" begin + v = ones(Float64, 11) + u = fill(1.) + + @inferred apply(op_l, v) + @inferred apply(op_r, v) + + @inferred apply_transpose(op_l, u, 4) + @inferred apply_transpose(op_l, u, Index(1,Lower)) + @inferred apply_transpose(op_l, u, Index(2,Lower)) + @inferred apply_transpose(op_l, u, Index(6,Interior)) + @inferred apply_transpose(op_l, u, Index(10,Upper)) + @inferred apply_transpose(op_l, u, Index(11,Upper)) + + @inferred apply_transpose(op_r, u, 4) + @inferred apply_transpose(op_r, u, Index(1,Lower)) + @inferred apply_transpose(op_r, u, Index(2,Lower)) + @inferred apply_transpose(op_r, u, Index(6,Interior)) + @inferred apply_transpose(op_r, u, Index(10,Upper)) + @inferred apply_transpose(op_r, u, Index(11,Upper)) + end + +end + +@testset "BoundaryRestriction" begin + op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") + g_1D = EquidistantGrid(11, 0.0, 1.0) + g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) + + @testset "Constructors" begin + @testset "1D" begin + e_l = BoundaryRestriction(g_1D,op.eClosure,Lower()) + @test e_l == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}()) + @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower()) + @test e_l isa BoundaryOperator{T,Lower} where T + @test e_l isa TensorMapping{T,0,1} where T + + e_r = BoundaryRestriction(g_1D,op.eClosure,Upper()) + @test e_r == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}()) + @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper()) + @test e_r isa BoundaryOperator{T,Upper} where T + @test e_r isa TensorMapping{T,0,1} where T + end + + @testset "2D" begin + e_w = BoundaryRestriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}()) + @test e_w isa InflatedTensorMapping + @test e_w isa TensorMapping{T,1,2} where T + end + end + + @testset "Application" begin + @testset "1D" begin + e_l = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}()) + e_r = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}()) + + v = evalOn(g_1D,x->1+x^2) + u = fill(3.124) + @test (e_l*v)[] == v[1] @test (e_r*v)[] == v[end] @test (e_r*v)[1] == v[end] @@ -249,6 +369,11 @@ end @testset "2D" begin + e_w = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}()) + e_e = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}()) + e_s = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}()) + e_n = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}()) + v = rand(11, 15) u = fill(3.124) @@ -278,45 +403,7 @@ @test e_s'*g_x == G_s @test e_n'*g_x == G_n end - - @testset "Regions" begin - u = fill(3.124) - @test (e_l'*u)[Index(1,Lower)] == 3.124 - @test (e_l'*u)[Index(2,Lower)] == 0 - @test (e_l'*u)[Index(6,Interior)] == 0 - @test (e_l'*u)[Index(10,Upper)] == 0 - @test (e_l'*u)[Index(11,Upper)] == 0 - - @test (e_r'*u)[Index(1,Lower)] == 0 - @test (e_r'*u)[Index(2,Lower)] == 0 - @test (e_r'*u)[Index(6,Interior)] == 0 - @test (e_r'*u)[Index(10,Upper)] == 0 - @test (e_r'*u)[Index(11,Upper)] == 3.124 - end end - - @testset "Inferred" begin - v = ones(Float64, 11) - u = fill(1.) - - @inferred apply(e_l, v) - @inferred apply(e_r, v) - - @inferred apply_transpose(e_l, u, 4) - @inferred apply_transpose(e_l, u, Index(1,Lower)) - @inferred apply_transpose(e_l, u, Index(2,Lower)) - @inferred apply_transpose(e_l, u, Index(6,Interior)) - @inferred apply_transpose(e_l, u, Index(10,Upper)) - @inferred apply_transpose(e_l, u, Index(11,Upper)) - - @inferred apply_transpose(e_r, u, 4) - @inferred apply_transpose(e_r, u, Index(1,Lower)) - @inferred apply_transpose(e_r, u, Index(2,Lower)) - @inferred apply_transpose(e_r, u, Index(6,Interior)) - @inferred apply_transpose(e_r, u, Index(10,Upper)) - @inferred apply_transpose(e_r, u, Index(11,Upper)) - end - end # # @testset "NormalDerivative" begin