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

Commit 2606a28

Browse files
xtalaxxtalax
authored
Merge pull request #12 from JuliaDiffEq/master
rebase
2 parents b6adfa2 + d28a0f5 commit 2606a28

10 files changed

+149
-33
lines changed

‎src/derivative_operators/concretization.jl‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -6,18 +6,18 @@ function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T
66
coeff = A.coefficients
77
stl_2 = div(stl,2)
88
for i in 1:A.boundary_point_count
9-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
9+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
1010
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
14-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
14+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
1515
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
1616
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
20-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
20+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
2121
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
@@ -32,18 +32,18 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) wh
3232
bstl = A.boundary_stencil_length
3333
coeff = A.coefficients
3434
for i in 1:A.boundary_point_count
35-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
35+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
3636
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
40-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
40+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
4141
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
4242
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
46-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
46+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
4747
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
@@ -62,22 +62,22 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) whe
6262
stl_2 = div(stl,2)
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
65-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
65+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
6666
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
70-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
70+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
7171
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
7272
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
76-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
76+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
7777
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
80-
return L
80+
return L
8181
end
8282

8383
function Base.convert(::Type{Array},A::DerivativeOperator{T}) where T

‎src/derivative_operators/convolutions.jl‎

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -15,8 +15,8 @@ function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
1515
mid = div(A.stencil_length,2)
1616
for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)
1717
xtempi = zero(T)
18-
cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i] : stencil
19-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
18+
cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil
19+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
2020
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
2121
for idx in 1:A.stencil_length
2222
xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx]
@@ -29,10 +29,10 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::D
2929
stencil = A.low_boundary_coefs
3030
coeff = A.coefficients
3131
for i in 1 : A.boundary_point_count
32-
xtempi = stencil[i][1]*x[1]
3332
cur_stencil = stencil[i]
34-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
33+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
3534
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
35+
xtempi = cur_coeff*stencil[i][1]*x[1]
3636
for idx in 2:A.boundary_stencil_length
3737
xtempi += cur_coeff * cur_stencil[idx] * x[idx]
3838
end
@@ -44,10 +44,10 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
4444
stencil = A.high_boundary_coefs
4545
coeff = A.coefficients
4646
for i in 1 : A.boundary_point_count
47-
xtempi = stencil[i][end]*x[end]
4847
cur_stencil = stencil[i]
49-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
48+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
5049
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
50+
xtempi = cur_coeff*stencil[i][end]*x[end]
5151
for idx in (A.boundary_stencil_length-1):-1:1
5252
xtempi += cur_coeff * cur_stencil[end-idx] * x[end-idx]
5353
end
@@ -67,7 +67,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
6767
for i in (2+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)-1
6868
xtempi = zero(T)
6969
cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil
70-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : true
70+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true
7171
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
7272
@inbounds for idx in 1:A.stencil_length
7373
xtempi += cur_coeff * cur_stencil[idx] * x[(i-1) - (mid-idx) + 1]
@@ -81,7 +81,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
8181
coeff = A.coefficients
8282
for i in 1 : A.boundary_point_count
8383
cur_stencil = stencil[i]
84-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : true
84+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
8585
xtempi = cur_coeff*cur_stencil[1]*_x.l
8686
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
8787
@inbounds for idx in 2:A.boundary_stencil_length
@@ -95,7 +95,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
9595
i = 1 + A.boundary_point_count
9696
xtempi = zero(T)
9797
cur_stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-A.boundary_point_count] : A.stencil_coefs
98-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : true
98+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true
9999
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
100100
xtempi = cur_coeff*cur_stencil[1]*_x.l
101101
@inbounds for idx in 2:A.stencil_length
@@ -114,7 +114,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
114114
i = length(x_temp)-A.boundary_point_count
115115
xtempi = zero(T)
116116
cur_stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-A.boundary_point_count] : A.stencil_coefs
117-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : true
117+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true
118118
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
119119
xtempi = cur_coeff*cur_stencil[end]*_x.r
120120
@inbounds for idx in 1:A.stencil_length-1
@@ -123,7 +123,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
123123
x_temp[i] = xtempi
124124
for i in 1 : A.boundary_point_count
125125
cur_stencil = stencil[i]
126-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[bc_start + i] : true
126+
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[bc_start + i] : coeff isa Number ? coeff : true
127127
xtempi = cur_coeff*cur_stencil[end]*_x.r
128128
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(cur_stencil) : cur_stencil
129129
@inbounds for idx in A.stencil_length:-1:1

‎src/derivative_operators/derivative_operator.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ struct CenteredDifference{N} end
2020
function CenteredDifference{N}(derivative_order::Int,
2121
approximation_order::Int, dx::T,
2222
len::Int, coeff_func=nothing) where {T<:Real,N}
23-
23+
@assert approximation_order>1 "approximation_order must be greater than 1."
2424
stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
2525
boundary_stencil_length = derivative_order + approximation_order
2626
dummy_x = -div(stencil_length,2) : div(stencil_length,2)

