Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/constantlaplace.jl @ 286:7247e85dc1e8 tensor_mappings

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 22 Jun 2020 19:47:20 +0200
parents
children dd621017b695
comparison
equal deleted inserted replaced
285:e21dcda55163 286:7247e85dc1e8
1 #TODO: Naming?! What is this? It is a 1D tensor operator but what is then the
2 # potentially multi-D laplace tensor mapping then?
3 # Ideally I would like the below to be the laplace operator in 1D, while the
4 # multi-D operator is a a tuple of the 1D-operator. Possible via recursive
5 # definitions? Or just bad design?
6 """
7 ConstantLaplaceOperator{T<:Real,N,M,K} <: TensorOperator{T,1}
8 Implements the Laplace tensor operator `L` with constant grid spacing and coefficients
9 in 1D dimension
10 """
11 struct ConstantLaplaceOperator{T<:Real,N,M,K} <: TensorOperator{T,1}
12 h_inv::T # The grid spacing could be included in the stencil already. Preferable?
13 a::T # TODO: Better name?
14 innerStencil::Stencil{T,N}
15 closureStencils::NTuple{M,Stencil{T,K}}
16 parity::Parity
17 end
18
19 @enum Parity begin
20 odd = -1
21 even = 1
22 end
23
24 LazyTensors.domain_size(L::ConstantLaplaceOperator, range_size::NTuple{1,Integer}) = range_size
25
26 function LazyTensors.apply(L::ConstantLaplaceOperator{T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
27 return L.a*apply_2nd_derivative(L, L.h_inv, v, I[1])
28 end
29
30 # Apply for different regions Lower/Interior/Upper or Unknown region
31 @inline function apply_2nd_derivative(L::ConstantLaplaceOperator, h_inv::Real, v::AbstractVector, i::Index{Lower})
32 return @inbounds h_inv*h_inv*apply_stencil(L.closureStencils[Int(i)], v, Int(i))
33 end
34
35 @inline function apply_2nd_derivative(L::ConstantLaplaceOperator, h_inv::Real, v::AbstractVector, i::Index{Interior})
36 return @inbounds h_inv*h_inv*apply_stencil(L.innerStencil, v, Int(i))
37 end
38
39 @inline function apply_2nd_derivative(L::ConstantLaplaceOperator, h_inv::Real, v::AbstractVector, i::Index{Upper})
40 N = length(v) # Can we use range_size here instead?
41 return @inbounds h_inv*h_inv*Int(L.parity)*apply_stencil_backwards(L.closureStencils[N-Int(i)+1], v, Int(i))
42 end
43
44 @inline function apply_2nd_derivative(L::ConstantLaplaceOperator, h_inv::Real, v::AbstractVector, index::Index{Unknown})
45 N = length(v) # Can we use range_size here instead?
46 r = getregion(Int(index), closuresize(L), N)
47 i = Index(Int(index), r)
48 return apply_2nd_derivative(op, h_inv, v, i)
49 end
50
51 function closuresize(L::ConstantLaplaceOperator{T<:Real,N,M,K})::Int
52 return M
53 end