Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 124:631eb9b35d72 cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make grid spacing a property of EquidistantGrid. Create plotting module for sbplib.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 13 Feb 2019 10:37:52 +0100
parents 4c0c02a80cd4
children 1aaeb46ba5f4
comparison
equal deleted inserted replaced
122:6c6979ff17f4 124:631eb9b35d72
57 return nothing 57 return nothing
58 end 58 end
59 59
60 # Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable 60 # Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable
61 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T 61 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
62 N = D.grid.numberOfPointsPerDim 62 for I ∈ regionindices(D.grid.size, closureSize(D.op), (r1,r2))
63 closuresize = closureSize(D.op)
64 for I ∈ regionindices(N, closuresize, (r1,r2))
65 @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2])) 63 @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
66 @inbounds u[I] = apply(D, v, indextuple) 64 @inbounds u[I] = apply(D, v, indextuple)
67 end 65 end
68 return nothing 66 return nothing
69 end 67 end
81 return nothing 79 return nothing
82 end 80 end
83 81
84 using TiledIteration 82 using TiledIteration
85 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T 83 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
86 N = D.grid.numberOfPointsPerDim 84 ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2))
87 closuresize = closureSize(D.op)
88 ri = regionindices(N, closuresize, (r1,r2))
89
90 for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2)) # TBD: Is this the right way, the right size? 85 for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2)) # TBD: Is this the right way, the right size?
91 for j ∈ tileaxs[2], i ∈ tileaxs[1] 86 for j ∈ tileaxs[2], i ∈ tileaxs[1]
92 I = ri[i,j] 87 I = ri[i,j]
93 u[i,j] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2]))) 88 u[i,j] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
94 end 89 end
112 error("not implemented") 107 error("not implemented")
113 end 108 end
114 109
115 # u = L*v 110 # u = L*v
116 function apply(L::Laplace{1}, v::AbstractVector, i::Int) 111 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
117 h = Grid.spacings(L.grid)[1] 112 uᵢ = L.a * apply(L.op, L.grid.spacing[1], v, i)
118 uᵢ = L.a * apply(L.op, h, v, i)
119 return uᵢ 113 return uᵢ
120 end 114 end
121 115
122 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2} 116 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
123 h = Grid.spacings(L.grid)
124 # 2nd x-derivative 117 # 2nd x-derivative
125 @inbounds vx = view(v, :, Int(I[2])) 118 @inbounds vx = view(v, :, Int(I[2]))
126 @inbounds uᵢ = L.a*apply(L.op, h[1], vx , I[1]) 119 @inbounds uᵢ = L.a*apply(L.op, L.grid.spacing[1], vx , I[1])
127 # 2nd y-derivative 120 # 2nd y-derivative
128 @inbounds vy = view(v, Int(I[1]), :) 121 @inbounds vy = view(v, Int(I[1]), :)
129 @inbounds uᵢ += L.a*apply(L.op, h[2], vy, I[2]) 122 @inbounds uᵢ += L.a*apply(L.op, L.grid.spacing[2], vy, I[2])
130 return uᵢ 123 return uᵢ
131 end 124 end
132 125
133 # Slow but maybe convenient? 126 # Slow but maybe convenient?
134 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2}) 127 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})