Mercurial > repos > public > sbplib_julia
changeset 103:a274d6384e91 cell_based_test
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 | 1862e901febb |
| files | diffOp.jl index.jl |
| diffstat | 2 files changed, 44 insertions(+), 5 deletions(-) [+] |
line wrap: on
line diff
--- a/diffOp.jl Wed Feb 06 09:09:55 2019 +0100 +++ b/diffOp.jl Thu Feb 07 17:15:26 2019 +0100 @@ -67,18 +67,35 @@ return uᵢ end +function apply!(L::Laplace{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T + regions = (Lower, Interior, Upper) + for r1 ∈ regions + for r2 ∈ regions + apply!(L, u, v, r1, r2) + end + end + return nothing +end + +function apply!(L::Laplace{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::R1, r2::R2) where {T, R1, R2} + N = L.grid.numberOfPointsPerDim; + closuresize = closureSize(L.op); + for I ∈ regionindices(N, closuresize, (r1,r2)) + @inbounds indextuple = (Index(I[1], r1), Index(I[2], r2)) + @inbounds u[I] = apply(L, v, indextuple) + end + return nothing +end + using UnsafeArrays - function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2} h = Grid.spacings(L.grid) - # 2nd x-derivative @inbounds vx = uview(v, :, Int(I[2])) - @inbounds uᵢ = apply(L.op, h[1], vx , I[1]) + @inbounds uᵢ = L.a*apply(L.op, h[1], vx , I[1]) # 2nd y-derivative @inbounds vy = uview(v, Int(I[1]), :) - @inbounds uᵢ += apply(L.op, h[2], vy, I[2]) - + @inbounds uᵢ += L.a*apply(L.op, h[2], vy, I[2]) return uᵢ end
--- a/index.jl Wed Feb 06 09:09:55 2019 +0100 +++ b/index.jl Thu Feb 07 17:15:26 2019 +0100 @@ -38,3 +38,25 @@ end IndexTuple(t::Vararg{Tuple{T, DataType}}) where T<:Integer = Index.(t) + +function regionindices(gridsize::NTuple{Dim,Integer}, closuresize::Integer, region::NTuple{Dim,DataType}) where Dim + return regionindices(gridsize, ntuple(x->closuresize,Dim), region) +end + +function regionindices(gridsize::NTuple{Dim,Integer}, closuresize::NTuple{Dim,Integer}, region::NTuple{Dim,DataType}) where Dim + regions = map(getunitrange,gridsize,closuresize,region) + return CartesianIndices(regions) +end + +function getunitrange(gridsize::Integer, closuresize::Integer, region::R) where R + if region == Lower + r = 1:closuresize + elseif region == Interior + r = (closuresize+1):(gridsize - closuresize) + elseif region == Upper + r = (gridsize - closuresize + 1):gridsize + else + error("Unspecified region") + end + return r +end