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

Commit 3a57c2c

Browse files
xtalaxxtalax
committed
tests
1 parent 6319af2 commit 3a57c2c

File tree

2 files changed

+235
-192
lines changed

2 files changed

+235
-192
lines changed

‎src/derivative_operators/derivative_operator.jl‎

Lines changed: 40 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -206,42 +206,43 @@ function CompleteCenteredDifference(derivative_order::Int,
206206
)
207207
end
208208

209-
function CompleteCenteredDifference{N}(derivative_order::Int,
210-
approximation_order::Int, dx::AbstractVector{T}, len::Integer) where {T<:Real,N}
211-
stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
209+
function CompleteCenteredDifference(derivative_order::Int,
210+
approximation_order::Int, x::AbstractVector{T}) where {T<:Real,N}
211+
stencil_length = derivative_order + approximation_order - 1 + (derivative_order + approximation_order) % 2
212212
boundary_stencil_length = derivative_order + approximation_order
213-
stencil_x = zeros(T, stencil_length)
214-
boundary_point_count = div(stencil_length,2)
215-
216-
interior_x = boundary_point_count+1:len+2-boundary_point_count
217-
dummy_x = -div(stencil_length,2) : div(stencil_length,2)-1
218-
low_boundary_x = [zero(T); cumsum(dx[1:boundary_stencil_length-1])]
219-
high_boundary_x = cumsum(dx[end-boundary_stencil_length+1:end])
213+
stencil_x = zeros(T, stencil_length)
214+
boundary_point_count = endpoint = div(stencil_length, 2)
215+
len = length(x)
216+
dx = [x[i+1] - x[i] for i in 1:length(x)-1]
217+
interior_x = boundary_point_count+1:len-boundary_point_count
218+
dummy_x = -div(stencil_length, 2):div(stencil_length, 2)-1
219+
low_boundary_x = [zero(T); cumsum(dx[1:boundary_stencil_length-1])]
220+
high_boundary_x = cumsum(dx[end-boundary_stencil_length+1:end])
220221
# 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.
221-
deriv_spots = (-div(stencil_length,2)+1) : -1
222+
deriv_spots = (-div(stencil_length, 2)+1):-1
223+
224+
stencil_coefs = [convert(SVector{stencil_length,T}, calculate_weights(derivative_order, x[i], x[i-endpoint:i+endpoint])) for i in interior_x]
225+
low_boundary_coefs = SVector{boundary_stencil_length,T}[convert(SVector{boundary_stencil_length,T},
226+
calculate_weights(derivative_order, low_boundary_x[i+1], low_boundary_x)) for i in 0:boundary_point_count-1]
222227

223-
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]
224-
_low_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T},
225-
calculate_weights(derivative_order, low_boundary_x[i+1], low_boundary_x)) for i in 0:boundary_point_count-1]
226-
low_boundary_coefs = convert(SVector{boundary_point_count},_low_boundary_coefs)
227-
_high_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T},
228-
calculate_weights(derivative_order, high_boundary_x[end-i], high_boundary_x)) for i in boundary_point_count+1:-1:0]
229-
high_boundary_coefs = convert(SVector{boundary_point_count},_high_boundary_coefs)
228+
229+
high_boundary_coefs = SVector{boundary_stencil_length,T}[convert(SVector{boundary_stencil_length,T},
230+
calculate_weights(derivative_order, high_boundary_x[end-i], high_boundary_x)) for i in boundary_point_count+1:-1:0]
230231

231232
offside = 0
232233
coefficients = nothing
233234

234-
DerivativeOperator{T,N,false,typeof(dx),typeof(stencil_coefs),
235+
DerivativeOperator{eltype(x),Nothing,false,typeof(dx),typeof(stencil_coefs),
235236
typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),
236-
typeof(coeff_func)}(
237+
Nothing}(
237238
derivative_order, approximation_order, dx,
238239
len, stencil_length,
239240
stencil_coefs,
240241
boundary_stencil_length,
241242
boundary_point_count,
242243
low_boundary_coefs,
243-
high_boundary_coefs,offside,coefficients,coeff_func
244-
)
244+
high_boundary_coefs, offside, coefficients, nothing
245+
)
245246
end
246247

247248

@@ -293,13 +294,13 @@ function CompleteHalfCenteredDifference(derivative_order::Int,
293294
end
294295

295296
function CompleteHalfCenteredDifference(derivative_order::Int,
296-
approximation_order::Int, dx::T) where {T<:AbstractVector{<:Real},N}
297+
approximation_order::Int, x::T) where {T<:AbstractVector{<:Real},N}
297298
@assert approximation_order > 1 "approximation_order must be greater than 1."
298299
centered_stencil_length = approximation_order + 2 * Int(floor(derivative_order / 2)) + (approximation_order % 2)
299300
boundary_stencil_length = derivative_order + approximation_order
300301
endpoint = div(centered_stencil_length, 2)
301-
x = [0.0, cumsum(dx)...]
302302
hx = [(x[i] + x[i+1]) / 2 for i in 1:length(x)-1]
303+
dx = [x[i+1] - x[i] for i in 1:length(x)-1]
303304

304305

