Mercurial > repos > public > sbplib_julia
changeset 793:1a4b3aecf9e5 operator_storage_array_of_table
Add ConstantInteriorScalingOperator which will be used to implement diagonal inner products
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Thu, 22 Jul 2021 23:33:58 +0200 |
parents | 5b6628a288a7 |
children | e0c5eeb3cbbb |
files | src/SbpOperators/SbpOperators.jl src/SbpOperators/volumeops/constant_interior_scaling_operator.jl test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl |
diffstat | 3 files changed, 80 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/SbpOperators.jl Mon Jul 19 22:27:52 2021 +0200 +++ b/src/SbpOperators/SbpOperators.jl Thu Jul 22 23:33:58 2021 +0200 @@ -8,6 +8,7 @@ include("d2.jl") include("readoperator.jl") include("volumeops/volume_operator.jl") +include("volumeops/constant_interior_scaling_operator.jl") include("volumeops/derivatives/second_derivative.jl") include("volumeops/laplace/laplace.jl") include("volumeops/inner_products/inner_product.jl")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Thu Jul 22 23:33:58 2021 +0200 @@ -0,0 +1,50 @@ +export ConstantInteriorScalingOperator + +""" + ConstantInteriorScalingOperator{T,N} <: TensorMapping{T,1,1} + +A one-dimensional operator scaling a vector. The first and last `N` points are +scaled with individual weights while all interior points are scaled the same. +""" +struct ConstantInteriorScalingOperator{T,N} <: TensorMapping{T,1,1} + interior_weight::T + closure_weights::NTuple{N,T} + size::Int + + function ConstantInteriorScalingOperator(interior_weight::T, closure_weights::NTuple{N,T}, size::Int) where {T,N} + if size < 2length(closure_weights) + throw(DomainError(size, "size must be larger that two times the closure size.")) + end + + return new{T,N}(interior_weight, closure_weights, size) + end +end + +function ConstantInteriorScalingOperator(grid::EquidistantGrid{1}, interior_weight, closure_weights) + return ConstantInteriorScalingOperator(interior_weight, Tuple(closure_weights), size(grid)) +end + +closure_size(::ConstantInteriorScalingOperator{T,N}) where {T,N} = N + +LazyTensors.range_size(op::ConstantInteriorScalingOperator) = (op.size,) +LazyTensors.domain_size(op::ConstantInteriorScalingOperator) = (op.size,) + +# TBD: @inbounds in apply methods? +function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T + return op.closure_weights[Int(i)]*v[Int(i)] +end + +function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T + return op.interior_weight*v[Int(i)] +end + +function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T + return op.closure_weights[op.size[1]-Int(i)+1]*v[Int(i)] +end + +function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i) where T + r = getregion(i, closure_size(op), op.size[1]) + return LazyTensors.apply(op, v, Index(i, r)) +end + +LazyTensors.apply_transpose(op::ConstantInteriorScalingOperator, v, i) = apply(op, v, i)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl Thu Jul 22 23:33:58 2021 +0200 @@ -0,0 +1,29 @@ +using Test + +using Sbplib.LazyTensors +using Sbplib.SbpOperators + +@testset "ConstantInteriorScalingOperator" begin + @test ConstantInteriorScalingOperator(1, (2,3), 10) isa ConstantInteriorScalingOperator{Int,2} + @test ConstantInteriorScalingOperator(1., (2.,3.), 10) isa ConstantInteriorScalingOperator{Float64,2} + + a = ConstantInteriorScalingOperator(4, (2,3), 10) + v = ones(Int, 10) + @test a*v == [2,3,4,4,4,4,4,4,3,2] + @test a'*v == [2,3,4,4,4,4,4,4,3,2] + + @test range_size(a) == (10,) + @test domain_size(a) == (10,) + + + a = ConstantInteriorScalingOperator(.5, (.1,.2), 7) + v = ones(7) + + @test a*v == [.1,.2,.5,.5,.5,.2,.1] + @test a'*v == [.1,.2,.5,.5,.5,.2,.1] + + @test range_size(a) == (7,) + @test domain_size(a) == (7,) + + @test_throws DomainError ConstantInteriorScalingOperator(4,(2,3), 3) +end