Mercurial > repos > public > sbplib_julia

annotate diffOp.jl @ 93:93df72e2b135 stencil_index

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement apply for 2D-Laplace which takes an StencilIndex as input
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 04 Feb 2019 09:13:48 +0100
parents 8d505e9bc715
children 84b1ad5a3755
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
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
3 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
4 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
5 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
6
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
7 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
8 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
9 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
10
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
11 function matrixRepresentation(D::DiffOp)
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
12 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
13 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
14
55
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
15 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
16 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
17 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
18
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
19 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
20 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
21 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
22
55
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
23 abstract type Closure end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
24
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
25 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
26 error("not implemented")
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
27 end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
28
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
29 abstract type Penalty end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
30
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
31 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
32 error("not implemented")
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
33 end
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
34
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
35 # DiffOp must have a grid!!!
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
36 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
37 for i ∈ 1:Grid.numberOfPoints(D.grid)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
38 u[i] = apply(D, v, i)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
39 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
40 return nothing
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
41 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
42
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
43 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
44 u = zeros(eltype(v), size(v))
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
45 apply!(D,v,u)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
46 return u
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
47 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
48
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
49 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOp
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
50 grid::Grid.EquidistantGrid{Dim,T}
76
81d9510cb2d0 Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents: 73
diff changeset
51 a::T
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
52 op::D2{Float64,N,M,K}
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
53 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
54
6
cb8e50ca9e15 Add attempt att apply methods for Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 2
diff changeset
55 # u = L*v
65
7054230b639c Make dimension a type parameter in Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
56 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
57 h = Grid.spacings(L.grid)[1]
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
58 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
59 return uᵢ
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
60 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
61
82
45dece5e4928 Use unsafe views
Jonatan Werpers <jonatan@werpers.com>
parents: 70
diff changeset
62 using UnsafeArrays
45dece5e4928 Use unsafe views
Jonatan Werpers <jonatan@werpers.com>
parents: 70
diff changeset
63
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
64 # u = L*v
93
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
65 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
66 h = Grid.spacings(L.grid)[1]
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
67 uᵢ = L.a * apply(L.op, h, v, i)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
68 return uᵢ
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
69 end
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
70
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
71 function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
72 lowerind, upperind, interiorind = stencilindices(L.grid)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
73 for si ∈ lowerind
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
74 u[si.globalindex] = apply(L, v, si)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
75 end
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
76 for si ∈ upperind
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
77 u[si.globalindex] = apply(L, v, si)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
78 end
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
79 for si ∈ interiorind
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
80 u[si.globalindex] = apply(L, v, si)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
81 end
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
82 return nothing
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
83 end
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
84
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
85 @inline function apply(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
86 h = Grid.spacings(L.grid)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
87 li = LinearIndices(L.grid.numberOfPointsPerDim)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
88 # 2nd x-derivative
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
89 @inbounds vx = uview(v, uview(li,:,si.gridindex[2]))
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
90 uᵢ = apply(L.op, h[1], vx , si.gridindex[1], si)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
91 # 2nd y-derivative
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
92 @inbounds vy = uview(v, uview(li,si.gridindex[1],:))
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
93 uᵢ += apply(L.op, h[2], vy, si.gridindex[2], si)
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
94 return uᵢ
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
95 end
93df72e2b135 Implement apply for 2D-Laplace which takes an StencilIndex as input
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 85
diff changeset
96
65
7054230b639c Make dimension a type parameter in Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
97 function apply(L::Laplace{2}, v::AbstractVector, i::Int)
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
98 h = Grid.spacings(L.grid)
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
99
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
100 li = LinearIndices(L.grid.numberOfPointsPerDim)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
101 ci = CartesianIndices(L.grid.numberOfPointsPerDim)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
102 I = ci[i]
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
103
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
104 # 2nd x-derivative
83
Jonatan Werpers <jonatan@werpers.com>
parents: 82 78
diff changeset
105 @inbounds vx = uview(v, uview(li,:,I[2]))
70
e4fa13137d12 Use double view
Jonatan Werpers <jonatan@werpers.com>
parents: 66
diff changeset
106 uᵢ = 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
107 # 2nd y-derivative
83
Jonatan Werpers <jonatan@werpers.com>
parents: 82 78
diff changeset
108 @inbounds vy = uview(v, uview(li,I[1],:))
70
e4fa13137d12 Use double view
Jonatan Werpers <jonatan@werpers.com>
parents: 66
diff changeset
109 uᵢ += apply(L.op, h[2], vy, I[2])
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
110
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
111 return uᵢ
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
112 end