changeset 1097:95464a1af340 feature/dissipation_operators

Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 10 May 2022 21:14:38 +0200
parents 9f0121e465a5 (current diff) 74c54996de6a (diff)
children 6f51160c7ca7
files src/LazyTensors/lazy_tensor_operations.jl test/LazyTensors/lazy_tensor_operations_test.jl
diffstat 8 files changed, 108 insertions(+), 2 deletions(-) [+]
line wrap: on
line diff
--- a/src/LazyTensors/LazyTensors.jl	Tue May 10 21:10:31 2022 +0200
+++ b/src/LazyTensors/LazyTensors.jl	Tue May 10 21:14:38 2022 +0200
@@ -3,9 +3,10 @@
 export TensorApplication
 export TensorTranspose
 export TensorComposition
-export DenseTensor
 export IdentityTensor
 export ScalingTensor
+export DiagonalTensor
+export DenseTensor
 export InflatedTensor
 export LazyOuterProduct
 export ⊗
--- a/src/LazyTensors/lazy_tensor_operations.jl	Tue May 10 21:10:31 2022 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Tue May 10 21:14:38 2022 +0200
@@ -98,7 +98,6 @@
     apply_transpose(c.t2, c.t1'*v, I...)
 end
 
-
 """
     TensorComposition(tm, tmi::IdentityTensor)
     TensorComposition(tmi::IdentityTensor, tm)
@@ -120,6 +119,8 @@
     return tmi
 end
 
+Base.:*(a::T, tm::LazyTensor{T}) where T = TensorComposition(ScalingTensor{T,range_dim(tm)}(a,range_size(tm)), tm)
+Base.:*(tm::LazyTensor{T}, a::T) where T = a*tm
 
 """
     InflatedTensor{T,R,D} <: LazyTensor{T,R,D}
--- a/src/LazyTensors/tensor_types.jl	Tue May 10 21:10:31 2022 +0200
+++ b/src/LazyTensors/tensor_types.jl	Tue May 10 21:14:38 2022 +0200
@@ -37,6 +37,24 @@
 
 
 """
+    DiagonalTensor{T,D,...} <: LazyTensor{T,D,D}
+    DiagonalTensor(a::AbstractArray)
+
+A lazy tensor with diagonal `a`.
+"""
+struct DiagonalTensor{T,D,AT<:AbstractArray{T,D}} <: LazyTensor{T,D,D}
+    diagonal::AT
+end
+
+range_size(tm::DiagonalTensor) = size(tm.diagonal)
+domain_size(tm::DiagonalTensor) = size(tm.diagonal)
+
+
+LazyTensors.apply(tm::DiagonalTensor{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = tm.diagonal[I...]*v[I...]
+LazyTensors.apply_transpose(tm::DiagonalTensor{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = tm.diagonal[I...]*v[I...]
+
+
+"""
     DenseTensor{T,R,D,...}(A, range_indicies, domain_indicies)
 
 LazyTensor defined by the AbstractArray A. `range_indicies` and `domain_indicies` define which indicies of A should
--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Tue May 10 21:10:31 2022 +0200
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Tue May 10 21:14:38 2022 +0200
@@ -181,6 +181,14 @@
 
     @test (Ã∘B̃*ComplexF64[1.,2.,3.,4.])[1] isa ComplexF64
     @test ((Ã∘B̃)'*ComplexF64[1.,2.])[1] isa ComplexF64
+
+    a = 2.
+    v = rand(3)
+    @test a*Ã isa TensorComposition
+    @test a*Ã == Ã*a
+    @test range_size(a*Ã) == range_size(Ã)
+    @test domain_size(a*Ã) == domain_size(Ã)
+    @test a*Ã*v == a.*A*v
 end
 
 
--- a/test/LazyTensors/tensor_types_test.jl	Tue May 10 21:10:31 2022 +0200
+++ b/test/LazyTensors/tensor_types_test.jl	Tue May 10 21:14:38 2022 +0200
@@ -1,5 +1,6 @@
 using Test
 using Sbplib.LazyTensors
+using BenchmarkTools
 
 @testset "IdentityTensor" begin
     @test IdentityTensor{Float64}((4,5)) isa IdentityTensor{T,2} where T
@@ -58,6 +59,47 @@
     @inferred (st'*v)[2,2]
 end
 
+@testset "DiagonalTensor" begin
+    @test DiagonalTensor([1,2,3,4]) isa LazyTensor{Int,1,1}
+    @test DiagonalTensor([1 2 3; 4 5 6]) isa LazyTensor{Int,2,2}
+    @test DiagonalTensor([1. 2. 3.; 4. 5. 6.]) isa LazyTensor{Float64,2,2}
+
+    @test range_size(DiagonalTensor([1,2,3,4])) == (4,)
+    @test domain_size(DiagonalTensor([1,2,3,4])) == (4,)
+
+    @test range_size(DiagonalTensor([1 2 3; 4 5 6])) == (2,3)
+    @test domain_size(DiagonalTensor([1 2 3; 4 5 6])) == (2,3)
+
+    @testset "apply size=$sz" for sz ∈ [(4,),(3,2),(3,4,2)]
+        diag = rand(sz...)
+        tm = DiagonalTensor(diag)
+
+        v = rand(sz...)
+
+        @test tm*v == diag.*v
+        @test tm'*v == diag.*v
+    end
+
+
+    @testset "allocations size=$sz" for sz ∈ [(4,),(3,2),(3,4,2)]
+        diag = rand(sz...)
+        tm = DiagonalTensor(diag)
+
+        v = rand(sz...)
+
+        @test tm*v == diag.*v
+        @test tm'*v == diag.*v
+    end
+
+    sz = (3,2)
+    diag = rand(sz...)
+    tm = DiagonalTensor(diag)
+
+    v = rand(sz...)
+    LazyTensors.apply(tm,v, 2,1)
+    @test (@ballocated LazyTensors.apply($tm,$v, 2,1)) == 0
+end
+
 
 @testset "DenseTensor" begin
     # Test a standard matrix-vector product
--- a/test/Manifest.toml	Tue May 10 21:10:31 2022 +0200
+++ b/test/Manifest.toml	Tue May 10 21:14:38 2022 +0200
@@ -12,6 +12,12 @@
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[deps.BenchmarkTools]]
+deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "1.3.1"
+
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
 git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
@@ -93,6 +99,12 @@
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
 version = "1.4.1"
 
+[[deps.JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.3"
+
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
@@ -161,6 +173,12 @@
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
+[[deps.Parsers]]
+deps = ["Dates"]
+git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "2.2.3"
+
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -175,6 +193,10 @@
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
+[[deps.Profile]]
+deps = ["Printf"]
+uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
+
 [[deps.REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
--- a/test/Project.toml	Tue May 10 21:10:31 2022 +0200
+++ b/test/Project.toml	Tue May 10 21:14:38 2022 +0200
@@ -1,4 +1,5 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue May 10 21:10:31 2022 +0200
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue May 10 21:14:38 2022 +0200
@@ -71,6 +71,7 @@
     end
 
     @testset "Accuracy on function" begin
+        # 1D
         g = EquidistantGrid(30, 0.,1.)
         v = evalOn(g, x->exp(x))
         @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)]
@@ -79,6 +80,18 @@
 
             @test D₁*v ≈ v rtol=tol
         end
+
+        # 2D
+        g = EquidistantGrid((30,60), (0.,0.),(1.,2.))
+        v = evalOn(g, (x,y)->exp(0.8x+1.2*y))
+        @testset for (order, tol) ∈ [(2, 6e-3),(4, 3e-4)]
+            stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
+            Dx = first_derivative(g, stencil_set, 1)
+            Dy = first_derivative(g, stencil_set, 2)
+
+            @test Dx*v ≈ 0.8v rtol=tol
+            @test Dy*v ≈ 1.2v rtol=tol
+        end
     end
 end