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