Skip to content

Commit

Permalink
Merge pull request #112 from ranocha/fix_111
Browse files Browse the repository at this point in the history 
fix fourth order variable coefficient operator of Mattsson2012
  • Loading branch information
@ranocha
ranocha committed Aug 4, 2021
2 parents f6e90a6 + 97160b1 commit d372a2d
Show file tree
Hide file tree
Showing 2 changed files with 93 additions and 76 deletions.
56 changes: 28 additions & 28 deletions src/SBP_coefficients/Mattsson2012.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -516,30 +516,30 @@ function convolve_boundary_coefficients!(dest::AbstractVector, cache::Mattsson20
+ (d363*b[3] + d364*b[4] + d365*b[5]) * u[6]
)
dest[ 4] = α * (
(d411*b[1] + d413*b[3] + d414) * u[1]
+ (d421*b[1] + d423*b[3] + d424) * u[2]
(d411*b[1] + d413*b[3] + d414*b[4]) * u[1]
+ (d421*b[1] + d423*b[3] + d424*b[4]) * u[2]
+ (d431*b[1] + d433*b[3] + d434*b[4] + d435*b[5]) * u[3]
+ (d441*b[1] + d443*b[3] + d444*b[4] + d445*b[5] + d446*b[6]) * u[4]
+ (d453*b[3] + d454*b[4] + d455*b[5] + d456*b[6]) * u[5]
+ (d463*b[3] + d464*b[4] + d465*b[5] + d466*b[6]) * u[6]
)
dest[ 5] = α * (
(d513*b[3] + d514) * u[1]
+ (d523*b[3] + d524) * u[2]
(d513*b[3] + d514*b[4]) * u[1]
+ (d523*b[3] + d524*b[4]) * u[2]
+ (d533*b[3] + d534*b[4] + d535*b[5]) * u[3]
+ (d543*b[3] + d544*b[4] + d545*b[5] + d546*b[6]) * u[4]
+ (d553*b[3] + d554*b[4] + d555*b[5] + d556*b[6] + d557*b[7]) * u[5]
+ (d563*b[3] + d564*b[4] + d565*b[5] + d566*b[6] + d567*b[7]) * u[6]
+ (d575*b[5] + d576*b[6] + d577*b[7]) * u[7]
)
dest[ 6] = α * (
(d613*b[3] + d614) * u[1]
+ (d623*b[3] + d624) * u[2]
(d613*b[3] + d614*b[4]) * u[1]
+ (d623*b[3] + d624*b[4]) * u[2]
+ (d633*b[3] + d634*b[4] + d635*b[5]) * u[3]
+ (d643*b[3] + d644*b[4] + d645*b[5] + d646*b[6]) * u[4]
+ (d653*b[3] + d654*b[4] + d655*b[5] + d656*b[6] + d657*b[7]) * u[5]
+ (d663*b[3] + d664*b[4] + d665*b[5] + d666*b[6] + d667*b[7] + d668*b[8]) * u[6]
+ (d675*b[5] + d676*b[6] + d677*b[7] + d678) * u[7]
+ (d675*b[5] + d676*b[6] + d677*b[7] + d678*b[8]) * u[7]
+ (d686*b[6] + d687*b[7] + d688*b[8]) * u[8]
)

