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

Commit af22cbe

Browse files
ChrisRackauckasChrisRackauckas
committed
condense operator implementations
1 parent 98d6c10 commit af22cbe

File tree

7 files changed

+175
-586
lines changed

7 files changed

+175
-586
lines changed

‎src/DiffEqOperators.jl‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -25,8 +25,6 @@ include("derivative_operators/BC_operators.jl")
2525

2626
### Derivative Operators
2727
include("derivative_operators/fornberg.jl")
28-
include("derivative_operators/upwind_operator.jl")
29-
include("derivative_operators/derivative_irreg_operator.jl")
3028
include("derivative_operators/derivative_operator.jl")
3129
include("derivative_operators/abstract_operator_functions.jl")
3230
include("derivative_operators/convolutions.jl")
@@ -44,6 +42,7 @@ end
4442

4543
export MatrixFreeOperator
4644
export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, JacVecOperator, getops
47-
export AbstractDerivativeOperator, DerivativeOperator, CenteredDifference, UpwindOperator, FiniteDifference
45+
export AbstractDerivativeOperator, DerivativeOperator,
46+
CenteredDifference, UpwindDifference
4847
export RobinBC, GeneralBC
4948
end # module

‎src/derivative_operators/abstract_operator_functions.jl‎

Lines changed: 18 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Colon) =
2020
@inline function getindex(A::AbstractDerivativeOperator, i::Int, j::Int)
2121
@boundscheck checkbounds(A, i, j)
2222
bpc = A.boundary_point_count
23-
N = A.dimension
23+
N = A.len
2424
bsl = A.boundary_stencil_length
2525
slen = A.stencil_length
2626
if bpc > 0 && 1<=i<=bpc
@@ -104,16 +104,16 @@ end
104104

105105
# Base.length(A::AbstractDerivativeOperator) = A.stencil_length
106106
Base.ndims(A::AbstractDerivativeOperator) = 2
107-
Base.size(A::AbstractDerivativeOperator) = (A.dimension, A.dimension + 2)
107+
Base.size(A::AbstractDerivativeOperator) = (A.len, A.len + 2)
108108
Base.size(A::AbstractDerivativeOperator,i::Integer) = size(A)[i]
109109
Base.length(A::AbstractDerivativeOperator) = reduce(*, size(A))
110110

111111
#=
112112
For the evenly spaced grid we have a symmetric matrix
113113
=#
114-
Base.transpose(A::Union{DerivativeOperator,UpwindOperator}) = A
115-
Base.adjoint(A::Union{DerivativeOperator,UpwindOperator}) = A
116-
LinearAlgebra.issymmetric(::Union{DerivativeOperator,UpwindOperator}) = true
114+
Base.transpose(A::DerivativeOperator) = A
115+
Base.adjoint(A::DerivativeOperator) = A
116+
LinearAlgebra.issymmetric(::DerivativeOperator) = true
117117

118118
#=
119119
Fallback methods that use the full representation of the operator
@@ -164,3 +164,16 @@ end
164164
function *(A::AbstractDerivativeOperator,B::AbstractDerivativeOperator)
165165
return BandedMatrix(A)*BandedMatrix(B)
166166
end
167+
168+
################################################################################
169+
170+
function DiffEqBase.update_coefficients!(A::AbstractDerivativeOperator,u,p,t)
171+
if A.coeff_func !== nothing
172+
A.coeff_func(A.coefficients,u,p,t)
173+
end
174+
end
175+
176+
################################################################################
177+
178+
(L::DerivativeOperator)(u,p,t) = L*u
179+
(L::DerivativeOperator)(du,u,p,t) = mul!(du,L,u)

‎src/derivative_operators/concretization.jl‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.dimension) where T
1+
function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T
22
L = zeros(T, N, N+2)
33
bl = A.boundary_point_count
44
stl = A.stencil_length
@@ -16,7 +16,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.dimension) where
1616
return L / A.dx^A.derivative_order
1717
end
1818

19-
function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.dimension) where T
19+
function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) where T
2020
L = spzeros(T, N, N+2)
2121
bl = A.boundary_point_count
2222
stl = A.stencil_length
@@ -34,12 +34,11 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.dimensi
3434
return L / A.dx^A.derivative_order
3535
end
3636

37-
function SparseArrays.sparse(A::AbstractDerivativeOperator{T}, N::Int=A.dimension) where T
37+
function SparseArrays.sparse(A::AbstractDerivativeOperator{T}, N::Int=A.len) where T
3838
SparseMatrixCSC(A,N)
3939
end
4040

41-
function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.dimension) where T
42-
N = A.dimension
41+
function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) where T
4342
bl = A.boundary_point_count
4443
stl = A.stencil_length
4544
bstl = A.boundary_stencil_length

‎src/derivative_operators/convolutions.jl‎

Lines changed: 47 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
# mul! done by convolutions
2-
function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::Union{DerivativeOperator{T},UpwindOperator{T}}, x::AbstractVector{T}) where T<:Real
2+
function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::DerivativeOperator, x::AbstractVector{T}) where T<:Real
33
convolve_BC_left!(x_temp, x, A)
44
convolve_interior!(x_temp, x, A)
55
convolve_BC_right!(x_temp, x, A)
@@ -9,102 +9,101 @@ end
99
################################################
1010

