Skip to content

Commit

Permalink
Ensure second derivatives give exactly periodic results
Browse files Browse the repository at this point in the history 
...without even rounding errors, as these might drift and accumulate
into a non-negligible error after many timesteps.
  • Loading branch information
@johnomotani
johnomotani committed Sep 11, 2024
1 parent e999753 commit ad9af46
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 8 deletions.
21 changes: 20 additions & 1 deletion moment_kinetics/src/calculus.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,15 @@ function second_derivative!(d2f, f, coord, spectral)
error("Distributed memory MPI not yet supported here")
end
d2f[1] -= C * coord.scratch2_2d[end,end]
d2f[end] += C * coord.scratch2_2d[1,1]
# With the first derivative from the opposite end of the grid, d2f[end] here
# should be equal to d2f[1] up to rounding errors...
@boundscheck isapprox(d2f[end] + C * coord.scratch2_2d[1,1], d2f[1]; atol=1.0e-14)
# ...but because arithmetic operations were in a different order, there may be
# rounding errors, so set the two ends exactly equal to ensure consistency for the
# rest of the code - we assume that duplicate versions of the 'same point' on
# element boundaries (due to periodic bc or distributed-MPI block boundaries) are
# exactly equal.
d2f[end] = d2f[1]
else
error("Unsupported bc '$(coord.bc)'")
end
Expand Down Expand Up @@ -220,6 +228,17 @@ function second_derivative!(d2f, f, coord, spectral::weak_discretization_info)
* "distributed coordinate")
end
mass_matrix_solve!(d2f, coord.scratch3, spectral)

if coord.bc == "periodic"
# d2f[end] here should be equal to d2f[1] up to rounding errors...
@boundscheck isapprox(d2f[end], d2f[1]; atol=1.0e-14)
# ...but in the matrix operations arithmetic operations are not necessarily in
# exactly the same order, there may be rounding errors, so set the two ends
# exactly equal to ensure consistency for the rest of the code - we assume that
# duplicate versions of the 'same point' on element boundaries (due to periodic bc
# or distributed-MPI block boundaries) are exactly equal.
d2f[end] = d2f[1]
end
end

function laplacian_derivative!(d2f, f, coord, spectral::weak_discretization_info)
Expand Down
24 changes: 17 additions & 7 deletions moment_kinetics/test/calculus_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ function runtests()
@testset "calculus" verbose=use_verbose begin
println("calculus tests")
@testset "fundamental theorem of calculus" begin
@testset "$discretization $ngrid $nelement" for
@testset "$discretization $ngrid $nelement $cheb_option" for
(discretization, element_spacing_option, etol, cheb_option) ∈ (("finite_difference", "uniform", 1.0e-15, ""), ("chebyshev_pseudospectral", "uniform", 1.0e-15, "FFT"), ("chebyshev_pseudospectral", "uniform", 2.0e-15, "matrix"), ("chebyshev_pseudospectral", "sqrt", 1.0e-2, "FFT"), ("gausslegendre_pseudospectral", "uniform", 1.0e-14, "")),
ngrid ∈ (5,6,7,8,9,10), nelement ∈ (1, 2, 3, 4, 5)

Expand Down Expand Up @@ -111,6 +111,7 @@ function runtests()
rtol = 1.e2 / (nelement*(ngrid-1))^4
@test isapprox(df, expected_df, rtol=rtol, atol=1.e-15,
norm=maxabs_norm)
@test df[1] == df[end]
end
end

Expand Down Expand Up @@ -164,6 +165,7 @@ function runtests()
rtol = 1.e2 / (nelement*(ngrid-1))^order
@test isapprox(df, expected_df, rtol=rtol, atol=1.e-15,
norm=maxabs_norm)
df[1] == df[end]
end
end
end
Expand Down Expand Up @@ -278,7 +280,7 @@ function runtests()
end

@testset "Chebyshev pseudospectral derivatives (4 argument), periodic" verbose=false begin
@testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈
@testset "$nelement $ngrid $cheb_option" for (nelement, ngrid, rtol) ∈
(
(1, 5, 8.e-1),
(1, 6, 2.e-1),
Expand Down Expand Up @@ -469,11 +471,12 @@ function runtests()

@test isapprox(df, expected_df, rtol=rtol, atol=1.e-14,
norm=maxabs_norm)
@test df[1] == df[end]
end
end

@testset "Chebyshev pseudospectral derivatives upwinding (5 argument), periodic" verbose=false begin
@testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈
@testset "$nelement $ngrid $cheb_option" for (nelement, ngrid, rtol) ∈
(
(1, 5, 8.e-1),
(1, 6, 2.e-1),
Expand Down Expand Up @@ -668,12 +671,14 @@ function runtests()

@test isapprox(df, expected_df, rtol=rtol, atol=1.e-12,
norm=maxabs_norm)
@test df[1] == df[end]
end
end
end

@testset "Chebyshev pseudospectral derivatives (4 argument), polynomials" verbose=false begin
@testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"),
@testset "$nelement $ngrid $bc $element_spacing_option $cheb_option" for
bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"),
nelement ∈ (1:5), ngrid ∈ (3:33), cheb_option in ("FFT","matrix")

# define inputs needed for the test
Expand Down Expand Up @@ -721,7 +726,8 @@ function runtests()
end

@testset "Chebyshev pseudospectral derivatives upwinding (5 argument), polynomials" verbose=false begin
@testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"),
@testset "$nelement $ngrid $bc $element_spacing_option $cheb_option" for
bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"),
nelement ∈ (1:5), ngrid ∈ (3:33), cheb_option in ("FFT","matrix")

# define inputs needed for the test
Expand Down Expand Up @@ -884,6 +890,7 @@ function runtests()

@test isapprox(df, expected_df, rtol=rtol, atol=1.e-14,
norm=maxabs_norm)
@test df[1] == df[end]
end
end

Expand Down Expand Up @@ -1002,6 +1009,7 @@ function runtests()

@test isapprox(df, expected_df, rtol=rtol, atol=1.e-12,
norm=maxabs_norm)
@test df[1] == df[end]
end
end
end
Expand Down Expand Up @@ -1107,7 +1115,7 @@ function runtests()
end

@testset "Chebyshev pseudospectral second derivatives (4 argument), periodic" verbose=false begin
@testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈
@testset "$nelement $ngrid $cheb_option" for (nelement, ngrid, rtol) ∈
(
(1, 5, 8.e-1),
(1, 6, 2.e-1),
Expand Down Expand Up @@ -1298,11 +1306,12 @@ function runtests()

@test isapprox(d2f, expected_d2f, rtol=rtol, atol=1.e-10,
norm=maxabs_norm)
@test d2f[1] == d2f[end]
end
end

@testset "Chebyshev pseudospectral cylindrical laplacian derivatives (4 argument), zero" verbose=false begin
@testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈
@testset "$nelement $ngrid $cheb_option" for (nelement, ngrid, rtol) ∈
(
(4, 7, 2.e-1),
(4, 8, 2.e-1),
Expand Down Expand Up @@ -1497,6 +1506,7 @@ function runtests()

@test isapprox(d2f, expected_d2f, rtol=rtol, atol=1.e-10,
norm=maxabs_norm)
@test d2f[1] == d2f[end]
end
end

Expand Down

0 comments on commit ad9af46

Please sign in to comment.