Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/laplace.jl @ 228:5acef2d5db2e boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Move Laplace operator and related structs/functions to separate file.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 26 Jun 2019 14:38:01 +0200
parents
children a5fdc00d5070
comparison
equal deleted inserted replaced
225:3ab0c61f1367 228:5acef2d5db2e
1 struct NormalDerivative{N,M,K}
2 op::D2{Float64,N,M,K}
3 grid::EquidistantGrid
4 bId::CartesianBoundary
5 end
6
7 function apply_transpose(d::NormalDerivative, v::AbstractArray, I::Integer)
8 u = selectdim(v,3-dim(d.bId),I)
9 return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
10 end
11
12 # Not correct abstraction level
13 # TODO: Not type stable D:<
14 function apply(d::NormalDerivative, v::AbstractArray, I::Tuple{Integer,Integer})
15 i = I[dim(d.bId)]
16 j = I[3-dim(d.bId)]
17 N_i = d.grid.size[dim(d.bId)]
18
19 r = getregion(i, closureSize(d.op), N_i)
20
21 if r != region(d.bId)
22 return 0
23 end
24
25 if r == Lower
26 # Note, closures are indexed by offset. Fix this D:<
27 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j]
28 elseif r == Upper
29 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-j]*v[j]
30 end
31 end
32
33 struct BoundaryValue{N,M,K}
34 op::D2{Float64,N,M,K}
35 grid::EquidistantGrid
36 bId::CartesianBoundary
37 end
38
39 function apply(e::BoundaryValue, v::AbstractArray, I::Tuple{Integer,Integer})
40 i = I[dim(e.bId)]
41 j = I[3-dim(e.bId)]
42 N_i = e.grid.size[dim(e.bId)]
43
44 r = getregion(i, closureSize(e.op), N_i)
45
46 if r != region(e.bId)
47 return 0
48 end
49
50 if r == Lower
51 # Note, closures are indexed by offset. Fix this D:<
52 return e.op.eClosure[i-1]*v[j]
53 elseif r == Upper
54 return e.op.eClosure[N_i-j]*v[j]
55 end
56 end
57
58 function apply_transpose(e::BoundaryValue, v::AbstractArray, I::Integer)
59 u = selectdim(v,3-dim(e.bId),I)
60 return apply_e(e.op, u, region(e.bId))
61 end
62
63 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
64 grid::EquidistantGrid{Dim,T}
65 a::T
66 op::D2{Float64,N,M,K}
67 # e::BoundaryValue
68 # d::NormalDerivative
69 end
70
71 function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
72 error("not implemented")
73 end
74
75 # u = L*v
76 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
77 uᵢ = L.a * SbpOperators.apply(L.op, L.grid.spacing[1], v, i)
78 return uᵢ
79 end
80
81 @inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
82 # 2nd x-derivative
83 @inbounds vx = view(v, :, Int(I[2]))
84 @inbounds uᵢ = L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[1], vx , I[1])
85 # 2nd y-derivative
86 @inbounds vy = view(v, Int(I[1]), :)
87 @inbounds uᵢ += L.a*SbpOperators.apply(L.op, L.grid.inverse_spacing[2], vy, I[2])
88 # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors
89 return uᵢ
90 end
91
92 # Slow but maybe convenient?
93 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})
94 I = Index{Unknown}.(Tuple(i))
95 apply(L, v, I)
96 end
97
98
99 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
100
101 function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}
102 e = BoundaryValue(L.op, L.grid, Bid())
103 d = NormalDerivative(L.op, L.grid, Bid())
104 Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid())
105 # TODO: Implement BoundaryQuadrature method
106
107 return -L.Hi*e*Hᵧ*(d'*v - g)
108 # Need to handle d'*v - g so that it is an AbstractArray that TensorMappings can act on
109 end
110
111 struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
112 tau::Float64
113 end
114
115 function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid}
116 e = BoundaryValue(L.op, L.grid, Bid())
117 d = NormalDerivative(L.op, L.grid, Bid())
118 Hᵧ = BoundaryQuadrature(L.op, L.grid, Bid())
119 # TODO: Implement BoundaryQuadrature method
120
121 return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g)
122 # Need to handle scalar multiplication and addition of TensorMapping
123 end
124
125 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
126 # return apply(s.L, v, i) +
127 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) +
128 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) +
129 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) +
130 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i)
131 # end