Skip to content
This repository was archived by the owner on Jul 19, 2023. It is now read-only.

Commit 3ee4441

Browse files
ChrisRackauckasChrisRackauckas
committed
add a dimension parameter to derivative operators
1 parent b169937 commit 3ee4441

File tree

4 files changed

+50
-43
lines changed

4 files changed

+50
-43
lines changed

‎src/derivative_operators/abstract_operator_functions.jl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -167,9 +167,9 @@ end
167167

168168
################################################################################
169169

170-
function *(coeff_func::Function, A::DerivativeOperator{T}) where T
170+
function *(coeff_func::Function, A::DerivativeOperator{T,N,Wind}) where {T,N,Wind}
171171
coefficients = A.coefficients === nothing ? Vector{T}(undef,A.len) : A.coefficients
172-
DerivativeOperator{T,typeof(A.dx),typeof(A.stencil_coefs),
172+
DerivativeOperator{T,N,Wind,typeof(A.dx),typeof(A.stencil_coefs),
173173
typeof(A.low_boundary_coefs),typeof(coefficients),
174174
typeof(coeff_func)}(
175175
A.derivative_order, A.approximation_order,
@@ -178,7 +178,7 @@ function *(coeff_func::Function, A::DerivativeOperator{T}) where T
178178
A.boundary_stencil_length,
179179
A.boundary_point_count,
180180
A.low_boundary_coefs,
181-
A.high_boundary_coefs,coefficients,coeff_func,A.use_winding
181+
A.high_boundary_coefs,coefficients,coeff_func
182182
)
183183
end
184184

‎src/derivative_operators/concretization.jl‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -7,18 +7,18 @@ function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T
77
stl_2 = div(stl,2)
88
for i in 1:A.boundary_point_count
99
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
10-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
10+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
1111
L[i,1:bstl] = cur_coeff * cur_stencil
1212
end
1313
for i in bl+1:N-bl
1414
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
1515
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
16-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(stencil) : stencil
16+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
1717
L[i,i+1-stl_2:i+1+stl_2] = cur_coeff * cur_stencil
1818
end
1919
for i in N-bl+1:N
2020
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
21-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
21+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
2222
L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
2323
end
2424
return L / A.dx^A.derivative_order
@@ -33,18 +33,18 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) wh
3333
coeff = A.coefficients
3434
for i in 1:A.boundary_point_count
3535
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
36-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
36+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
3737
L[i,1:bstl] = cur_coeff * cur_stencil
3838
end
3939
for i in bl+1:N-bl
4040
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
4141
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
42-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(stencil) : stencil
42+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
4343
L[i,i+1-stl_2:i+1+stl_2] = cur_coeff * cur_stencil
4444
end
4545
for i in N-bl+1:N
4646
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
47-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
47+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
4848
L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
4949
end
5050
return L / A.dx^A.derivative_order
@@ -63,18 +63,18 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) whe
6363
L = BandedMatrix{T}(Zeros(N, N+2), (max(stl-3,0,bstl),max(stl-1,0,bstl)))
6464
for i in 1:A.boundary_point_count
6565
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
66-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
66+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
6767
L[i,1:bstl] = cur_coeff * cur_stencil
6868
end
6969
for i in bl+1:N-bl
7070
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
7171
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
72-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(stencil) : stencil
72+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
7373
L[i,i+1-stl_2:i+1+stl_2] = cur_coeff * cur_stencil
7474
end
7575
for i in N-bl+1:N
7676
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
77-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
77+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
7878
L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
7979
end
8080
return L / A.dx^A.derivative_order

