Mercurial > repos > public > sbplib_julia

annotate diffOp.jl @ 154:3193bac1c086 boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
More sketching of how things might work
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 23 Apr 2019 09:53:34 +0200
parents 754c36796ac8
children 89b63bdf1ea8
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
112
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
47 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
48 apply_region!(D, u, v, Lower, Lower)
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
49 apply_region!(D, u, v, Lower, Interior)
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
50 apply_region!(D, u, v, Lower, Upper)
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
51 apply_region!(D, u, v, Interior, Lower)
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
52 apply_region!(D, u, v, Interior, Interior)
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
53 apply_region!(D, u, v, Interior, Upper)
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
54 apply_region!(D, u, v, Upper, Lower)
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
55 apply_region!(D, u, v, Upper, Interior)
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
56 apply_region!(D, u, v, Upper, Upper)
108
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
57 return nothing
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
58 end
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
59
114
d24497780ebd Add comment about maybe splitting apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 113
diff changeset
60 # Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable
112
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
61 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
124
631eb9b35d72 Make grid spacing a property of EquidistantGrid. Create plotting module for sbplib.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 118
diff changeset
62 for I ∈ regionindices(D.grid.size, closureSize(D.op), (r1,r2))
110
ee071b8ed58c Use deafult constructor for Index in apply!
Jonatan Werpers <jonatan@werpers.com>
parents: 109
diff changeset
63 @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
108
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
64 @inbounds u[I] = apply(D, v, indextuple)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
65 end
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
66 return nothing
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
67 end
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
68
113
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
69 function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
70 apply_region_tiled!(D, u, v, Lower, Lower)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
71 apply_region_tiled!(D, u, v, Lower, Interior)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
72 apply_region_tiled!(D, u, v, Lower, Upper)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
73 apply_region_tiled!(D, u, v, Interior, Lower)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
74 apply_region_tiled!(D, u, v, Interior, Interior)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
75 apply_region_tiled!(D, u, v, Interior, Upper)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
76 apply_region_tiled!(D, u, v, Upper, Lower)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
77 apply_region_tiled!(D, u, v, Upper, Interior)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
78 apply_region_tiled!(D, u, v, Upper, Upper)
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
79 return nothing
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
80 end
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
81
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
82 using TiledIteration
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
83 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
124
631eb9b35d72 Make grid spacing a property of EquidistantGrid. Create plotting module for sbplib.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 118
diff changeset
84 ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2))
136
c6aaf061c0a9 Fix incorrect indexing of solution vector in apply_region_tiled
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 129
diff changeset
85 # TODO: Pass Tilesize to function
c6aaf061c0a9 Fix incorrect indexing of solution vector in apply_region_tiled
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 129
diff changeset
86 for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
117
ff7f377433b4 Change loop order to follow memory layout
Jonatan Werpers <jonatan@werpers.com>
parents: 114
diff changeset
87 for j ∈ tileaxs[2], i ∈ tileaxs[1]
113
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
88 I = ri[i,j]
136
c6aaf061c0a9 Fix incorrect indexing of solution vector in apply_region_tiled
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 129
diff changeset
89 u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
113
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
90 end
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
91 end
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
92 return nothing
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
93 end
3b89aa6dc7f2 Add apply_tiled! that tiles the iteration to optimize cache usage. Doesn't improve runtime at all at the moment
Jonatan Werpers <jonatan@werpers.com>
parents: 112
diff changeset
94
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
95 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
96 u = zeros(eltype(v), size(v))
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
97 apply!(D,v,u)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
98 return u
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
99 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
100
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
101 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
102 grid::Grid.EquidistantGrid{Dim,T}
76
81d9510cb2d0 Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents: 73
diff changeset
103 a::T
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
104 op::D2{Float64,N,M,K}
154
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
105 e::BoundaryValue
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
106 d::NormalDerivative
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
107 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
108
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
109 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
110 error("not implemented")
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
111 end
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
112
6
cb8e50ca9e15 Add attempt att apply methods for Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 2
diff changeset
113 # u = L*v
65
7054230b639c Make dimension a type parameter in Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
114 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
124
631eb9b35d72 Make grid spacing a property of EquidistantGrid. Create plotting module for sbplib.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 118
diff changeset
115 uᵢ = L.a * apply(L.op, L.grid.spacing[1], v, i)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
116 return uᵢ
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
117 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
118
129
1aaeb46ba5f4 Improve efficiency of apply by the following:
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 124
diff changeset
119 @inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
120 # 2nd x-derivative
118
4c0c02a80cd4 Change uview to view. It seems the compiler is now able to remove the allocation
Jonatan Werpers <jonatan@werpers.com>
parents: 117
diff changeset
121 @inbounds vx = view(v, :, Int(I[2]))
129
1aaeb46ba5f4 Improve efficiency of apply by the following:
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 124
diff changeset
122 @inbounds uᵢ = L.a*apply(L.op, L.grid.inverse_spacing[1], vx , I[1])
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
123 # 2nd y-derivative
118
4c0c02a80cd4 Change uview to view. It seems the compiler is now able to remove the allocation
Jonatan Werpers <jonatan@werpers.com>
parents: 117
diff changeset
124 @inbounds vy = view(v, Int(I[1]), :)
129
1aaeb46ba5f4 Improve efficiency of apply by the following:
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 124
diff changeset
125 @inbounds uᵢ += L.a*apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
126 return uᵢ
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
127 end
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
128
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
129 # Slow but maybe convenient?
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
130 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
131 I = Index{Unknown}.(Tuple(i))
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
132 apply(L, v, I)
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
133 end
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
134
154
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
135 struct BoundaryOperator
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
136
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
137 end
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
138
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
139 struct BoundaryValue
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
140 op::D2{Float64,N,M,K}
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
141 end
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
142
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
143 function apply(e::BoundaryValue)
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
144
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
145 end
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
146
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
147 function apply_adjoint(e::BoundaryValue)
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
148
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
149 end
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
150
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
151 struct NormalDerivative
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
152 op::D2{Float64,N,M,K}
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
153 end
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
154
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
155 function apply(e::NormalDerivative)
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
156
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
157 end
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
158
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
159 function apply_adjoint(e::NormalDerivative)
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
160
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
161 end
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
162
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
163
143
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
164 # Boundary operators
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
165
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
166 function apply_e(L::Laplace{2}, v::AbstractArray{T,2} where T, ::CartesianBoundary{1,R}, j::Int) where R
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
167 @inbounds vy = view(v, :, j)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
168 return apply_e(L.op,vy, R)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
169 end
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
170
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
171 function apply_e(L::Laplace{2}, v::AbstractArray{T,2} where T, ::CartesianBoundary{2,R}, i::Int) where R
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
172 @inbounds vx = view(v, i, :)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
173 return apply_e(L.op, vy, R)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
174 end
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
175
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
176 function apply_d(L::Laplace{2}, v::AbstractArray{T,2} where T, ::CartesianBoundary{1,R}, j::Int) where R
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
177 @inbounds vy = view(v, :, j)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
178 return apply_d(L.op,vy, R)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
179 end
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
180
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
181 function apply_d(L::Laplace{2}, v::AbstractArray{T,2} where T, ::CartesianBoundary{2,R}, i::Int) where R
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
182 @inbounds vx = view(v, i, :)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
183 return apply_d(L.op, vy, R)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
184 end
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
185
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
186
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
187 function apply_e_T(L::Laplace{2}, v::AbstractArray{T,2} where T, boundaryId, i::Int)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
188
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
189 end
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
190
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
191 function apply_d_T(L::Laplace{2}, v::AbstractArray{T,2} where T, boundaryId, i::Int)
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
192
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
193 end
755246142200 More sketching of the sat interface
Jonatan Werpers <jonatan@werpers.com>
parents: 135
diff changeset
194
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
195 """
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
196 A BoundaryCondition should implement the method
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
197 sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...)
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
198 """
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
199 abstract type BoundaryCondition end
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
200
145
e0c8f5cf3a3f Rename boundary indentifiers and use them inte the sat functions
Jonatan Werpers <jonatan@werpers.com>
parents: 143
diff changeset
201 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
202 tau::Float64
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
203 end
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
204
145
e0c8f5cf3a3f Rename boundary indentifiers and use them inte the sat functions
Jonatan Werpers <jonatan@werpers.com>
parents: 143
diff changeset
205 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
206 end
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
207
145
e0c8f5cf3a3f Rename boundary indentifiers and use them inte the sat functions
Jonatan Werpers <jonatan@werpers.com>
parents: 143
diff changeset
208 function sat(L::Laplace{2}, bc::Neumann{CartesianBoundary{1,R}}, v::AbstractArray{T,2} where T, g::AbstractVector{T}, i::CartesianIndex{2}) where R
e0c8f5cf3a3f Rename boundary indentifiers and use them inte the sat functions
Jonatan Werpers <jonatan@werpers.com>
parents: 143
diff changeset
209
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
210 # Hi * e * H_gamma * (d'*v - g)
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
211 # e, d, H_gamma applied based on bc.boundaryId
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
212 end
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
213
153
754c36796ac8 Add missing where statement
Jonatan Werpers <jonatan@werpers.com>
parents: 152
diff changeset
214 function sat(L::Laplace{2}, bc::Dirichlet{CartesianBoundary{1,R}}, v::AbstractArray{T,2} where T, g::AbstractVector{T}, i::CartesianIndex{2}) where R
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
215 # Hi * (tau/h*e + sig*d) * H_gamma * (e'*v - g)
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
216 # e, d, H_gamma applied based on bc.boundaryId
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
217 end
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
218
154
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
219 function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
220 return apply(s.L, v, i) +
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
221 sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) +
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
222 sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) +
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
223 sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) +
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
224 sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) +
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
225 end