Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 103:a274d6384e91 cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 07 Feb 2019 17:15:26 +0100
parents 6b6d680f2e25
children 44cd6b4371de
comparison
equal deleted inserted replaced
99:6b6d680f2e25 103:a274d6384e91
65 h = Grid.spacings(L.grid)[1] 65 h = Grid.spacings(L.grid)[1]
66 uᵢ = L.a * apply(L.op, h, v, i) 66 uᵢ = L.a * apply(L.op, h, v, i)
67 return uᵢ 67 return uᵢ
68 end 68 end
69 69
70 function apply!(L::Laplace{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
71 regions = (Lower, Interior, Upper)
72 for r1 ∈ regions
73 for r2 ∈ regions
74 apply!(L, u, v, r1, r2)
75 end
76 end
77 return nothing
78 end
79
80 function apply!(L::Laplace{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::R1, r2::R2) where {T, R1, R2}
81 N = L.grid.numberOfPointsPerDim;
82 closuresize = closureSize(L.op);
83 for I ∈ regionindices(N, closuresize, (r1,r2))
84 @inbounds indextuple = (Index(I[1], r1), Index(I[2], r2))
85 @inbounds u[I] = apply(L, v, indextuple)
86 end
87 return nothing
88 end
89
70 using UnsafeArrays 90 using UnsafeArrays
71
72 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2} 91 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
73 h = Grid.spacings(L.grid) 92 h = Grid.spacings(L.grid)
74
75 # 2nd x-derivative 93 # 2nd x-derivative
76 @inbounds vx = uview(v, :, Int(I[2])) 94 @inbounds vx = uview(v, :, Int(I[2]))
77 @inbounds uᵢ = apply(L.op, h[1], vx , I[1]) 95 @inbounds uᵢ = L.a*apply(L.op, h[1], vx , I[1])
78 # 2nd y-derivative 96 # 2nd y-derivative
79 @inbounds vy = uview(v, Int(I[1]), :) 97 @inbounds vy = uview(v, Int(I[1]), :)
80 @inbounds uᵢ += apply(L.op, h[2], vy, I[2]) 98 @inbounds uᵢ += L.a*apply(L.op, h[2], vy, I[2])
81
82 return uᵢ 99 return uᵢ
83 end 100 end
84 101
85 # Slow but maybe convenient? 102 # Slow but maybe convenient?
86 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2}) 103 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})