‎src/derivative_operators/convolutions.jl‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
1818
xtempi = zero(T)
1919
cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i] : stencil
2020
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
21-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
21+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
2222
for idx in 1:A.stencil_length
2323
xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx]
2424
end
@@ -33,7 +33,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::D
3333
xtempi = stencil[i][1]*x[1]
3434
cur_stencil = stencil[i]
3535
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
36-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
36+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
3737
for idx in 2:A.boundary_stencil_length
3838
xtempi += cur_coeff * cur_stencil[idx] * x[idx]
3939
end
@@ -48,7 +48,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
4848
xtempi = stencil[i][end]*x[end]
4949
cur_stencil = stencil[i]
5050
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
51-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
51+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
5252
for idx in (A.boundary_stencil_length-1):-1:1
5353
xtempi += cur_coeff * cur_stencil[end-idx] * x[end-idx]
5454
end
@@ -69,7 +69,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
6969
xtempi = zero(T)
7070
cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_length] : stencil
7171
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true
72-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
72+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
7373
@inbounds for idx in 1:A.stencil_length
7474
xtempi += cur_coeff * cur_stencil[idx] * x[i - (mid-idx) + 1]
7575
end
@@ -84,7 +84,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
8484
xtempi = stencil[i][1]*x.l
8585
cur_stencil = stencil[i]
8686
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true
87-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
87+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
8888
@inbounds for idx in 2:A.stencil_length
8989
xtempi += cur_coeff * cur_stencil[idx] * x.u[idx-1]
9090
end
@@ -100,7 +100,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
100100
xtempi = stencil[i][end]*x.r
101101
cur_stencil = stencil[i]
102102
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true
103-
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
103+
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
104104
@inbounds for idx in A.stencil_length:-1:2
105105
xtempi += cur_coeff * cur_stencil[end-idx] * x.u[end-idx+1]
106106
end

‎src/derivative_operators/derivative_operator.jl‎

Lines changed: 32 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
struct DerivativeOperator{T<:Real,T2,S1,S2<:SVector,T3,F} <: AbstractDerivativeOperator{T}
1+
struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SVector,T3,F} <: AbstractDerivativeOperator{T}
22
derivative_order :: Int
33
approximation_order :: Int
44
dx :: T2
@@ -11,12 +11,13 @@ struct DerivativeOperator{T<:Real,T2,S1,S2<:SVector,T3,F} <: AbstractDerivativeO
1111
high_boundary_coefs :: S2
1212
coefficients :: T3
1313
coeff_func :: F
14-
use_winding :: Bool
1514
end
1615

17-
function CenteredDifference(derivative_order::Int,
16+
struct CenteredDifference{N} end
17+
18+
function CenteredDifference{N}(derivative_order::Int,
1819
approximation_order::Int, dx::T,
19-
len::Int, coeff_func=nothing) where {T<:Real}
20+
len::Int, coeff_func=nothing) where {T<:Real,N}
2021

