Mercurial > repos > public > sbplib_julia

annotate diffOp.jl @ 390:074862bcd626 patch_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Close branch
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 01 Oct 2020 08:28:03 +0200
parents 2882e1318cc3
children
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
33
e8d7137b3f07 Raname apply(DiffOp) to apply!(DiffOp)
Jonatan Werpers <jonatan@werpers.com>
parents: 31
diff changeset
3 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector)
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
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
35 # Differential operator for a*d^2/dx^2
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
36 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
37 grid::Grid.EquidistantGrid{Dim,T}
76
81d9510cb2d0 Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents: 73
diff changeset
38 a::T
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
39 op::D2{Float64,N,M,K}
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
40 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
41
6
cb8e50ca9e15 Add attempt att apply methods for Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 2
diff changeset
42 # u = L*v
76
81d9510cb2d0 Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents: 73
diff changeset
43 function apply!(L::Laplace{1}, u::AbstractVector, v::AbstractVector)
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
44 h = Grid.spacings(L.grid)[1]
43
ef060ab3b035 remove stride and remove some bugs
Ylva Rydin <ylva.rydin@telia.com>
parents: 34
diff changeset
45 apply!(L.op, u, v, h)
48
be53b58581a6 Fix bug with call by reference in Laplace1D apply
Jonatan Werpers <jonatan@werpers.com>
parents: 43
diff changeset
46 u .= L.a * u
9
aafc4bdbe675 Add h in stencil application
Jonatan Werpers <jonatan@werpers.com>
parents: 8
diff changeset
47 return nothing
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
48 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
50 # u = L*v
87
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
51 using UnsafeArrays
81
7f72e7e14659 Add benchmarktest and mark all apply functions with @inline and @inbounds
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 80
diff changeset
52 @inline function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
88
170e5447bc19 Reduce allocations
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 87
diff changeset
53 fill!(u,0)
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
54 h = Grid.spacings(L.grid)
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
55
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
56 li = LinearIndices(L.grid.numberOfPointsPerDim)
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
57 n_x, n_y = L.grid.numberOfPointsPerDim
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
58
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
59
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
60 # For each x
90
2882e1318cc3 Add TODO:s on how to remove allocations in apply! for Laplace2d
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 89
diff changeset
61 #TODO: Can remove temporary if we pass u and indices to apply!
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
62 temp = zeros(eltype(u), n_y)
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
63 for i ∈ 1:n_x
87
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
64 @inbounds indices = uview(li,i,:)
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
65 @inbounds apply!(L.op, temp, uview(v, indices), h[2])
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
66 for i ∈ eachindex(indices)
89
c0729ade65da Merge with default
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 88 84
diff changeset
67 @inbounds u[indices[i]] = temp[i]
87
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
68 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
69 end
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
70
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
71 # For each y
90
2882e1318cc3 Add TODO:s on how to remove allocations in apply! for Laplace2d
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 89
diff changeset
72 #TODO: Can remove temporary if we pass u and indices to apply!
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
73 temp = zeros(eltype(u), n_x)
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
74 for i ∈ 1:n_y
87
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
75 @inbounds indices = uview(li,:,i)
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
76 @inbounds apply!(L.op, temp, uview(v, indices), h[1])
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
77 for i ∈ eachindex(indices)
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
78 @inbounds u[indices[i]] += temp[i]
38733e84ef1a Clean up bounds checking in stencil. change to using uview in DiffOp
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 81
diff changeset
79 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
80 end
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
81
88
170e5447bc19 Reduce allocations
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 87
diff changeset
82 for i ∈ eachindex(u)
170e5447bc19 Reduce allocations
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 87
diff changeset
83 @inbounds u[i] = L.a*u[i]
170e5447bc19 Reduce allocations
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 87
diff changeset
84 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
85
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
86 return nothing
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
87 end