Expand Down Expand Up @@ -569,30 +569,30 @@ function convolve_boundary_coefficients!(dest::AbstractVector, cache::Mattsson20
+ (d363*b[end-2] + d364*b[end-3] + d365*b[end-4]) * u[end-5]
)
dest[end-3] = α * (
(d411*b[end] + d413*b[end-2] + d414) * u[end]
+ (d421*b[end] + d423*b[end-2] + d424) * u[end-1]
(d411*b[end] + d413*b[end-2] + d414*b[end-3]) * u[end]
+ (d421*b[end] + d423*b[end-2] + d424*b[end-3]) * u[end-1]
+ (d431*b[end] + d433*b[end-2] + d434*b[end-3] + d435*b[end-4]) * u[end-2]
+ (d441*b[end] + d443*b[end-2] + d444*b[end-3] + d445*b[end-4] + d446*b[end-5]) * u[end-3]
+ (d453*b[end-2] + d454*b[end-3] + d455*b[end-4] + d456*b[end-5]) * u[end-4]
+ (d463*b[end-2] + d464*b[end-3] + d465*b[end-4] + d466*b[end-5]) * u[end-5]
)
dest[end-4] = α * (
(d513*b[end-2] + d514) * u[end]
+ (d523*b[end-2] + d524) * u[end-1]
(d513*b[end-2] + d514*b[end-3]) * u[end]
+ (d523*b[end-2] + d524*b[end-3]) * u[end-1]
+ (d533*b[end-2] + d534*b[end-3] + d535*b[end-4]) * u[end-2]
+ (d543*b[end-2] + d544*b[end-3] + d545*b[end-4] + d546*b[end-5]) * u[end-3]
+ (d553*b[end-2] + d554*b[end-3] + d555*b[end-4] + d556*b[end-5] + d557*b[end-6]) * u[end-4]
+ (d563*b[end-2] + d564*b[end-3] + d565*b[end-4] + d566*b[end-5] + d567*b[end-6]) * u[end-5]
+ (d575*b[end-4] + d576*b[end-5] + d577*b[end-6]) * u[end-6]
)
dest[end-5] = α * (
(d613*b[end-2] + d614) * u[end]
+ (d623*b[end-2] + d624) * u[end-1]
(d613*b[end-2] + d614*b[end-3]) * u[end]
+ (d623*b[end-2] + d624*b[end-3]) * u[end-1]
+ (d633*b[end-2] + d634*b[end-3] + d635*b[end-4]) * u[end-2]
+ (d643*b[end-2] + d644*b[end-3] + d645*b[end-4] + d646*b[end-5]) * u[end-3]
+ (d653*b[end-2] + d654*b[end-3] + d655*b[end-4] + d656*b[end-5] + d657*b[end-6]) * u[end-4]
+ (d663*b[end-2] + d664*b[end-3] + d665*b[end-4] + d666*b[end-5] + d667*b[end-6] + d668*b[end-7]) * u[end-5]
+ (d675*b[end-4] + d676*b[end-5] + d677*b[end-6] + d678) * u[end-6]
+ (d675*b[end-4] + d676*b[end-5] + d677*b[end-6] + d678*b[end-7]) * u[end-6]
+ (d686*b[end-5] + d687*b[end-6] + d688*b[end-7]) * u[end-7]
)
end
Expand Down Expand Up @@ -650,30 +650,30 @@ function convolve_boundary_coefficients!(dest::AbstractVector, cache::Mattsson20
+ (d363*b[3] + d364*b[4] + d365*b[5]) * u[6]
) + β*dest[3]
dest[ 4] = α * (
(d411*b[1] + d413*b[3] + d414) * u[1]
+ (d421*b[1] + d423*b[3] + d424) * u[2]
(d411*b[1] + d413*b[3] + d414*b[4]) * u[1]
+ (d421*b[1] + d423*b[3] + d424*b[4]) * u[2]
+ (d431*b[1] + d433*b[3] + d434*b[4] + d435*b[5]) * u[3]
+ (d441*b[1] + d443*b[3] + d444*b[4] + d445*b[5] + d446*b[6]) * u[4]
+ (d453*b[3] + d454*b[4] + d455*b[5] + d456*b[6]) * u[5]
+ (d463*b[3] + d464*b[4] + d465*b[5] + d466*b[6]) * u[6]
) + β*dest[4]
dest[ 5] = α * (
(d513*b[3] + d514) * u[1]
+ (d523*b[3] + d524) * u[2]
(d513*b[3] + d514*b[4]) * u[1]
+ (d523*b[3] + d524*b[4]) * u[2]
+ (d533*b[3] + d534*b[4] + d535*b[5]) * u[3]
+ (d543*b[3] + d544*b[4] + d545*b[5] + d546*b[6]) * u[4]
+ (d553*b[3] + d554*b[4] + d555*b[5] + d556*b[6] + d557*b[7]) * u[5]
+ (d563*b[3] + d564*b[4] + d565*b[5] + d566*b[6] + d567*b[7]) * u[6]
+ (d575*b[5] + d576*b[6] + d577*b[7]) * u[7]
) + β*dest[5]
dest[ 6] = α * (
(d613*b[3] + d614) * u[1]
+ (d623*b[3] + d624) * u[2]
(d613*b[3] + d614*b[4]) * u[1]
+ (d623*b[3] + d624*b[4]) * u[2]
+ (d633*b[3] + d634*b[4] + d635*b[5]) * u[3]
+ (d643*b[3] + d644*b[4] + d645*b[5] + d646*b[6]) * u[4]
+ (d653*b[3] + d654*b[4] + d655*b[5] + d656*b[6] + d657*b[7]) * u[5]
+ (d663*b[3] + d664*b[4] + d665*b[5] + d666*b[6] + d667*b[7] + d668*b[8]) * u[6]
+ (d675*b[5] + d676*b[6] + d677*b[7] + d678) * u[7]
+ (d675*b[5] + d676*b[6] + d677*b[7] + d678*b[8]) * u[7]
+ (d686*b[6] + d687*b[7] + d688*b[8]) * u[8]
) + β*dest[6]

