Mercurial > repos > public > sbplib_julia

changeset 513:547639572208 feature/boundary_ops

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Get some kind of tested working implementation.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 23 Nov 2020 20:22:14 +0100
parents 5a8cfcc0765d
children 14e722e8607d
files src/SbpOperators/boundaryops/boundary_restriction.jl test/testSbpOperators.jl
diffstat 2 files changed, 100 insertions(+), 66 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Nov 23 13:11:19 2020 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Nov 23 20:22:14 2020 +0100
@@ -2,11 +2,11 @@
     r = region(boundary)
     d = dim(boundary)
     d_orth = 3-d # orthogonal dimension
-    e = BoundaryRestriction(restrict(grid, d), closureStencil, r)
-    I = IdentityMapping(size(restrict(g,d_orth)))
-    if r == Lower
+    e = BoundaryRestriction(restrict(grid, d), closureStencil, r())
+    I = IdentityMapping{T}(size(restrict(grid,d_orth)))
+    if d == 1
         return e⊗I
-    elseif r == Upper
+    elseif d == 2
         return I⊗e
     else
         # throw error
@@ -26,7 +26,7 @@
 export BoundaryRestriction
 
 function BoundaryRestriction(grid::EquidistantGrid{1,T}, closureStencil::Stencil{T,M}, region::Region) where {T,M}
-    return BoundaryRestriction{T,M,typeof(region)}(stencil,size(grid))
+    return BoundaryRestriction{T,M,typeof(region)}(closureStencil,size(grid))
 end
 
 LazyTensors.range_size(e::BoundaryRestriction) = (1,)
@@ -56,3 +56,12 @@
     end
     return e.stencil[Int(i)-1]*v[1]
 end
+
+" Transpose of a restriction is an inflation or prolongation.
+  Inflates the scalar (1-element) vector to a vector of size of the grid"
+function LazyTensors.apply_transpose(e::BoundaryRestriction{T,M,Upper}, v::AbstractVector{T}, i) where {T,M}
+    @boundscheck if !(0 < Int(i) <= e.size[1])
+        throw(BoundsError())
+    end
+    return e.stencil[e.size[1] - Int(i)]*v[1]
+end
--- a/test/testSbpOperators.jl	Mon Nov 23 13:11:19 2020 +0100
+++ b/test/testSbpOperators.jl	Mon Nov 23 20:22:14 2020 +0100
@@ -172,67 +172,92 @@
     @test_broken Qinv*(Q*v) ≈ v
     @test Qinv*v == Qinv'*v
 end
