Mercurial > repos > public > sbplib_julia

annotate diffOp.jl @ 108:d0a28888528a cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change input type of apply(::Laplace) to ::DiffOpCartesian
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 08 Feb 2019 12:47:21 +0100
parents e3e32299fe72
children b3fbef345810
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
108
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
47 function apply!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
48 apply!(D, u, v, Lower, Lower)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
49 apply!(D, u, v, Lower, Interior)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
50 apply!(D, u, v, Lower, Upper)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
51 apply!(D, u, v, Interior, Lower)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
52 apply!(D, u, v, Interior, Interior)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
53 apply!(D, u, v, Interior, Upper)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
54 apply!(D, u, v, Upper, Lower)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
55 apply!(D, u, v, Upper, Interior)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
56 apply!(D, u, v, Upper, Upper)
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
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
60 function apply!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
61 N = D.grid.numberOfPointsPerDim
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
62 closuresize = closureSize(D.op)
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
63 for I ∈ regionindices(N, closuresize, (r1,r2))
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
64 @inbounds indextuple = (Index(I[1], r1), Index(I[2], r2))
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
65 @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
66 end
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
67 return nothing
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
68 end
d0a28888528a Change input type of apply(::Laplace) to ::DiffOpCartesian
Jonatan Werpers <jonatan@werpers.com>
parents: 107
diff changeset
69
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
70 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
71 u = zeros(eltype(v), size(v))
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
72 apply!(D,v,u)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
73 return u
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
74 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
75
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
76 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
77 grid::Grid.EquidistantGrid{Dim,T}
76
81d9510cb2d0 Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents: 73
diff changeset
78 a::T
84
48079bd39969 Change to using tuples in stencils and ops
Jonatan Werpers <jonatan@werpers.com>
parents: 76
diff changeset
79 op::D2{Float64,N,M,K}
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
80 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
81
95
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
82 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
83 error("not implemented")
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
84 end
9d53ecca34f7 Switch to using cartesian indicies
Jonatan Werpers <jonatan@werpers.com>
parents: 85
diff changeset
85
6
cb8e50ca9e15 Add attempt att apply methods for Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 2
diff changeset
86 # u = L*v
65
7054230b639c Make dimension a type parameter in Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
87 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
88 h = Grid.spacings(L.grid)[1]
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
89 uᵢ = L.a * apply(L.op, h, v, i)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
90 return uᵢ
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
91 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
92
82
45dece5e4928 Use unsafe views
Jonatan Werpers <jonatan@werpers.com>
parents: 70
diff changeset
93 using UnsafeArrays
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
94 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
95 h = Grid.spacings(L.grid)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
96 # 2nd x-derivative
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
97 @inbounds vx = uview(v, :, Int(I[2]))
103
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
98 @inbounds uᵢ = L.a*apply(L.op, h[1], vx , I[1])
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
99 # 2nd y-derivative
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
100 @inbounds vy = uview(v, Int(I[1]), :)
103
a274d6384e91 Apply 2d Laplace one region at a time (removing the need for branching in the innermost loop)
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 99
diff changeset
101 @inbounds uᵢ += L.a*apply(L.op, h[2], vy, I[2])
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
102 return uᵢ
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
103 end
99
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
104
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
105 # Slow but maybe convenient?
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
106 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
107 I = Index{Unknown}.(Tuple(i))
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
108 apply(L, v, I)
6b6d680f2e25 Allow apply(::Laplace) to take Index types
Jonatan Werpers <jonatan@werpers.com>
parents: 95
diff changeset
109 end