Mercurial > repos > public > sbplib_julia

changeset 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 c4cc2a1a68b9
files StencilIndex.jl diffOp.jl sbpD2.jl stencil.jl
diffstat 4 files changed, 87 insertions(+), 51 deletions(-) [+]
line wrap: on
line diff
--- a/StencilIndex.jl	Mon Feb 04 09:13:48 2019 +0100
+++ b/StencilIndex.jl	Mon Feb 04 16:09:07 2019 +0100
@@ -1,45 +1,55 @@
-abstract type StencilIndex end
+abstract type Region end
+struct Interior <: Region end
+struct Lower    <: Region end
+struct Upper    <: Region end
 
-function Base.getindex(si::StencilIndex, i::Int)
-    return si.gridindex[i]
+struct StencilIndex{R<:Region, T<:Integer}
+    localindex::CartesianIndex
+    globalindex::T
+
+    StencilIndex{R}(li::CartesianIndex, gi::T) where {R<:Region,T<:Integer} = new{R, T}(li, gi)
+    StencilIndex(li::CartesianIndex, gi::T, ::Type{R}) where {R<:Region,T<:Integer} = StencilIndex{R}(li, gi)
+
+    # Index(t::Tuple{T, Type{R}}) where {R<:Region,T<:Integer} = Index{t[2]}(t[1])
+    # Above doesn't work, below does but is less type strict
+    #Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2]}(t[1])
 end
 
-struct LowerClosureIndex <: StencilIndex
-    globalindex::Integer
-    gridindex::CartesianIndex
-end
-
-struct UpperClosureIndex <: StencilIndex
-    globalindex::Integer
-    gridindex::CartesianIndex
-end
-
-struct InteriorIndex <: StencilIndex
-    globalindex::Integer
-    gridindex::CartesianIndex
+function Base.getindex(si::StencilIndex, i::Int)
+     return si.localindex[i]
 end
 
-# TODO: The design of StencilIndex is wrong. Use Jonatans design instead.
-# TODO: This should take a Stencil or DiffOp so that we can extract all the
-# indices in the closures.
+#Index(t::Vararg{Tuple{T, DataType}}) where T = Index.(t)
 # TODO: Where to place this function?
-function stencilindices(grid::Grid.EquidistantGrid)
-    lowerclosure = Vector{LowerClosureIndex}(undef, 0)
-    upperclosure = Vector{UpperClosureIndex}(undef, 0)
-    interior = Vector{InteriorIndex}(undef, 0)
-    # TODO: Fix such that the indices of the entire closure width is included.
-    islower = x -> (x == 1)
-    isupper = x -> (x in grid.numberOfPointsPerDim)
-    ci = CartesianIndices(grid.numberOfPointsPerDim)
-    for i ∈ 1:Grid.numberOfPoints(grid)
-        I = Tuple(ci[i])
-        if any(islower, I)
-            push!(lowerclosure, LowerClosureIndex(i,ci[i]))
-        elseif any(isupper, I)
-            push!(upperclosure, UpperClosureIndex(i,ci[i]))
-        else
-            push!(interior, InteriorIndex(i,ci[i]))
+
+function stencilindices(diffOp)
+    N = diffOp.grid.numberOfPointsPerDim
+
+    lowerclosure = Vector{Vector{StencilIndex{Lower, Int64}}}(undef,0)
+    upperclosure = Vector{Vector{StencilIndex{Upper, Int64}}}(undef,0)
+    interior = Vector{Vector{StencilIndex{Interior, Int64}}}(undef,0)
+    cSize = closureSize(diffOp.op)
+    ci = CartesianIndices(diffOp.grid.numberOfPointsPerDim)
+
+    # TODO: Loop over all points or one loop for each region?
+    for j = 1:Grid.numberOfDimensions(diffOp.grid)
+        templ = Vector{StencilIndex{Lower,Int64}}(undef, 0)
+        tempu = Vector{StencilIndex{Upper,Int64}}(undef, 0)
+        tempi = Vector{StencilIndex{Interior,Int64}}(undef, 0)
+        for i ∈ 1:Grid.numberOfPoints(diffOp.grid)
+            val = ci[i][j]
+            if val ∈ range(1; length=cSize)
+                push!(templ, StencilIndex{Lower}(ci[i],i))
+            elseif val ∈ range(N[j] - cSize+1, length=cSize)
+                push!(tempu, StencilIndex{Upper}(ci[i],i))
+            else
+                push!(tempi, StencilIndex{Interior}(ci[i],i))
+            end
         end
