From f7a6ded45e3e89748f0601c0167187a1f7c331b8 Mon Sep 17 00:00:00 2001 From: Omar Ashour Date: Sun, 14 Feb 2021 14:29:04 -0800 Subject: [PATCH] Add support for Kahan & Li's s17 8th order int. --- docs/src/man/integrators.md | 2 + src/CubicSolvers.jl | 90 +++++++++++++++++++++++++++++++++++++ src/NonlinearSchrodinger.jl | 2 +- src/Simulation.jl | 4 +- test/runtests.jl | 3 +- 5 files changed, 97 insertions(+), 4 deletions(-) diff --git a/docs/src/man/integrators.md b/docs/src/man/integrators.md index 62da33f..ad2d972 100644 --- 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 index ccd3514..94071a6 100644 --- 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 \ No newline at end of file diff --git a/src/NonlinearSchrodinger.jl b/src/NonlinearSchrodinger.jl index 19fbe89..0aa2a09 100644 --- 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 index 2b7545e..376332b 100644 --- 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 index e2be2ae..8b84caf 100644 --- 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)