305306
low_boundary_x = @view(x[1:boundary_stencil_length])
@@ -314,21 +315,20 @@ function CompleteHalfCenteredDifference(derivative_order::Int,
314315

315316
R_boundary_deriv_spots = hx[length(hx):-1:length(hx)-div(centered_stencil_length, 2)+1]
316317

317-
stencil_coefs = [convert(SVector{centered_stencil_length,T}, calculate_weights(derivative_order, hx[i], x[i-endpoint+1:i+endpoint])) for i in endpoint+1:length(hx)-endpoint]
318+
stencil_coefs = [convert(SVector{centered_stencil_length,eltype(x)}, calculate_weights(derivative_order, hx[i], x[i-endpoint+1:i+endpoint])) for i in endpoint+1:length(hx)-endpoint]
318319
# For each boundary point, for each tappoint in the half offset central difference stencil, we need to calculate the coefficients for the stencil.
319320

320321

321-
_low_boundary_coefs = [convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, offset, low_boundary_x)) for offset in L_boundary_deriv_spots]
322-
low_boundary_coefs = convert(SVector{boundary_point_count}, _low_boundary_coefs)
322+
low_boundary_coefs = [convert(SVector{boundary_stencil_length,eltype(x)}, calculate_weights(derivative_order, offset, low_boundary_x)) for offset in L_boundary_deriv_spots]
323+
324+
high_boundary_coefs = [convert(SVector{boundary_stencil_length,eltype(x)}, calculate_weights(derivative_order, offset, high_boundary_x)) for offset in R_boundary_deriv_spots]
323325

324-
_high_boundary_coefs = [convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, offset, high_boundary_x)) for offset in R_boundary_deriv_spots]
325-
high_boundary_coefs = convert(SVector{boundary_point_count}, _high_boundary_coefs)
326326
# _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]
327327

328328
offside = 0
329329
coefficients = nothing
330330

331-
DerivativeOperator{T,Nothing,false,T,typeof(stencil_coefs),
331+
DerivativeOperator{eltype(x),Nothing,false,typeof(dx),typeof(stencil_coefs),
332332
typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),Nothing}(
333333
derivative_order, approximation_order, dx, 1, centered_stencil_length,
334334
stencil_coefs,
@@ -510,9 +510,8 @@ function CompleteUpwindDifference(derivative_order::Int,
510510
end
511511

512512
# TODO implement the non-uniform grid
513-
function CompleteUpwindDifference{N}(derivative_order::Int,
514-
approximation_order::Int, dx::AbstractVector{T},
515-
len::Int, coeff_func=1; offside::Int=0) where {T<:Real,N}
513+
function CompleteUpwindDifference(derivative_order::Int,
514+
approximation_order::Int, x::AbstractVector{T}, offside::Int=0) where {T<:Real,N}
516515

517516
@assert offside > -1 "Number of offside points should be non-negative"
518517
@assert offside <= div(derivative_order + approximation_order - 1, 2) "Number of offside points should not exceed the primary wind points"
@@ -521,11 +520,7 @@ function CompleteUpwindDifference{N}(derivative_order::Int,
521520
boundary_stencil_length = derivative_order + approximation_order
522521
boundary_point_count = boundary_stencil_length - 1 - offside
523522

524-
endpoint = div(centered_stencil_length, 2)
525-
526-
x = [0.0, cumsum(dx)...]
527-
hx = [(x[i] + x[i+1]) / 2 for i in 1:length(x)-1]
528-
523+
dx = [x[i+1] - x[i] for i in 1:length(x)-1]
529524

530525
low_boundary_x = @view(x[1:boundary_stencil_length])
531526
high_boundary_x = @view(x[end-boundary_stencil_length+1:end])
@@ -536,23 +531,23 @@ function CompleteUpwindDifference{N}(derivative_order::Int,
536531

537532
R_boundary_deriv_spots = x[end-boundary_stencil_length+1:end]
538533

539-
stencil_coefs = [convert(SVector{centered_stencil_length,T}, calculate_weights(derivative_order, x[i], x[i-offside:i+stencil_length-1-offside])) for i in boundary_stencil_length-offside:length(x)-boundary_stencil_length]
534+
stencil_coefs = [convert(SVector{stencil_length,eltype(x)}, calculate_weights(derivative_order, x[i], @view(x[i-offside:i+stencil_length-1-offside]))) for i in boundary_stencil_length-offside:length(x)-boundary_stencil_length]
540535
# For each boundary point, for each tappoint in the half offset central difference stencil, we need to calculate the coefficients for the stencil.
541536

542537

543-
_low_boundary_coefs = [convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, offset, low_boundary_x)) for offset in L_boundary_deriv_spots]
544-
low_boundary_coefs = convert(SVector{boundary_point_count}, _low_boundary_coefs)
538+
low_boundary_coefs = [convert(SVector{boundary_stencil_length,eltype(x)}, calculate_weights(derivative_order, offset, low_boundary_x)) for offset in L_boundary_deriv_spots]
539+
540+
541+
high_boundary_coefs = [convert(SVector{boundary_stencil_length,eltype(x)}, calculate_weights(derivative_order, offset, high_boundary_x)) for offset in R_boundary_deriv_spots]
545542

546-
_high_boundary_coefs = [convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, offset, high_boundary_x)) for offset in R_boundary_deriv_spots]
547-
high_boundary_coefs = convert(SVector{boundary_point_count}, _high_boundary_coefs)
548543
# _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]
549544

550545
offside = 0
551546
coefficients = nothing
552547

553-
DerivativeOperator{T,Nothing,false,T,typeof(stencil_coefs),
548+
DerivativeOperator{eltype(x),Nothing,false,typeof(dx),typeof(stencil_coefs),
554549
typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),Nothing}(
555-
derivative_order, approximation_order, dx, 1, centered_stencil_length,
550+
derivative_order, approximation_order, dx, 1, stencil_length,
556551
stencil_coefs,
557552
boundary_stencil_length,
558553
boundary_point_count,

0 commit comments

Comments
 (0)