Mercurial > repos > public > sbplib_julia
changeset 225:3ab0c61f1367 boundary_conditions
Merge in package_refactor
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 26 Jun 2019 14:02:28 +0200 |
| parents | 2aa33d0eef90 (current diff) 31bd875eedf5 (diff) |
| children | eb8525066f9b 5acef2d5db2e |
| files | AbstractGrid.jl EquidistantGrid.jl d2_2nd.txt d2_4th.txt diffOp.jl grid.jl h_2nd.txt h_4th.txt index.jl sbpD2.jl stencil.jl |
| diffstat | 36 files changed, 875 insertions(+), 607 deletions(-) [+] |
line wrap: on
line diff
--- a/AbstractGrid.jl Tue Jun 25 17:26:39 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -abstract type AbstractGrid end - -function dimension(grid::AbstractGrid) - error("Not implemented for abstact type AbstractGrid") -end - -function points(grid::AbstractGrid) - error("Not implemented for abstact type AbstractGrid") -end - -# Evaluate function f on the grid g -function evalOn(g::AbstractGrid, f::Function) - F(x) = f(x...) - return F.(points(g)) -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DiffOps/Manifest.toml Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,61 @@ +# This file is machine-generated - editing it directly is not advised + +[[Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[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" + +[[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]] +deps = ["RegionIndices"] +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"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DiffOps/Project.toml Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,16 @@ +name = "DiffOps" +uuid = "39474f48-97ec-11e9-01fc-6ddcbe5918df" +authors = ["Jonatan Werpers <jonatan.werpers@it.uu.se>"] +version = "0.1.0" + +[deps] +Grids = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a" +RegionIndices = "5d527584-97f1-11e9-084c-4540c7ecf219" +SbpOperators = "204357d8-97fd-11e9-05e7-010897a14cd0" +TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"]
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DiffOps/src/DiffOps.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,236 @@ +module DiffOps + +using RegionIndices +using SbpOperators +using Grids + +export Laplace + +abstract type DiffOp end + +# TBD: The "error("not implemented")" thing seems to be hiding good error information. How to fix that? Different way of saying that these should be implemented? +function apply(D::DiffOp, v::AbstractVector, i::Int) + error("not implemented") +end + +function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real + error("not implemented") +end + +function matrixRepresentation(D::DiffOp) + error("not implemented") +end + +abstract type DiffOpCartesian{Dim} <: DiffOp end + +# DiffOp must have a grid of dimension Dim!!! +function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim} + for I ∈ eachindex(D.grid) + u[I] = apply(D, v, I) + end + + return nothing +end +export apply! + +function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T + apply_region!(D, u, v, Lower, Lower) + apply_region!(D, u, v, Lower, Interior) + apply_region!(D, u, v, Lower, Upper) + apply_region!(D, u, v, Interior, Lower) + apply_region!(D, u, v, Interior, Interior) + apply_region!(D, u, v, Interior, Upper) + apply_region!(D, u, v, Upper, Lower) + apply_region!(D, u, v, Upper, Interior) + apply_region!(D, u, v, Upper, Upper) + return nothing +end + +# Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable +function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T + for I ∈ regionindices(D.grid.size, closureSize(D.op), (r1,r2)) + @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2])) + @inbounds u[I] = apply(D, v, indextuple) + end + return nothing +end +export apply_region! + +function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T + apply_region_tiled!(D, u, v, Lower, Lower) + apply_region_tiled!(D, u, v, Lower, Interior) + apply_region_tiled!(D, u, v, Lower, Upper) + apply_region_tiled!(D, u, v, Interior, Lower) + apply_region_tiled!(D, u, v, Interior, Interior) + apply_region_tiled!(D, u, v, Interior, Upper) + apply_region_tiled!(D, u, v, Upper, Lower) + apply_region_tiled!(D, u, v, Upper, Interior) + apply_region_tiled!(D, u, v, Upper, Upper) + return nothing +end + +using TiledIteration +function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T + ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2)) + # TODO: Pass Tilesize to function + for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2)) + for j ∈ tileaxs[2], i ∈ tileaxs[1] + I = ri[i,j] + u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2]))) + end + end + return nothing +end +export apply_region_tiled! + +function apply(D::DiffOp, v::AbstractVector)::AbstractVector + u = zeros(eltype(v), size(v)) + apply!(D,v,u) + return u +end + +export apply + +struct NormalDerivative{N,M,K} + op::D2{Float64,N,M,K} + grid::EquidistantGrid + bId::CartesianBoundary +end + +function apply_transpose(d::NormalDerivative, v::AbstractArray, I::Integer) + u = selectdim(v,3-dim(d.bId),I) + return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId)) +end + +# Not correct abstraction level +# TODO: Not type stable D:< +function apply(d::NormalDerivative, v::AbstractArray, I::Tuple{Integer,Integer}) + i = I[dim(d.bId)] + j = I[3-dim(d.bId)] + N_i = d.grid.size[dim(d.bId)] + + r = getregion(i, closureSize(d.op), N_i) + + if r != region(d.bId) + return 0 + end + + if r == Lower + # Note, closures are indexed by offset. Fix this D:< + return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j] + elseif r == Upper + return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-j]*v[j] + end +end + +struct BoundaryValue{N,M,K} + op::D2{Float64,N,M,K} + grid::EquidistantGrid + bId::CartesianBoundary +end + +function apply(e::BoundaryValue, v::AbstractArray, I::Tuple{Integer,Integer}) + i = I[dim(e.bId)] + j = I[3-dim(e.bId)] + N_i = e.grid.size[dim(e.bId)] + + r = getregion(i, closureSize(e.op), N_i) + + if r != region(e.bId) + return 0 + end + + if r == Lower + # Note, closures are indexed by offset. Fix this D:< + return e.op.eClosure[i-1]*v[j] + elseif r == Upper + return e.op.eClosure[N_i-j]*v[j] + end +end + +function apply_transpose(e::BoundaryValue, v::AbstractArray, I::Integer) + u = selectdim(v,3-dim(e.bId),I) + return apply_e(e.op, u, region(e.bId)) +end + +struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim} + grid::EquidistantGrid{Dim,T} + a::T + op::D2{Float64,N,M,K} + # e::BoundaryValue + # d::NormalDerivative +end + +function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim + error("not implemented") +end + +# u = L*v +function apply(L::Laplace{1}, v::AbstractVector, i::Int) + uᵢ = L.a * SbpOperators.apply(L.op, L.grid.spacing[1], v, i) + return uᵢ +end + +@inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2} + # 2nd x-derivative + @inbounds vx = view(v, :, Int(I[2])) + @inbounds uᵢ = L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[1], vx , I[1]) + # 2nd y-derivative + @inbounds vy = view(v, Int(I[1]), :) + @inbounds uᵢ += L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[2], vy, I[2]) + # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors + return uᵢ +end + +# Slow but maybe convenient? +function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2}) + I = Index{Unknown}.(Tuple(i)) + apply(L, v, I) +end + +struct BoundaryOperator + +end + + +""" +A BoundaryCondition should implement the method + sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...) +""" +abstract type BoundaryCondition end + +struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end + +function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} + e = BoundaryValue(L.op, L.grid, Bid()) + d = NormalDerivative(L.op, L.grid, Bid()) + Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid()) + # TODO: Implement BoundaryQuadrature method + + return -L.Hi*e*Hᵧ*(d'*v - g) + # Need to handle d'*v - g so that it is an AbstractArray that TensorMappings can act on +end + +struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition + tau::Float64 +end + +function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid} + e = BoundaryValue(L.op, L.grid, Bid()) + d = NormalDerivative(L.op, L.grid, Bid()) + Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid()) + # TODO: Implement BoundaryQuadrature method + + return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g) + # Need to handle scalar multiplication and addition of TensorMapping +end + +# function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D +# return apply(s.L, v, i) + +# sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + +# sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + +# sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) + +# sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) +# end + +end # module \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DiffOps/test/runtests.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,4 @@ +using DiffOps +using Test + +@test_broken false
--- a/EquidistantGrid.jl Tue Jun 25 17:26:39 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,56 +0,0 @@ -# EquidistantGrid is a grid with equidistant grid spacing per coordinat -# direction. The domain is defined through the two points P1 = x̄₁, P2 = x̄₂ -# by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto -# the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2) -# the domain is defined as (-1,1)x(0,2). - -struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid - size::NTuple{Dim, Int} # First coordinate direction stored first - limit_lower::NTuple{Dim, T} - limit_upper::NTuple{Dim, T} - inverse_spacing::NTuple{Dim, T} # The reciprocal of the grid spacing - - # General constructor - function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T - @assert all(size.>0) - @assert all(limit_upper.-limit_lower .!= 0) - inverse_spacing = (size.-1)./abs.(limit_upper.-limit_lower) - return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing) - end -end - -function Base.eachindex(grid::EquidistantGrid) - CartesianIndices(grid.size) -end - -# Returns the number of dimensions of an EquidistantGrid. -# -# @Input: grid - an EquidistantGrid -# @Return: dimension - The dimension of the grid -function dimension(grid::EquidistantGrid) - return length(grid.size) -end - -# Returns the spacing of the grid -# -function spacing(grid::EquidistantGrid) - return 1.0./grid.inverse_spacing -end - -# Computes the points of an EquidistantGrid as an array of tuples with -# the same dimension as the grid. -# -# @Input: grid - an EquidistantGrid -# @Return: points - the points of the grid. -function points(grid::EquidistantGrid) - # TODO: Make this return an abstract array? - indices = Tuple.(CartesianIndices(grid.size)) - h = spacing(grid) - return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices) -end - -function pointsalongdim(grid::EquidistantGrid, dim::Integer) - @assert dim<=dimension(grid) - @assert dim>0 - points = collect(range(grid.limit_lower[dim],stop=grid.limit_upper[dim],length=grid.size[dim])) -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Grids/Manifest.toml Wed Jun 26 14:02:28 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/Grids/Project.toml Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,13 @@ +name = "Grids" +uuid = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a" +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/Grids/src/AbstractGrid.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,24 @@ +""" + AbstractGrid + +Should implement + dimension(grid::AbstractGrid) + points(grid::AbstractGrid) + +""" +abstract type AbstractGrid end + +function dimension end +function points end +export dimension, points + +""" + evalOn(g::AbstractGrid, f::Function) + +Evaluate function f on the grid g +""" +function evalOn(g::AbstractGrid, f::Function) + F(x) = f(x...) + return F.(points(g)) +end +export evalOn
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Grids/src/EquidistantGrid.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,58 @@ +# EquidistantGrid is a grid with equidistant grid spacing per coordinat +# direction. The domain is defined through the two points P1 = x̄₁, P2 = x̄₂ +# by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto +# the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2) +# the domain is defined as (-1,1)x(0,2). + +export EquidistantGrid + +struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid + size::NTuple{Dim, Int} # First coordinate direction stored first + limit_lower::NTuple{Dim, T} + limit_upper::NTuple{Dim, T} + inverse_spacing::NTuple{Dim, T} # The reciprocal of the grid spacing + + # General constructor + function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T + @assert all(size.>0) + @assert all(limit_upper.-limit_lower .!= 0) + inverse_spacing = (size.-1)./abs.(limit_upper.-limit_lower) + return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing) + end +end + +function Base.eachindex(grid::EquidistantGrid) + CartesianIndices(grid.size) +end + +# Returns the number of dimensions of an EquidistantGrid. +# +# @Input: grid - an EquidistantGrid +# @Return: dimension - The dimension of the grid +function dimension(grid::EquidistantGrid) + return length(grid.size) +end + +# Returns the spacing of the grid +# +function spacing(grid::EquidistantGrid) + return 1.0./grid.inverse_spacing +end + +# Computes the points of an EquidistantGrid as an array of tuples with +# the same dimension as the grid. +# +# @Input: grid - an EquidistantGrid +# @Return: points - the points of the grid. +function points(grid::EquidistantGrid) + # TODO: Make this return an abstract array? + indices = Tuple.(CartesianIndices(grid.size)) + h = spacing(grid) + return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices) +end + +function pointsalongdim(grid::EquidistantGrid, dim::Integer) + @assert dim<=dimension(grid) + @assert dim>0 + points = collect(range(grid.limit_lower[dim],stop=grid.limit_upper[dim],length=grid.size[dim])) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Grids/src/Grids.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,15 @@ +module Grids + +using RegionIndices + +export BoundaryIdentifier, CartesianBoundary + +abstract type BoundaryIdentifier end +struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end +dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim +region(::CartesianBoundary{Dim, R}) where {Dim, R} = R + +include("AbstractGrid.jl") +include("EquidistantGrid.jl") + +end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Grids/test/runtests.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,4 @@ +using Grids +using Test + +@test_broken false
--- a/Manifest.toml Tue Jun 25 17:26:39 2019 +0200 +++ b/Manifest.toml Wed Jun 26 14:02:28 2019 +0200 @@ -1,6 +1,72 @@ # This file is machine-generated - editing it directly is not advised +[[Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[DiffOps]] +deps = ["Grids", "RegionIndices", "SbpOperators", "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]] +deps = ["RegionIndices"] +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 Tue Jun 25 17:26:39 2019 +0200 +++ b/Project.toml Wed Jun 26 14:02:28 2019 +0200 @@ -1,2 +1,6 @@ [deps] +DiffOps = "39474f48-97ec-11e9-01fc-6ddcbe5918df" +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/RegionIndices/Project.toml Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,10 @@ +name = "RegionIndices" +uuid = "5d527584-97f1-11e9-084c-4540c7ecf219" +authors = ["Jonatan Werpers <jonatan.werpers@it.uu.se>"] +version = "0.1.0" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"]
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/RegionIndices/src/RegionIndices.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,81 @@ +module RegionIndices + +abstract type Region end +struct Interior <: Region end +struct Lower <: Region end +struct Upper <: Region end +struct Unknown <: Region end + +export Region, Interior, Lower, Upper, Unknown + +struct Index{R<:Region, T<:Integer} + i::T + + Index{R,T}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i) + Index{R}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i) + Index(i::T, ::Type{R}) where {R<:Region,T<:Integer} = Index{R,T}(i) + Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed? +end + +export Index + +# Index(R::Type{<:Region}) = Index{R} + +## Vill kunna skriva +## IndexTupleType(Int, (Lower, Interior)) +Index(R::Type{<:Region}, T::Type{<:Integer}) = Index{R,T} +IndexTupleType(T::Type{<:Integer},R::NTuple{N, DataType} where N) = Tuple{Index.(R, T)...} + +Base.convert(::Type{T}, i::Index{R,T} where R) where T = i.i +Base.convert(::Type{CartesianIndex}, I::NTuple{N,Index} where N) = CartesianIndex(convert.(Int, I)) + +Base.Int(I::Index) = I.i + +function Index(i::Integer, boundary_width::Integer, dim_size::Integer) + return Index{getregion(i,boundary_width,dim_size)}(i) +end + +IndexTuple(t::Vararg{Tuple{T, DataType}}) where T<:Integer = Index.(t) +export IndexTuple + +# TODO: Use the values of the region structs, e.g. Lower(), for the region parameter instead of the types. +# For example the following works: +# (Lower(),Upper()) isa NTuple{2, Region} -> true +# typeof((Lower(),Upper())) -> Tuple{Lower,Upper} +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(getrange,gridsize,closuresize,region) + return CartesianIndices(regions) +end + +export regionindices + +function getregion(i::Integer, boundary_width::Integer, dim_size::Integer) + if 0 < i <= boundary_width + return Lower + elseif boundary_width < i <= dim_size-boundary_width + return Interior + elseif dim_size-boundary_width < i <= dim_size + return Upper + else + error("Bounds error") # TODO: Make this more standard + end +end + +function getrange(gridsize::Integer, closuresize::Integer, region::DataType) + 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 + +end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/RegionIndices/test/runtests.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,4 @@ +using RegionIndices +using Test + +@test_broken false
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SbpOperators/Manifest.toml Wed Jun 26 14:02:28 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 14:02:28 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/operators/d2_2nd.txt Wed Jun 26 14:02:28 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/operators/d2_4th.txt Wed Jun 26 14:02:28 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/operators/h_2nd.txt Wed Jun 26 14:02:28 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/operators/h_4th.txt Wed Jun 26 14:02:28 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/SbpOperators.jl Wed Jun 26 14:02:28 2019 +0200 @@ -0,0 +1,163 @@ +module SbpOperators + +using RegionIndices + +include("stencil.jl") + +export D2, closureSize, apply, readOperator, apply_e, apply_d + +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/stencil.jl Wed Jun 26 14:02:28 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 14:02:28 2019 +0200 @@ -0,0 +1,4 @@ +using SbpOperators +using Test + +@test_broken false
--- a/d2_2nd.txt Tue Jun 25 17:26:39 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 Tue Jun 25 17:26:39 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/diffOp.jl Tue Jun 25 17:26:39 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,220 +0,0 @@ -abstract type DiffOp end - -# TBD: The "error("not implemented")" thing seems to be hiding good error information. How to fix that? Different way of saying that these should be implemented? -function apply(D::DiffOp, v::AbstractVector, i::Int) - error("not implemented") -end - -function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real - error("not implemented") -end - -function matrixRepresentation(D::DiffOp) - error("not implemented") -end - -abstract type DiffOpCartesian{Dim} <: DiffOp end - -# DiffOp must have a grid of dimension Dim!!! -function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim} - for I ∈ eachindex(D.grid) - u[I] = apply(D, v, I) - end - - return nothing -end - -function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T - apply_region!(D, u, v, Lower, Lower) - apply_region!(D, u, v, Lower, Interior) - apply_region!(D, u, v, Lower, Upper) - apply_region!(D, u, v, Interior, Lower) - apply_region!(D, u, v, Interior, Interior) - apply_region!(D, u, v, Interior, Upper) - apply_region!(D, u, v, Upper, Lower) - apply_region!(D, u, v, Upper, Interior) - apply_region!(D, u, v, Upper, Upper) - return nothing -end - -# Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable -function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T - for I ∈ regionindices(D.grid.size, closureSize(D.op), (r1,r2)) - @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2])) - @inbounds u[I] = apply(D, v, indextuple) - end - return nothing -end - -function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T - apply_region_tiled!(D, u, v, Lower, Lower) - apply_region_tiled!(D, u, v, Lower, Interior) - apply_region_tiled!(D, u, v, Lower, Upper) - apply_region_tiled!(D, u, v, Interior, Lower) - apply_region_tiled!(D, u, v, Interior, Interior) - apply_region_tiled!(D, u, v, Interior, Upper) - apply_region_tiled!(D, u, v, Upper, Lower) - apply_region_tiled!(D, u, v, Upper, Interior) - apply_region_tiled!(D, u, v, Upper, Upper) - return nothing -end - -using TiledIteration -function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T - ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2)) - # TODO: Pass Tilesize to function - for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2)) - for j ∈ tileaxs[2], i ∈ tileaxs[1] - I = ri[i,j] - u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2]))) - end - end - return nothing -end - -function apply(D::DiffOp, v::AbstractVector)::AbstractVector - u = zeros(eltype(v), size(v)) - apply!(D,v,u) - return u -end - -struct NormalDerivative{N,M,K} - op::D2{Float64,N,M,K} - grid::EquidistantGrid - bId::CartesianBoundary -end - -function apply_transpose(d::NormalDerivative, v::AbstractArray, I::Integer) - u = selectdim(v,3-dim(d.bId),I) - return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId)) -end - -# Not correct abstraction level -# TODO: Not type stable D:< -function apply(d::NormalDerivative, v::AbstractArray, I::Tuple{Integer,Integer}) - i = I[dim(d.bId)] - j = I[3-dim(d.bId)] - N_i = d.grid.size[dim(d.bId)] - - r = getregion(i, closureSize(d.op), N_i) - - if r != region(d.bId) - return 0 - end - - if r == Lower - # Note, closures are indexed by offset. Fix this D:< - return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j] - elseif r == Upper - return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-j]*v[j] - end -end - -struct BoundaryValue{N,M,K} - op::D2{Float64,N,M,K} - grid::EquidistantGrid - bId::CartesianBoundary -end - -function apply(e::BoundaryValue, v::AbstractArray, I::Tuple{Integer,Integer}) - i = I[dim(e.bId)] - j = I[3-dim(e.bId)] - N_i = e.grid.size[dim(e.bId)] - - r = getregion(i, closureSize(e.op), N_i) - - if r != region(e.bId) - return 0 - end - - if r == Lower - # Note, closures are indexed by offset. Fix this D:< - return e.op.eClosure[i-1]*v[j] - elseif r == Upper - return e.op.eClosure[N_i-j]*v[j] - end -end - -function apply_transpose(e::BoundaryValue, v::AbstractArray, I::Integer) - u = selectdim(v,3-dim(e.bId),I) - return apply_e(e.op, u, region(e.bId)) -end - -struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim} - grid::EquidistantGrid{Dim,T} - a::T - op::D2{Float64,N,M,K} - e::BoundaryValue - d::NormalDerivative -end - -function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim - error("not implemented") -end - -# u = L*v -function apply(L::Laplace{1}, v::AbstractVector, i::Int) - uᵢ = L.a * apply(L.op, L.grid.spacing[1], v, i) - return uᵢ -end - -@inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2} - # 2nd x-derivative - @inbounds vx = view(v, :, Int(I[2])) - @inbounds uᵢ = L.a*apply(L.op, L.grid.inverse_spacing[1], vx , I[1]) - # 2nd y-derivative - @inbounds vy = view(v, Int(I[1]), :) - @inbounds uᵢ += L.a*apply(L.op, L.grid.inverse_spacing[2], vy, I[2]) - return uᵢ -end - -# Slow but maybe convenient? -function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2}) - I = Index{Unknown}.(Tuple(i)) - apply(L, v, I) -end - -struct BoundaryOperator - -end - - -""" -A BoundaryCondition should implement the method - sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...) -""" -abstract type BoundaryCondition end - -struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end - -function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid} - e = BoundaryValue(L.op, L.grid, Bid()) - d = NormalDerivative(L.op, L.grid, Bid()) - Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid()) - # TODO: Implement BoundaryQuadrature method - - return -L.Hi*e*Hᵧ*(d'*v - g) - # Need to handle d'*v - g so that it is an AbstractArray that TensorMappings can act on -end - -struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition - tau::Float64 -end - -function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid} - e = BoundaryValue(L.op, L.grid, Bid()) - d = NormalDerivative(L.op, L.grid, Bid()) - Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid()) - # TODO: Implement BoundaryQuadrature method - - return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g) - # Need to handle scalar multiplication and addition of TensorMapping -end - -# function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D -# return apply(s.L, v, i) + -# sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + -# sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + -# sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) + -# sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) -# end
--- a/grid.jl Tue Jun 25 17:26:39 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -abstract type BoundaryIdentifier end -struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end -dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim -region(::CartesianBoundary{Dim, R}) where {Dim, R} = R - -include("AbstractGrid.jl") -include("EquidistantGrid.jl")
--- a/h_2nd.txt Tue Jun 25 17:26:39 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 Tue Jun 25 17:26:39 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/index.jl Tue Jun 25 17:26:39 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ -abstract type Region end -struct Interior <: Region end -struct Lower <: Region end -struct Upper <: Region end -struct Unknown <: Region end - -struct Index{R<:Region, T<:Integer} - i::T - - Index{R,T}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i) - Index{R}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i) - Index(i::T, ::Type{R}) where {R<:Region,T<:Integer} = Index{R,T}(i) - Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed? -end - -# Index(R::Type{<:Region}) = Index{R} - -## Vill kunna skriva -## IndexTupleType(Int, (Lower, Interior)) -Index(R::Type{<:Region}, T::Type{<:Integer}) = Index{R,T} -IndexTupleType(T::Type{<:Integer},R::NTuple{N, DataType} where N) = Tuple{Index.(R, T)...} - -Base.convert(::Type{T}, i::Index{R,T} where R) where T = i.i -Base.convert(::Type{CartesianIndex}, I::NTuple{N,Index} where N) = CartesianIndex(convert.(Int, I)) - -Base.Int(I::Index) = I.i - -function Index(i::Integer, boundary_width::Integer, dim_size::Integer) - return Index{getregion(i,boundary_width,dim_size)}(i) -end - -IndexTuple(t::Vararg{Tuple{T, DataType}}) where T<:Integer = Index.(t) - -# TODO: Use the values of the region structs, e.g. Lower(), for the region parameter instead of the types. -# For example the following works: -# (Lower(),Upper()) isa NTuple{2, Region} -> true -# typeof((Lower(),Upper())) -> Tuple{Lower,Upper} -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(getrange,gridsize,closuresize,region) - return CartesianIndices(regions) -end - -function getregion(i::Integer, boundary_width::Integer, dim_size::Integer) - if 0 < i <= boundary_width - return Lower - elseif boundary_width < i <= dim_size-boundary_width - return Interior - elseif dim_size-boundary_width < i <= dim_size - return Upper - else - error("Bounds error") # TODO: Make this more standard - end -end - -function getrange(gridsize::Integer, closuresize::Integer, region::DataType) - 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
--- a/sbp.jl Tue Jun 25 17:26:39 2019 +0200 +++ b/sbp.jl Wed Jun 26 14:02:28 2019 +0200 @@ -1,8 +1,9 @@ module sbp -include("index.jl") -include("grid.jl") -include("stencil.jl") -include("sbpD2.jl") -include("diffOp.jl") + +using Grids +using RegionIndices +using SbpOperators +using DiffOps + include("TimeStepper.jl") end # module
--- a/sbpD2.jl Tue Jun 25 17:26:39 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 Tue Jun 25 17:26:39 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