Mercurial > repos > public > sbplib_julia

changeset 219:69a6049e14d9 package_refactor

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create package SbpOperators
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2019 13:12:38 +0200
parents 03375aa30edd
children 5bba704a89d8
files Manifest.toml Project.toml SbpOperators/Manifest.toml SbpOperators/Project.toml SbpOperators/src/SbpOperators.jl SbpOperators/src/d2_2nd.txt SbpOperators/src/d2_4th.txt SbpOperators/src/h_2nd.txt SbpOperators/src/h_4th.txt SbpOperators/src/stencil.jl SbpOperators/test/runtests.jl d2_2nd.txt d2_4th.txt h_2nd.txt h_4th.txt sbpD2.jl stencil.jl
diffstat 17 files changed, 316 insertions(+), 234 deletions(-) [+]
line wrap: on
line diff
--- a/Manifest.toml	Wed Jun 26 12:57:20 2019 +0200
+++ b/Manifest.toml	Wed Jun 26 13:12:38 2019 +0200
@@ -1,21 +1,71 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
 [[DiffOps]]
+deps = ["RegionIndices", "TiledIteration"]
 path = "DiffOps"
 uuid = "39474f48-97ec-11e9-01fc-6ddcbe5918df"
 version = "0.1.0"
 
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
 [[Grids]]
+deps = ["RegionIndices"]
 path = "Grids"
 uuid = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
 version = "0.1.0"
 
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
 [[LazyTensors]]
 path = "LazyTensors"
 uuid = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
 version = "0.1.0"
 
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[OffsetArrays]]
+git-tree-sha1 = "1af2f79c7eaac3e019a0de41ef63335ff26a0a57"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "0.11.1"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
 [[RegionIndices]]
 path = "RegionIndices"
 uuid = "5d527584-97f1-11e9-084c-4540c7ecf219"
 version = "0.1.0"
+
+[[SbpOperators]]
+path = "SbpOperators"
+uuid = "204357d8-97fd-11e9-05e7-010897a14cd0"
+version = "0.1.0"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[TiledIteration]]
+deps = ["OffsetArrays", "Test"]
+git-tree-sha1 = "58f6f07d3b54a363ec283a8f5fc9fb4ecebde656"
+uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
+version = "0.2.3"
--- a/Project.toml	Wed Jun 26 12:57:20 2019 +0200
+++ b/Project.toml	Wed Jun 26 13:12:38 2019 +0200
@@ -3,3 +3,4 @@
 Grids = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
 LazyTensors = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
 RegionIndices = "5d527584-97f1-11e9-084c-4540c7ecf219"
