Mercurial > repos > public > sbplib_julia

annotate diffOp.jl @ 116:cfe7d091aca2 cell_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Close head
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 08 Feb 2019 23:57:27 +0100
parents 17b971a0b852
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
1 abstract type DiffOp end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
2
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
3 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
4 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
5 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
6
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
7 function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
8 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
9 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
10
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
11 function matrixRepresentation(D::DiffOp)
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
12 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
13 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
14
55
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
15 function boundaryCondition(D::DiffOp,b::Grid.BoundaryId,type)::(Closure, Penalty)
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
16 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
17 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
18
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
19 function interface(Du::DiffOp, Dv::DiffOp, b::Grid.BoundaryId; type)
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
20 error("not implemented")
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
21 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
22
55
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
23 abstract type Closure end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
24
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
25 function apply(c::Closure, v::AbstractVector, i::Int)
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
26 error("not implemented")
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
27 end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
28
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
29 abstract type Penalty end
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
30
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
31 function apply(c::Penalty, g, i::Int)
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
32 error("not implemented")
c62ea0112d4d Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents: 54
diff changeset
33 end
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
34
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
35 # DiffOp must have a grid!!!
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
36 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector)
60
7ab0043aca53 Only do inner stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
37 li = LinearIndices(D.grid.numberOfPointsPerDim)
7ab0043aca53 Only do inner stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
38
7ab0043aca53 Only do inner stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
39 Nx, Ny = D.grid.numberOfPointsPerDim
7ab0043aca53 Only do inner stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
40
62
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
41 h = Grid.spacings(D.grid)
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
42 li = LinearIndices(D.grid.numberOfPointsPerDim)
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
43 ci = CartesianIndices(D.grid.numberOfPointsPerDim)
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
44 is = D.op.innerStencil
60
7ab0043aca53 Only do inner stencil
Jonatan Werpers <jonatan@werpers.com>
parents: 57
diff changeset
45 for i ∈ view(li, 5:Nx-4, 5:Ny-4)
62
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
46 I = ci[i]
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
47
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
48
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
49 u[i] = zero(eltype(v))
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
50 for j ∈ is.range[1]:is.range[2]
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
51 u[i] += is[j]*v[(I[2]-1)*Nx - 1 + I[1]+j]/h[1]^2
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
52 end
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
53
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
54 for j ∈ is.range[1]:is.range[2]
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
55 u[i] += is[j]*v[(I[2]-1)*Nx - 1 + I[1]+ j*Nx]/h[2]^2
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
56 end
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
57 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
58
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
59 return nothing
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
60 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
61
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
62 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
63 u = zeros(eltype(v), size(v))
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
64 apply!(D,v,u)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
65 return u
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
66 end
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
67
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
68 # Differential operator for a*d^2/dx^2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
69 struct Laplace1D <: DiffOp
57
178a203f3e6d Add types to Laplace operators
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
70 grid::Grid.EquidistantGrid
178a203f3e6d Add types to Laplace operators
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
71 a::Real
178a203f3e6d Add types to Laplace operators
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
72 op::D2{Float64}
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
73 end
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
74
6
cb8e50ca9e15 Add attempt att apply methods for Laplace
Jonatan Werpers <jonatan@werpers.com>
parents: 2
diff changeset
75 # u = L*v
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
76 function apply(L::Laplace1D, v::AbstractVector, i::Int)
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
77 h = Grid.spacings(L.grid)[1]
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
78 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
79 return uᵢ
2
43be32298ae2 Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
80 end
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
81
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
82
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
83 # Differential operator for a*d^2/dx^2 + a*d^2/dy^2
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
84 struct Laplace2D <: DiffOp
57
178a203f3e6d Add types to Laplace operators
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
85 grid::Grid.EquidistantGrid
178a203f3e6d Add types to Laplace operators
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
86 a::Real
178a203f3e6d Add types to Laplace operators
Jonatan Werpers <jonatan@werpers.com>
parents: 56
diff changeset
87 op::D2{Float64}
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
88 end
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
89
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
90 # u = L*v
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
91 function apply(L::Laplace2D, v::AbstractVector, i::Int)
54
4300a3fbd818 switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents: 49
diff changeset
92 h = Grid.spacings(L.grid)
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
93
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
94 li = LinearIndices(L.grid.numberOfPointsPerDim)
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
95 ci = CartesianIndices(L.grid.numberOfPointsPerDim)
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
96 I = ci[i]
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
97
62
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
98 uᵢ = apply(L.op.innerStencil, view(v, li[:,I[2]]), I[1])/h[1]^2
17b971a0b852 Inline all the things
Jonatan Werpers <jonatan@werpers.com>
parents: 60
diff changeset
99 uᵢ += apply(L.op.innerStencil, view(v, li[I[1],:]), I[2])/h[2]^2
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
100
56
27a8d3021a1c Convert apply functions to cell-based
Ylva Rydin <ylva.rydin@telia.com>
parents: 55
diff changeset
101 return uᵢ
49
947f7579ba9c Add Laplace2d
Jonatan Werpers <jonatan@werpers.com>
parents: 48
diff changeset
102 end