Mercurial > repos > public > sbplib_julia

annotate DiffOps/src/DiffOps.jl @ 226:f8b452336b19 package_refactor

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 Wed, 26 Jun 2019 14:21:37 +0200
parents b3506cfbb9d8
children eb8525066f9b 5acef2d5db2e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
215
3a93d8a799ce Add missing module in DiffOps
Jonatan Werpers <jonatan@werpers.com>
parents: 211
diff changeset
1 module DiffOps
3a93d8a799ce Add missing module in DiffOps
Jonatan Werpers <jonatan@werpers.com>
parents: 211
diff changeset
2
3a93d8a799ce Add missing module in DiffOps
Jonatan Werpers <jonatan@werpers.com>
parents: 211
diff changeset
3 using RegionIndices
221
235f0a771c8f Make all packages load properly
Jonatan Werpers <jonatan@werpers.com>
parents: 215
diff changeset
4 using SbpOperators
235f0a771c8f Make all packages load properly
Jonatan Werpers <jonatan@werpers.com>
parents: 215
diff changeset
5 using Grids
215
3a93d8a799ce Add missing module in DiffOps
Jonatan Werpers <jonatan@werpers.com>
parents: 211
diff changeset
6
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
7 export Laplace
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
8
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
9 abstract type DiffOp end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
10
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
11 # 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
12 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
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
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
16 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
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
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
20 function matrixRepresentation(D::DiffOp)
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
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
24 abstract type DiffOpCartesian{Dim} <: DiffOp end
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
25
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
26 # DiffOp must have a grid of dimension Dim!!!
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
27 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
28 for I ∈ eachindex(D.grid)
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
29 u[I] = apply(D, v, I)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
30 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
31
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
32 return nothing
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
33 end
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
34 export apply!
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
35
112
98c788cba9bf Rename the apply! that applies in regions to apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 111
diff changeset
36 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
37 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
38 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
39 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
40 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
41 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
42 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
43 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
44 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
45 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
46 return nothing
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
47 end
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
48
114
d24497780ebd Add comment about maybe splitting apply_region!
Jonatan Werpers <jonatan@werpers.com>
parents: 113
diff changeset
49 # 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
50 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
51 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
52 @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
53 @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
54 end
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
55 return nothing
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
56 end
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
57 export apply_region!
108
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
58
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
59 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
60 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
61 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
62 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
63 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
64 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
65 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
66 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
67 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
68 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
69 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
70 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
71
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 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
73 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
74 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
75 # 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
76 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
77 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
78 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
79 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
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 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
82 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
83 end
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
84 export apply_region_tiled!
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
85
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
86 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
87 u = zeros(eltype(v), size(v))
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
88 apply!(D,v,u)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
89 return u
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
90 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
91
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
92 export apply
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
93
168
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
94 struct NormalDerivative{N,M,K}
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
95 op::D2{Float64,N,M,K}
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
96 grid::EquidistantGrid
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
97 bId::CartesianBoundary
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
98 end
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
99
169
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
100 function apply_transpose(d::NormalDerivative, v::AbstractArray, I::Integer)
172
a8bc71608588 Fix bug in dimension selection in apply_tanspose of boundary operators
Jonatan Werpers <jonatan@werpers.com>
parents: 171
diff changeset
101 u = selectdim(v,3-dim(d.bId),I)
171
d407611ed71a Add BondaryValue operator
Jonatan Werpers <jonatan@werpers.com>
parents: 169
diff changeset
102 return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
169
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
103 end
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
104
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
105 # Not correct abstraction level
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
106 # TODO: Not type stable D:<
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
107 function apply(d::NormalDerivative, v::AbstractArray, I::Tuple{Integer,Integer})
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
108 i = I[dim(d.bId)]
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
109 j = I[3-dim(d.bId)]
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
110 N_i = d.grid.size[dim(d.bId)]
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
111
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
112 r = getregion(i, closureSize(d.op), N_i)
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
113
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
114 if r != region(d.bId)
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
115 return 0
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
116 end
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
117
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
118 if r == Lower
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
119 # Note, closures are indexed by offset. Fix this D:<
171
d407611ed71a Add BondaryValue operator
Jonatan Werpers <jonatan@werpers.com>
parents: 169
diff changeset
120 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j]
169
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
121 elseif r == Upper
171
d407611ed71a Add BondaryValue operator
Jonatan Werpers <jonatan@werpers.com>
parents: 169
diff changeset
122 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-j]*v[j]
169
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
123 end
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
124 end
24ee4def7ffb Move NormalDerivative methods to type definition
Jonatan Werpers <jonatan@werpers.com>
parents: 168
diff changeset
125
173
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
126 struct BoundaryValue{N,M,K}
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
127 op::D2{Float64,N,M,K}
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
128 grid::EquidistantGrid
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
129 bId::CartesianBoundary
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
130 end
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
131
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
132 function apply(e::BoundaryValue, v::AbstractArray, I::Tuple{Integer,Integer})
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
133 i = I[dim(e.bId)]
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
134 j = I[3-dim(e.bId)]
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
135 N_i = e.grid.size[dim(e.bId)]
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
136
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
137 r = getregion(i, closureSize(e.op), N_i)
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
138
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
139 if r != region(e.bId)
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
140 return 0
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
141 end
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
142
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
143 if r == Lower
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
144 # Note, closures are indexed by offset. Fix this D:<
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
145 return e.op.eClosure[i-1]*v[j]
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
146 elseif r == Upper
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
147 return e.op.eClosure[N_i-j]*v[j]
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
148 end
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
149 end
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
150
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
151 function apply_transpose(e::BoundaryValue, v::AbstractArray, I::Integer)
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
152 u = selectdim(v,3-dim(e.bId),I)
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
153 return apply_e(e.op, u, region(e.bId))
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
154 end
fabd475bb258 Move definition of BoundaryValue before definition of Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 172
diff changeset
155
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
156 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
168
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
157 grid::EquidistantGrid{Dim,T}
76
81d9510cb2d0 Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents: 73
diff changeset
158 a::T
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
159 op::D2{Float64,N,M,K}
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
160 # e::BoundaryValue
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
161 # d::NormalDerivative
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
162 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
163
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
164 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
165 error("not implemented")
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
166 end
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
167
6
cb8e50ca9e15 Add attempt att apply methods for Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 2
diff changeset
168 # u = L*v
65
7054230b639c Make dimension a type parameter in Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
169 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
170 uᵢ = L.a * SbpOperators.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
171 return uᵢ
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
172 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
173
129
1aaeb46ba5f4 Improve efficiency of apply by the following:
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 124
diff changeset
174 @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
175 # 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
176 @inbounds vx = view(v, :, Int(I[2]))
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
177 @inbounds uᵢ = L.a*SbpOperators.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
178 # 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
179 @inbounds vy = view(v, Int(I[1]), :)
223
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
180 @inbounds uᵢ += L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
b3506cfbb9d8 Add some missing exports
Jonatan Werpers <jonatan@werpers.com>
parents: 221
diff changeset
181 # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
182 return uᵢ
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
183 end
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
184
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
185 # Slow but maybe convenient?
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
186 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
187 I = Index{Unknown}.(Tuple(i))
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
188 apply(L, v, I)
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
189 end
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
190
154
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
191 struct BoundaryOperator
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
192
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
193 end
3193bac1c086 More sketching of how things might work
Jonatan Werpers <jonatan@werpers.com>
parents: 153
diff changeset
194
168
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
195
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
196 """
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
197 A BoundaryCondition should implement the method
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
198 sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...)
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
199 """
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
200 abstract type BoundaryCondition end
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
201
174
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
202 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
203
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
204 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
205 e = BoundaryValue(L.op, L.grid, Bid())
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
206 d = NormalDerivative(L.op, L.grid, Bid())
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
207 Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid())
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
208 # TODO: Implement BoundaryQuadrature method
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
209
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
210 return -L.Hi*e*Hᵧ*(d'*v - g)
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
211 # Need to handle d'*v - g so that it is an AbstractArray that TensorMappings can act on
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
212 end
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
213
145
e0c8f5cf3a3f Rename boundary indentifiers and use them inte the sat functions
Jonatan Werpers <jonatan@werpers.com>
parents: 143
diff changeset
214 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
215 tau::Float64
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
216 end
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
217
174
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
218 function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid}
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
219 e = BoundaryValue(L.op, L.grid, Bid())
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
220 d = NormalDerivative(L.op, L.grid, Bid())
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
221 Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid())
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
222 # TODO: Implement BoundaryQuadrature method
145
e0c8f5cf3a3f Rename boundary indentifiers and use them inte the sat functions
Jonatan Werpers <jonatan@werpers.com>
parents: 143
diff changeset
223
174
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
224 return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g)
187295479984 Sketch implementation of sat methods for Neumann and Dirichlet
Jonatan Werpers <jonatan@werpers.com>
parents: 173
diff changeset
225 # Need to handle scalar multiplication and addition of TensorMapping
135
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
226 end
bb1cc9c7877c Add outline of idea for implemenation of sats
Jonatan Werpers <jonatan@werpers.com>
parents: 129
diff changeset
227
168
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
228 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
229 # return apply(s.L, v, i) +
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
230 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) +
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
231 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) +
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
232 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) +
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
233 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i)
45840a8127d6 Fix typos and add NormalDerivative
Jonatan Werpers <jonatan@werpers.com>
parents: 166
diff changeset
234 # end
215
3a93d8a799ce Add missing module in DiffOps
Jonatan Werpers <jonatan@werpers.com>
parents: 211
diff changeset
235
3a93d8a799ce Add missing module in DiffOps
Jonatan Werpers <jonatan@werpers.com>
parents: 211
diff changeset
236 end # module