+SbpOperators = "204357d8-97fd-11e9-05e7-010897a14cd0"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/Manifest.toml	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,6 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[RegionIndices]]
+path = "../RegionIndices"
+uuid = "5d527584-97f1-11e9-084c-4540c7ecf219"
+version = "0.1.0"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/Project.toml	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,13 @@
+name = "SbpOperators"
+uuid = "204357d8-97fd-11e9-05e7-010897a14cd0"
+authors = ["Jonatan Werpers <jonatan.werpers@it.uu.se>"]
+version = "0.1.0"
+
+[deps]
+RegionIndices = "5d527584-97f1-11e9-084c-4540c7ecf219"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/SbpOperators.jl	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,161 @@
+module SbpOperators
+
+using RegionIndices
+
+include("stencil.jl")
+
+abstract type ConstantStencilOperator end
+
+# Apply for different regions Lower/Interior/Upper or Unknown region
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
+    return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i))
+end
+
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
+    return @inbounds h*h*apply(op.innerStencil, v, Int(i))
+end
+
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
+    N = length(v)
+    return @inbounds h*h*Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
+end
+
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})
+    cSize = closureSize(op)
+    N = length(v)
+
+    i = Int(index)
+
+    if 0 < i <= cSize
+        return apply(op, h, v, Index{Lower}(i))
+    elseif cSize < i <= N-cSize
+        return apply(op, h, v, Index{Interior}(i))
+    elseif N-cSize < i <= N
+        return apply(op, h, v, Index{Upper}(i))
+    else
+        error("Bounds error") # TODO: Make this more standard
+    end
+end
+
+
+# Wrapper functions for using regular indecies without specifying regions
+@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
+    return apply(op, h, v, Index{Unknown}(i))
+end
+
+@enum Parity begin
+    odd = -1
+    even = 1
+end
+
+struct D2{T,N,M,K} <: ConstantStencilOperator
+    quadratureClosure::NTuple{M,T}
+    innerStencil::Stencil{T,N}
+    closureStencils::NTuple{M,Stencil{T,K}}
+    eClosure::Stencil{T,M}
+    dClosure::Stencil{T,M}
+    parity::Parity
+end
+
+function closureSize(D::D2)::Int
+    return length(D.quadratureClosure)
+end
+
+function readOperator(D2fn, Hfn)
+    d = readSectionedFile(D2fn)
+    h = readSectionedFile(Hfn)
+
+    # Create inner stencil
+    innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
+    width = length(innerStencilWeights)
+    r = (-div(width,2), div(width,2))
+
+    innerStencil = Stencil(r, innerStencilWeights)
+
+    # Create boundary stencils
+    boundarySize = length(d["boundary_stencils"])
+    closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
+
+    for i ∈ 1:boundarySize
+        stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
+        width = length(stencilWeights)
+        r = (1-i,width-i)
+        closureStencils = (closureStencils..., Stencil(r, stencilWeights))
+    end
+
+    quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
+    eClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize))
+    dClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize))
+
+    d2 = D2(
+        quadratureClosure,
+        innerStencil,
+        closureStencils,
+        eClosure,
+        dClosure,
+        even
+    )
+
+    return d2
+end
+
+
+function apply_e(op::D2, v::AbstractVector, ::Type{Lower})
+    apply(op.eClosure,v,1)
+end
+
+function apply_e(op::D2, v::AbstractVector, ::Type{Upper})
+    apply(flip(op.eClosure),v,length(v))
+end
+
+
+function apply_d(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower})
+    -h_inv*apply(op.dClosure,v,1)
+end
+
+function apply_d(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper})
+    -h_inv*apply(flip(op.dClosure),v,length(v))
+end
+
+function readSectionedFile(filename)::Dict{String, Vector{String}}
+    f = open(filename)
+    sections = Dict{String, Vector{String}}()
+    currentKey = ""
+
+    for ln ∈ eachline(f)
+        if ln == "" || ln[1] == '#' # Skip comments and empty lines
+            continue
+        end
+
+        if isletter(ln[1]) # Found start of new section
+            if ~haskey(sections, ln)
+                sections[ln] =  Vector{String}()
+            end
+            currentKey = ln
+            continue
+        end
+
+        push!(sections[currentKey], ln)
+    end
+
+    return sections
+end
+
+function stringToTuple(T::DataType, s::String)
+    return Tuple(stringToVector(T,s))
+end
+
+function stringToVector(T::DataType, s::String)
+    return T.(eval.(Meta.parse.(split(s))))
+end
+
+
+function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
+    if N >= n
+        return t
+    else
+        return pad_tuple((t..., zero(T)), n)
+    end
+end
+
+end # module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/d2_2nd.txt	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,13 @@
+# D2 order 2
+
+boundary_stencils
+1 -2 1
+
+inner_stencil
+1 -2 1
+
+e
+1
+
+d1
+-3/2 2 -1/2
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/d2_4th.txt	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,16 @@
+# D2 order 4
+
+boundary_stencils
+     2    -5       4      -1     0     0
+     1    -2       1       0     0     0
+ -4/43 59/43 -110/43   59/43 -4/43     0
+ -1/49     0   59/49 -118/49 64/49 -4/49
+
+inner_stencil
+-1/12     4/3  -5/2   4/3 -1/12
+
+e
+1
+
+d1
+-11/6 3 -3/2 1/3
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/h_2nd.txt	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,7 @@
+# H order 2
+
+closure
+1/2
+
+inner_stencil
+1
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/h_4th.txt	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,7 @@
+# H order 4
+
+closure
+17/48 59/48 43/48 49/48
+
+inner_stencil
+1
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/src/stencil.jl	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,38 @@
+struct Stencil{T<:Real,N}
+    range::Tuple{Int,Int}
+    weights::NTuple{N,T}
+
+    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T <: Real, N}
+        @assert range[2]-range[1]+1 == N
+        new{T,N}(range,weights)
+    end
+end
+
+function flip(s::Stencil)
+    range = (-s.range[2], -s.range[1])
+    return Stencil(range, reverse(s.weights))
+end
+
+# Provides index into the Stencil based on offset for the root element
+@inline function Base.getindex(s::Stencil, i::Int)
+    @boundscheck if i < s.range[1] || s.range[2] < i
+        return eltype(s.weights)(0)
+    end
+    return s.weights[1 + i - s.range[1]]
+end
+
+Base.@propagate_inbounds @inline function apply(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
+    w = s.weights[1]*v[i + s.range[1]]
+    @simd for k ∈ 2:N
+        w += s.weights[k]*v[i + s.range[1] + k-1]
+    end
+    return w
+end
+
+Base.@propagate_inbounds @inline function apply_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
+    w = s.weights[N]*v[i - s.range[2]]
+    @simd for k ∈ N-1:-1:1
+        w += s.weights[k]*v[i - s.range[1] - k + 1]
+    end
+    return w
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SbpOperators/test/runtests.jl	Wed Jun 26 13:12:38 2019 +0200
@@ -0,0 +1,4 @@
+using SbpOperators
+using Test
+
+@test_broken false
--- a/d2_2nd.txt	Wed Jun 26 12:57:20 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,13 +0,0 @@
-# D2 order 2
-
-boundary_stencils
-1 -2 1
-
-inner_stencil
-1 -2 1
-
-e
-1
-
-d1
--3/2 2 -1/2
--- a/d2_4th.txt	Wed Jun 26 12:57:20 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,16 +0,0 @@
-# D2 order 4
-
-boundary_stencils
-     2    -5       4      -1     0     0
-     1    -2       1       0     0     0
- -4/43 59/43 -110/43   59/43 -4/43     0
- -1/49     0   59/49 -118/49 64/49 -4/49
-
-inner_stencil
--1/12     4/3  -5/2   4/3 -1/12
-
-e
-1
-
-d1
--11/6 3 -3/2 1/3
\ No newline at end of file
--- a/h_2nd.txt	Wed Jun 26 12:57:20 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-# H order 2
-
-closure
-1/2
-
-inner_stencil
-1
--- a/h_4th.txt	Wed Jun 26 12:57:20 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-# H order 4
-
-closure
-17/48 59/48 43/48 49/48
-
-inner_stencil
-1
\ No newline at end of file
--- a/sbpD2.jl	Wed Jun 26 12:57:20 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,153 +0,0 @@
-abstract type ConstantStencilOperator end
-
-# Apply for different regions Lower/Interior/Upper or Unknown region
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Lower})
-    return @inbounds h*h*apply(op.closureStencils[Int(i)], v, Int(i))
-end
-
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Interior})
-    return @inbounds h*h*apply(op.innerStencil, v, Int(i))
-end
-
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Index{Upper})
-    N = length(v)
-    return @inbounds h*h*Int(op.parity)*apply_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
-end
-
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, index::Index{Unknown})
-    cSize = closureSize(op)
-    N = length(v)
-
-    i = Int(index)
-
-    if 0 < i <= cSize
-        return apply(op, h, v, Index{Lower}(i))
-    elseif cSize < i <= N-cSize
-        return apply(op, h, v, Index{Interior}(i))
-    elseif N-cSize < i <= N
-        return apply(op, h, v, Index{Upper}(i))
-    else
-        error("Bounds error") # TODO: Make this more standard
-    end
-end
-
-
-# Wrapper functions for using regular indecies without specifying regions
-@inline function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int)
-    return apply(op, h, v, Index{Unknown}(i))
-end
-
-@enum Parity begin
-    odd = -1
-    even = 1
-end
-
-struct D2{T,N,M,K} <: ConstantStencilOperator
-    quadratureClosure::NTuple{M,T}
-    innerStencil::Stencil{T,N}
-    closureStencils::NTuple{M,Stencil{T,K}}
-    eClosure::Stencil{T,M}
-    dClosure::Stencil{T,M}
-    parity::Parity
-end
-
-function closureSize(D::D2)::Int
-    return length(D.quadratureClosure)
-end
-
-function readOperator(D2fn, Hfn)
-    d = readSectionedFile(D2fn)
-    h = readSectionedFile(Hfn)
-
-    # Create inner stencil
-    innerStencilWeights = stringToTuple(Float64, d["inner_stencil"][1])
-    width = length(innerStencilWeights)
-    r = (-div(width,2), div(width,2))
-
-    innerStencil = Stencil(r, innerStencilWeights)
-
-    # Create boundary stencils
-    boundarySize = length(d["boundary_stencils"])
-    closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
-
-    for i ∈ 1:boundarySize
-        stencilWeights = stringToTuple(Float64, d["boundary_stencils"][i])
-        width = length(stencilWeights)
-        r = (1-i,width-i)
-        closureStencils = (closureStencils..., Stencil(r, stencilWeights))
-    end
-
-    quadratureClosure = pad_tuple(stringToTuple(Float64, h["closure"][1]), boundarySize)
-    eClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["e"][1]), boundarySize))
-    dClosure = Stencil((0,boundarySize-1), pad_tuple(stringToTuple(Float64, d["d1"][1]), boundarySize))
-
-    d2 = D2(
-        quadratureClosure,
-        innerStencil,
-        closureStencils,
-        eClosure,
-        dClosure,
-        even
-    )
-
-    return d2
-end
-
-
-function apply_e(op::D2, v::AbstractVector, ::Type{Lower})
-    apply(op.eClosure,v,1)
-end
-
-function apply_e(op::D2, v::AbstractVector, ::Type{Upper})
-    apply(flip(op.eClosure),v,length(v))
-end
-
-
-function apply_d(op::D2, h_inv::Real, v::AbstractVector, ::Type{Lower})
-    -h_inv*apply(op.dClosure,v,1)
-end
-
-function apply_d(op::D2, h_inv::Real, v::AbstractVector, ::Type{Upper})
-    -h_inv*apply(flip(op.dClosure),v,length(v))
-end
-
-function readSectionedFile(filename)::Dict{String, Vector{String}}
-    f = open(filename)
-    sections = Dict{String, Vector{String}}()
-    currentKey = ""
-
-    for ln ∈ eachline(f)
-        if ln == "" || ln[1] == '#' # Skip comments and empty lines
-            continue
-        end
-
-        if isletter(ln[1]) # Found start of new section
-            if ~haskey(sections, ln)
-                sections[ln] =  Vector{String}()
-            end
-            currentKey = ln
-            continue
-        end
-
-        push!(sections[currentKey], ln)
-    end
-
-    return sections
-end
-
-function stringToTuple(T::DataType, s::String)
-    return Tuple(stringToVector(T,s))
-end
-
-function stringToVector(T::DataType, s::String)
-    return T.(eval.(Meta.parse.(split(s))))
-end
-
-
-function pad_tuple(t::NTuple{N, T}, n::Integer) where {N,T}
-    if N >= n
-        return t
-    else
-        return pad_tuple((t..., zero(T)), n)
-    end
-end
--- a/stencil.jl	Wed Jun 26 12:57:20 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,38 +0,0 @@
-struct Stencil{T<:Real,N}
-    range::Tuple{Int,Int}
-    weights::NTuple{N,T}
-
-    function Stencil(range::Tuple{Int,Int},weights::NTuple{N,T}) where {T <: Real, N}
-        @assert range[2]-range[1]+1 == N
-        new{T,N}(range,weights)
-    end
-end
-
-function flip(s::Stencil)
-    range = (-s.range[2], -s.range[1])
-    return Stencil(range, reverse(s.weights))
-end
-
-# Provides index into the Stencil based on offset for the root element
-@inline function Base.getindex(s::Stencil, i::Int)
-    @boundscheck if i < s.range[1] || s.range[2] < i
-        return eltype(s.weights)(0)
-    end
-    return s.weights[1 + i - s.range[1]]
-end
-
-Base.@propagate_inbounds @inline function apply(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
-    w = s.weights[1]*v[i + s.range[1]]
-    @simd for k ∈ 2:N
-        w += s.weights[k]*v[i + s.range[1] + k-1]
-    end
-    return w
-end
-
-Base.@propagate_inbounds @inline function apply_backwards(s::Stencil{T,N}, v::AbstractVector, i::Int) where {T,N}
-    w = s.weights[N]*v[i - s.range[2]]
-    @simd for k ∈ N-1:-1:1
-        w += s.weights[k]*v[i - s.range[1] - k + 1]
-    end
-    return w
-end