‎src/derivative_operators/derivative_operator_functions.jl‎

Lines changed: 19 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function LinearAlgebra.mul!(x_temp::AbstractArray{T}, A::DerivativeOperator{T,N}, M::AbstractArray{T}) where {T<:Real,N}
1+
function LinearAlgebra.mul!(x_temp::AbstractArray{T}, A::DerivativeOperator{T,N}, M::AbstractArray{T}) where {T,N}
22

33
# Check that x_temp has correct dimensions
44
v = zeros(ndims(x_temp))
@@ -27,8 +27,8 @@ end
2727

2828
for MT in [2,3]
2929
@eval begin
30-
function LinearAlgebra.mul!(x_temp::AbstractArray{T,$MT}, A::DerivativeOperator{T,N,false,T2,S1}, M::AbstractArray{T,$MT}) where {T<:Real,N,T2,SL,S1<:SArray{Tuple{SL},T,1,SL}}
31-
30+
function LinearAlgebra.mul!(x_temp::AbstractArray{T,$MT}, A::DerivativeOperator{T,N,false,T2,S1,S2,T3}, M::AbstractArray{T,$MT}) where
31+
{T,N,T2,SL,S1<:SArray{Tuple{SL},T,1,SL},S2,T3<:Union{Nothing,Number}}
3232
# Check that x_temp has correct dimensions
3333
v = zeros(ndims(x_temp))
3434
v[N] = 2
@@ -57,7 +57,8 @@ for MT in [2,3]
5757
W = zeros(Wdims...)
5858
Widx = Any[Wdims...]
5959
setindex!(Widx,:,N)
60-
W[Widx...] = s
60+
coeff = A.coefficients === nothing ? true : A.coefficients
61+
W[Widx...] = coeff*s
6162

6263
cv = DenseConvDims(_M, W, padding=pad,flipkernel=true)
6364
conv!(_x_temp, _M, W, cv)
@@ -91,3 +92,17 @@ function *(A::DerivativeOperator{T,N},M::AbstractArray{T}) where {T<:Real,N}
9192
LinearAlgebra.mul!(x_temp, A, M)
9293
return x_temp
9394
end
95+
96+
function *(c::Number, A::DerivativeOperator{T,N,Wind}) where {T,N,Wind}
97+
coefficients = A.coefficients === nothing ? one(T)*c : c*A.coefficients
98+
DerivativeOperator{T,N,Wind,typeof(A.dx),typeof(A.stencil_coefs),
99+
typeof(A.low_boundary_coefs),typeof(coefficients),
100+
typeof(A.coeff_func)}(
101+
A.derivative_order, A.approximation_order,
102+
A.dx, A.len, A.stencil_length,
103+
A.stencil_coefs,
104+
A.boundary_stencil_length,
105+
A.boundary_point_count,
106+
A.low_boundary_coefs,
107+
A.high_boundary_coefs,coefficients,A.coeff_func)
108+
end