-#
-# @testset "BoundaryValue" begin
-#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-#     g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
-#
-#     e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
-#     e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
-#     e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
-#     e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
-#
-#     v = zeros(Float64, 4, 5)
-#     v[:,5] = [1, 2, 3,4]
-#     v[:,4] = [1, 2, 3,4]
-#     v[:,3] = [4, 5, 6, 7]
-#     v[:,2] = [7, 8, 9, 10]
-#     v[:,1] = [10, 11, 12, 13]
-#
-#     @test e_w  isa TensorMapping{T,2,1} where T
-#     @test e_w' isa TensorMapping{T,1,2} where T
-#
-#     @test domain_size(e_w, (3,2)) == (2,)
-#     @test domain_size(e_e, (3,2)) == (2,)
-#     @test domain_size(e_s, (3,2)) == (3,)
-#     @test domain_size(e_n, (3,2)) == (3,)
-#
-#     @test size(e_w'*v) == (5,)
-#     @test size(e_e'*v) == (5,)
-#     @test size(e_s'*v) == (4,)
-#     @test size(e_n'*v) == (4,)
-#
-#     @test e_w'*v == [10,7,4,1.0,1]
-#     @test e_e'*v == [13,10,7,4,4.0]
-#     @test e_s'*v == [10,11,12,13.0]
-#     @test e_n'*v == [1,2,3,4.0]
-#
-#     g_x = [1,2,3,4.0]
-#     g_y = [5,4,3,2,1.0]
-#
-#     G_w = zeros(Float64, (4,5))
-#     G_w[1,:] = g_y
-#
-#     G_e = zeros(Float64, (4,5))
-#     G_e[4,:] = g_y
-#
-#     G_s = zeros(Float64, (4,5))
-#     G_s[:,1] = g_x
-#
-#     G_n = zeros(Float64, (4,5))
-#     G_n[:,5] = g_x
-#
-#     @test size(e_w*g_y) == (UnknownDim,5)
-#     @test size(e_e*g_y) == (UnknownDim,5)
-#     @test size(e_s*g_x) == (4,UnknownDim)
-#     @test size(e_n*g_x) == (4,UnknownDim)
-#
-#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
-#     @test_broken e_w*g_y == G_w
-#     @test_broken e_e*g_y == G_e
-#     @test_broken e_s*g_x == G_s
-#     @test_broken e_n*g_x == G_n
-# end
+
+@testset "BoundaryRestrictrion" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+
+    g = EquidistantGrid(4, 0.0, 1.0)
+
+    e_l = BoundaryRestriction(g,op.eClosure,Lower())
+    e_r = BoundaryRestriction(g,op.eClosure,Upper())
+
+    v = evalOn(g,x->1+x^2)
+    u = [3.124] #How to handle scalars having to be arrays? It's kind of ugly.
+
+    e_l*v isa LazyTensorMappingApplication
+    @test (e_l*v)[Index{Lower}(1)] == v[1]
+    @test (e_r*v)[Index{Upper}(4)] == v[end]
+    @test e_l'*u == [u[1], 0, 0, 0]
+    @test e_r'*u == [0, 0, 0, u[1]]
+    @test_throws BoundsError (e_l*v)[Index{Lower}(3)]
+    @test_throws BoundsError (e_r*v)[Index{Upper}(3)]
+
+
+
+
+
+    g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
+
+    e_w = boundary_restriction(g, op.eClosure, CartesianBoundary{1,Lower}())
+    e_e = boundary_restriction(g, op.eClosure, CartesianBoundary{1,Upper}())
+    e_s = boundary_restriction(g, op.eClosure, CartesianBoundary{2,Lower}())
+    e_n = boundary_restriction(g, op.eClosure, CartesianBoundary{2,Upper}())
+
+    v = zeros(Float64, 4, 5)
+    v[:,5] = [1, 2, 3,4]
+    v[:,4] = [1, 2, 3,4]
+    v[:,3] = [4, 5, 6, 7]
+    v[:,2] = [7, 8, 9, 10]
+    v[:,1] = [10, 11, 12, 13]
+
+    @test_broken e_w isa TensorMapping{T,1,2} where T
+    @test_broken e_w'  isa TensorMapping{T,2,1} where T
+
+
+
+    @test domain_size(e_w) == (4,5)
+    @test domain_size(e_e) == (4,5)
+    @test domain_size(e_s) == (4,5)
+    @test domain_size(e_n) == (4,5)
+
+    @test range_size(e_w) == (1,5)
+    @test range_size(e_e) == (1,5)
+    @test range_size(e_s) == (4,1)
+    @test range_size(e_n) == (4,1)
+
+    e_w*v isa LazyTensorMappingApplication
+
+    @test_broken e_w'*v == [10,7,4,1.0,1]
+    @test_broken e_e'*v == [13,10,7,4,4.0]
+    @test_broken e_s'*v == [10,11,12,13.0]
+    @test_broken e_n'*v == [1,2,3,4.0]
+
+    g_x = [1,2,3,4.0]
+    g_y = [5,4,3,2,1.0]
+
+    G_w = zeros(Float64, (4,5))
+    G_w[1,:] = g_y
+
+    G_e = zeros(Float64, (4,5))
+    G_e[4,:] = g_y
+
+    G_s = zeros(Float64, (4,5))
+    G_s[:,1] = g_x
+
+    G_n = zeros(Float64, (4,5))
+    G_n[:,5] = g_x
+
+    @test_broken size(e_w*g_y) == (UnknownDim,5)
+    @test_broken size(e_e*g_y) == (UnknownDim,5)
+    @test_broken size(e_s*g_x) == (4,UnknownDim)
+    @test_broken size(e_n*g_x) == (4,UnknownDim)
+
+    # These tests should be moved to where they are possible (i.e we know what the grid should be)
+    @test_broken e_w*g_y == G_w
+    @test_broken e_e*g_y == G_e
+    @test_broken e_s*g_x == G_s
+    @test_broken e_n*g_x == G_n
+end
 #
 # @testset "NormalDerivative" begin
 #     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")