Update for GeometricIntegrators v0.7.

Project.toml

@@ -20,7 +20,7 @@ ApproxFun = "0.10, 0.11, 0.12"
 Documenter = "0.23, 0.24, 0.25, 0.26"
 DoubleFloats = "1"
 FFTW = "1"
-GeometricIntegrators = "0.5, 0.6"
+GeometricIntegrators = "0.7"
 HDF5 = "0.10, 0.12, 0.13, 0.14"
 OffsetArrays = "1"
 ProgressMeter = "1"

src/poincare_invariant_1st.jl

@@ -27,11 +27,7 @@ function PoincareInvariant1st(f_equ::Function, f_loop::Function, Θ::ΘT, Δt::T
     end
 
     # compute initial conditions
-    q₀ = zeros(DT, (d, nloop))
-
-    for i in 1:nloop
-        q₀[:,i] .= f_loop(i/nloop)
-    end
+    q₀ = [f_loop(i/nloop) for i in 1:nloop]
 
     # initialise euation
     equ = f_equ(q₀)
@@ -48,27 +44,27 @@ end
 
 
 function evaluate_poincare_invariant(pinv::PoincareInvariant1st, sol::Solution)
-    p = zero(sol.q)
-    g = zero(sol.q)
-    v = zeros(size(sol.q,1), size(sol.q,3))
-    γ = zeros(size(sol.q,1), size(sol.q,3))
+    p = zeros(size(sol.q[begin],1), size(sol.q,1), size(sol.q,2))
+    g = zeros(size(sol.q[begin],1), size(sol.q,1), size(sol.q,2))
+    v = zeros(size(sol.q[begin],1), size(sol.q,2))
+    γ = zeros(size(sol.q[begin],1), size(sol.q,2))
 
     compute_one_form(sol.t, sol.q, p, pinv.Θ)
 
     if isdefined(sol, :λ)
         compute_correction(sol.t, sol.q, sol.λ, g, pinv.equ.g)
     end
 
-    for i in axes(sol.q,2)
-        compute_velocity(sol.q[:,i,:], v)
-        pinv.I[i] = compute_loop_integral(p[:,i,:], v)
+    for i in axes(sol.q,1)
+        compute_velocity(hcat(sol.q[i,:]...), v)
+        pinv.I[i] = compute_loop_integral(hcat(sol.p[i,:]...), v)
 
         if isdefined(sol, :p)
-            pinv.J[i] = compute_loop_integral(sol.d[:,i,:], v)
+            pinv.J[i] = compute_loop_integral(hcat(sol.d[i,:]...), v)
         end
 
         if isdefined(sol, :λ)
-            compute_velocity(sol.λ[:,i,:], γ)
+            compute_velocity(hcat(sol.λ[i,:]...), γ)
             pinv.L[i] = compute_loop_integral(p[:,i,:] .- pinv.Δt .* g[:,i,:], v .- pinv.Δt .* γ)
         end
     end

src/poincare_invariant_1st_canonical.jl

@@ -24,13 +24,8 @@ function PoincareInvariant1stCanonical(f_equ::Function, f_loop_q::Function, f_lo
     end
 
     # compute initial conditions
-    q₀ = zeros(DT, (d, nloop))
-    p₀ = zeros(DT, (d, nloop))
-
-    for i in 1:nloop
-        q₀[:,i] .= f_loop_q(i/nloop)
-        p₀[:,i] .= f_loop_p(i/nloop)
-    end
+    q₀ = [f_loop_q(i/nloop) for i in 1:nloop]
+    p₀ = [f_loop_p(i/nloop) for i in 1:nloop]
 
     # initialise euation
     equ = f_equ(q₀, p₀)
@@ -45,11 +40,11 @@ end
 
 
 function evaluate_poincare_invariant(pinv::PoincareInvariant1stCanonical, sol::SolutionPODE)
-    v = zeros(size(sol.q,1), size(sol.q,3))
+    v = zeros(size(sol.q[begin],1), size(sol.q,2))
 
-    for i in axes(sol.q,2)
-        compute_velocity(sol.q[:,i,:], v)
-        pinv.I[i] = compute_loop_integral(sol.p[:,i,:], v)
+    for i in axes(sol.q,1)
+        compute_velocity(hcat(sol.q[i,:]...), v)
+        pinv.I[i] = compute_loop_integral(hcat(sol.p[i,:]...), v)
     end
 
     pinv.I

src/poincare_invariant_2nd_approxfun.jl

@@ -28,11 +28,7 @@ function PoincareInvariant2ndApproxFun(f_equ::Function, f_surface::Function, ω:
     c = ApproxFun.points(ApproxFun.Chebyshev(0..1)^2, nx*ny)
 
     # compute initial conditions
-    q₀ = zeros(DT, (nd, length(c)))
-
-    for i in eachindex(c)
-        q₀[:,i] .= f_surface(c[i][1], c[i][2])
-    end
+    q₀ = [f_surface(c[i][1], c[i][2]) for i in eachindex(c)]
 
     # initialise euation
     equ = f_equ(q₀)
@@ -53,7 +49,7 @@ function PoincareInvariant2ndApproxFun(f_equ::Function, f_surface::Function, ω:
     Dx = ApproxFun.Derivative(SC, [1,0])
     Dy = ApproxFun.Derivative(SC, [0,1])
 
-    fqc = Fun(SC, ApproxFun.transform(SC, q₀[1,:]))
+    fqc = Fun(SC, ApproxFun.transform(SC, hcat(q₀...)[1,:]))
     fqu = Fun(Dx * fqc, SU)
     nc = ApproxFun.ncoefficients(fqu)
     nv = length(ApproxFun.values(fqu))
@@ -98,18 +94,18 @@ function evaluate_poincare_invariant(pinv::PoincareInvariant2ndApproxFun{DT}, so
     #     local λ = permutedims(sol.λ, [3,1,2])
     # end
 
-    verbosity ≤ 1 ? prog = Progress(size(sol.q,2), 5) : nothing
+    verbosity ≤ 1 ? prog = Progress(size(sol.q,1), 5) : nothing
 
-    for i in axes(sol.q,2)
+    for i in axes(sol.q,1)
         verbosity > 1 ? println("      it = ", i) : nothing
-        pinv.I[i] = compute_noncanonical_invariant(pinv, sol.t[i], sol.q[:,i,:])
+        pinv.I[i] = compute_noncanonical_invariant(pinv, sol.t[i], hcat(sol.q[i,:]...))
         pinv.ΔI[i] = abs(pinv.I[0]) < sqrt(eps()) ? pinv.I[i] : (pinv.I[i] .- pinv.I[0]) ./ pinv.I[0]
         verbosity > 1 ? println("           I_q = ", pinv.I[i], ",   ε_q = ", pinv.ΔI[i]) : nothing
     end
 
     if isdefined(sol, :p)
-        for i in axes(sol.q,2)
-            pinv.J[i] = compute_canonical_invariant(pinv, sol.q[:,i,:], sol.p[:,i,:])
+        for i in axes(sol.q,1)
+            pinv.J[i] = compute_canonical_invariant(pinv, hcat(sol.q[i,:]...), hcat(sol.p[i,:]...))
             pinv.ΔJ[i] = abs(pinv.J[0]) < sqrt(eps()) ? pinv.J[i] : (pinv.J[i] .- pinv.J[0]) ./ pinv.J[0]
             verbosity > 1 ? println("           I_p = ", pinv.J[i], ",   ε_p = ", pinv.ΔJ[i]) : nothing
         end
@@ -124,8 +120,8 @@ function evaluate_poincare_invariant_correction(pinv::PoincareInvariant2ndApprox
     local verbosity = get_config(:verbosity)
 
     if isdefined(sol, :λ)
-        for i in axes(sol.q,2)
-            pinv.K[i] = compute_noncanonical_correction(pinv, sol.t[i], sol.q[:,i,:], sol.λ[:,i,:])
+        for i in axes(sol.q,1)
+            pinv.K[i] = compute_noncanonical_correction(pinv, sol.t[i], hcat(sol.q[i,:]...), hcat(sol.λ[i,:]...))
             pinv.L[i] = pinv.I[i] - pinv.Δt^2 * pinv.K[i]
             verbosity > 1 ? println("           I_λ = ", pinv.L[i], ",   ε_λ = ", (pinv.L[i]-pinv.L[0])/pinv.L[0]) : nothing
             verbosity > 1 ? println("           K_λ = ", pinv.K[i]) : nothing
@@ -185,13 +181,8 @@ function PoincareInvariant2ndApproxFunCanonical(f_equ::Function, f_surface_q::Fu
     c = ApproxFun.points(ApproxFun.Chebyshev(0..1)^2, nx*ny)
 
     # compute initial conditions
-    q₀ = zeros(DT, (nd, length(c)))
-    p₀ = zeros(DT, (nd, length(c)))
-
-    for i in eachindex(c)
-        q₀[:,i] .= f_surface_q(c[i][1], c[i][2])
-        p₀[:,i] .= f_surface_p(c[i][1], c[i][2])
-    end
+    q₀ = [f_surface_q(c[i][1], c[i][2]) for i in eachindex(c)]
+    p₀ = [f_surface_p(c[i][1], c[i][2]) for i in eachindex(c)]
 
     # initialise euation
     equ = f_equ(q₀, p₀)
@@ -208,7 +199,7 @@ function PoincareInvariant2ndApproxFunCanonical(f_equ::Function, f_surface_q::Fu
     Dx = ApproxFun.Derivative(SC, [1,0])
     Dy = ApproxFun.Derivative(SC, [0,1])
 
-    fq = Fun(Dx * Fun(SC, ApproxFun.transform(SC, q₀[1,:])), SU)
+    fq = Fun(Dx * Fun(SC, ApproxFun.transform(SC, hcat(q₀...)[1,:])), SU)
     nc = ApproxFun.ncoefficients(fq)
     nv = length(ApproxFun.values(fq))
 
@@ -220,11 +211,11 @@ end
 function evaluate_poincare_invariant(pinv::PoincareInvariant2ndApproxFunCanonical{DT}, sol::Solution) where {DT}
     local verbosity = get_config(:verbosity)
 
-    verbosity ≤ 1 ? prog = Progress(size(sol.q,2), 5) : nothing
+    verbosity ≤ 1 ? prog = Progress(size(sol.q,1), 5) : nothing
 
-    for i in axes(sol.q,2)
+    for i in axes(sol.q,1)
         verbosity > 1 ? println("      it = ", i) : nothing
-        pinv.I[i] = compute_canonical_invariant(pinv, sol.q[:,i,:], sol.p[:,i,:])
+        pinv.I[i] = compute_canonical_invariant(pinv, hcat(sol.q[i,:]...), hcat(sol.p[i,:]...))
         verbosity > 1 ? println("           I_p = ", pinv.I[i], ",   ε_p = ", (pinv.I[i]-pinv.I[0])/pinv.I[0]) : nothing
 
         verbosity ≤ 1 ? next!(prog) : nothing

src/poincare_invariant_2nd_opq.jl

@@ -47,15 +47,9 @@ function PoincareInvariant2ndOPQ(f_equ::Function, f_surface::Function, ω::ΩT,
     xGrid = OPQ.grid(OPQ.ChebyshevT()[:,1:nx])
     yGrid = OPQ.grid(OPQ.ChebyshevT()[:,1:ny])
 
-    q₀ = zeros(DT, (d, nx, ny))
+    q₀ = reshape([f_surface((xGrid[i] + 1) / 2, (yGrid[j] + 1) / 2) for i in 1:nx, j in 1:ny], nx*ny)
 
-    for j in 1:ny
-        for i in 1:nx
-            q₀[:,i,j] = f_surface((xGrid[i] + 1) / 2, (yGrid[j] + 1) / 2)
-        end
-    end
-
-    equ = f_equ(reshape(q₀, (d, nx*ny)))
+    equ = f_equ(q₀)
 
 
     # compute initial conditions
@@ -452,12 +446,12 @@ function evaluate_poincare_invariant(pinv::PoincareInvariant2ndOPQ{DT}, sol::Sol
 
     for i in axes(sol.q,2)
         verbosity > 1 ? println("      it = ", i) : nothing
-        @views pinv.I[i]  = surface_integral(pinv, sol.t[i], sol.q[:,i,:])
+        @views pinv.I[i]  = surface_integral(pinv, sol.t[i], hcat(sol.q[i,:]...))
         pinv.ΔI[i] = abs(pinv.I[0]) < sqrt(eps()) ? pinv.I[i] : (pinv.I[i] .- pinv.I[0]) ./ pinv.I[0]
         verbosity > 1 ? println("           I_q = ", pinv.I[i], ",   ε_q = ", pinv.ΔI[i]) : nothing
 
         # if hasproperty(sol, :p)
-        #     pinv.J[i] = surface_integral_canonical(pinv, sol.q[:,i,:], sol.p[:,i,:])
+        #     pinv.J[i] = surface_integral_canonical(pinv, hcat(sol.q[i,:]...), hcat(sol.p[i,:]...))
         #     pinv.ΔJ[i] = abs(pinv.J[0]) < sqrt(eps()) ? pinv.J[i] : (pinv.J[i] .- pinv.J[0]) ./ pinv.J[0]
         #     verbosity > 1 ? println("           I_p = ", pinv.J[i], ",   ε_p = ", pinv.ΔJ[i]) : nothing
         # end

src/poincare_invariant_2nd_trapezoidal.jl

@@ -30,15 +30,9 @@ function PoincareInvariant2ndTrapezoidal(f_equ::Function, f_surface::Function, 
     end
 
     # compute initial conditions
-    q₀ = zeros(DT, (d, nx, ny))
+    q₀ = reshape([f_surface(i/nx, j/ny) for i in 1:nx, j in 1:ny], nx*ny)
 
-    for j in 1:ny
-        for i in 1:nx
-            q₀[:,i,j] = f_surface(i/nx, j/ny)
-        end
-    end
-
-    equ = f_equ(reshape(q₀, (d, nx*ny)))
+    equ = f_equ(q₀)
 
     # create arrays for results
     nt = div(ntime, nsave)
@@ -211,14 +205,14 @@ end
 function evaluate_poincare_invariant(pinv::PoincareInvariant2ndTrapezoidal{DT}, sol::Solution) where {DT}
     local verbosity = get_config(:verbosity)
 
-    for i in axes(sol.q,2)
+    for i in axes(sol.q,1)
         verbosity > 1 ? println("      it = ", i) : nothing
-        @views pinv.I[i]  = surface_integral(sol.t[i], sol.q[:,i,:], pinv.ω, pinv.nx, pinv.ny)
+        @views pinv.I[i]  = surface_integral(sol.t[i], hcat(sol.q[i,:]...), pinv.ω, pinv.nx, pinv.ny)
         pinv.ΔI[i] = abs(pinv.I[0]) < sqrt(eps()) ? pinv.I[i] : (pinv.I[i] .- pinv.I[0]) ./ pinv.I[0]
         verbosity > 1 ? println("           I_q = ", pinv.I[i], ",   ε_q = ", pinv.ΔI[i]) : nothing
 
         if hasproperty(sol, :p)
-            @views pinv.J[i] = surface_integral_canonical(sol.q[:,i,:], sol.p[:,i,:], pinv.nx, pinv.ny)
+            @views pinv.J[i] = surface_integral_canonical(hcat(sol.q[i,:]...), hcat(sol.p[i,:]...), pinv.nx, pinv.ny)
             pinv.ΔJ[i] = abs(pinv.J[0]) < sqrt(eps()) ? pinv.J[i] : (pinv.J[i] .- pinv.J[0]) ./ pinv.J[0]
             verbosity > 1 ? println("           I_p = ", pinv.J[i], ",   ε_p = ", pinv.ΔJ[i]) : nothing
         end

test/poincare_invariant_2nd_module.jl

@@ -28,6 +28,12 @@ module PoincareInvariant2ndTest
         nothing
     end
 
+    function gcϑ(t, q)
+        ϑ = zero(q)
+        gcϑ(t, q, ϑ)
+        return ϑ
+    end
+
     function gcω(t, q, Ω)
         Ω[1,1] = 0
         Ω[1,2] =-B(t,q)
@@ -62,10 +68,7 @@ module PoincareInvariant2ndTest
     end
 
     function gc_surface_p(s,t)
-        q = gc_surface_q(s,t)
-        p = zero(q)
-        gcϑ(zero(eltype(q)), q, p)
-        p
+        gcϑ(0, gc_surface_q(s,t))
     end
 
     function gc_dummy(t, a, b, c)
@@ -77,21 +80,7 @@ module PoincareInvariant2ndTest
     end
 
     function gc_dummy_iode(q₀)
-        p₀ = zero(q₀)
-
-        if ndims(q₀) == 1
-            gcϑ(zero(eltype(q₀)), q₀, p₀)
-        else
-            tq = zeros(eltype(q₀), size(q₀,1))
-            tp = zeros(eltype(p₀), size(p₀,1))
-
-            for i in axes(q₀,2)
-                tq .= q₀[:,i]
-                gcϑ(zero(eltype(q₀)), tq, tp)
-                p₀[:,i] .= tp
-            end
-        end
-
+        p₀ = [gcϑ(0,q) for q in q₀]
         IODE(gc_dummy, gc_dummy, gc_dummy, q₀, p₀; v̄=gc_dummy)
     end