Expand Down Expand Up @@ -703,30 +703,30 @@ function convolve_boundary_coefficients!(dest::AbstractVector, cache::Mattsson20
+ (d363*b[end-2] + d364*b[end-3] + d365*b[end-4]) * u[end-5]
) + β*dest[end-2]
dest[end-3] = α * (
(d411*b[end] + d413*b[end-2] + d414) * u[end]
+ (d421*b[end] + d423*b[end-2] + d424) * u[end-1]
(d411*b[end] + d413*b[end-2] + d414*b[end-3]) * u[end]
+ (d421*b[end] + d423*b[end-2] + d424*b[end-3]) * u[end-1]
+ (d431*b[end] + d433*b[end-2] + d434*b[end-3] + d435*b[end-4]) * u[end-2]
+ (d441*b[end] + d443*b[end-2] + d444*b[end-3] + d445*b[end-4] + d446*b[end-5]) * u[end-3]
+ (d453*b[end-2] + d454*b[end-3] + d455*b[end-4] + d456*b[end-5]) * u[end-4]
+ (d463*b[end-2] + d464*b[end-3] + d465*b[end-4] + d466*b[end-5]) * u[end-5]
) + β*dest[end-3]
dest[end-4] = α * (
(d513*b[end-2] + d514) * u[end]
+ (d523*b[end-2] + d524) * u[end-1]
(d513*b[end-2] + d514*b[end-3]) * u[end]
+ (d523*b[end-2] + d524*b[end-3]) * u[end-1]
+ (d533*b[end-2] + d534*b[end-3] + d535*b[end-4]) * u[end-2]
+ (d543*b[end-2] + d544*b[end-3] + d545*b[end-4] + d546*b[end-5]) * u[end-3]
+ (d553*b[end-2] + d554*b[end-3] + d555*b[end-4] + d556*b[end-5] + d557*b[end-6]) * u[end-4]
+ (d563*b[end-2] + d564*b[end-3] + d565*b[end-4] + d566*b[end-5] + d567*b[end-6]) * u[end-5]
+ (d575*b[end-4] + d576*b[end-5] + d577*b[end-6]) * u[end-6]
) + β*dest[end-4]
dest[end-5] = α * (
(d613*b[end-2] + d614) * u[end]
+ (d623*b[end-2] + d624) * u[end-1]
(d613*b[end-2] + d614*b[end-3]) * u[end]
+ (d623*b[end-2] + d624*b[end-3]) * u[end-1]
+ (d633*b[end-2] + d634*b[end-3] + d635*b[end-4]) * u[end-2]
+ (d643*b[end-2] + d644*b[end-3] + d645*b[end-4] + d646*b[end-5]) * u[end-3]
+ (d653*b[end-2] + d654*b[end-3] + d655*b[end-4] + d656*b[end-5] + d657*b[end-6]) * u[end-4]
+ (d663*b[end-2] + d664*b[end-3] + d665*b[end-4] + d666*b[end-5] + d667*b[end-6] + d668*b[end-7]) * u[end-5]
+ (d675*b[end-4] + d676*b[end-5] + d677*b[end-6] + d678) * u[end-6]
+ (d675*b[end-4] + d676*b[end-5] + d677*b[end-6] + d678*b[end-7]) * u[end-6]
+ (d686*b[end-5] + d687*b[end-6] + d688*b[end-7]) * u[end-7]
) + β*dest[end-5]
end
Expand Down
113 changes: 65 additions & 48 deletions test/var_coef_operators_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,58 +4,75 @@ using SummationByPartsOperators

