Mercurial > repos > public > sbplib_julia

annotate diffOp.jl @ 106:8010aadc4ff1 cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up in Laplace 2d apply!
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 08 Feb 2019 10:59:01 +0100
parents 44cd6b4371de
children e3e32299fe72
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
1 abstract type DiffOp end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
2
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
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?
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
4 function apply(D::DiffOp, v::AbstractVector, i::Int)
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
5 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
6 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
7
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
8 function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
9 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
10 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
11
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
12 function matrixRepresentation(D::DiffOp)
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
13 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
14 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
15
55
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
16 function boundaryCondition(D::DiffOp,b::Grid.BoundaryId,type)::(Closure, Penalty)
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
17 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
18 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
19
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
20 function interface(Du::DiffOp, Dv::DiffOp, b::Grid.BoundaryId; type)
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
21 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
22 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
23
55
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
24 abstract type Closure end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
25
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
26 function apply(c::Closure, v::AbstractVector, i::Int)
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
27 error("not implemented")
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
28 end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
29
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
30 abstract type Penalty end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
31
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
32 function apply(c::Penalty, g, i::Int)
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
33 error("not implemented")
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
34 end
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
35
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
36 abstract type DiffOpCartesian{Dim} <: DiffOp end
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
37
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
38 # DiffOp must have a grid of dimension Dim!!!
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
39 function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim}
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
40 for I ∈ eachindex(D.grid)
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
41 u[I] = apply(D, v, I)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
42 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
43
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
44 return nothing
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
45 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
46
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
47 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
48 u = zeros(eltype(v), size(v))
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
49 apply!(D,v,u)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
50 return u
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
51 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
52
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
53 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
54 grid::Grid.EquidistantGrid{Dim,T}
76
81d9510cb2d0 Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents: 73
diff changeset
55 a::T
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
56 op::D2{Float64,N,M,K}
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
57 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
58
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
59 function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
60 error("not implemented")
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
61 end
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
62
6
cb8e50ca9e15 Add attempt att apply methods for Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 2
diff changeset
63 # u = L*v
65
7054230b639c Make dimension a type parameter in Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
64 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
65 h = Grid.spacings(L.grid)[1]
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
66 uᵢ = L.a * apply(L.op, h, v, i)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
67 return uᵢ
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
68 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
69
103
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
70 function apply!(L::Laplace{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
71 regions = (Lower, Interior, Upper)
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
72 for r1 ∈ regions
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
73 for r2 ∈ regions
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
74 apply!(L, u, v, r1, r2)
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
75 end
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
76 end
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
77 return nothing
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
78 end
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
79
105
44cd6b4371de Specify region types in Laplace 2d apply!
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 103
diff changeset
80 function apply!(L::Laplace{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
44cd6b4371de Specify region types in Laplace 2d apply!
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 103
diff changeset
81 N = L.grid.numberOfPointsPerDim
44cd6b4371de Specify region types in Laplace 2d apply!
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 103
diff changeset
82 closuresize = closureSize(L.op)
106
8010aadc4ff1 Clean up in Laplace 2d apply!
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 105
diff changeset
83 for I ∈ regionindices(N, closuresize, (r1,r2))
103
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
84 @inbounds indextuple = (Index(I[1], r1), Index(I[2], r2))
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
85 @inbounds u[I] = apply(L, v, indextuple)
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
86 end
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
87 return nothing
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
88 end
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
89
82
45dece5e4928 Use unsafe views
Jonatan Werpers <jonatan@werpers.com>
parents: 70
diff changeset
90 using UnsafeArrays
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
91 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
92 h = Grid.spacings(L.grid)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
93 # 2nd x-derivative
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
94 @inbounds vx = uview(v, :, Int(I[2]))
103
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
95 @inbounds uᵢ = L.a*apply(L.op, h[1], vx , I[1])
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
96 # 2nd y-derivative
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
97 @inbounds vy = uview(v, Int(I[1]), :)
103
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
98 @inbounds uᵢ += L.a*apply(L.op, h[2], vy, I[2])
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
99 return uᵢ
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
100 end
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
101
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
102 # Slow but maybe convenient?
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
103 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
104 I = Index{Unknown}.(Tuple(i))
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
105 apply(L, v, I)
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
106 end