+        push!(lowerclosure,templ)
+        push!(upperclosure,tempu)
+        push!(interior,tempi)
     end
     return lowerclosure, upperclosure, interior
 end
+
--- a/diffOp.jl	Mon Feb 04 09:13:48 2019 +0100
+++ b/diffOp.jl	Mon Feb 04 16:09:07 2019 +0100
@@ -69,28 +69,55 @@
 end
 
 function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
-    lowerind, upperind, interiorind = stencilindices(L.grid)
-    for si ∈ lowerind
-       u[si.globalindex] = apply(L, v, si)
+    lowerind, upperind, interiorind = stencilindices(L)
+    # First x!
+    for si ∈ lowerind[1]
+       u[si.globalindex] = applyx(L, v, si)
+    end
+    for si ∈ upperind[1]
+       u[si.globalindex] = applyx(L, v, si)
     end
-    for si ∈ upperind
-       u[si.globalindex] = apply(L, v, si)
+    for si ∈ interiorind[1]
+        u[si.globalindex] = applyx(L, v, si)
     end
-    for si ∈ interiorind
-        u[si.globalindex] = apply(L, v, si)
+    # NOW y!
+    for si ∈ lowerind[2]
+       u[si.globalindex] += applyy(L, v, si)
+    end
+    for si ∈ upperind[2]
+       u[si.globalindex] += applyy(L, v, si)
+    end
+    for si ∈ interiorind[2]
+        u[si.globalindex] += applyy(L, v, si)
     end
     return nothing
 end
 
-@inline function apply(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
+@inline function applyx(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
     h = Grid.spacings(L.grid)
     li = LinearIndices(L.grid.numberOfPointsPerDim)
     # 2nd x-derivative
-    @inbounds vx = uview(v, uview(li,:,si.gridindex[2]))
-    uᵢ  = apply(L.op, h[1], vx , si.gridindex[1], si)
+    @inbounds vx = uview(v, uview(li,:,si.localindex[2]))
+    uᵢ  = apply(L.op, h[1], vx , si.localindex[1], si)
+    return uᵢ
+end
+
+@inline function applyy(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
+    h = Grid.spacings(L.grid)
+    li = LinearIndices(L.grid.numberOfPointsPerDim)
     # 2nd y-derivative
-    @inbounds vy = uview(v, uview(li,si.gridindex[1],:))
-    uᵢ += apply(L.op, h[2], vy, si.gridindex[2], si)
+    # @show typeof(si.localindex)
+    # @show si.localindex
+
+    # @show typeof(si.localindex[1])
+    # @show typeof(si.localindex[2])
+
+    # @show typeof(uview(li, si.localindex[1],:))
+    # @show uview(li,si.localindex[1],:)
+
+
+    @inbounds vy = uview(v, uview(li,si.localindex[1],:))
+    uᵢ = apply(L.op, h[2], vy, si.localindex[2], si)
     return uᵢ
 end
 
--- a/sbpD2.jl	Mon Feb 04 09:13:48 2019 +0100
+++ b/sbpD2.jl	Mon Feb 04 16:09:07 2019 +0100
@@ -15,15 +15,15 @@
     return uᵢ
 end
 
-Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::InteriorIndex)
+Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::StencilIndex{Interior})
     return apply(op.innerStencil, v,  i)/h^2
 end
 
-Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::LowerClosureIndex)
+Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::StencilIndex{Lower})
     return apply(op.closureStencils[i], v, i)/h^2
 end
 
-Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::UpperClosureIndex)
+Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::StencilIndex{Upper})
     N = length(v)
     return Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2 #TODO: Write an applybackwards instead?
 end
--- a/stencil.jl	Mon Feb 04 09:13:48 2019 +0100
+++ b/stencil.jl	Mon Feb 04 16:09:07 2019 +0100
@@ -19,7 +19,6 @@
 
 Base.@propagate_inbounds function apply(s::Stencil, v::AbstractVector, i::Int)
     w = zero(eltype(v))
-    @show s.range[1]:s.range[2]
     for j ∈ s.range[1]:s.range[2]
         @inbounds weight = s[j]
         w += weight*v[i+j]