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