Add support for Kahan & Li's s17 8th order int.

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

Hide file tree

Showing 5 changed files with 97 additions and 4 deletions.
diff --git a/docs/src/man/integrators.md b/docs/src/man/integrators.md
@@ -36,6 +36,8 @@ T6A_KLs9!
 T6B_KLs9!
 T8A_Ss15!
 T8B_Ss15!
+T8A_KLs17!
+T8B_KLs17!
 T1A_H!
 T2A_H!
 T1A_SS!

diff --git a/src/CubicSolvers.jl b/src/CubicSolvers.jl
@@ -726,4 +726,94 @@ function T8B_Ss15!(ψₒ, ψᵢ, dx, ops)
     T2B!(ψₒ, ψₒ, γ₁₃*dx, ops)
     T2B!(ψₒ, ψₒ, γ₁₄*dx, ops)
     T2B!(ψₒ, ψₒ, γ₁₅*dx, ops)
+end
+
+"""
+    T8A_KLs17!(ψₒ, ψᵢ, dx, ops)
+
+Compute `ψₒ`, i.e. `ψᵢ` advanced a step `dx` forward using Kahan & Li's s17 Symplectic Eighth 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 T8A_KLs17!(ψₒ, ψᵢ, dx, ops)
+    γ₁ = 0.13020248308889008087881763
+    γ₂ = 0.56116298177510838456196441 
+    γ₃ =  -0.38947496264484728640807860 
+    γ₄ = 0.15884190655515560089621075 
+    γ₅ = -0.39590389413323757733623154 
+    γ₆ = 0.18453964097831570709183254
+    γ₇ = 0.25837438768632204729397911
+    γ₈ = 0.29501172360931029887096624 
+    γ₉ = -0.60550853383003451169892108
+    γ₁₇ = γ₁
+    γ₁₆ = γ₂
+    γ₁₅ = γ₃
+    γ₁₄ = γ₄
+    γ₁₃ = γ₅
+    γ₁₂ = γ₆
+    γ₁₁ = γ₇
+    γ₁₀ = γ₈
+
+    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)
+    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
+
+"""
+    T8B_KLs17!(ψₒ, ψᵢ, dx, ops)
+
+Compute `ψₒ`, i.e. `ψᵢ` advanced a step `dx` forward using Kahan & Li's s17 Symplectic Eighth 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 T8B_KLs17!(ψₒ, ψᵢ, dx, ops)
+    γ₁ = 0.13020248308889008087881763
+    γ₂ = 0.56116298177510838456196441 
+    γ₃ =  -0.38947496264484728640807860 
+    γ₄ = 0.15884190655515560089621075 
+    γ₅ = -0.39590389413323757733623154 
+    γ₆ = 0.18453964097831570709183254
+    γ₇ = 0.25837438768632204729397911
+    γ₈ = 0.29501172360931029887096624 
+    γ₉ = -0.60550853383003451169892108
+    γ₁₇ = γ₁
+    γ₁₆ = γ₂
+    γ₁₅ = γ₃
+    γ₁₄ = γ₄
+    γ₁₃ = γ₅
+    γ₁₂ = γ₆
+    γ₁₁ = γ₇
+    γ₁₀ = γ₈
+
+    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)
+    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
diff --git a/src/NonlinearSchrodinger.jl b/src/NonlinearSchrodinger.jl
@@ -26,7 +26,7 @@ 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!, T6A_KLs9!, T6B_KLs9!
-export T8A_Ss15!, T8B_Ss15!
+export T8A_Ss15!, T8B_Ss15!, T8A_KLs17!, T8B_KLs17!
 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!, T6A_KLs9!, 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!, T8A_KLs17!, 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!, T6B_KLs9!, 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!, T8B_KLs17!)
         soln_loop_B(sim, ops, ind_p)
     else
         throw(ArgumentError("Unknown algorithm $(sim.T̂)"))

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -151,8 +151,9 @@ using Test
         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!, T6A_KLs9!, T8A_Ss15!,
+                T8A_KLs17!,
                 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!, T6B_KLs9!,T8B_Ss15!)
+                T4B_CMP!, T6B_CMP!, T8B_CMP!, T6B_Ss14!, T6B_Ys7!, T6B_KLs9!, T8B_Ss15!, T8B_KLs17!)
             sim = Sim(λ, box, ψ₀, algo)
             @info "Testing Algorithm" algo
             solve!(sim)