Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 1367:71e89507dd9a feature/variable_derivatives

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove size field from SecondDerivativeVariable since it duplicates info from the coefficient field
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 26 May 2023 14:53:19 +0200
parents 4684c7f1c4cb
children 26ad90b42efd
comparison
equal deleted inserted replaced
1365:4684c7f1c4cb 1367:71e89507dd9a
7 A second derivative operator in direction `Dir` with a variable coefficient. 7 A second derivative operator in direction `Dir` with a variable coefficient.
8 """ 8 """
9 struct SecondDerivativeVariable{Dir,T,D,M,IStencil<:NestedStencil{T},CStencil<:NestedStencil{T},TArray<:AbstractArray} <: LazyTensor{T,D,D} 9 struct SecondDerivativeVariable{Dir,T,D,M,IStencil<:NestedStencil{T},CStencil<:NestedStencil{T},TArray<:AbstractArray} <: LazyTensor{T,D,D}
10 inner_stencil::IStencil 10 inner_stencil::IStencil
11 closure_stencils::NTuple{M,CStencil} 11 closure_stencils::NTuple{M,CStencil}
12 size::NTuple{D,Int}
13 coefficient::TArray 12 coefficient::TArray
14 13
15 function SecondDerivativeVariable{Dir, D}(inner_stencil::NestedStencil{T}, closure_stencils::NTuple{M,NestedStencil{T}}, size::NTuple{D,Int}, coefficient::AbstractArray) where {Dir,T,D,M} 14 function SecondDerivativeVariable{Dir, D}(inner_stencil::NestedStencil{T}, closure_stencils::NTuple{M,NestedStencil{T}}, coefficient::AbstractArray) where {Dir,T,D,M}
16 IStencil = typeof(inner_stencil) 15 IStencil = typeof(inner_stencil)
17 CStencil = eltype(closure_stencils) 16 CStencil = eltype(closure_stencils)
18 TArray = typeof(coefficient) 17 TArray = typeof(coefficient)
19 return new{Dir,T,D,M,IStencil,CStencil,TArray}(inner_stencil,closure_stencils,size, coefficient) 18 return new{Dir,T,D,M,IStencil,CStencil,TArray}(inner_stencil, closure_stencils, coefficient)
20 end 19 end
21 end 20 end
22 21
23 function SecondDerivativeVariable(g::TensorGrid, coeff::AbstractArray, inner_stencil::NestedStencil, closure_stencils, dir) 22 function SecondDerivativeVariable(g::TensorGrid, coeff::AbstractArray, inner_stencil::NestedStencil, closure_stencils, dir)
24 check_coefficient(g, coeff) 23 check_coefficient(g, coeff)
25 24
26 Δxᵢ = spacing(g.grids[dir]) 25 Δxᵢ = spacing(g.grids[dir])
27 scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2) 26 scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2)
28 scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2) 27 scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2)
29 return SecondDerivativeVariable{dir, ndims(g)}(scaled_inner_stencil, scaled_closure_stencils, size(g), coeff) 28 return SecondDerivativeVariable{dir, ndims(g)}(scaled_inner_stencil, scaled_closure_stencils, coeff)
30 end 29 end
31 30
32 function SecondDerivativeVariable(g::EquidistantGrid, coeff::AbstractVector, inner_stencil::NestedStencil, closure_stencils) 31 function SecondDerivativeVariable(g::EquidistantGrid, coeff::AbstractVector, inner_stencil::NestedStencil, closure_stencils)
33 return SecondDerivativeVariable(TensorGrid(g), coeff, inner_stencil, closure_stencils, 1) 32 return SecondDerivativeVariable(TensorGrid(g), coeff, inner_stencil, closure_stencils, 1)
34 end 33 end
76 75
77 derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir 76 derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir
78 77
79 closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils) 78 closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils)
80 79
81 LazyTensors.range_size(op::SecondDerivativeVariable) = op.size 80 LazyTensors.range_size(op::SecondDerivativeVariable) = size(op.coefficient)
82 LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size 81 LazyTensors.domain_size(op::SecondDerivativeVariable) = size(op.coefficient)
83 82
84 83
85 function derivative_view(op, a, I) 84 function derivative_view(op, a, I)
86 d = derivative_direction(op) 85 d = derivative_direction(op)
87 86
108 function apply_upper(op::SecondDerivativeVariable, v, I...) 107 function apply_upper(op::SecondDerivativeVariable, v, I...)
109 ṽ = derivative_view(op, v, I) 108 ṽ = derivative_view(op, v, I)
110 c̃ = derivative_view(op, op.coefficient, I) 109 c̃ = derivative_view(op, op.coefficient, I)
111 110
112 i = I[derivative_direction(op)] 111 i = I[derivative_direction(op)]
113 stencil = op.closure_stencils[op.size[derivative_direction(op)]-i+1] 112 sz = domain_size(op)[derivative_direction(op)]
113 stencil = op.closure_stencils[sz-i+1]
114 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i) 114 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i)
115 end 115 end
116 116
117 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I::Vararg{Index}) 117 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I::Vararg{Index})
118 if I[derivative_direction(op)] isa Index{Lower} 118 if I[derivative_direction(op)] isa Index{Lower}
126 end 126 end
127 end 127 end
128 128
129 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...) 129 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...)
130 dir = derivative_direction(op) 130 dir = derivative_direction(op)
131 sz = domain_size(op)[dir]
131 132
132 i = I[dir] 133 i = I[dir]
133 134
134 I = map(i->Index(i, Interior), I) 135 I = map(i->Index(i, Interior), I)
135 if 0 < i <= closure_size(op) 136 if 0 < i <= closure_size(op)
136 I = Base.setindex(I, Index(i, Lower), dir) 137 I = Base.setindex(I, Index(i, Lower), dir)
137 return LazyTensors.apply(op, v, I...) 138 return LazyTensors.apply(op, v, I...)
138 elseif closure_size(op) < i <= op.size[dir]-closure_size(op) 139 elseif closure_size(op) < i <= sz-closure_size(op)
139 I = Base.setindex(I, Index(i, Interior), dir) 140 I = Base.setindex(I, Index(i, Interior), dir)
140 return LazyTensors.apply(op, v, I...) 141 return LazyTensors.apply(op, v, I...)
141 elseif op.size[dir]-closure_size(op) < i <= op.size[dir] 142 elseif sz-closure_size(op) < i <= sz
142 I = Base.setindex(I, Index(i, Upper), dir) 143 I = Base.setindex(I, Index(i, Upper), dir)
143 return LazyTensors.apply(op, v, I...) 144 return LazyTensors.apply(op, v, I...)
144 else 145 else
145 error("Bounds error") # TODO: Make this more standard 146 error("Bounds error") # TODO: Make this more standard
146 end 147 end
168 169
169 function apply_upper(op::SecondDerivativeVariable{1}, v, i, j) 170 function apply_upper(op::SecondDerivativeVariable{1}, v, i, j)
170 ṽ = @view v[:,j] 171 ṽ = @view v[:,j]
171 c̃ = @view op.coefficient[:,j] 172 c̃ = @view op.coefficient[:,j]
172 173
173 stencil = op.closure_stencils[op.size[derivative_direction(op)]-i+1] 174 sz = domain_size(op)[derivative_direction(op)]
175 stencil = op.closure_stencils[sz-i+1]
174 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i) 176 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i)
175 end 177 end
176 178
177 179
178 ## y-direction 180 ## y-direction
192 194
193 function apply_upper(op::SecondDerivativeVariable{2}, v, i, j) 195 function apply_upper(op::SecondDerivativeVariable{2}, v, i, j)
194 ṽ = @view v[i,:] 196 ṽ = @view v[i,:]
195 c̃ = @view op.coefficient[i,:] 197 c̃ = @view op.coefficient[i,:]
196 198
197 stencil = op.closure_stencils[op.size[derivative_direction(op)]-j+1] 199 sz = domain_size(op)[derivative_direction(op)]
200 stencil = op.closure_stencils[sz-j+1]
198 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, j) 201 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, j)
199 end 202 end