From 4ce32b752643f1352c507234c041bfe1ed2fea0b Mon Sep 17 00:00:00 2001 From: xtalax Date: Tue, 4 Jan 2022 16:40:44 +0000 Subject: [PATCH 01/14] Added CompleteCenteredDifference for MOL --- .../derivative_operator.jl | 43 +++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 82f08cbf5..f369bf85a 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -105,6 +105,49 @@ function CenteredDifference{N}(derivative_order::Int, ) end +function CompleteCenteredDifference{N}(derivative_order::Int, + approximation_order::Int, dx::T, + len::Int, coeff_func=1) 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 +dummy_x = -div(stencil_length,2) : div(stencil_length,2) +left_boundary_x = -1:(boundary_stencil_length-1) +right_boundary_x = reverse(-boundary_stencil_length+1:1) + +boundary_point_count = div(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[2:div(stencil_length,2)] +#R_boundary_deriv_spots = right_boundary_x[2:div(stencil_length,2)] + +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, left_boundary_x)) for x0 in L_boundary_deriv_spots] +low_boundary_coefs = convert(SVector{boundary_point_count},_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))) + +offside = 0 + +coefficients = fill!(Vector{T}(undef,len),0) + +compute_coeffs!(coeff_func, coefficients) + + + +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, +stencil_coefs, +boundary_stencil_length, +boundary_point_count, +low_boundary_coefs, +high_boundary_coefs,offside,coefficients,coeff_func +) +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) From 85b8be8959354362ce8c1d8981dc2a677e89dee0 Mon Sep 17 00:00:00 2001 From: xtalax Date: Tue, 4 Jan 2022 18:35:31 +0000 Subject: [PATCH 02/14] fixed complete differences --- .../derivative_operator.jl | 131 ++++++++++++------ 1 file changed, 90 insertions(+), 41 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index f369bf85a..966d94c38 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -105,49 +105,8 @@ function CenteredDifference{N}(derivative_order::Int, ) end -function CompleteCenteredDifference{N}(derivative_order::Int, - approximation_order::Int, dx::T, - len::Int, coeff_func=1) 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 -dummy_x = -div(stencil_length,2) : div(stencil_length,2) -left_boundary_x = -1:(boundary_stencil_length-1) -right_boundary_x = reverse(-boundary_stencil_length+1:1) - -boundary_point_count = div(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[2:div(stencil_length,2)] -#R_boundary_deriv_spots = right_boundary_x[2:div(stencil_length,2)] - -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, left_boundary_x)) for x0 in L_boundary_deriv_spots] -low_boundary_coefs = convert(SVector{boundary_point_count},_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))) - -offside = 0 - -coefficients = fill!(Vector{T}(undef,len),0) - -compute_coeffs!(coeff_func, coefficients) - -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, -stencil_coefs, -boundary_stencil_length, -boundary_point_count, -low_boundary_coefs, -high_boundary_coefs,offside,coefficients,coeff_func -) -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) @@ -204,6 +163,52 @@ function CenteredDifference{N}(derivative_order::Int, ) end +""" +A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme. +""" +function CompleteCenteredDifference{N}(derivative_order::Int, + approximation_order::Int, dx::T, + len::Int, coeff_func=1) 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 + dummy_x = -div(stencil_length,2) : div(stencil_length,2) + left_boundary_x = 0:(boundary_stencil_length-1) + right_boundary_x = reverse(-boundary_stencil_length+1:0) + + boundary_point_count = div(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[0:div(stencil_length,2)] + R_boundary_deriv_spots = right_boundary_x[0:div(stencil_length,2)] + + 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, left_boundary_x)) for x0 in L_boundary_deriv_spots] + 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))) + + offside = 0 + + coefficients = fill!(Vector{T}(undef,len),0) + + compute_coeffs!(coeff_func, coefficients) + + + + 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, + stencil_coefs, + boundary_stencil_length, + boundary_point_count, + low_boundary_coefs, + high_boundary_coefs,offside,coefficients,coeff_func + ) +end + struct UpwindDifference{N} end """ @@ -330,6 +335,50 @@ function UpwindDifference{N}(derivative_order::Int, ) end +""" +A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the upwind scheme. +""" +function CompleteUpwindDifference{N}(derivative_order::Int, + approximation_order::Int, dx::T, + len::Int, coeff_func=1; offside::Int=0) where {T<:Real,N} + + @assert offside > -1 "Number of offside points should be non-negative" + @assert offside <= div(derivative_order + approximation_order - 1,2) "Number of offside points should not exceed the primary wind points" + + stencil_length = derivative_order + approximation_order + boundary_stencil_length = derivative_order + approximation_order + boundary_point_count = boundary_stencil_length - 1 - offside + + # TODO: Clean up the implementation here so that it is more readable and easier to extend in the future + dummy_x = (0.0 - offside) : stencil_length - 1.0 - offside + stencil_coefs = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, 0.0, dummy_x)) + + low_boundary_x = 0.0:(boundary_stencil_length-1) + L_boundary_deriv_spots = 0.0:boundary_stencil_length - 2.0 - offside + _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, low_boundary_x)) for x0 in L_boundary_deriv_spots] + low_boundary_coefs = convert(SVector{boundary_point_count},_low_boundary_coefs) + + high_boundary_x = 0.0:-1.0:-(boundary_stencil_length-1.0) + R_boundary_deriv_spots = 0.0:-1.0:-(boundary_stencil_length-2.0) + _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, high_boundary_x)) for x0 in R_boundary_deriv_spots] + high_boundary_coefs = convert(SVector{boundary_point_count + offside},_high_boundary_coefs) + + coefficients = zeros(T,len) + compute_coeffs!(coeff_func, coefficients) + + DerivativeOperator{T,N,true,T,typeof(stencil_coefs), + typeof(low_boundary_coefs),typeof(high_boundary_coefs),Vector{T}, + typeof(coeff_func)}( + derivative_order, approximation_order, dx, len, stencil_length, + stencil_coefs, + boundary_stencil_length, + boundary_point_count, + low_boundary_coefs, + high_boundary_coefs,offside,coefficients,coeff_func + ) +end + + CenteredDifference(args...) = CenteredDifference{1}(args...) UpwindDifference(args...;kwargs...) = UpwindDifference{1}(args...;kwargs...) nonlinear_diffusion(args...) = nonlinear_diffusion{1}(args...) From 5a5792ce375bc95a71630c69182482fed1711ffc Mon Sep 17 00:00:00 2001 From: xtalax Date: Wed, 5 Jan 2022 14:36:37 +0000 Subject: [PATCH 03/14] a bracket --- src/derivative_operators/derivative_operator.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 966d94c38..18e2b109f 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -184,7 +184,7 @@ function CompleteCenteredDifference{N}(derivative_order::Int, 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, left_boundary_x)) 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_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))) From 3ded861c6ccf26f2b2895d746957725f97730e84 Mon Sep 17 00:00:00 2001 From: xtalax Date: Wed, 5 Jan 2022 14:40:10 +0000 Subject: [PATCH 04/14] Removed coeff computation for complete --- src/derivative_operators/derivative_operator.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 18e2b109f..b4cb6955e 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -193,9 +193,6 @@ function CompleteCenteredDifference{N}(derivative_order::Int, coefficients = fill!(Vector{T}(undef,len),0) - compute_coeffs!(coeff_func, coefficients) - - DerivativeOperator{T,N,false,T,typeof(stencil_coefs), typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients), @@ -363,8 +360,8 @@ function CompleteUpwindDifference{N}(derivative_order::Int, _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, high_boundary_x)) for x0 in R_boundary_deriv_spots] high_boundary_coefs = convert(SVector{boundary_point_count + offside},_high_boundary_coefs) - coefficients = zeros(T,len) - compute_coeffs!(coeff_func, coefficients) + coefficients = fill!(Vector{T}(undef,len),0) + DerivativeOperator{T,N,true,T,typeof(stencil_coefs), typeof(low_boundary_coefs),typeof(high_boundary_coefs),Vector{T}, From c16e2a5f9074a4f9a514cf824af05bd82039e0f4 Mon Sep 17 00:00:00 2001 From: xtalax Date: Wed, 5 Jan 2022 14:49:07 +0000 Subject: [PATCH 05/14] Exports --- src/DiffEqOperators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index 61bd400d0..40bb384d8 100644 --- a/src/DiffEqOperators.jl +++ b/src/DiffEqOperators.jl @@ -76,7 +76,7 @@ end export MatrixFreeOperator export AnalyticalJacVecOperator, JacVecOperator, getops export AbstractDerivativeOperator, DerivativeOperator, - CenteredDifference, UpwindDifference, nonlinear_diffusion, nonlinear_diffusion!, + CenteredDifference, UpwindDifference, CompleteCenteredDifference, CompleteUpwindDifference, nonlinear_diffusion, nonlinear_diffusion!, GradientOperator, Gradient, CurlOperator, Curl, DivergenceOperator, Divergence export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC, MultiDimDirectionalBC, ComposedMultiDimBC From e9cf75d6009940334214c153adfdeea7a961e182 Mon Sep 17 00:00:00 2001 From: xtalax Date: Wed, 5 Jan 2022 17:17:32 +0000 Subject: [PATCH 06/14] Structs --- src/derivative_operators/derivative_operator.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index b4cb6955e..01e01ad16 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -163,6 +163,8 @@ function CenteredDifference{N}(derivative_order::Int, ) end +struct CompleteCenteredDifference{N} end + """ A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme. """ @@ -332,6 +334,8 @@ function UpwindDifference{N}(derivative_order::Int, ) end +struct CompleteUpwindDifference{N} end + """ A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the upwind scheme. """ From d1a2fcebaafa842f46abad00fd79585479c957cd Mon Sep 17 00:00:00 2001 From: xtalax Date: Wed, 5 Jan 2022 19:56:38 +0000 Subject: [PATCH 07/14] Added forward and reverse half difference for nonlinear laplacian --- .../derivative_operator.jl | 85 +++++++++++++++++++ 1 file changed, 85 insertions(+) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 01e01ad16..ea91d0272 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -208,6 +208,91 @@ function CompleteCenteredDifference{N}(derivative_order::Int, ) end +struct CompleteCenteredDifference{N} end + +""" +A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme. +""" +function FoorwardHalfDifference(derivative_order::Int, + approximation_order::Int, dx::T, + len::Int, coeff_func=1) 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 + dummy_x = -div(stencil_length,2) : div(stencil_length,2) + left_boundary_x = 0:(boundary_stencil_length-1) + right_boundary_x = reverse(-boundary_stencil_length+1:0) + + boundary_point_count = div(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[0:div(stencil_length,2)] + R_boundary_deriv_spots = right_boundary_x[0:div(stencil_length,2)] + + stencil_coefs = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, convert(T, 0.5), 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 + convert(T, 0.5), left_boundary_x)) for x0 in L_boundary_deriv_spots] + 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))) + + offside = 0 + + coefficients = fill!(Vector{T}(undef,len),0) + + + DerivativeOperator{T,1,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, + stencil_coefs, + boundary_stencil_length, + boundary_point_count, + low_boundary_coefs, + high_boundary_coefs,offside,coefficients,coeff_func + ) +end + +function ReverseHalfDifference(derivative_order::Int, + approximation_order::Int, dx::T, + len::Int, coeff_func=1) 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 + dummy_x = -div(stencil_length,2) : div(stencil_length,2) + left_boundary_x = 0:(boundary_stencil_length-1) + right_boundary_x = reverse(-boundary_stencil_length+1:0) + + boundary_point_count = div(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[0:div(stencil_length,2)] + R_boundary_deriv_spots = right_boundary_x[0:div(stencil_length,2)] + + stencil_coefs = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, convert(T, -0.5), 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 + convert(T, -0.5), left_boundary_x)) for x0 in L_boundary_deriv_spots] + 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))) + + offside = 0 + + coefficients = fill!(Vector{T}(undef,len),0) + + + DerivativeOperator{T,1,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, + stencil_coefs, + boundary_stencil_length, + boundary_point_count, + low_boundary_coefs, + high_boundary_coefs,offside,coefficients,coeff_func + ) +end + struct UpwindDifference{N} end """ From aaa06456b207bc20707081df7fcf5a358e8490c6 Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 6 Jan 2022 18:48:15 +0000 Subject: [PATCH 08/14] Fixed ccd and chalfcd --- .../derivative_operator.jl | 74 +++++-------------- 1 file changed, 19 insertions(+), 55 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index ea91d0272..4f9dea949 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/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}, @@ -181,8 +181,8 @@ function CompleteCenteredDifference{N}(derivative_order::Int, boundary_point_count = div(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[0:div(stencil_length,2)] - R_boundary_deriv_spots = right_boundary_x[0:div(stencil_length,2)] + L_boundary_deriv_spots = left_boundary_x[1:div(stencil_length,2)] + R_boundary_deriv_spots = right_boundary_x[1:div(stencil_length,2)] 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, left_boundary_x)) for x0 in L_boundary_deriv_spots] @@ -208,69 +208,33 @@ function CompleteCenteredDifference{N}(derivative_order::Int, ) end -struct CompleteCenteredDifference{N} end """ -A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme. +A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the half centered scheme. See table 2 in https://web.njit.edu/~jiang/math712/fornberg.pdf """ -function FoorwardHalfDifference(derivative_order::Int, +function CompleteHalfCenteredDifference(derivative_order::Int, approximation_order::Int, dx::T, - len::Int, coeff_func=1) where {T<:Real,N} + len::Int, isforward=false) 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 + stencil_length = approximation_order + 2*floor(derivative_order/2) + 2*(approximation_order%2) boundary_stencil_length = derivative_order + approximation_order - dummy_x = -div(stencil_length,2) : div(stencil_length,2) - left_boundary_x = 0:(boundary_stencil_length-1) - right_boundary_x = reverse(-boundary_stencil_length+1:0) + endpoint = (stencil_length-1)/2 + dummy_x = -endpoint : endpoint + left_boundary_x = -1:(boundary_stencil_length-1) + right_boundary_x = reverse(-boundary_stencil_length+1:1) - boundary_point_count = div(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[0:div(stencil_length,2)] - R_boundary_deriv_spots = right_boundary_x[0:div(stencil_length,2)] + # ? Is fornberg valid when taking an x0 outside of the stencil i.e at the boundary? + xoffset = isforward ? convert(T, 0.5) : convert(T, -0.5) - stencil_coefs = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, convert(T, 0.5), 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 + convert(T, 0.5), left_boundary_x)) for x0 in L_boundary_deriv_spots] - 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))) - - offside = 0 - - coefficients = fill!(Vector{T}(undef,len),0) - - - DerivativeOperator{T,1,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, - stencil_coefs, - boundary_stencil_length, - boundary_point_count, - low_boundary_coefs, - high_boundary_coefs,offside,coefficients,coeff_func - ) -end - -function ReverseHalfDifference(derivative_order::Int, - approximation_order::Int, dx::T, - len::Int, coeff_func=1) 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 - dummy_x = -div(stencil_length,2) : div(stencil_length,2) - left_boundary_x = 0:(boundary_stencil_length-1) - right_boundary_x = reverse(-boundary_stencil_length+1:0) - - boundary_point_count = div(stencil_length,2)# -1 due to the ghost point + boundary_point_count = div(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[0:div(stencil_length,2)] - R_boundary_deriv_spots = right_boundary_x[0:div(stencil_length,2)] + L_boundary_deriv_spots = left_boundary_x[1:div(stencil_length,2)] + R_boundary_deriv_spots = right_boundary_x[1:div(stencil_length,2)] - stencil_coefs = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, convert(T, -0.5), 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 + convert(T, -0.5), left_boundary_x)) for x0 in L_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+xoffset, left_boundary_x)) for x0 in L_boundary_deriv_spots] 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] @@ -281,7 +245,7 @@ function ReverseHalfDifference(derivative_order::Int, coefficients = fill!(Vector{T}(undef,len),0) - DerivativeOperator{T,1,false,T,typeof(stencil_coefs), + 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, From 53ffe772ef3b82895cec0244313b4b4d0c5a6142 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 7 Jan 2022 22:18:40 +0000 Subject: [PATCH 09/14] Finished up half offset weight generation --- src/derivative_operators/derivative_operator.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 4f9dea949..73e3fafc2 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -182,7 +182,6 @@ function CompleteCenteredDifference{N}(derivative_order::Int, # 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)] - R_boundary_deriv_spots = right_boundary_x[1:div(stencil_length,2)] 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, left_boundary_x)) for x0 in L_boundary_deriv_spots] @@ -217,28 +216,30 @@ function CompleteHalfCenteredDifference(derivative_order::Int, len::Int, isforward=false) 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 boundary_stencil_length = derivative_order + approximation_order endpoint = (stencil_length-1)/2 dummy_x = -endpoint : endpoint - left_boundary_x = -1:(boundary_stencil_length-1) - right_boundary_x = reverse(-boundary_stencil_length+1:1) + left_boundary_x = 0:(boundary_stencil_length-1) + right_boundary_x = reverse(-boundary_stencil_length+1:0) # ? Is fornberg valid when taking an x0 outside of the stencil i.e at the boundary? - xoffset = isforward ? convert(T, 0.5) : convert(T, -0.5) + xoffset = convert(T, 0.5).*(-div(centered_stencil_length,2) : div(centered_stencil_length,2)) boundary_point_count = div(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)] - R_boundary_deriv_spots = right_boundary_x[1:div(stencil_length,2)] 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+xoffset, left_boundary_x)) for x0 in L_boundary_deriv_spots] + # 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 = [convert(SVector{centered_stencil_length}, SVector{boundary_stencil_length, T}[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)) # _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))) + 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 From 3690385673da63a0e9986de6dbb9f2c182861b0f Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 13 Jan 2022 16:35:16 +0000 Subject: [PATCH 10/14] . --- src/derivative_operators/derivative_operator.jl | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 73e3fafc2..fdcb6fe05 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -163,12 +163,10 @@ function CenteredDifference{N}(derivative_order::Int, ) end -struct CompleteCenteredDifference{N} end - """ A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme. """ -function CompleteCenteredDifference{N}(derivative_order::Int, +function CompleteCenteredDifference(derivative_order::Int, approximation_order::Int, dx::T, len::Int, coeff_func=1) where {T<:Real,N} @assert approximation_order>1 "approximation_order must be greater than 1." @@ -209,11 +207,10 @@ end """ -A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the half centered scheme. See table 2 in https://web.njit.edu/~jiang/math712/fornberg.pdf +A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the half offset centered scheme. See table 2 in https://web.njit.edu/~jiang/math712/fornberg.pdf """ function CompleteHalfCenteredDifference(derivative_order::Int, - approximation_order::Int, dx::T, - len::Int, isforward=false) 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 = approximation_order + 2*floor(derivative_order/2) + 2*(approximation_order%2) centered_stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2 @@ -234,7 +231,7 @@ function CompleteHalfCenteredDifference(derivative_order::Int, stencil_coefs = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, zero(T), 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 = [convert(SVector{centered_stencil_length}, SVector{boundary_stencil_length, T}[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 = [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)) # _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] @@ -384,14 +381,12 @@ function UpwindDifference{N}(derivative_order::Int, ) end -struct CompleteUpwindDifference{N} end - """ A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the upwind scheme. """ -function CompleteUpwindDifference{N}(derivative_order::Int, +function CompleteUpwindDifference(derivative_order::Int, approximation_order::Int, dx::T, - len::Int, coeff_func=1; offside::Int=0) where {T<:Real,N} + offside::Int=0) where {T<:Real,N} @assert offside > -1 "Number of offside points should be non-negative" @assert offside <= div(derivative_order + approximation_order - 1,2) "Number of offside points should not exceed the primary wind points" From e833891d7a9fb615aa5a22af53160c8c32c67393 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 21 Jan 2022 19:05:13 +0000 Subject: [PATCH 11/14] Added half offset --- src/DiffEqOperators.jl | 2 +- .../derivative_operator.jl | 60 +++++++++---------- 2 files changed, 31 insertions(+), 31 deletions(-) diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index ebf1c0720..6cad3788b 100644 --- a/src/DiffEqOperators.jl +++ b/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 diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index fdcb6fe05..0f2323685 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/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 From bfed2d8cb83d3727c7c08da41195056aba72a141 Mon Sep 17 00:00:00 2001 From: xtalax Date: Wed, 26 Jan 2022 16:07:27 +0000 Subject: [PATCH 12/14] fixed fornberg for order 0 and weird dispatch --- src/composite_operators.jl | 49 ++++++++++++++++------------ src/derivative_operators/fornberg.jl | 4 ++- 2 files changed, 31 insertions(+), 22 deletions(-) diff --git a/src/composite_operators.jl b/src/composite_operators.jl index ce5560cf7..2ae1b2840 100644 --- a/src/composite_operators.jl +++ b/src/composite_operators.jl @@ -19,19 +19,19 @@ struct DiffEqOperatorCombination{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{ new{T,typeof(ops),typeof(cache)}(ops, cache) end end -+(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorCombination(ops) -+(op::AbstractDiffEqLinearOperator, A::AbstractMatrix) = op + DiffEqArrayOperator(A) -+(op::AbstractDiffEqLinearOperator{T}, α::UniformScaling) where T = op + UniformScaling(T(α.λ))(size(op,1)) -+(A::AbstractMatrix, op::AbstractDiffEqLinearOperator) = op + A -+(α::UniformScaling, op::AbstractDiffEqLinearOperator) = op + α -+(L1::DiffEqOperatorCombination, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorCombination((L1.ops..., L2)) -+(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1, L2.ops...)) -+(L1::DiffEqOperatorCombination, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1.ops..., L2.ops...)) --(L1::AbstractDiffEqLinearOperator, L2::AbstractDiffEqLinearOperator) = L1 + (-L2) --(L::AbstractDiffEqLinearOperator, A::AbstractMatrix) = L + (-A) --(A::AbstractMatrix, L::AbstractDiffEqLinearOperator) = A + (-L) --(L::AbstractDiffEqLinearOperator, α::UniformScaling) = L + (-α) --(α::UniformScaling, L::AbstractDiffEqLinearOperator) = α + (-L) +Base.:+(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorCombination(ops) +Base.:+(op::AbstractDiffEqLinearOperator, A::AbstractMatrix) = op + DiffEqArrayOperator(A) +Base.:+(op::AbstractDiffEqLinearOperator{T}, α::UniformScaling) where T = op + UniformScaling(T(α.λ))(size(op,1)) +Base.:+(A::AbstractMatrix, op::AbstractDiffEqLinearOperator) = op + A +Base.:+(α::UniformScaling, op::AbstractDiffEqLinearOperator) = op + α +Base.:+(L1::DiffEqOperatorCombination, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorCombination((L1.ops..., L2)) +Base.:+(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1, L2.ops...)) +Base.:+(L1::DiffEqOperatorCombination, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1.ops..., L2.ops...)) +Base.:-(L1::AbstractDiffEqLinearOperator, L2::AbstractDiffEqLinearOperator) = L1 + (-L2) +Base.:-(L::AbstractDiffEqLinearOperator, A::AbstractMatrix) = L + (-A) +Base.:-(A::AbstractMatrix, L::AbstractDiffEqLinearOperator) = A + (-L) +Base.:-(L::AbstractDiffEqLinearOperator, α::UniformScaling) = L + (-α) +Base.:-(α::UniformScaling, L::AbstractDiffEqLinearOperator) = α + (-L) getops(L::DiffEqOperatorCombination) = L.ops Matrix(L::DiffEqOperatorCombination) = sum(Matrix, L.ops) convert(::Type{AbstractMatrix}, L::DiffEqOperatorCombination) = @@ -40,8 +40,8 @@ convert(::Type{AbstractMatrix}, L::DiffEqOperatorCombination) = size(L::DiffEqOperatorCombination, args...) = size(L.ops[1], args...) getindex(L::DiffEqOperatorCombination, i::Int) = sum(op -> op[i], L.ops) getindex(L::DiffEqOperatorCombination, I::Vararg{Int, N}) where {N} = sum(op -> op[I...], L.ops) -*(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op * x, L.ops) -*(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops) +Base.:*(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op * x, L.ops) +Base.:*(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops) /(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op / x, L.ops) \(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops) function mul!(y::AbstractVector, L::DiffEqOperatorCombination, b::AbstractVector) @@ -77,13 +77,20 @@ struct DiffEqOperatorComposition{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{ new{T,typeof(ops),typeof(caches)}(ops, caches) end end -*(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorComposition(reverse(ops)) +# this is needed to not break dispatch in MethodOfLines +function Base.:*(ops::AbstractDiffEqLinearOperator...) + try + return DiffEqOperatorComposition(reverse(ops)) + catch e + return 1 + end +end ∘(L1::AbstractDiffEqLinearOperator, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1)) -*(L1::DiffEqOperatorComposition, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1.ops...)) +Base.:*(L1::DiffEqOperatorComposition, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1.ops...)) ∘(L1::DiffEqOperatorComposition, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1.ops...)) -*(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1)) +Base.:*(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1)) ∘(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1)) -*(L1::DiffEqOperatorComposition, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1.ops...)) +Base.:*(L1::DiffEqOperatorComposition, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1.ops...)) ∘(L1::DiffEqOperatorComposition, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1.ops...)) getops(L::DiffEqOperatorComposition) = L.ops Matrix(L::DiffEqOperatorComposition) = prod(Matrix, reverse(L.ops)) @@ -93,8 +100,8 @@ convert(::Type{AbstractMatrix}, L::DiffEqOperatorComposition) = size(L::DiffEqOperatorComposition) = (size(L.ops[end], 1), size(L.ops[1], 2)) size(L::DiffEqOperatorComposition, m::Integer) = size(L)[m] opnorm(L::DiffEqOperatorComposition) = prod(opnorm, L.ops) -*(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op*acc, L.ops; init=x) -*(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x) +Base.:*(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op*acc, L.ops; init=x) +Base.:*(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x) /(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op/acc, L.ops; init=x) /(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc/op, L.ops; init=x) \(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op\acc, reverse(L.ops); init=x) diff --git a/src/derivative_operators/fornberg.jl b/src/derivative_operators/fornberg.jl index 0b24ae706..1406d8bbf 100644 --- a/src/derivative_operators/fornberg.jl +++ b/src/derivative_operators/fornberg.jl @@ -55,6 +55,8 @@ function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507 - Modified Fornberg Algorithm =# _C = C[:,end] - _C[div(N,2)+1] -= sum(_C) + if order != 0 + _C[div(N,2)+1] -= sum(_C) + end return _C end From e5a3da28f9010b0540adb01c486db6f0d8ccd8eb Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 27 Jan 2022 15:09:45 +0000 Subject: [PATCH 13/14] removed unnessecary boundary point --- src/derivative_operators/derivative_operator.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 0f2323685..9dd7a9a1a 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -221,7 +221,7 @@ function CompleteHalfCenteredDifference(derivative_order::Int, 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 = 1.5:boundary_stencil_length+0.5 + xoffset = 1.5:boundary_stencil_length-0.5 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. From 473d512639ea6f988ea3160593204ebaf7f135b1 Mon Sep 17 00:00:00 2001 From: xtalax Date: Sat, 29 Jan 2022 18:20:45 +0000 Subject: [PATCH 14/14] tests passing --- .../derivative_operator.jl | 4 +-- .../derivative_operators_interface.jl | 28 +++++++++++++++++++ 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 9dd7a9a1a..83f308b1e 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -220,10 +220,10 @@ function CompleteHalfCenteredDifference(derivative_order::Int, left_boundary_x = 1:(boundary_stencil_length) right_boundary_x = reverse(-boundary_stencil_length:-1) + boundary_point_count = div(centered_stencil_length,2) # -1 due to the ghost point # ? Is fornberg valid when taking an x0 outside of the stencil i.e at the boundary? - xoffset = 1.5:boundary_stencil_length-0.5 + xoffset = range(1.5, length=boundary_point_count, step = 1.0) - 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 = xoffset[1:div(centered_stencil_length,2)] diff --git a/test/DerivativeOperators/derivative_operators_interface.jl b/test/DerivativeOperators/derivative_operators_interface.jl index 63df393fa..952ea9914 100644 --- a/test/DerivativeOperators/derivative_operators_interface.jl +++ b/test/DerivativeOperators/derivative_operators_interface.jl @@ -144,6 +144,34 @@ end @test BandedMatrix(L) ≈ correct end +@testset "Correctness of Uniform Stencils, Complete" begin + weights = [] + + push!(weights, ([-0.5,0,0.5], [1.,-2.,1.], [-1/2,1.,0.,-1.,1/2])) + push!(weights, ([1/12, -2/3,0,2/3,-1/12], [-1/12,4/3,-5/2,4/3,-1/12], [1/8,-1.,13/8,0.,-13/8,1.,-1/8])) + + for d in 1:3 + for (i,a) in enumerate([2,4]) + D = CompleteCenteredDifference(d,a,1.0) + + @test all(isapprox.(D.stencil_coefs, weights[i][d], atol=1e-10)) + end + end +end + +@testset "Correctness of Uniform Stencils, Complete Half" begin + weights = (([.5, .5], [-1/16, 9/16, 9/16, -1/16]), + ([-1., 1.], [1/24, -9/8, 9/8, -1/24])) + for (i,a) in enumerate([2,4]) + for d in 0:1 + D = CompleteHalfCenteredDifference(d,a,1.0) + @test all(D.stencil_coefs .≈ weights[d+1][i]) + end + end + +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)