Skip to content

Commit

Permalink
Add difficult test cases
Browse files Browse the repository at this point in the history
  • Loading branch information
@Luapulu
Luapulu committed Mar 23, 2022
1 parent 61eb395 commit c714417
Showing 1 changed file with 77 additions and 63 deletions.
140 changes: 77 additions & 63 deletions test/SecondPoincareInvariants/test_SecondPoincareInvariants.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,8 @@ end

@safetestset "Consistency Between Implementations" begin
using PoincareInvariants
using PoincareInvariants.SecondPoincareInvariants.Chebyshev: ChebyshevPlan
using PoincareInvariants.SecondPoincareInvariants.FiniteDifferences: FiniteDiffPlan
using PoincareInvariants.SecondPoincareInvariants: ChebyshevPlan
using PoincareInvariants.SecondPoincareInvariants: FiniteDiffPlan

D = 6
Ω(z, t, p) = [
Expand All @@ -52,91 +52,105 @@ end
x + y
)

Idefault = let pinv = PI2{Float64}(Ω, D, 10_000)
Idefault = let pinv = PI2{Float64}(Ω, D, 1_000)
compute!(pinv, getpoints(f, pinv), 0, nothing)
end

Icheb = let pinv = PI2{Float64}(Ω, D, 10_000, ChebyshevPlan)
Icheb = let pinv = PI2{Float64}(Ω, D, 1_000, ChebyshevPlan)
compute!(pinv, getpoints(f, pinv), 0, nothing)
end

Ifindiff = let pinv = PI2{Float64}(Ω, D, (100, 100), FiniteDiffPlan)
Ifindiff = let pinv = PI2{Float64}(Ω, D, (31, 33), FiniteDiffPlan)
compute!(pinv, getpoints(f, pinv), 0, nothing)
end

@test Icheb ≈ Idefault rtol=10eps()
@test Icheb ≈ Ifindiff rtol=1e-3
end

@safetestset "Free Particle" begin
using PoincareInvariants

function free_particle!(points, t)
mid = length(points) ÷ 2
for i in 1:mid
points[i] .+= points[mid+i] .* t
end
end

@testset "Square" begin
D = 8
N = 1_000
Ω(z, t, p) = CanonicalSymplecticMatrix(D)
pinv = SecondPoincareInvariant{Float64}(Ω, D, N)

phasepoints = getpoints(pinv) do x, y
(x, 0, 0, 0, y, 0, 0, 0)
end
@safetestset "Linear Symplectic Maps" begin
@safetestset "In 2D" begin
using PoincareInvariants
using PoincareInvariants.SecondPoincareInvariants: FiniteDiffPlan

@test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 10
Ω(z, t, p) = CanonicalSymplecticMatrix(2)
pinv1 = PI2{Float64}(Ω, 2, 1_000)
pinv2 = PI2{Float64}(Ω, 2, 1_000, FiniteDiffPlan)

free_particle!(phasepoints, 10)
@test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 20
A = [1 5;
0 1]
B = [2 0;
0 0.5]
f(r, θ) = A * B * [r*cospi(θ), r*sinpi(θ)]
# half circle with radius one

free_particle!(phasepoints, 100)
@test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 100
@test abs(π/2 - compute!(pinv1, getpoints(f, pinv1), 0, nothing)) / eps() < 10
@test abs(π/2 - compute!(pinv2, getpoints(f, pinv2), 0, nothing)) / π/2 < 1e-3
end

@testset "Quarter Circle" begin
D = 4
N = 2_000
Ω(z, t, p) = CanonicalSymplecticMatrix(D)
pinv = SecondPoincareInvariant{Float64}(Ω, D, N)
@safetestset "In 8D" begin
using PoincareInvariants
using PoincareInvariants.SecondPoincareInvariants: FiniteDiffPlan

Ω(z, t, p) = CanonicalSymplecticMatrix(8)
pinv1 = PI2{Float64}(Ω, 8, 1_000)
pinv2 = PI2{Float64}(Ω, 8, 1_000, FiniteDiffPlan)

