Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 94:84b1ad5a3755 stencil_index

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Made everything work(?) but also go really slow. And also not type-stable.
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 04 Feb 2019 16:09:07 +0100
parents 93df72e2b135
children
comparison
equal deleted inserted replaced
93:93df72e2b135 94:84b1ad5a3755
67 uᵢ = L.a * apply(L.op, h, v, i) 67 uᵢ = L.a * apply(L.op, h, v, i)
68 return uᵢ 68 return uᵢ
69 end 69 end
70 70
71 function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector) 71 function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
72 lowerind, upperind, interiorind = stencilindices(L.grid) 72 lowerind, upperind, interiorind = stencilindices(L)
73 for si ∈ lowerind 73 # First x!
74 u[si.globalindex] = apply(L, v, si) 74 for si ∈ lowerind[1]
75 u[si.globalindex] = applyx(L, v, si)
75 end 76 end
76 for si ∈ upperind 77 for si ∈ upperind[1]
77 u[si.globalindex] = apply(L, v, si) 78 u[si.globalindex] = applyx(L, v, si)
78 end 79 end
79 for si ∈ interiorind 80 for si ∈ interiorind[1]
80 u[si.globalindex] = apply(L, v, si) 81 u[si.globalindex] = applyx(L, v, si)
82 end
83 # NOW y!
84 for si ∈ lowerind[2]
85 u[si.globalindex] += applyy(L, v, si)
86 end
87 for si ∈ upperind[2]
88 u[si.globalindex] += applyy(L, v, si)
89 end
90 for si ∈ interiorind[2]
91 u[si.globalindex] += applyy(L, v, si)
81 end 92 end
82 return nothing 93 return nothing
83 end 94 end
84 95
85 @inline function apply(L::Laplace{2}, v::AbstractVector, si::StencilIndex) 96 @inline function applyx(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
86 h = Grid.spacings(L.grid) 97 h = Grid.spacings(L.grid)
87 li = LinearIndices(L.grid.numberOfPointsPerDim) 98 li = LinearIndices(L.grid.numberOfPointsPerDim)
88 # 2nd x-derivative 99 # 2nd x-derivative
89 @inbounds vx = uview(v, uview(li,:,si.gridindex[2])) 100 @inbounds vx = uview(v, uview(li,:,si.localindex[2]))
90 uᵢ = apply(L.op, h[1], vx , si.gridindex[1], si) 101 uᵢ = apply(L.op, h[1], vx , si.localindex[1], si)
102 return uᵢ
103 end
104
105 @inline function applyy(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
106 h = Grid.spacings(L.grid)
107 li = LinearIndices(L.grid.numberOfPointsPerDim)
91 # 2nd y-derivative 108 # 2nd y-derivative
92 @inbounds vy = uview(v, uview(li,si.gridindex[1],:)) 109 # @show typeof(si.localindex)
93 uᵢ += apply(L.op, h[2], vy, si.gridindex[2], si) 110 # @show si.localindex
111
112 # @show typeof(si.localindex[1])
113 # @show typeof(si.localindex[2])
114
115 # @show typeof(uview(li, si.localindex[1],:))
116 # @show uview(li,si.localindex[1],:)
117
118
119 @inbounds vy = uview(v, uview(li,si.localindex[1],:))
120 uᵢ = apply(L.op, h[2], vy, si.localindex[2], si)
94 return uᵢ 121 return uᵢ
95 end 122 end
96 123
97 function apply(L::Laplace{2}, v::AbstractVector, i::Int) 124 function apply(L::Laplace{2}, v::AbstractVector, i::Int)
98 h = Grid.spacings(L.grid) 125 h = Grid.spacings(L.grid)