Mercurial > repos > public > sbplib_julia

annotate diffOp.jl @ 87:38733e84ef1a patch_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up bounds checking in stencil. change to using uview in DiffOp
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 25 Jan 2019 13:40:15 +0100
parents 7f72e7e14659
children 170e5447bc19
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
80
700a74c41b26 Improve type stability
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 76
diff changeset
36 struct Laplace{Dim, T<:Real} <: DiffOp
700a74c41b26 Improve type stability
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
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
73
d04569696918 Add types to Laplace operators
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
39 op::D2{Float64}
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)
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
53 u .= 0*u # Fix this?
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
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
61 temp = zeros(eltype(u), n_y)
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
62 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
63 @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
64 @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
65 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
66 @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
67 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
68 end
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
69
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
70 # For each y
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
71 temp = zeros(eltype(u), n_x)
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
72 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
73 @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
74 @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
75 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
76 @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
77 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
78 end
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
79
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
80 u .= L.a*u
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
81
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
82 return nothing
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
83 end