1111
# Against a standard vector, assume already padded and just apply the stencil
12-
function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
12+
function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real}
1313
@assert length(x_temp)+2 == length(x)
14-
coeffs = A.stencil_coefs
14+
stencil = A.stencil_coefs
15+
coeff = A.coefficients
1516
mid = div(A.stencil_length,2)
1617
for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)
1718
xtempi = zero(T)
19+
cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i] : stencil
20+
cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i] : true
21+
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
1822
for idx in 1:A.stencil_length
19-
xtempi += coeffs[idx] * x[i - mid + idx]
23+
xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx]
2024
end
2125
x_temp[i] = xtempi
2226
end
2327
end
2428

25-
function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
26-
coeffs = A.low_boundary_coefs
29+
function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real}
30+
stencil = A.low_boundary_coefs
31+
coeff = A.coefficients
2732
for i in 1 : A.boundary_point_count
28-
xtempi = coeffs[i][1]*x[1]
33+
xtempi = stencil[i][1]*x[1]
34+
cur_stencil = stencil[i]
35+
cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i] : true
36+
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
2937
for idx in 2:A.boundary_stencil_length
30-
xtempi += coeffs[i][idx] * x[idx]
38+
xtempi += cur_coeff * cur_stencil[idx] * x[idx]
3139
end
3240
x_temp[i] = xtempi
3341
end
3442
end
3543

36-
function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
37-
coeffs = A.high_boundary_coefs
44+
function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real}
45+
stencil = A.high_boundary_coefs
46+
coeff = A.coefficients
3847
for i in 1 : A.boundary_point_count
39-
xtempi = coeffs[i][end]*x[end]
48+
xtempi = stencil[i][end]*x[end]
49+
cur_stencil = stencil[i]
50+
cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i] : true
51+
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
4052
for idx in (A.boundary_stencil_length-1):-1:1
41-
xtempi += coeffs[i][end-idx] * x[end-idx]
53+
xtempi += cur_coeff * cur_stencil[end-idx] * x[end-idx]
4254
end
4355
x_temp[end-A.boundary_point_count+i] = xtempi
4456
end
4557
end
4658

4759
###########################################
48-
# NOT NEEDED ANYMORE
4960

5061
# Against A BC-padded vector, specialize the computation to explicitly use the left, right, and middle parts
51-
function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
52-
coeffs = A.stencil_coefs
62+
function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real}
63+
stencil = A.stencil_coefs
64+
coeff = A.coefficients
5365
x = _x.u
5466
mid = div(A.stencil_length,2) + 1
5567
# Just do the middle parts
5668
for i in (1+A.boundary_length) : (length(x_temp)-A.boundary_length)
5769
xtempi = zero(T)
70+
cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_length] : stencil
71+
cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true
72+
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
5873
@inbounds for idx in 1:A.stencil_length
59-
xtempi += coeffs[idx] * x[i - (mid-idx) + 1]
74+
xtempi += cur_coeff * cur_stencil[idx] * x[i - (mid-idx) + 1]
6075
end
6176
x_temp[i] = xtempi
6277
end
6378
end
6479

65-
function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
66-
coeffs = A.low_boundary_coefs
80+
function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real}
81+
stencil = A.low_boundary_coefs
82+
coeff = A.coefficients
6783
for i in 1 : A.boundary_length
68-
xtempi = coeffs[i][1]*x.l
84+
xtempi = stencil[i][1]*x.l
85+
cur_stencil = stencil[i]
86+
cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true
87+
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
6988
@inbounds for idx in 2:A.stencil_length
70-
xtempi += coeffs[i][idx] * x.u[idx-1]
89+
xtempi += cur_coeff * cur_stencil[idx] * x.u[idx-1]
7190
end
7291
x_temp[i] = xtempi
7392
end
7493
end
7594

76-
function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
77-
coeffs = A.low_boundary_coefs
95+
function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real}
96+
stencil = A.low_boundary_coefs
97+
coeff = A.coefficients
7898
bc_start = length(x.u) - A.stencil_length
7999
for i in 1 : A.boundary_length
80-
xtempi = coeffs[i][end]*x.r
100+
xtempi = stencil[i][end]*x.r
101+
cur_stencil = stencil[i]
102+
cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true
103+
cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
81104
@inbounds for idx in A.stencil_length:-1:2
82-
xtempi += coeffs[i][end-idx] * x.u[end-idx+1]
105+
xtempi += cur_coeff * cur_stencil[end-idx] * x.u[end-idx+1]
83106
end
84107
x_temp[i] = xtempi
85108
end
86109
end
87-
88-
##########################################
89-
90-
#= INTERIOR CONVOLUTION =#
91-
function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::UpwindOperator{T,S,LBC,RBC}) where {T<:Real,S<:SVector,LBC,RBC}
92-
N = length(x)
93-
stencil_length = length(A.up_stencil_coefs)
94-
stencil_rem = 1 - stencil_length%2
95-
Threads.@threads for i in A.boundary_point_count[1]+1 : N-A.boundary_point_count[2]
96-
xtempi = zero(T)
97-
if A.directions[][i] == false
98-
@inbounds for j in 1 : length(A.up_stencil_coefs)
99-
xtempi += A.up_stencil_coefs[j] * x[i+j-1-stencil_rem]
100-
end
101-
else
102-
@inbounds for j in -length(A.down_stencil_coefs)+1 : 0
103-
xtempi += A.down_stencil_coefs[j+stencil_length] * x[i+j+stencil_rem]
104-
# println("i = $i, j = $j, s_idx = $(j+stencil_length), x_idx = $(i+j+stencil_rem), $(A.down_stencil_coefs[j+stencil_length]) * $(x[i+j+stencil_rem]), xtempi = $xtempi")
105-
end
106-
end
107-
108-
x_temp[i] = xtempi
109-
end
110-
end

0 commit comments

Comments
 (0)