‎src/matrixfree_operators.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,8 @@ function LinearAlgebra.opnorm(M::MatrixFreeOperator, p::Real)
3030
argument. E.g. `(p::Real) -> p == Inf ? 100 : error("only Inf norm is
3131
defined")`
3232
""")
33-
return M.opnorm(p)
33+
opn = M.opnorm
34+
return opn isa Number ? opn : M.opnorm(p)
3435
end
3536

3637
# Interface

‎test/derivative_operators_interface.jl‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -155,7 +155,7 @@ end
155155
@test A[60:100,500:600] == M[60:100,500:600]
156156
end
157157

158-
@testset begin "Operations on matrices"
158+
@testset "Operations on matrices" begin
159159
N = 51
160160
M = 101
161161
d_order = 2
@@ -182,6 +182,9 @@ end
182182
A*G*B
183183
end
184184

185+
# These tests are broken due to the implementation 2.2*LD creating a DerivativeOperator
186+
# rather than an Array
187+
#=
185188
@testset "Linear combinations of operators" begin
186189
N = 10
187190
Random.seed!(0); LA = DiffEqArrayOperator(rand(N,N+2))
@@ -195,8 +198,9 @@ end
195198
fullL[:,i] = L*v
196199
v[i] = 0.0
197200
end
198-
@test convert(AbstractMatrix,L) ≈ fullL
201+
@test_broken convert(AbstractMatrix,L) ≈ fullL
199202
for p in [1,2,Inf]
200-
@test opnorm(L,p) ≈ opnorm(fullL,p)
203+
@test_broken opnorm(L,p) ≈ opnorm(fullL,p)
201204
end
202205
end
206+
=#

‎test/differentiation_dimension.jl‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -299,3 +299,71 @@ end
299299
end
300300
end
301301
end
302+
303+
@testset "Differentiating with coefficients" begin
304+
305+
# Three dimensions
306+
N = 1
307+
µ = 5.5
308+
L = CenteredDifference{N}(3,4,0.1,30)
309+
310+
µL = µ*L
311+
M = zeros(32,32,32)
312+
for i in 1:32
313+
for j in 1:32
314+
for k in 1:32
315+
M[i,j,k] = cos(0.1i)
316+
end
317+
end
318+
end
319+
320+
correct_slice = 5.5*(L*M[:,1,1])
321+
M_temp = zeros(30,32,32)
322+
mul!(M_temp, µL, M)
323+
324+
for i in 1:32
325+
for j in 1:32
326+
@test M_temp[:,i,j] ≈ correct_slice
327+
end
328+
end
329+
330+
N = 2
331+
µL = µ*CenteredDifference{N}(3,4,0.1,30)
332+
M = zeros(32,32,32)
333+
for i in 1:32
334+
for j in 1:32
335+
for k in 1:32
336+
M[i,j,k] = cos(0.1j)
337+
end
338+
end
339+
end
340+
341+
M_temp = zeros(32,30,32)
342+
mul!(M_temp, µL, M)
343+
344+
for i in 1:32
345+
for j in 1:32
346+
@test M_temp[i,:,j] ≈ correct_slice
347+
end
348+
end
349+
350+
N = 3
351+
µL = µ*CenteredDifference{N}(3,4,0.1,30)
352+
M = zeros(32,32,32)
353+
for i in 1:32
354+
for j in 1:32
355+
for k in 1:32
356+
M[i,j,k] = cos(0.1k)
357+
end
358+
end
359+
end
360+
361+
M_temp = zeros(32,32,30)
362+
mul!(M_temp, µL, M)
363+
364+
for i in 1:32
365+
for j in 1:32
366+
@test M_temp[i,j,:] ≈ correct_slice
367+
end
368+
end
369+
end

‎test/error_benchmarking.jl‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
using DiffEqOperators, Plots
2+
3+
n = 10000
4+
@show n
5+
x = range(0.0; length = n, stop = 2π)
6+
dx = x[2]-x[1]
7+
y = exp.(π*x)
8+
y_ = y[2:(end-1)]
9+
10+
for dor in 1:4, aor in 2:6
11+
12+
D1 = CenteredDifference(dor,aor,dx,length(x))
13+
14+
#test result
15+
@show dor
16+
@show aor
17+
#take derivatives
18+
err_abs = abs.(D1*y .- (π^dor)*y_)
19+
err_percent = 100*err_abs./abs.(((π^dor)*y_))
20+
max_error = maximum(err_percent) # test operator with known derivative of exp(kx)
21+
avg_error = sum(err_percent)/length(err_percent)
22+
23+
@show max_error
24+
@show avg_error
25+
plot(x[2:(end-1)], err_percent, title="Percentage Error, n=$n aor = $aor, dor = $dor")
26+
27+
#test result
28+
end

‎test/jacvec_operators.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,7 @@ ff2 = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0,autodif
5757
for ff in [ff1, ff2]
5858
prob = ODEProblem(ff,u0,tspan)
5959
@test solve(prob,TRBDF2()).retcode == :Success
60-
@test solve(prob,TRBDF2(linsolve=LinSolveGMRES(tol=1e-10))).retcode == :Success
60+
@test solve(prob,TRBDF2(linsolve=LinSolveGMRES())).retcode == :Success
6161
@test solve(prob,Exprb32()).retcode == :Success
6262
@test_broken sol = solve(prob,Rosenbrock23())
6363
@test_broken sol = solve(prob,Rosenbrock23(linsolve=LinSolveGMRES(tol=1e-10)))

‎test/matrixfree.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ using DiffEqOperators, OrdinaryDiffEq
77
@. du = p*du
88
end
99
p = 1.
10-
A = MatrixFreeOperator(f, (p,), size=(5,5), opnorm=(p)->5)
10+
A = MatrixFreeOperator(f, (p,), size=(5,5), opnorm=5)
1111
b = rand(5)
1212
@test is_constant(A)
1313
prob = ODEProblem(A, b, (0,1.), p)
@@ -34,7 +34,7 @@ using DiffEqOperators, OrdinaryDiffEq
3434
mul!(df, A, f)
3535
end
3636
p = (A1,A2)
37-
O = MatrixFreeOperator(Q!, (p, 0.), size=(2,2), opnorm=(p)->10)
37+
O = MatrixFreeOperator(Q!, (p, 0.), size=(2,2))
3838
# solve DE numerically
3939
T = 2
4040
f_0 = [2.0; 1/2]

0 commit comments

Comments
 (0)