Merge pull request #227 from fgerick/irregular_grid_fix

Irregular grid fix

Author: ChrisRackauckas

src/derivative_operators/concretization.jl

@@ -28,7 +28,7 @@ function Base.copyto!(L::AbstractMatrix{T}, A::DerivativeOperator{T}, N::Int) wh
 
     for i in bl+1:N-bl
         cur_coeff   = get_coeff(i)
-        stencil     = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
+        stencil     = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-bl] : A.stencil_coefs
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
         L[i,i+1-stencil_pivot:i-stencil_pivot+stencil_length] = cur_coeff * cur_stencil
     end

src/derivative_operators/derivative_operator.jl

@@ -67,6 +67,18 @@ function CenteredDifference{N}(derivative_order::Int,
         )
 end
 
+function generate_coordinates(i::Int, stencil_x, dummy_x, dx::AbstractVector{T}) where T <: Real
+    len = length(stencil_x)
+    stencil_x .= stencil_x.*zero(T)
+    for idx in 1:div(len,2)
+        shifted_idx1 = index(idx, len)
+        shifted_idx2 = index(-idx, len)
+        stencil_x[shifted_idx1] = stencil_x[shifted_idx1-1] + dx[i+idx-1]
+        stencil_x[shifted_idx2] = stencil_x[shifted_idx2+1] - dx[i-idx]
+    end
+    return stencil_x
+end
+
 function CenteredDifference{N}(derivative_order::Int,
                             approximation_order::Int, dx::AbstractVector{T},
                             len::Int, coeff_func=nothing) where {T<:Real,N}
@@ -78,28 +90,18 @@ function CenteredDifference{N}(derivative_order::Int,
 
     interior_x              = boundary_point_count+2:len+1-boundary_point_count
     dummy_x                 = -div(stencil_length,2) : div(stencil_length,2)-1
-    boundary_x              = -boundary_stencil_length+1:0
-
+    low_boundary_x          = [zero(T); cumsum(dx[1:boundary_stencil_length-1])]
+    high_boundary_x         = cumsum(dx[end-boundary_stencil_length+1:end])
     # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point.
     deriv_spots             = (-div(stencil_length,2)+1) : -1
-    boundary_deriv_spots    = boundary_x[2:div(stencil_length,2)]
-
-    function generate_coordinates(i, stencil_x, dummy_x, dx)
-        len = length(stencil_x)
-        stencil_x .= stencil_x.*zero(T)
-        for idx in 1:div(len,2)
-            shifted_idx1 = index(idx, len)
-            shifted_idx2 = index(-idx, len)
-            stencil_x[shifted_idx1] = stencil_x[shifted_idx1-1] + dx[i+idx-1]
-            stencil_x[shifted_idx2] = stencil_x[shifted_idx2+1] - dx[i-idx]
-        end
-        return stencil_x
-    end
 
-    stencil_coefs           = convert(SVector{length(interior_x)}, [convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), generate_coordinates(i, stencil_x, dummy_x, dx))) for i in interior_x])
-    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots]
+    stencil_coefs           = [convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), generate_coordinates(i, stencil_x, dummy_x, dx))) for i in interior_x]
+    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T},
+                                                                  calculate_weights(derivative_order, low_boundary_x[i+1], low_boundary_x)) for i in 1:boundary_point_count]
     low_boundary_coefs      = convert(SVector{boundary_point_count},_low_boundary_coefs)
-    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]))
+    _high_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T},
+                                                                  calculate_weights(derivative_order, high_boundary_x[end-i], high_boundary_x)) for i in boundary_point_count:-1:1]
+    high_boundary_coefs      = convert(SVector{boundary_point_count},_high_boundary_coefs)
 
     coefficients            = coeff_func isa Nothing ? nothing : Vector{T}(undef,len)
 

test/derivative_operators_interface.jl

@@ -146,55 +146,56 @@ end
 
 @testset "Correctness of Non-Uniform Stencils" begin
     x = [0., 0.08, 0.1, 0.15, 0.19, 0.26, 0.29]
+    nx = length(x)
     dx = diff(x)
 
     # Second-Order First Derivative
-    L = CenteredDifference(1, 2, dx, 5)
+    L = CenteredDifference(1, 2, dx, nx-2)
     correct = analyticCtrOneTwoIrr()
 
     # Check that stencils agree with correct
     for (i,coefs) in enumerate(L.stencil_coefs)
         @test Array(coefs) ≈ correct[i,correct[i,:] .!= 0.]
     end
-    @test_broken Array(L) ≈ correct # All of these concretizations
-    @test_broken sparse(L) ≈ correct # only give the first three
-    @test_broken BandedMatrix(L) ≈ correct # rows of the computed stencil coefficients
+    @test Array(L) ≈ correct
+    @test sparse(L) ≈ correct
+    @test BandedMatrix(L) ≈ correct
 
     # Second-Order Second Derivative
-    L = CenteredDifference(2, 2, dx, 5)
+    L = CenteredDifference(2, 2, dx, nx-2)
     correct = analyticCtrTwoTwoIrr()
 
     # Check that stencils agree with correct
     for (i,coefs) in enumerate(L.stencil_coefs)
         @test Array(coefs) ≈ correct[i,correct[i,:] .!= 0.]
     end
