Merge pull request #137 from JuliaDiffEq/stencil_scaling

Stencil scaling for regular grids.

Author: ChrisRackauckas

src/derivative_operators/concretization.jl

@@ -21,7 +21,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
         L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
     end
-    return L / A.dx^A.derivative_order
+    return L
 end
 
 function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) where T
@@ -47,7 +47,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) wh
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
         L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
     end
-    return L / A.dx^A.derivative_order
+    return L
 end
 
 function SparseArrays.sparse(A::AbstractDerivativeOperator{T}, N::Int=A.len) where T
@@ -77,7 +77,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) whe
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
         L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
     end
-    return L / A.dx^A.derivative_order
+    return L 
 end
 
 function Base.convert(::Type{Array},A::DerivativeOperator{T}) where T

src/derivative_operators/convolutions.jl

@@ -3,7 +3,6 @@ function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::DerivativeOperator, x:
     convolve_BC_left!(x_temp, x, A)
     convolve_interior!(x_temp, x, A)
     convolve_BC_right!(x_temp, x, A)
-    rmul!(x_temp, @.(1/(A.dx^A.derivative_order)))
 end
 
 ################################################

src/derivative_operators/derivative_operator.jl

@@ -30,8 +30,8 @@ function CenteredDifference{N}(derivative_order::Int,
     deriv_spots             = (-div(stencil_length,2)+1) : -1
     boundary_deriv_spots    = boundary_x[2:div(stencil_length,2)]
 
-    stencil_coefs           = convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_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}, (1/dx^derivative_order) * calculate_weights(derivative_order, zero(T), dummy_x))
+    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots]
     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]))
 

src/derivative_operators/derivative_operator_functions.jl

@@ -80,7 +80,6 @@ for MT in [2,3]
                     convolve_BC_right!(view(x_temp, idx...), view(M, idx...), A)
                 end
             end
-            mul!(x_temp,x_temp,1/A.dx^A.derivative_order)
         end
     end
 end

test/bc_coeff_compositions.jl

@@ -197,19 +197,19 @@ end
 @testset "Test Left Division L4 (fourth order)" begin
 
     # Test \ homogenous and inhomogenous BC
-    dx = 0.00001
-    x = 0.0001:dx:0.01
+    dx = 0.01
+    x = 0.01:dx:0.2
     N = length(x)
     u = sin.(x)
 
     L = CenteredDifference(4, 4, dx, N)
-    Q = RobinBC(1.0, 0.0, sin(0), dx, 1.0, 0.0, sin(0.01+dx), dx)
+    Q = RobinBC(1.0, 0.0, sin(0.0), dx, 1.0, 0.0, sin(0.2+dx), dx)
     A = L*Q
 
     analytic_L = fourth_deriv_approx_stencil(N) ./ dx^4
     analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
     analytic_AL = analytic_L*analytic_QL
-    analytic_Qb = [zeros(N+1); sin(0.01+dx)]
+    analytic_Qb = [zeros(N+1); sin(0.2+dx)]
     analytic_Ab = analytic_L*analytic_Qb
 
     analytic_u = analytic_AL \ (u - analytic_Ab)