Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/laplace/secondderivative.jl @ 311:f2d6ec89dfc5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make apply dispatch on Index instead of Tuples of Index
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 09 Sep 2020 21:16:13 +0200
parents 27a0bca5e1f2
children 9cc5d1498b2d
comparison
equal deleted inserted replaced
310:ece3f6f8a1d4 311:f2d6ec89dfc5
1 """ 1 """
2 SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1} 2 SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1}
3 Implements the Laplace tensor operator `L` with constant grid spacing and coefficients 3 Implements the Laplace tensor operator `L` with constant grid spacing and coefficients
4 in 1D dimension 4 in 1D dimension
5 """ 5 """
6 struct SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1} 6
7 struct SecondDerivative{T,N,M,K} <: TensorOperator{T,1}
7 h_inv::T # The grid spacing could be included in the stencil already. Preferable? 8 h_inv::T # The grid spacing could be included in the stencil already. Preferable?
8 innerStencil::Stencil{T,N} 9 innerStencil::Stencil{T,N}
9 closureStencils::NTuple{M,Stencil{T,K}} 10 closureStencils::NTuple{M,Stencil{T,K}}
10 parity::Parity 11 parity::Parity
11 #TODO: Write a nice constructor 12 #TODO: Write a nice constructor
12 end 13 end
13 14 export SecondDerivative
14 @enum Parity begin
15 odd = -1
16 even = 1
17 end
18 15
19 LazyTensors.domain_size(D2::SecondDerivative, range_size::NTuple{1,Integer}) = range_size 16 LazyTensors.domain_size(D2::SecondDerivative, range_size::NTuple{1,Integer}) = range_size
20 17
21 @inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T 18 #TODO: The 1D tensor mappings should not have to dispatch on 1D tuples if we write LazyTensor.apply for vararg right?!?!
22 return @inbounds apply(D2, v, I[1]) 19 # Currently have to index the Tuple{Index} in each method in order to call the stencil methods which is ugly.
20 # I thought I::Vararg{Index,R} fell back to just Index for R = 1
21
22 # Apply for different regions Lower/Interior/Upper or Unknown region
23 @inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Lower}) where T
24 return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(I)], v, Int(I))
23 end 25 end
24 26
25 function LazyTensors.apply_transpose(D2::SecondDerivative{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T = LazyTensors.apply(D2, v, I) 27 @inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Interior}) where T
26 28 return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(I))
27 # Apply for different regions Lower/Interior/Upper or Unknown region
28 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Lower})
29 return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(i)], v, Int(i))
30 end 29 end
31 30
32 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Interior}) 31 @inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Upper}) where T
33 return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(i)) 32 N = length(v) # TODO: Use domain_size here instead? N = domain_size(D2,size(v))
33 return @inbounds D2.h_inv*D2.h_inv*Int(D2.parity)*apply_stencil_backwards(D2.closureStencils[N-Int(I)+1], v, Int(I))
34 end 34 end
35 35
36 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, i::Index{Upper}) 36 @inline function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, index::Index{Unknown}) where T
37 N = length(v) # TODO: Use domain_size here instead? 37 N = length(v) # TODO: Use domain_size here instead?
38 return @inbounds D2.h_inv*D2.h_inv*Int(D2.parity)*apply_stencil_backwards(D2.closureStencils[N-Int(i)+1], v, Int(i)) 38 r = getregion(Int(index), closuresize(D2), N)
39 I = Index(Int(index), r)
40 return LazyTensors.apply(D2, v, I)
39 end 41 end
40 42
41 @inline function LazyTensors.apply(D2::SecondDerivative, v::AbstractVector, index::Index{Unknown}) 43
42 N = length(v) # TODO: Use domain_size here instead? 44 @inline function LazyTensors.apply_transpose(D2::SecondDerivative, v::AbstractVector, I::Index)
43 r = getregion(Int(index), closuresize(L), N) 45 return LazyTensors.apply(D2, v, I)
44 i = Index(Int(index), r)
45 return apply(D2, v, i)
46 end 46 end
47 47
48 function closuresize(D2::SecondDerivative{T<:Real,N,M,K}) where {T,N,M,K} 48
49 function closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K}
49 return M 50 return M
50 end 51 end