Mercurial > repos > public > sbplib_julia
changeset 592:4781e759d92f refactor/toml_operator_format
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 02 Dec 2020 14:20:24 +0100 |
| parents | 089d4cb65146 (current diff) 0c6f8331c190 (diff) |
| children | fa03dae0ff0b |
| files | |
| diffstat | 5 files changed, 221 insertions(+), 189 deletions(-) [+] |
line wrap: on
line diff
--- a/src/Grids/Grids.jl Wed Dec 02 14:19:37 2020 +0100 +++ b/src/Grids/Grids.jl Wed Dec 02 14:20:24 2020 +0100 @@ -7,7 +7,7 @@ abstract type BoundaryIdentifier end struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim -region(::CartesianBoundary{Dim, R}) where {Dim, R} = R #TODO: Should return R() +region(::CartesianBoundary{Dim, R}) where {Dim, R} = R() export dim, region
--- a/src/SbpOperators/BoundaryValue.jl Wed Dec 02 14:19:37 2020 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,127 +0,0 @@ -""" - BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} - -Implements the boundary operator `e` as a TensorMapping -""" -struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} - eClosure::Stencil{T,M} - bId::CartesianBoundary -end -export BoundaryValue - -# TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue? -# Can we give special treatment to TensorMappings that go to a higher dim? -function LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T - if dim(e.bId) == 1 - return (UnknownDim, domain_size[1]) - elseif dim(e.bId) == 2 - return (domain_size[1], UnknownDim) - end -end -LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],) -# TODO: Make a nicer solution for 3-dim(e.bId) - -# TODO: Make this independent of dimension -function LazyTensors.apply(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T - i = I[dim(e.bId)] - j = I[3-dim(e.bId)] - N_i = size(e.grid)[dim(e.bId)] - return apply_boundary_value(e.op, v[j], i, N_i, region(e.bId)) -end - -function LazyTensors.apply_transpose(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T - u = selectdim(v,3-dim(e.bId),Int(I[1])) - return apply_boundary_value_transpose(e.op, u, region(e.bId)) -end - -function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower}) - @boundscheck if length(v) < closuresize(op) - throw(BoundsError()) - end - apply_stencil(op.eClosure,v,1) -end - -function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Upper}) - @boundscheck if length(v) < closuresize(op) - throw(BoundsError()) - end - apply_stencil_backwards(op.eClosure,v,length(v)) -end -export apply_boundary_value_transpose - -function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index, N::Integer, ::Type{Lower}) - @boundscheck if !(0<length(Int(i)) <= N) - throw(BoundsError()) - end - op.eClosure[Int(i)-1]*v -end - -function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index, N::Integer, ::Type{Upper}) - @boundscheck if !(0<length(Int(i)) <= N) - throw(BoundsError()) - end - op.eClosure[N-Int(i)]*v -end -export apply_boundary_value - - -""" - BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} - -Implements the boundary operator `e` as a TensorMapping -""" -struct BoundaryValue{D,T,M,R} <: TensorMapping{T,D,1} - e:BoundaryOperator{T,M,R} - bId::CartesianBoundary -end - -function LazyTensors.apply_transpose(bv::BoundaryValue{T,M,Lower}, v::AbstractVector{T}, i::Index) where T - u = selectdim(v,3-dim(bv.bId),Int(I[1])) - return apply_transpose(bv.e, u, I) -end - - -""" - BoundaryOperator{T,N,R} <: TensorMapping{T,1,1} - -Implements the boundary operator `e` as a TensorMapping -""" -export BoundaryOperator -struct BoundaryOperator{T,M,R<:Region} <: TensorMapping{T,1,1} - closure::Stencil{T,M} -end - -function LazyTensors.range_size(e::BoundaryOperator, domain_size::NTuple{1,Integer}) - return UnknownDim -end - -LazyTensors.domain_size(e::BoundaryOperator{T}, range_size::NTuple{1,Integer}) where T = range_size - -function LazyTensors.apply_transpose(e::BoundaryOperator{T,M,Lower}, v::AbstractVector{T}, i::Index{Lower}) where T - @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here? - throw(BoundsError()) - end - apply_stencil(e.closure,v,Int(i)) -end - -function LazyTensors.apply_transpose(e::BoundaryOperator{T,M,Upper}}, v::AbstractVector{T}, i::Index{Upper}) where T - @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here? - throw(BoundsError()) - end - apply_stencil_backwards(e.closure,v,Int(i)) -end - -function LazyTensors.apply_transpose(e::BoundaryOperator{T}, v::AbstractVector{T}, i::Index) where T - @boundscheck if length(v) < closuresize(e) #TODO: Use domain_size here? - throw(BoundsError()) - end - return eltype(v)(0) -end - -#TODO: Implement apply in a meaningful way. Should it return a vector or a single value (perferable?) Should fit into the -function LazyTensors.apply(e::BoundaryOperator, v::AbstractVector, i::Index) - @boundscheck if !(0<length(Int(i)) <= length(v)) - throw(BoundsError()) - end - return e.closure[Int(i)].*v -end
--- a/src/SbpOperators/SbpOperators.jl Wed Dec 02 14:19:37 2020 +0100 +++ b/src/SbpOperators/SbpOperators.jl Wed Dec 02 14:20:24 2020 +0100 @@ -14,5 +14,6 @@ include("quadrature/quadrature.jl") include("quadrature/inverse_diagonal_inner_product.jl") include("quadrature/inverse_quadrature.jl") +include("boundaryops/boundary_restriction.jl") end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Wed Dec 02 14:20:24 2020 +0100 @@ -0,0 +1,81 @@ +""" + 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) + +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 + +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 + +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
--- a/test/testSbpOperators.jl Wed Dec 02 14:19:37 2020 +0100 +++ b/test/testSbpOperators.jl Wed Dec 02 14:20:24 2020 +0100 @@ -180,67 +180,144 @@ @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_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 + + 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 + end + + @testset "2D" begin + e_w = boundary_restriction(g_2D,op.eClosure,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}()) + + 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}()) + + @testset "Sizes" begin + @testset "1D" begin + @test domain_size(e_l) == (11,) + @test domain_size(e_r) == (11,) + + @test range_size(e_l) == () + @test range_size(e_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 range_size(e_w) == (15,) + @test range_size(e_e) == (15,) + @test range_size(e_s) == (11,) + @test range_size(e_n) == (11,) + end + end + + + @testset "Application" begin + @testset "1D" begin + 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] + @test e_l'*u == [u[]; zeros(10)] + @test e_r'*u == [zeros(10); u[]] + end + + @testset "2D" begin + v = rand(11, 15) + u = fill(3.124) + + @test e_w*v == v[1,:] + @test e_e*v == v[end,:] + @test e_s*v == v[:,1] + @test e_n*v == v[:,end] + + + g_x = rand(11) + g_y = rand(15) + + G_w = zeros(Float64, (11,15)) + G_w[1,:] = g_y + + G_e = zeros(Float64, (11,15)) + G_e[end,:] = g_y + + G_s = zeros(Float64, (11,15)) + G_s[:,1] = g_x + + G_n = zeros(Float64, (11,15)) + G_n[:,end] = g_x + + @test e_w'*g_y == G_w + @test e_e'*g_y == G_e + @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 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")