Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 134:79699dda29be

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in cell_based_test
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 21 Feb 2019 16:27:28 +0100
parents 1aaeb46ba5f4
children bb1cc9c7877c c6aaf061c0a9
comparison
equal deleted inserted replaced
84:48079bd39969 134:79699dda29be
1 abstract type DiffOp end 1 abstract type DiffOp end
2 2
3 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector) 3 # 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?
4 function apply(D::DiffOp, v::AbstractVector, i::Int)
4 error("not implemented") 5 error("not implemented")
5 end 6 end
6 7
7 function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real 8 function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real
8 error("not implemented") 9 error("not implemented")
30 31
31 function apply(c::Penalty, g, i::Int) 32 function apply(c::Penalty, g, i::Int)
32 error("not implemented") 33 error("not implemented")
33 end 34 end
34 35
35 # Differential operator for a*d^2/dx^2 36 abstract type DiffOpCartesian{Dim} <: DiffOp end
36 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOp 37
38 # DiffOp must have a grid of dimension Dim!!!
39 function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim}
40 for I ∈ eachindex(D.grid)
41 u[I] = apply(D, v, I)
42 end
43
44 return nothing
45 end
46
47 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
48 apply_region!(D, u, v, Lower, Lower)
49 apply_region!(D, u, v, Lower, Interior)
50 apply_region!(D, u, v, Lower, Upper)
51 apply_region!(D, u, v, Interior, Lower)
52 apply_region!(D, u, v, Interior, Interior)
53 apply_region!(D, u, v, Interior, Upper)
54 apply_region!(D, u, v, Upper, Lower)
55 apply_region!(D, u, v, Upper, Interior)
56 apply_region!(D, u, v, Upper, Upper)
57 return nothing
58 end
59
60 # Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable
61 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
62 for I ∈ regionindices(D.grid.size, closureSize(D.op), (r1,r2))
63 @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
64 @inbounds u[I] = apply(D, v, indextuple)
65 end
66 return nothing
67 end
68
69 function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
70 apply_region_tiled!(D, u, v, Lower, Lower)
71 apply_region_tiled!(D, u, v, Lower, Interior)
72 apply_region_tiled!(D, u, v, Lower, Upper)
73 apply_region_tiled!(D, u, v, Interior, Lower)
74 apply_region_tiled!(D, u, v, Interior, Interior)
75 apply_region_tiled!(D, u, v, Interior, Upper)
76 apply_region_tiled!(D, u, v, Upper, Lower)
77 apply_region_tiled!(D, u, v, Upper, Interior)
78 apply_region_tiled!(D, u, v, Upper, Upper)
79 return nothing
80 end
81
82 using TiledIteration
83 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
84 ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2))
85 for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2)) # TBD: Is this the right way, the right size?
86 for j ∈ tileaxs[2], i ∈ tileaxs[1]
87 I = ri[i,j]
88 u[i,j] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
89 end
90 end
91 return nothing
92 end
93
94 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
95 u = zeros(eltype(v), size(v))
96 apply!(D,v,u)
97 return u
98 end
99
100 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
37 grid::Grid.EquidistantGrid{Dim,T} 101 grid::Grid.EquidistantGrid{Dim,T}
38 a::T 102 a::T
39 op::D2{Float64,N,M,K} 103 op::D2{Float64,N,M,K}
40 end 104 end
41 105
42 # u = L*v 106 function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
43 function apply!(L::Laplace{1}, u::AbstractVector, v::AbstractVector) 107 error("not implemented")
44 h = Grid.spacings(L.grid)[1]
45 apply!(L.op, u, v, h)
46 u .= L.a * u
47 return nothing
48 end 108 end
49 109
50 # u = L*v 110 # u = L*v
51 function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector) 111 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
52 u .= 0*u 112 uᵢ = L.a * apply(L.op, L.grid.spacing[1], v, i)
53 h = Grid.spacings(L.grid) 113 return uᵢ
114 end
54 115
55 li = LinearIndices(L.grid.numberOfPointsPerDim) 116 @inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
56 n_x, n_y = L.grid.numberOfPointsPerDim 117 # 2nd x-derivative
118 @inbounds vx = view(v, :, Int(I[2]))
119 @inbounds uᵢ = L.a*apply(L.op, L.grid.inverse_spacing[1], vx , I[1])
120 # 2nd y-derivative
121 @inbounds vy = view(v, Int(I[1]), :)
122 @inbounds uᵢ += L.a*apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
123 return uᵢ
124 end
57 125
58 126 # Slow but maybe convenient?
59 # For each x 127 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})
60 temp = zeros(eltype(u), n_y) 128 I = Index{Unknown}.(Tuple(i))
61 for i ∈ 1:n_x 129 apply(L, v, I)
62
63 v_i = view(v, li[i,:])
64 apply!(L.op, temp, v_i, h[2])
65
66 u[li[i,:]] += temp
67 end
68
69 # For each y
70 temp = zeros(eltype(u), n_x)
71 for i ∈ 1:n_y
72 v_i = view(v, li[:,i])
73 apply!(L.op, temp, v_i, h[1])
74
75 u[li[:,i]] += temp
76 end
77
78 u .= L.a*u
79
80 return nothing
81 end 130 end