Add support for Kahan & Li's s9 6th order int.

oashour · Feb 14, 2021 · 2b9a36a · 2b9a36a
1 parent b156e09
commit 2b9a36a
Show file tree

Hide file tree

Showing 5 changed files with 66 additions and 5 deletions.
diff --git a/docs/src/man/integrators.md b/docs/src/man/integrators.md
@@ -32,6 +32,8 @@ T6A_Ss14!
 T6B_Ss14!
 T6A_Ys7!
 T6B_Ys7!
+T6A_KLs9!
+T6B_KLs9!
 T8A_Ss15!
 T8B_Ss15!
 T1A_H!

diff --git a/src/CubicSolvers.jl b/src/CubicSolvers.jl
@@ -510,6 +510,64 @@ function T6B_Ys7!(ψₒ, ψᵢ, dx, ops)
     T2B!(ψₒ, ψₒ, γ₇*dx, ops)
 end
 
+"""
+    T6A_KLs9!(ψₒ, ψᵢ, dx, ops)
+
+Compute `ψₒ`, i.e. `ψᵢ` advanced a step `dx` forward using Kahan & Li's s9 Symplectic Sixth order integrator of type A. The structure `ops::Operators` contains the FFT plans and the kinetic energy operators.  
+
+See also: [`solve!`](@ref), [`Operators`](@ref)
+"""
+function T6A_KLs9!(ψₒ, ψᵢ, dx, ops)
+    γ₁ = 0.39216144400731413927925056
+    γ₂ = 0.33259913678935943859974864
+    γ₃ = -0.70624617255763935980996482
+    γ₄ =  0.08221359629355080023149045
+    γ₅ = 0.79854399093482996339895035
+    γ₉ = γ₁
+    γ₈ = γ₂
+    γ₇ = γ₃
+    γ₆ = γ₄
+
+    T2A!(ψₒ, ψᵢ, γ₁*dx, ops)
+    T2A!(ψₒ, ψₒ, γ₂*dx, ops)
+    T2A!(ψₒ, ψₒ, γ₃*dx, ops)
+    T2A!(ψₒ, ψₒ, γ₄*dx, ops)
+    T2A!(ψₒ, ψₒ, γ₅*dx, ops)
+    T2A!(ψₒ, ψₒ, γ₆*dx, ops)
+    T2A!(ψₒ, ψₒ, γ₇*dx, ops)
+    T2A!(ψₒ, ψₒ, γ₈*dx, ops)
+    T2A!(ψₒ, ψₒ, γ₉*dx, ops)
+end
+
+"""
+    T6B_KLs9!(ψₒ, ψᵢ, dx, ops)
+
+Compute `ψₒ`, i.e. `ψᵢ` advanced a step `dx` forward using Kahan & Li's s9 Symplectic Sixth order integrator of type B. The structure `ops::Operators` contains the FFT plans and the kinetic energy operators.  
+
+See also: [`solve!`](@ref), [`Operators`](@ref)
+"""
+function T6B_KLs9!(ψₒ, ψᵢ, dx, ops)
+    γ₁ = 0.39216144400731413927925056
+    γ₂ = 0.33259913678935943859974864
+    γ₃ = -0.70624617255763935980996482
+    γ₄ =  0.08221359629355080023149045
+    γ₅ = 0.79854399093482996339895035
+    γ₉ = γ₁
+    γ₈ = γ₂
+    γ₇ = γ₃
+    γ₆ = γ₄
+
+    T2B!(ψₒ, ψᵢ, γ₁*dx, ops)
+    T2B!(ψₒ, ψₒ, γ₂*dx, ops)
+    T2B!(ψₒ, ψₒ, γ₃*dx, ops)
+    T2B!(ψₒ, ψₒ, γ₄*dx, ops)
+    T2B!(ψₒ, ψₒ, γ₅*dx, ops)
+    T2B!(ψₒ, ψₒ, γ₆*dx, ops)
+    T2B!(ψₒ, ψₒ, γ₇*dx, ops)
+    T2B!(ψₒ, ψₒ, γ₈*dx, ops)
+    T2B!(ψₒ, ψₒ, γ₉*dx, ops)
+end
+
 """
     T6A_Ss14!(ψₒ, ψᵢ, dx, ops)
 

diff --git a/src/NonlinearSchrodinger.jl b/src/NonlinearSchrodinger.jl
@@ -25,7 +25,8 @@ export T1A!, T1B!, T2A!, T2B!
 export T4A_TJ!, T4B_TJ!, T6A_TJ!, T6B_TJ!, T8A_TJ!, T8B_TJ!
 export T4A_SF!, T4B_SF!, T6A_SF!, T6B_SF!, T8A_SF!, T8B_SF!
 export T4A_CMP!, T4B_CMP!, T6A_CMP!, T6B_CMP!, T8A_CMP!, T8B_CMP!
-export T6A_Ss14!, T6B_Ss14!, T6A_Ys7!, T6B_Ys7!, T8A_Ss15!, T8B_Ss15!
+export T6A_Ss14!, T6B_Ss14!, T6A_Ys7!, T6B_Ys7!, T6A_KLs9!, T6B_KLs9!
+export T8A_Ss15!, T8B_Ss15!
 export T1A_H!, T2A_H!
 export T1A_SS!
 

diff --git a/src/Simulation.jl b/src/Simulation.jl
@@ -24,9 +24,9 @@ function solve!(sim::Sim)
 
     @info "Starting evolution"
     # Define these tuples as global consts
-    if sim.T̂ ∈ (T1A!, T2A!, T4A_TJ!, T6A_TJ!, T8A_TJ!, T4A_SF!, T4A_SF!, T6A_SF!, T8A_SF!, T4A_CMP!, T6A_CMP!, T8A_CMP!, T6A_Ss14!, T6A_Ys7!, T8A_Ss15!, T1A_H!, T2A_H!, T1A_SS!)
+    if sim.T̂ ∈ (T1A!, T2A!, T4A_TJ!, T6A_TJ!, T8A_TJ!, T4A_SF!, T4A_SF!, T6A_SF!, T8A_SF!, T4A_CMP!, T6A_CMP!, T8A_CMP!, T6A_Ss14!, T6A_Ys7!, T6A_KLs9!, T8A_Ss15!, T1A_H!, T2A_H!, T1A_SS!)
         soln_loop_A(sim, ops, ind_p)
-    elseif sim.T̂ ∈ (T1B!, T2B!, T4B_TJ!, T6B_TJ!, T8B_TJ!, T4B_SF!, T6B_SF!, T8B_SF!, T4B_CMP!, T6B_CMP!, T8B_CMP!, T6B_Ss14!, T6B_Ys7!, T8B_Ss15!)
+    elseif sim.T̂ ∈ (T1B!, T2B!, T4B_TJ!, T6B_TJ!, T8B_TJ!, T4B_SF!, T6B_SF!, T8B_SF!, T4B_CMP!, T6B_CMP!, T8B_CMP!, T6B_Ss14!, T6B_Ys7!, T6B_KLs9!, T8B_Ss15!)
         soln_loop_B(sim, ops, ind_p)
     else
         throw(ArgumentError("Unknown algorithm $(sim.T̂)"))

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -150,9 +150,9 @@ using Test
             @test abs.(sim.ψ[:, end]) ≈ abs.(ψ₀) atol=1e-2
         end
         for algo ∈ (T2A!, T4A_TJ!, T6A_TJ!, T8A_TJ!, T4A_SF!, T4A_SF!, T6A_SF!, T8A_SF!, 
-                T4A_CMP!, T6A_CMP!, T8A_CMP!, T6A_Ss14!, T6A_Ys7!, T8A_Ss15!,
+                T4A_CMP!, T6A_CMP!, T8A_CMP!, T6A_Ss14!, T6A_Ys7!, T6A_KLs9!, T8A_Ss15!,
                 T2B!, T4B_TJ!, T6B_TJ!, T8B_TJ!, T4B_SF!, T4B_SF!, T6B_SF!, T8B_SF!,
-                T4B_CMP!, T6B_CMP!, T8B_CMP!, T6B_Ss14!, T6B_Ys7!, T8B_Ss15!)
+                T4B_CMP!, T6B_CMP!, T8B_CMP!, T6B_Ss14!, T6B_Ys7!, T6B_KLs9!,T8B_Ss15!)
             sim = Sim(λ, box, ψ₀, algo)
             @info "Testing Algorithm" algo
             solve!(sim)