A = CanonicalSymplecticMatrix(8)

B = [1 0 0 0 1 5 8 1;
0 1 0 0 5 2 6 9;
0 0 1 0 8 6 3 7;
0 0 0 1 1 9 7 4;
0 0 0 0 1 0 0 0;
0 0 0 0 0 1 0 0;
0 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 1]

C = [1.0 0 2.0 0 0 0 0 0;
0 2.0 0 1.0 0 0 0 0;
0 1.0 0 -2.0 0 0 0 0;
2.0 0 -1.0 0 0 0 0 0;
0 0 0 0 0.2 0 0.4 0;
0 0 0 0 0 0.4 0 0.2;
0 0 0 0 0.2 0 0 -0.4;
0 0 0 0 0.4 0 -0.2 0]

f(r, θ) = A * B * C * [r*sinpi(θ), 0, 0, 0, r*cospi(θ), 0, 0, 0]
# half circle with radius one and negative orientation

@test abs(-π/2 - compute!(pinv1, getpoints(f, pinv1), 0, nothing)) / eps() < 10
@test abs(-π/2 - compute!(pinv2, getpoints(f, pinv2), 0, nothing)) / π/2 < 1e-3
end
end

phasepoints = getpoints(pinv) do r, θ
(0, r * cos(θ * π/2), 0, r * sin(θ * π/2))
end
@safetestset "Henon Map" begin
using PoincareInvariants
using PoincareInvariants.SecondPoincareInvariants: FiniteDiffPlan

@test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 10
Ω(z, t, p) = CanonicalSymplecticMatrix(4)
pinv1 = PI2{Float64}(Ω, 4, 1_000)
pinv2 = PI2{Float64}(Ω, 4, 1_000, FiniteDiffPlan)

free_particle!(phasepoints, 10)
@test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 20
init(x, y) = (x, -x, y, -y)
henon((q, p), a) = (p + a - q^2, -q)

free_particle!(phasepoints, 100)
@test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 200
# henon plus symplectic intermixing of q1/p1 with q2/p2
function f((q1, q2, p1, p2))
q1, p1 = henon((q1, p1), 2)
q2, p2 = henon((q2, p2), 3)
return (q1 - p2, q2 - p1, p1, p2)
end

@testset "Half Sphere in 4D" begin
D = 4
N = 1500
Ω(z, t, p) = CanonicalSymplecticMatrix(D)
pinv = SecondPoincareInvariant{Float64}(Ω, D, N)

phasepoints = getpoints(pinv) do θ, ϕ
sinθ, cosθ = sincospi(θ)
sinϕ, cosϕ = sincospi(ϕ)
return (sinθ * cosϕ, sinθ * sinϕ, cosθ, -1)
end
f1(x, y) = init(x, y) |> f
@test abs(2 - compute!(pinv1, getpoints(f1, pinv1), 0, nothing)) / eps() < 10
@test abs(2 - compute!(pinv2, getpoints(f1, pinv2), 0, nothing)) / eps() < 500

# invariant is -π
# TODO: do the calculation by hand and double check it
f2(x, y) = init(x, y) |> f |> f
@test abs(2 - compute!(pinv1, getpoints(f2, pinv1), 0, nothing)) / eps() < 50
@test abs(2 - compute!(pinv2, getpoints(f2, pinv2), 0, nothing)) < 0.01

@test abs(-π - compute!(pinv, phasepoints, 0, nothing)) / eps() < 15

free_particle!(phasepoints, 10)
@test abs(-π - compute!(pinv, phasepoints, 0, nothing)) / eps() < 15

free_particle!(phasepoints, 100)
@test abs(-π - compute!(pinv, phasepoints, 0, nothing)) / eps() < 150
end
f3(x, y) = init(x, y) |> f |> f |> f
@test abs(2 - compute!(pinv1, getpoints(f3, pinv1), 0, nothing)) / eps() < 50
@test abs(2 - compute!(pinv2, getpoints(f3, pinv2), 0, nothing)) < 0.5
end

0 comments on commit c714417

Please sign in to comment.