2122
stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
2223
boundary_stencil_length = derivative_order + approximation_order
@@ -33,21 +34,21 @@ function CenteredDifference(derivative_order::Int,
3334
high_boundary_coefs = convert(SVector{boundary_point_count},reverse(SVector{boundary_stencil_length, T}[reverse(low_boundary_coefs[i]) for i in 1:boundary_point_count]))
3435

3536
coefficients = coeff_func isa Nothing ? nothing : Vector{T}(undef,len)
36-
DerivativeOperator{T,T,typeof(stencil_coefs),
37-
typeof(low_boundary_coefs),typeof(coefficients),
38-
typeof(coeff_func)}(
37+
DerivativeOperator{T,N,false,T,typeof(stencil_coefs),
38+
typeof(low_boundary_coefs),typeof(coefficients),
39+
typeof(coeff_func)}(
3940
derivative_order, approximation_order, dx, len, stencil_length,
4041
stencil_coefs,
4142
boundary_stencil_length,
4243
boundary_point_count,
4344
low_boundary_coefs,
44-
high_boundary_coefs,coefficients,coeff_func,false
45+
high_boundary_coefs,coefficients,coeff_func
4546
)
4647
end
4748

48-
function CenteredDifference(derivative_order::Int,
49+
function CenteredDifference{N}(derivative_order::Int,
4950
approximation_order::Int, dx::AbstractVector{T},
50-
len::Int, coeff_func=nothing) where {T<:Real}
51+
len::Int, coeff_func=nothing) where {T<:Real,N}
5152

5253
stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
5354
boundary_stencil_length = derivative_order + approximation_order
@@ -65,9 +66,9 @@ function CenteredDifference(derivative_order::Int,
6566

6667
coefficients = coeff_func isa Nothing ? nothing : Vector{T}(undef,len)
6768

68-
DerivativeOperator{T,typeof(dx),typeof(stencil_coefs),
69-
typeof(low_boundary_coefs),typeof(coefficients),
70-
typeof(coeff_func)}(
69+
DerivativeOperator{T,N,false,typeof(dx),typeof(stencil_coefs),
70+
typeof(low_boundary_coefs),typeof(coefficients),
71+
typeof(coeff_func)}(
7172
derivative_order, approximation_order, dx,
7273
len, stencil_length,
7374
stencil_coefs,
@@ -78,9 +79,11 @@ function CenteredDifference(derivative_order::Int,
7879
)
7980
end
8081

81-
function UpwindDifference(derivative_order::Int,
82+
struct UpwindDifference{N} end
83+
84+
function UpwindDifference{N}(derivative_order::Int,
8285
approximation_order::Int, dx::T,
83-
len::Int, coeff_func=nothing) where {T<:Real}
86+
len::Int, coeff_func=nothing) where {T<:Real,N}
8487

8588
stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
8689
boundary_stencil_length = derivative_order + approximation_order
@@ -98,21 +101,21 @@ function UpwindDifference(derivative_order::Int,
98101

99102
coefficients = Vector{T}(undef,len)
100103

101-
DerivativeOperator{T,T,typeof(stencil_coefs),
102-
typeof(low_boundary_coefs),Vector{T},
103-
typeof(coeff_func)}(
104+
DerivativeOperator{T,N,true,T,typeof(stencil_coefs),
105+
typeof(low_boundary_coefs),Vector{T},
106+
typeof(coeff_func)}(
104107
derivative_order, approximation_order, dx, len, stencil_length,
105108
stencil_coefs,
106109
boundary_stencil_length,
107110
boundary_point_count,
108111
low_boundary_coefs,
109-
high_boundary_coefs,coefficients,coeff_func,true
112+
high_boundary_coefs,coefficients,coeff_func
110113
)
111114
end
112115

113-
function UpwindDifference(derivative_order::Int,
116+
function UpwindDifference{N}(derivative_order::Int,
114117
approximation_order::Int, dx::AbstractVector{T},
115-
len::Int, coeff_func=nothing) where {T<:Real}
118+
len::Int, coeff_func=nothing) where {T<:Real,N}
116119

117120
stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
118121
boundary_stencil_length = derivative_order + approximation_order
@@ -130,14 +133,18 @@ function UpwindDifference(derivative_order::Int,
130133

131134
coefficients = Vector{T}(undef,len)
132135

133-
DerivativeOperator{T,typeof(dx),typeof(stencil_coefs),
134-
typeof(low_boundary_coefs),Vector{T},
135-
typeof(coeff_func)}(
136+
DerivativeOperator{T,N,true,typeof(dx),typeof(stencil_coefs),
137+
typeof(low_boundary_coefs),Vector{T},
138+
typeof(coeff_func)}(
136139
derivative_order, approximation_order, dx, len, stencil_length,
137140
stencil_coefs,
138141
boundary_stencil_length,
139142
boundary_point_count,
140143
low_boundary_coefs,
141-
high_boundary_coefs,coefficients,coeff_func,true
144+
high_boundary_coefs,coefficients,coeff_func
142145
)
143146
end
147+
148+
CenteredDifference(args...) = CenteredDifference{1}(args...)
149+
UpwindDifference(args...) = UpwindDifference{1}(args...)
150+
use_winding(A::DerivativeOperator{T,N,Wind}) where {T,N,Wind} = Wind

0 commit comments

Comments
 (0)