test_list = (Mattsson2012(),)

# Test consistency with constant coefficient operators
for source in test_list, acc_order in (2,4,6), T in (Float32,Float64)
xmin = -one(T)
xmax = 2*one(T)
N = 101

D2 = derivative_operator(source, 2, acc_order, xmin, xmax, N)
D2var = try
var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one)
catch err
!isa(err, ArgumentError) && throw(err)
nothing
end
D2var === nothing && continue
@testset "Test consistency with constant coefficient operators" begin
for source in test_list, acc_order in (2,4,6), T in (Float32,Float64)
xmin = -one(T)
xmax = 2*one(T)
N = 101

D2 = derivative_operator(source, 2, acc_order, xmin, xmax, N)
D2var = try
var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one)
catch err
!isa(err, ArgumentError) && throw(err)
nothing
end
D2var === nothing && continue

for compact in (true, false)
show(IOContext(devnull, :compact=>false), D2var)
show(IOContext(devnull, :compact=>false), D2var.coefficients)
end

for compact in (true, false)
show(IOContext(devnull, :compact=>false), D2var)
show(IOContext(devnull, :compact=>false), D2var.coefficients)
@test maximum(abs, Matrix(D2) - Matrix(D2var)) < 10000*eps(T)
end
end


@testset "Test consistency for vanishing coefficients" begin
for source in test_list, T in (Float32, Float64), acc_order in (2,4,6)
xmin = zero(T)
xmax = 5*one(T)
N = 51
source = Mattsson2012()

@test maximum(abs, Matrix(D2) - Matrix(D2var)) < 10000*eps(T)
D2var = var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, zero)
x = grid(D2var)
u = rand(T, length(x))
@test maximum(abs, D2var * u) < 10 * eps(T)
end
end


# Compare mul! with β=0 and mul! without β.
for T in (Float32, Float64), acc_order in (2,4,6)
xmin = zero(T)
xmax = 5*one(T)
N = 51
source = Mattsson2012()

D2var_serial = var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one)
D2var_threads= var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one, ThreadedMode())
D2var_safe = var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one, SafeMode())
D2var_full = Matrix(D2var_serial)
D2var_sparse = sparse(D2var_serial)

x = grid(D2var_serial)
u = x.^5
dest1 = fill(zero(eltype(u)), length(u))
dest2 = fill(zero(eltype(u)), length(u))

mul!(dest1, D2var_serial, u, one(T), zero(T))
mul!(dest2, D2var_serial, u, one(T))
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
mul!(dest1, D2var_threads, u, one(T), zero(T))
mul!(dest2, D2var_threads, u, one(T))
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
mul!(dest2, D2var_serial, u)
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
mul!(dest2, D2var_full, u)
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
mul!(dest2, D2var_safe, u)
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
@testset "Compare mul! with β=0 and mul! without β" begin
for T in (Float32, Float64), acc_order in (2,4,6)
xmin = zero(T)
xmax = 5*one(T)
N = 51
source = Mattsson2012()

D2var_serial = var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one)
D2var_threads= var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one, ThreadedMode())
D2var_safe = var_coef_derivative_operator(source, 2, acc_order, xmin, xmax, N, one, SafeMode())
D2var_full = Matrix(D2var_serial)
D2var_sparse = sparse(D2var_serial)

x = grid(D2var_serial)
u = x.^5
dest1 = fill(zero(eltype(u)), length(u))
dest2 = fill(zero(eltype(u)), length(u))

mul!(dest1, D2var_serial, u, one(T), zero(T))
mul!(dest2, D2var_serial, u, one(T))
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
mul!(dest1, D2var_threads, u, one(T), zero(T))
mul!(dest2, D2var_threads, u, one(T))
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
mul!(dest2, D2var_serial, u)
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
mul!(dest2, D2var_full, u)
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
mul!(dest2, D2var_safe, u)
@test all(i->dest1[i] ≈ dest2[i], eachindex(u))
end
end

0 comments on commit d372a2d

Please sign in to comment.