Added half offset

Author: xtalax

src/DiffEqOperators.jl

@@ -73,7 +73,7 @@ end
 export MatrixFreeOperator
 export AnalyticalJacVecOperator, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator,
-       CenteredDifference, UpwindDifference, CompleteCenteredDifference, CompleteUpwindDifference, nonlinear_diffusion, nonlinear_diffusion!,
+       CenteredDifference, UpwindDifference, CompleteCenteredDifference, CompleteUpwindDifference, CompleteHalfCenteredDifference, nonlinear_diffusion, nonlinear_diffusion!,
        GradientOperator, Gradient, CurlOperator, Curl, DivergenceOperator, Divergence
 export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
        MultiDimDirectionalBC, ComposedMultiDimBC

src/derivative_operators/derivative_operator.jl

@@ -133,7 +133,7 @@ function CenteredDifference{N}(derivative_order::Int,
     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>
+    deriv_spots             = (-div(stencil_length,2)+1) : -1
 
     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},
@@ -163,12 +163,13 @@ function CenteredDifference{N}(derivative_order::Int,
         )
 end
 
+
+# TODO: Set weight to zero if its below a certain threshold
 """
 A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme.
 """
 function CompleteCenteredDifference(derivative_order::Int,
-    approximation_order::Int, dx::T,
-    len::Int, coeff_func=1) where {T<:Real,N}
+    approximation_order::Int, dx::T) where {T<:Real,N}
     @assert approximation_order>1 "approximation_order must be greater than 1."
     stencil_length          = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
     boundary_stencil_length = derivative_order + approximation_order
@@ -186,22 +187,22 @@ function CompleteCenteredDifference(derivative_order::Int,
     low_boundary_coefs      = convert(SVector{boundary_point_count},vcat(_low_boundary_coefs))
 
     # _high_boundary_coefs    = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, reverse(right_boundary_x))) for x0 in R_boundary_deriv_spots]
-    high_boundary_coefs      = convert(SVector{boundary_point_count},reverse(map(reverse, _low_boundary_coefs*(-1)^derivative_order)))
+    high_boundary_coefs      = convert(SVector{boundary_point_count},[reverse(v)*(-1)^derivative_order for v in _low_boundary_coefs])
 
     offside = 0
 
-    coefficients            = fill!(Vector{T}(undef,len),0)
+    coefficients            = nothing
 
 
-    DerivativeOperator{T,N,false,T,typeof(stencil_coefs),
+    DerivativeOperator{T,Nothing,false,T,typeof(stencil_coefs),
     typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),
-    typeof(coeff_func)}(
-    derivative_order, approximation_order, dx, len, stencil_length,
+    Nothing}(
+    derivative_order, approximation_order, dx, 1, stencil_length,
     stencil_coefs,
     boundary_stencil_length,
     boundary_point_count,
     low_boundary_coefs,
-    high_boundary_coefs,offside,coefficients,coeff_func
+    high_boundary_coefs,offside,coefficients,nothing
     )
 end
 
@@ -212,46 +213,45 @@ A helper function to compute the coefficients of a derivative operator including
 function CompleteHalfCenteredDifference(derivative_order::Int,
     approximation_order::Int, dx::T) where {T<:Real,N}
     @assert approximation_order>1 "approximation_order must be greater than 1."
-    stencil_length          = approximation_order + 2*floor(derivative_order/2) + 2*(approximation_order%2)
-    centered_stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
+    centered_stencil_length          = approximation_order + 2*Int(floor(derivative_order/2)) + (approximation_order%2)
     boundary_stencil_length = derivative_order + approximation_order
-    endpoint = (stencil_length-1)/2
-    dummy_x                 = -endpoint : endpoint
-    left_boundary_x         = 0:(boundary_stencil_length-1)
-    right_boundary_x        = reverse(-boundary_stencil_length+1:0)
+    endpoint = div(centered_stencil_length, 2)
+    dummy_x                 = 1-endpoint : endpoint
+    left_boundary_x         = 1:(boundary_stencil_length)
+    right_boundary_x        = reverse(-boundary_stencil_length:-1)
 
     # ? Is fornberg valid when taking an x0 outside of the stencil i.e at the boundary?
-    xoffset = convert(T, 0.5).*(-div(centered_stencil_length,2) : div(centered_stencil_length,2))
-
+    xoffset = 1.5:boundary_stencil_length+0.5
 
-    boundary_point_count    = div(stencil_length,2) # -1 due to the ghost point
+    boundary_point_count    = div(centered_stencil_length,2) # -1 due to the ghost point
     # 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  # unused
-    L_boundary_deriv_spots  = left_boundary_x[1:div(stencil_length,2)]
+    L_boundary_deriv_spots  = xoffset[1:div(centered_stencil_length,2)]
 
-    stencil_coefs           = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, zero(T), dummy_x))
+    stencil_coefs           = convert(SVector{centered_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, convert(T, 0.5), dummy_x))
     # For each boundary point, for each tappoint in the half offset central difference stencil, we need to calculate the coefficients for the stencil.
-    _low_boundary_coefs     = [Dict([offset => convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0+offset, left_boundary_x)) for offset in xoffset]) for x0 in L_boundary_deriv_spots]
-    low_boundary_coefs      = convert(SVector{boundary_point_count},vcat(_low_boundary_coefs))
 
+
+    _low_boundary_coefs     = [convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, offset, left_boundary_x)) for offset in L_boundary_deriv_spots]
+    low_boundary_coefs     = convert(SVector{boundary_point_count}, _low_boundary_coefs)
+
+    _high_boundary_coefs = [reverse(stencil)*(-1)^derivative_order for stencil in low_boundary_coefs]
+    high_boundary_coefs = convert(SVector{boundary_point_count}, _high_boundary_coefs)
     # _high_boundary_coefs    = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, reverse(right_boundary_x))) for x0 in R_boundary_deriv_spots]
-    inner_reverse(v::AbstractVector) = reverse(map.((x->reverse(x*(-1)^derivative_order),), v)) # reverse 3 ayers deep to get the high boundary cooeffs, but don't reverse in terms of the offset
-    high_boundary_coefs      = convert(SVector{boundary_point_count}, inner_reverse(_low_boundary_coefs))
 
     offside = 0
 
-    coefficients            = fill!(Vector{T}(undef,len),0)
+    coefficients            = nothing
 
 
-    DerivativeOperator{T,N,false,T,typeof(stencil_coefs),
-    typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),
-    typeof(coeff_func)}(
-    derivative_order, approximation_order, dx, len, stencil_length,
+    DerivativeOperator{T,Nothing,false,T,typeof(stencil_coefs),
+    typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),Nothing}(
+    derivative_order, approximation_order, dx, 1, centered_stencil_length,
     stencil_coefs,
     boundary_stencil_length,
     boundary_point_count,
     low_boundary_coefs,
-    high_boundary_coefs,offside,coefficients,coeff_func
+    high_boundary_coefs,offside,coefficients,nothing
     )
 end