Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/laplace.jl @ 247:ed29ee13e92e boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Restructure laplace.jl
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 27 Jun 2019 09:43:44 +0200
parents a827568fc251
children 05e7bbe0af97
comparison
equal deleted inserted replaced
244:a827568fc251 247:ed29ee13e92e
1 """ 1 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
2 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} 2 grid::EquidistantGrid{Dim,T}
3 3 a::T
4 Implements the boundary operator `d` as a TensorMapping 4 op::D2{Float64,N,M,K}
5 """ 5 # e::BoundaryValue
6 struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} 6 # d::NormalDerivative
7 op::D2{T,N,M,K}
8 grid::EquidistantGrid
9 bId::CartesianBoundary
10 end
11 export NormalDerivative
12
13 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
14 # Can we give special treatment to TensorMappings that go to a higher dim?
15 LazyTensors.range_size(e::NormalDerivative{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
16 LazyTensors.domain_size(e::NormalDerivative{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
17
18 # Not correct abstraction level
19 # TODO: Not type stable D:<
20 function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int})
21 i = I[dim(d.bId)]
22 j = I[3-dim(d.bId)]
23 N_i = size(d.grid)[dim(d.bId)]
24
25 if region(d.bId) == Lower
26 # Note, closures are indexed by offset. Fix this D:<
27 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j]
28 elseif region(d.bId) == Upper
29 return -d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-i]*v[j]
30 end
31 end 7 end
32 8
33 function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int}) 9 function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
34 u = selectdim(v,3-dim(d.bId),I[1]) 10 error("not implemented")
35 return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
36 end 11 end
12
13 # u = L*v
14 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
15 uᵢ = L.a * SbpOperators.apply(L.op, L.grid.spacing[1], v, i)
16 return uᵢ
17 end
18
19 @inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
20 # 2nd x-derivative
21 @inbounds vx = view(v, :, Int(I[2]))
22 @inbounds uᵢ = L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[1], vx , I[1])
23 # 2nd y-derivative
24 @inbounds vy = view(v, Int(I[1]), :)
25 @inbounds uᵢ += L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
26 # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors
27 return uᵢ
28 end
29
30 # Slow but maybe convenient?
31 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})
32 I = Index{Unknown}.(Tuple(i))
33 apply(L, v, I)
34 end
35
37 36
38 37
39 """ 38 """
40 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} 39 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
41 40
71 return apply_e(e.op, u, region(e.bId)) 70 return apply_e(e.op, u, region(e.bId))
72 end 71 end
73 72
74 73
75 74
76 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim} 75 """
77 grid::EquidistantGrid{Dim,T} 76 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
78 a::T 77
79 op::D2{Float64,N,M,K} 78 Implements the boundary operator `d` as a TensorMapping
80 # e::BoundaryValue 79 """
81 # d::NormalDerivative 80 struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
81 op::D2{T,N,M,K}
82 grid::EquidistantGrid
83 bId::CartesianBoundary
84 end
85 export NormalDerivative
86
87 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
88 # Can we give special treatment to TensorMappings that go to a higher dim?
89 LazyTensors.range_size(e::NormalDerivative{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
90 LazyTensors.domain_size(e::NormalDerivative{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
91
92 # Not correct abstraction level
93 # TODO: Not type stable D:<
94 function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int})
95 i = I[dim(d.bId)]
96 j = I[3-dim(d.bId)]
97 N_i = size(d.grid)[dim(d.bId)]
98
99 if region(d.bId) == Lower
100 # Note, closures are indexed by offset. Fix this D:<
101 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j]
102 elseif region(d.bId) == Upper
103 return -d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-i]*v[j]
104 end
82 end 105 end
83 106
84 function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim 107 function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
85 error("not implemented") 108 u = selectdim(v,3-dim(d.bId),I[1])
109 return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
86 end 110 end
87 111
88 # u = L*v
89 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
90 uᵢ = L.a * SbpOperators.apply(L.op, L.grid.spacing[1], v, i)
91 return uᵢ
92 end
93
94 @inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
95 # 2nd x-derivative
96 @inbounds vx = view(v, :, Int(I[2]))
97 @inbounds uᵢ = L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[1], vx , I[1])
98 # 2nd y-derivative
99 @inbounds vy = view(v, Int(I[1]), :)
100 @inbounds uᵢ += L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
101 # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors
102 return uᵢ
103 end
104
105 # Slow but maybe convenient?
106 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})
107 I = Index{Unknown}.(Tuple(i))
108 apply(L, v, I)
109 end
110 112
111 113
112 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end 114 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
113 115
114 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} 116 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}