-    @test_broken Array(L) ≈ correct # same issue as previous derivative
-    @test_broken sparse(L) ≈ correct
-    @test_broken BandedMatrix(L) ≈ correct
+    @test Array(L) ≈ correct
+    @test sparse(L) ≈ correct
+    @test BandedMatrix(L) ≈ correct
 
     # Fourth-Order Second Derivative
-    L = CenteredDifference(2, 4, dx, 5)
+    L = CenteredDifference(2, 4, dx, nx-2)
     correct = analyticCtrTwoFourIrr()
 
     # Check that stencils agree with correct
     for (i,coefs) in enumerate(L.stencil_coefs)
         @test Array(coefs) ≈ correct[i,correct[i,:] .!= 0.]
     end
-    @test_broken Array(L) ≈ correct # L.stencil_coefs is populated, but the concretization doesn't work. It appears to be an issue of improper calculation of indexing from the various lengths computed during construction (e.g. boundary_stencil_length, len) and potentially the fact that "len" doesn't seem to specify the number of grid points at which we compute finite differences but appears to specify the location of the last grid point at which we compute finite differences (so if X is a 5-length vector, entering len = 2 means computing FDs for X[2] and X[3])
-    @test_broken sparse(L) ≈ correct
-    @test_broken BandedMatrix(L) ≈ correct
+    @test Array(L)[2:end-1,:] ≈ correct
+    @test sparse(L)[2:end-1,:] ≈ correct
+    @test BandedMatrix(L)[2:end-1,:] ≈ correct
 
     # Second-Order Fourth Derivative
-    L = CenteredDifference(4, 2, dx, 5)
+    L = CenteredDifference(4, 2, dx, nx-2)
     correct = analyticCtrFourTwoIrr()
 
     # Check that stencils agree with correct
     for (i,coefs) in enumerate(L.stencil_coefs)
         @test Array(coefs) ≈ correct[i,correct[i,:] .!= 0.]
     end
-    @test_broken Array(L) ≈ correct
-    @test_broken sparse(L) ≈ correct
-    @test_broken BandedMatrix(L) ≈ correct
+    @test Array(L)[2:end-1,:] ≈ correct
+    @test sparse(L)[2:end-1,:] ≈ correct
+    @test BandedMatrix(L)[2:end-1,:] ≈ correct
 end
 
 # tests for full and sparse function

test/generic_operator_validation.jl

@@ -1,30 +1,66 @@
-using DiffEqOperators
-
-n = 100
-x=0.0:0.005:2π
-xprime = x[2:(end-1)]
-dx=diff(x)
-y = exp.(π*x)
-y_im = exp.(π*im*x)
-yim_ = y_im[2:(end-1)]
+using DiffEqOperators, SparseArrays, LinearAlgebra
+
+x = 0.0:0.01:π
+x_ = x[2:end-1]
+dx = diff(x)
+y = sin.(2x)
 y_ = y[2:(end-1)]
+dy = [2cos.(2x_),-4sin.(2x_),-8cos.(2x_),16y_,32cos.(2x_)]
+
+
+
+# test concretization against regular grid operator
+for dor in 1:4, aor in 2:2:6
+    Dr = CenteredDifference(dor,aor,dx[1],length(x)-2)
+    Dir = CenteredDifference(dor,aor,dx,length(x)-2)
+
+    @test sparse(Dr)≈sparse(Dir)
+    @test Array(Dr)≈Array(Dir)
 
-@test_broken for dor in 1:6, aor in 1:6
+end
+
+# test irregular grid operator with regular grid
+for dor in 1:4, aor in 2:6
 
-    D1 = CenteredDifference(dor,aor,dx,length(x))
+    D1 = CenteredDifference(dor,aor,dx,length(x)-2)
 
     #take derivative
     yprime1 = D1*y
 
     #test result
-    @test_broken yprime1 ≈ (π^dor)*y_ # test operator with known derivative of exp(kx)
+    @test yprime1 ≈ dy[dor] atol = 10.0^(1-aor)#2test operator with known derivative of exp(kx)
+
+    #TODO: implement specific tests for the left and right boundary regions, waiting until after update
+end
 
-    #take derivatives
-    y_imprime1 = D1*y_im
+
+# test irregular grid
+
+x = sin.(0.0:0.05:π)
+x = cumsum(x)
+x = x/x[end]*π
+x_ = x[2:end-1]
+dx = diff(x)
+y  = sin.(2x)
+y_ = y[2:(end-1)]
+dy = [2cos.(2x_), -4sin.(2x_), -8cos.(2x_), 16y_, 32cos.(2x_)]
+
+for dor in 1:4, aor in 4:10
+
+    D1 = CenteredDifference(dor,aor,dx,length(x)-2)
+
+    #take derivative
+    yprime1 = D1*y
 
     #test result
-    @test_broken y_imprime1 ≈ ((pi*im)^dor)*yim_ # test operator with known derivative of exp(jkx)
+    tol = 2*10.0^(2-aor)*maximum(dx)^(2-dor)
+    # error estimate is fairly difficult for high order derivatives and small dx.
+
+    # err = norm(yprime1.-dy[dor])
+    # @show err
+    # @show tol
 
+    @test yprime1 ≈ dy[dor] atol = tol #2test operator with known derivative of exp(kx)
 
     #TODO: implement specific tests for the left and right boundary regions, waiting until after update
 end