From 063c523390ddb0fbb497dba36f621eb97801ca79 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 10:24:42 -0400
Subject: [PATCH 01/17] First working version

---
 Project.toml                                  |   1 +
 lib/OrdinaryDiffEqBoris/Project.toml          |  25 ++
 .../src/OrdinaryDiffEqBoris.jl                |  26 ++
 lib/OrdinaryDiffEqBoris/src/alg_utils.jl      |   8 +
 lib/OrdinaryDiffEqBoris/src/algorithms.jl     |  32 ++
 lib/OrdinaryDiffEqBoris/src/boris_caches.jl   |  74 ++++
 .../src/boris_perform_step.jl                 | 326 ++++++++++++++++++
 lib/OrdinaryDiffEqBoris/test/runtests.jl      | 144 ++++++++
 8 files changed, 636 insertions(+)
 create mode 100644 lib/OrdinaryDiffEqBoris/Project.toml
 create mode 100644 lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
 create mode 100644 lib/OrdinaryDiffEqBoris/src/alg_utils.jl
 create mode 100644 lib/OrdinaryDiffEqBoris/src/algorithms.jl
 create mode 100644 lib/OrdinaryDiffEqBoris/src/boris_caches.jl
 create mode 100644 lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
 create mode 100644 lib/OrdinaryDiffEqBoris/test/runtests.jl

diff --git a/Project.toml b/Project.toml
index ab595de9c..6c88be42b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,6 +17,7 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+OrdinaryDiffEqBoris = "1e8b2ed2-ea1c-4b6e-b3f5-67c29d6bd197"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
diff --git a/lib/OrdinaryDiffEqBoris/Project.toml b/lib/OrdinaryDiffEqBoris/Project.toml
new file mode 100644
index 000000000..a5521972a
--- /dev/null
+++ b/lib/OrdinaryDiffEqBoris/Project.toml
@@ -0,0 +1,25 @@
+name = "OrdinaryDiffEqBoris"
+uuid = "1e8b2ed2-ea1c-4b6e-b3f5-67c29d6bd197"
+version = "0.1.0"
+
+[deps]
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[compat]
+OrdinaryDiffEqCore = "3"
+RecursiveArrayTools = "3.24, 4"
+SciMLBase = "2, 3"
+StaticArrays = "1"
+MuladdMacro = "0.2"
+julia = "1.10"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
diff --git a/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl b/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
new file mode 100644
index 000000000..e2f0944f9
--- /dev/null
+++ b/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
@@ -0,0 +1,26 @@
+module OrdinaryDiffEqBoris
+
+import OrdinaryDiffEqCore
+import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, OrdinaryDiffEqAdaptiveAlgorithm,
+    OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
+    alg_order, alg_cache, isfsal, initialize!, perform_step!
+
+
+using SciMLBase
+using RecursiveArrayTools
+using StaticArrays
+using MuladdMacro
+using LinearAlgebra
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("boris_caches.jl")
+include("boris_perform_step.jl")
+
+export Boris
+export MultistepBoris
+export AdaptiveBoris
+export boris_velocity_update # Exposing strictly for debugging if needed
+export update_velocity_multistep
+
+end
diff --git a/lib/OrdinaryDiffEqBoris/src/alg_utils.jl b/lib/OrdinaryDiffEqBoris/src/alg_utils.jl
new file mode 100644
index 000000000..a8b8e9f0e
--- /dev/null
+++ b/lib/OrdinaryDiffEqBoris/src/alg_utils.jl
@@ -0,0 +1,8 @@
+OrdinaryDiffEqCore.alg_order(alg::Boris) = 2
+OrdinaryDiffEqCore.isfsal(alg::Boris) = false
+
+OrdinaryDiffEqCore.alg_order(alg::MultistepBoris) = 2
+OrdinaryDiffEqCore.isfsal(alg::MultistepBoris) = false
+
+OrdinaryDiffEqCore.alg_order(alg::AdaptiveBoris) = 2
+OrdinaryDiffEqCore.isfsal(alg::AdaptiveBoris) = false
diff --git a/lib/OrdinaryDiffEqBoris/src/algorithms.jl b/lib/OrdinaryDiffEqBoris/src/algorithms.jl
new file mode 100644
index 000000000..1ec0a4271
--- /dev/null
+++ b/lib/OrdinaryDiffEqBoris/src/algorithms.jl
@@ -0,0 +1,32 @@
+"""
+    Boris()
+
+The standard Boris method for particle pushing in electric and magnetic fields.
+
+This solver expects a problem where `p` is structured as `(q2m, m, E, B, ...)`, which matches the signature used by TestParticle.jl's `TraceProblem`.
+"""
+struct Boris <: OrdinaryDiffEqAlgorithm end
+
+"""
+    MultistepBoris(; n=1, N=2)
+
+The Multistep/Hyper Boris method.
+`n` specifies the number of subcycles.
+`N` specifies the gyrophase correction order. `N=2` corresponds to the Multicycle solver, while `N=4` or `N=6` is the Hyper Boris solver.
+"""
+struct MultistepBoris <: OrdinaryDiffEqAlgorithm
+    n::Int
+    N::Int
+end
+MultistepBoris(; n::Int = 1, N::Int = 2) = MultistepBoris(n, N)
+
+"""
+    AdaptiveBoris(; safety=0.1)
+
+Adaptive Boris method with adaptive time stepping based on the local gyroperiod.
+The time step is evaluated as `dt = safety * 2π / |qB/m|`.
+"""
+struct AdaptiveBoris{T} <: OrdinaryDiffEqAlgorithm
+    safety::T
+end
+AdaptiveBoris(; safety = 0.1) = AdaptiveBoris(safety)
diff --git a/lib/OrdinaryDiffEqBoris/src/boris_caches.jl b/lib/OrdinaryDiffEqBoris/src/boris_caches.jl
new file mode 100644
index 000000000..2e0ec7758
--- /dev/null
+++ b/lib/OrdinaryDiffEqBoris/src/boris_caches.jl
@@ -0,0 +1,74 @@
+struct BorisConstantCache <: OrdinaryDiffEqConstantCache end
+
+struct BorisCache{uType, rateType} <: OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    tmp::uType
+    k::rateType
+end
+
+function OrdinaryDiffEqCore.alg_cache(
+        alg::Boris, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
+        ::Val{false}, args...; kwargs...
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    return BorisConstantCache()
+end
+
+function OrdinaryDiffEqCore.alg_cache(
+        alg::Boris, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
+        ::Val{true}, args...; kwargs...
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    return BorisCache(u, uprev, similar(u), similar(rate_prototype))
+end
+
+struct MultistepBorisConstantCache <: OrdinaryDiffEqConstantCache end
+
+struct MultistepBorisCache{uType, rateType} <: OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    tmp::uType
+    k::rateType
+end
+
+function OrdinaryDiffEqCore.alg_cache(
+        alg::MultistepBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
+        ::Val{false}, args...; kwargs...
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    return MultistepBorisConstantCache()
+end
+
+function OrdinaryDiffEqCore.alg_cache(
+        alg::MultistepBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
+        ::Val{true}, args...; kwargs...
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    return MultistepBorisCache(u, uprev, similar(u), similar(rate_prototype))
+end
+
+struct AdaptiveBorisConstantCache <: OrdinaryDiffEqConstantCache end
+
+struct AdaptiveBorisCache{uType, rateType} <: OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    tmp::uType
+    k::rateType
+end
+
+function OrdinaryDiffEqCore.alg_cache(
+        alg::AdaptiveBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
+        ::Val{false}, args...; kwargs...
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    return AdaptiveBorisConstantCache()
+end
+
+function OrdinaryDiffEqCore.alg_cache(
+        alg::AdaptiveBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
+        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
+        ::Val{true}, args...; kwargs...
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    return AdaptiveBorisCache(u, uprev, similar(u), similar(rate_prototype))
+end
diff --git a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
new file mode 100644
index 000000000..0a0d9b499
--- /dev/null
+++ b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
@@ -0,0 +1,326 @@
+const TN_MAG_THRESHOLD = 1.0e-4
+
+@inline @muladd function boris_velocity_update(v, E, B, qdt_2m)
+    t_rotate = qdt_2m * B
+    t_mag2 = sum(abs2, t_rotate)
+    s_rotate = 2 * t_rotate / (1 + t_mag2)
+
+    v⁻ = v + qdt_2m * E
+    v′ = v⁻ + (v⁻ × t_rotate)
+    v⁺ = v⁻ + (v′ × s_rotate)
+
+    v_new = v⁺ + qdt_2m * E
+
+    return v_new
+end
+
+function OrdinaryDiffEqCore.initialize!(integrator, cache::BorisConstantCache)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
+    integrator.stats.nf += 1
+    integrator.kshortsize = 0
+    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+end
+
+function OrdinaryDiffEqCore.perform_step!(integrator, cache::BorisConstantCache, repeat_step = false)
+    t = integrator.t
+    dt = integrator.dt
+    uprev = integrator.uprev
+    f = integrator.f
+    p = integrator.p
+
+    # Extract particles from uprev
+    # TestParticle.jl standard: u[1:3] is r, u[4:6] is v
+    r = uprev[SVector(1, 2, 3)]
+    v = uprev[SVector(4, 5, 6)]
+
+    # Expected p structure in TestParticle.jl: (q2m, m, Efunc, Bfunc, Ffunc)
+    q2m = p[1]
+    Efunc = p[3]
+    Bfunc = p[4]
+
+    # Evaluate fields at time t
+    # For a full time step, first we push velocity half step?
+    # Wait, the traditional Boris pushes v by half step, r by full step, v by half step.
+    # In TestParticle.jl, they do:
+    # v_new = velocity_updater(v, r, dt, t + dt, p) # actually it uses t + 0.5*dt somewhere
+
+    # We need to closely match what the Boris algorithm conceptually does.
+    # Actually, a standard SciML step goes from t to t+dt.
+    # If the user expects standard Boris, it is:
+    r_half = r + v * (dt / 2)
+    # Evaluate fields at t + dt/2, r_half
+    t_half = t + dt / 2
+    E = Efunc(r_half, t_half)
+    B = Bfunc(r_half, t_half)
+
+    qdt_2m = q2m * 0.5 * dt
+    v_new = boris_velocity_update(v, E, B, qdt_2m)
+
+    r_new = r_half + v_new * (dt / 2)
+
+    u_new = vcat(r_new, v_new)
+    return integrator.u = u_new
+end
+
+function OrdinaryDiffEqCore.initialize!(integrator, cache::BorisCache)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
+    integrator.stats.nf += 1
+    integrator.kshortsize = 0
+    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+end
+
+function OrdinaryDiffEqCore.perform_step!(integrator, cache::BorisCache, repeat_step = false)
+    # In-place version.
+    # Since TestParticle.jl mostly focuses on SVector for small state vectors,
+    # we can do a naive copy for the in-place cache.
+    t = integrator.t
+    dt = integrator.dt
+    uprev = integrator.uprev
+    f = integrator.f
+    p = integrator.p
+
+    r = SVector(uprev[1], uprev[2], uprev[3])
+    v = SVector(uprev[4], uprev[5], uprev[6])
+
+    q2m = p[1]
+    Efunc = p[3]
+    Bfunc = p[4]
+
+    r_half = r + v * (dt / 2)
+    t_half = t + dt / 2
+    E = Efunc(r_half, t_half)
+    B = Bfunc(r_half, t_half)
+
+    qdt_2m = q2m * 0.5 * dt
+    v_new = boris_velocity_update(v, E, B, qdt_2m)
+
+    r_new = r_half + v_new * (dt / 2)
+
+    integrator.u[1] = r_new[1]
+    integrator.u[2] = r_new[2]
+    integrator.u[3] = r_new[3]
+    integrator.u[4] = v_new[1]
+    integrator.u[5] = v_new[2]
+    return integrator.u[6] = v_new[3]
+end
+
+@muladd function update_velocity_multistep(v, r, dt, t, n::Int, N::Int, param)
+    q2m = param[1]
+    Efunc = param[3]
+    Bfunc = param[4]
+
+    E = Efunc(r, t)
+    B = Bfunc(r, t)
+
+    # t_n and e_n vectors
+    factor = q2m * dt / (2 * n)
+
+    t_n = factor * B # (q/m * dt/(2n)) * B
+    e_n = factor * E # (q/m * dt/(2n)) * E
+
+    # Hyper Boris N-th order gyrophase correction
+    if N != 2
+        t_mag2 = sum(abs2, t_n)
+        if N == 4
+            f_N = 1 + t_mag2 / 3
+            e_corr_factor = -1 / 3
+        else # N == 6
+            f_N = 1 + t_mag2 / 3 + 2 * t_mag2 * t_mag2 / 15
+            e_corr_factor = -1 / 3 - 2 * t_mag2 / 15
+        end
+
+        e_dot_t = e_n ⋅ t_n
+        e_n = f_N * e_n + (e_corr_factor * e_dot_t) * t_n
+        t_n = f_N * t_n
+    end
+
+    t_n_mag2 = sum(abs2, t_n)
+    t_n_mag = sqrt(t_n_mag2)
+
+    # Calculate coefficients
+    # Check for small t_n to avoid division by zero or precision loss
+    if t_n_mag < TN_MAG_THRESHOLD
+        # Taylor expansion limits as t_n -> 0
+        c_n1 = 1 - 2 * n * n * t_n_mag2
+
+        n_term1 = 2 * n
+        n_term3 = 4 * n * n * n
+
+        c_n2 = n_term1 - (n_term1 + n_term3) / 3 * t_n_mag2
+        c_n3 = 2 * n * n - (4 * n * n + 2 * n * n * n * n) / 3 * t_n_mag2
+        c_n6 = (n_term1 + n_term3) / 3
+    else
+        alpha_n = atan(t_n_mag)
+        n_alpha_n = n * alpha_n
+        sin_n_alpha, cos_n_alpha = sincos(n_alpha_n)
+        sin_2n_alpha = 2 * sin_n_alpha * cos_n_alpha
+        cos_2n_alpha = cos_n_alpha * cos_n_alpha - sin_n_alpha * sin_n_alpha
+
+        c_n1 = cos_2n_alpha
+        c_n2 = sin_2n_alpha / t_n_mag
+        c_n3 = 2 * sin_n_alpha * sin_n_alpha / t_n_mag2
+        c_n6 = (2 * n - c_n2) / t_n_mag2
+    end
+
+    c_n4 = c_n2
+    c_n5 = c_n3
+
+    v_dot_t = v ⋅ t_n
+    e_dot_t = e_n ⋅ t_n
+
+    v_cross_t = v × t_n
+    e_cross_t = e_n × t_n
+
+    # Update velocity
+    # Equation 39:
+    v_new = c_n1 * v +
+        c_n2 * v_cross_t +
+        c_n3 * v_dot_t * t_n +
+        c_n4 * e_n +
+        c_n5 * e_cross_t +
+        c_n6 * e_dot_t * t_n
+
+    return v_new
+end
+
+function OrdinaryDiffEqCore.initialize!(integrator, cache::MultistepBorisConstantCache)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
+    integrator.stats.nf += 1
+    integrator.kshortsize = 0
+    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+end
+
+function OrdinaryDiffEqCore.perform_step!(integrator, cache::MultistepBorisConstantCache, repeat_step = false)
+    t = integrator.t
+    dt = integrator.dt
+    uprev = integrator.uprev
+    p = integrator.p
+    alg = integrator.alg
+
+    r = uprev[SVector(1, 2, 3)]
+    v = uprev[SVector(4, 5, 6)]
+
+    # Half step update
+    r_half = r + v * (dt / 2)
+    t_half = t + dt / 2
+
+    # Update velocity using multistep
+    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, alg.N, p)
+
+    r_new = r_half + v_new * (dt / 2)
+
+    return integrator.u = vcat(r_new, v_new)
+end
+
+function OrdinaryDiffEqCore.initialize!(integrator, cache::MultistepBorisCache)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
+    integrator.stats.nf += 1
+    integrator.kshortsize = 0
+    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+end
+
+function OrdinaryDiffEqCore.perform_step!(integrator, cache::MultistepBorisCache, repeat_step = false)
+    t = integrator.t
+    dt = integrator.dt
+    uprev = integrator.uprev
+    p = integrator.p
+    alg = integrator.alg
+
+    r = SVector(uprev[1], uprev[2], uprev[3])
+    v = SVector(uprev[4], uprev[5], uprev[6])
+
+    r_half = r + v * (dt / 2)
+    t_half = t + dt / 2
+
+    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, alg.N, p)
+
+    r_new = r_half + v_new * (dt / 2)
+
+    integrator.u[1] = r_new[1]
+    integrator.u[2] = r_new[2]
+    integrator.u[3] = r_new[3]
+    integrator.u[4] = v_new[1]
+    integrator.u[5] = v_new[2]
+    return integrator.u[6] = v_new[3]
+end
+
+function OrdinaryDiffEqCore.initialize!(integrator, cache::AdaptiveBorisConstantCache)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
+    integrator.stats.nf += 1
+    integrator.kshortsize = 0
+    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+end
+
+function OrdinaryDiffEqCore.perform_step!(integrator, cache::AdaptiveBorisConstantCache, repeat_step = false)
+    t = integrator.t
+    dt = integrator.dt
+    uprev = integrator.uprev
+    p = integrator.p
+    alg = integrator.alg
+
+    r = uprev[SVector(1, 2, 3)]
+    v = uprev[SVector(4, 5, 6)]
+
+    q2m = p[1]
+    Efunc = p[3]
+    Bfunc = p[4]
+
+    r_half = r + v * (dt / 2)
+    t_half = t + dt / 2
+    E = Efunc(r_half, t_half)
+    B = Bfunc(r_half, t_half)
+
+    qdt_2m = q2m * 0.5 * dt
+    v_new = boris_velocity_update(v, E, B, qdt_2m)
+
+    r_new = r_half + v_new * (dt / 2)
+    integrator.u = vcat(r_new, v_new)
+
+    # Adaptive Step proposition based on local gyroperiod
+    Bmag = norm(Bfunc(r_new, t + dt))
+    dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
+    return set_proposed_dt!(integrator, dt_new)
+end
+
+function OrdinaryDiffEqCore.initialize!(integrator, cache::AdaptiveBorisCache)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
+    integrator.stats.nf += 1
+    integrator.kshortsize = 0
+    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+end
+
+function OrdinaryDiffEqCore.perform_step!(integrator, cache::AdaptiveBorisCache, repeat_step = false)
+    t = integrator.t
+    dt = integrator.dt
+    uprev = integrator.uprev
+    p = integrator.p
+    alg = integrator.alg
+
+    r = SVector(uprev[1], uprev[2], uprev[3])
+    v = SVector(uprev[4], uprev[5], uprev[6])
+
+    q2m = p[1]
+    Efunc = p[3]
+    Bfunc = p[4]
+
+    r_half = r + v * (dt / 2)
+    t_half = t + dt / 2
+    E = Efunc(r_half, t_half)
+    B = Bfunc(r_half, t_half)
+
+    qdt_2m = q2m * 0.5 * dt
+    v_new = boris_velocity_update(v, E, B, qdt_2m)
+
+    r_new = r_half + v_new * (dt / 2)
+
+    integrator.u[1] = r_new[1]
+    integrator.u[2] = r_new[2]
+    integrator.u[3] = r_new[3]
+    integrator.u[4] = v_new[1]
+    integrator.u[5] = v_new[2]
+    integrator.u[6] = v_new[3]
+
+    Bmag = norm(Bfunc(r_new, t + dt))
+    dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
+    return set_proposed_dt!(integrator, dt_new)
+end
diff --git a/lib/OrdinaryDiffEqBoris/test/runtests.jl b/lib/OrdinaryDiffEqBoris/test/runtests.jl
new file mode 100644
index 000000000..ec7e24175
--- /dev/null
+++ b/lib/OrdinaryDiffEqBoris/test/runtests.jl
@@ -0,0 +1,144 @@
+using Test
+using OrdinaryDiffEqBoris
+using SciMLBase
+using StaticArrays
+using LinearAlgebra: norm
+
+@testset "OrdinaryDiffEqBoris.jl" begin
+    # Definitions
+    constant_Ey(x, t) = SA[0.0, 1.0, 0.0]
+    constant_Bz(x, t) = SA[0.0, 0.0, 1.0]
+    ZeroField() = (x, t) -> SA[0.0, 0.0, 0.0]
+
+    dummy_f(u, p, t) = u # ODEProblem requires an f, but Boris uses p directly
+
+    @testset "E cross B drift - Standard Boris" begin
+        # E = (0, 1, 0), B = (0, 0, 1)
+        # Analytic solution: particle moves at constant velocity v_drift = (1, 0, 0)
+        param = (1.0, 1.0, constant_Ey, constant_Bz, ZeroField())
+        u0 = SA[0.0, 0.0, 0.0, 1.0, 0.0, 0.0]
+        tspan = (0.0, 10.0)
+        prob = ODEProblem(dummy_f, u0, tspan, param)
+
+        dt = 0.1
+        # Test standard Boris
+        sol_boris = solve(prob, Boris(); dt, adaptive = false)
+        @test sol_boris.u[end][1] ≈ 10.0 atol = 1.0e-6
+        @test sol_boris.u[end][4] ≈ 1.0 atol = 1.0e-6
+    end
+
+    @testset "E cross B drift - Multistep Boris" begin
+        param = (1.0, 1.0, constant_Ey, constant_Bz, ZeroField())
+        u0 = SA[0.0, 0.0, 0.0, 1.0, 0.0, 0.0]
+        tspan = (0.0, 10.0)
+        prob = ODEProblem(dummy_f, u0, tspan, param)
+
+        dt = 0.1
+        sol_multi_2 = solve(prob, MultistepBoris(n = 2, N = 2); dt, adaptive = false)
+        @test sol_multi_2.u[end][1] ≈ 10.0 atol = 1.0e-6
+        @test sol_multi_2.u[end][4] ≈ 1.0 atol = 1.0e-6
+
+        # Test Hyper Boris N=4
+        sol_hyper_4 = solve(prob, MultistepBoris(n = 2, N = 4); dt, adaptive = false)
+        @test sol_hyper_4.u[end][1] ≈ 10.0 atol = 1.0e-6
+    end
+
+    @testset "Gyrating particle" begin
+        param_gyro = (1.0, 1.0, ZeroField(), constant_Bz, ZeroField())
+        u0_gyro = SA[0.0, 0.0, 0.0, 1.0, 0.0, 0.0]
+        prob_gyro = ODEProblem(dummy_f, u0_gyro, (0.0, 2π), param_gyro)
+
+        # After one period, should return to origin
+        sol_1step_gyro = solve(prob_gyro, MultistepBoris(n = 1, N = 2); dt = 0.1, adaptive = false)
+        @test hypot(@views sol_1step_gyro.u[end][1:3]...) < 0.02
+
+        sol_hyper_6_gyro = solve(prob_gyro, MultistepBoris(n = 4, N = 6); dt = 0.1, adaptive = false)
+        @test hypot(@views sol_hyper_6_gyro.u[end][1:3]...) < 0.02
+    end
+
+    @testset "Adaptive Boris" begin
+        constant_E(x, t) = SA[0.0, 1.0e2, 0.0]
+        gradient_B(x, t) = SA[0.0, 0.0, 1.0 * (1.0 + x[1])]
+        param = (-1.0, 1.0, constant_E, gradient_B, ZeroField()) # q2m = -1
+
+        u0 = SA[0.0, 0.0, 0.0, 10.0, 0.0, 0.0]
+        tspan = (0.0, 10.0)
+        prob = ODEProblem(dummy_f, u0, tspan, param)
+
+        safety = 0.1
+        init_Bmag = norm(gradient_B(u0[1:3], 0.0))
+        dt_init = safety * 2π / (abs(param[1]) * init_Bmag)
+
+        sol = solve(prob, AdaptiveBoris(safety = safety); dt = dt_init, adaptive = false)
+
+        # Check if step sizes are not uniform (hence adapted)
+        dts = diff(sol.t)
+        @test !all(y -> isapprox(y, dts[1], rtol = 1.0e-5), dts)
+        @test sol.t[end] ≈ tspan[2]
+
+        # Energy conservation check
+        total_energy(u) = 0.5 * param[2] * norm(u[4:6])^2 + (param[1] * param[2]) * (-1.0e2 * u[2])
+        E_start = total_energy(sol.u[1])
+        E_end = total_energy(sol.u[end])
+        @test isapprox(E_end, E_start, rtol = 1.0e-3)
+    end
+
+    @testset "Nonzero tspan[1]" begin
+        # B field that is only present when t > 5
+        B_field(r, t) = t > 5 ? SA[0.0, 0.0, 0.01] : SA[0.0, 0.0, 0.0]
+        E_field(r, t) = SA[0.0, 0.0, 0.0]
+
+        # param = (q2m, m, E, B, F)
+        param = (-1.0e11, 1.0e-30, E_field, B_field, ZeroField())
+
+        # Start at t=10. If absolute time is used, B should be 0.01.
+        tspan = (10.0, 10.0 + 1.0e-7)
+        dt = 1.0e-9
+        u0 = SA[0.0, 0.0, 0.0, 1.0e5, 0.0, 0.0]
+        prob = ODEProblem(dummy_f, u0, tspan, param)
+
+        # Standard Boris
+        sol_boris = solve(prob, Boris(); dt = dt, adaptive = false)
+        # If B was 0.01, vx should have deviated from initial 1.0e5
+        @test abs(sol_boris.u[end][4] - 1.0e5) > 1.0e-4
+
+        # Adaptive Boris
+        sol_adaptive = solve(prob, AdaptiveBoris(safety = 0.1); dt = dt, adaptive = false)
+        @test abs(sol_adaptive.u[end][4] - 1.0e5) > 1.0e-4
+    end
+
+    @testset "SciML Output saving flags" begin
+        # Setup
+        x0 = [0.0, 0.0, 0.0]
+        v0 = [0.0, 1.0e5, 0.0]
+        stateinit = [x0..., v0...]
+        tspan = (0.0, 3.0e-8)
+        dt = 3.0e-11
+
+        zero_E(x, t) = SA[0.0, 0.0, 0.0]
+        uniform_B(x, t) = SA[0.0, 0.0, 0.01]
+        param = (-1.0e11, 1.0e-30, zero_E, uniform_B, zero_E)
+        prob = ODEProblem(dummy_f, stateinit, tspan, param)
+
+        # Baseline: save_everystep=true (default)
+        sol = solve(prob, Boris(); dt, adaptive = false)
+        @test length(sol.t) == 1001 # 3e-8 / 3e-11 = 1000 steps + start
+
+        # Scenario 2: Only final state
+        sol = solve(prob, Boris(); dt, save_everystep = false, save_start = false, save_on = false, adaptive = false)
+        @test length(sol.t) == 1
+        @test sol.t[1] ≈ tspan[2]
+
+        # Scenario 3: Start and End
+        sol = solve(prob, Boris(); dt, save_everystep = false, save_start = true, save_on = false, adaptive = false)
+        @test length(sol.t) == 2
+        @test sol.t[1] == tspan[1]
+        @test sol.t[end] ≈ tspan[2]
+
+        # Multistep Boris test flags
+        sol_ms = solve(prob, MultistepBoris(n = 2, N = 2); dt, save_everystep = false, save_start = true, save_on = false, adaptive = false)
+        @test length(sol_ms.t) == 2
+        @test sol_ms.t[1] == tspan[1]
+        @test sol_ms.t[end] ≈ tspan[2]
+    end
+end

From f484a9b3306dca30f546307ecb11291aae7f9c7f Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 10:43:56 -0400
Subject: [PATCH 02/17] Rexport SciMLBase; Follow the solver naming standard

---
 lib/OrdinaryDiffEqBoris/Project.toml          |  8 +++--
 .../src/OrdinaryDiffEqBoris.jl                | 10 ++----
 lib/OrdinaryDiffEqBoris/src/alg_utils.jl      | 12 +++----
 lib/OrdinaryDiffEqBoris/src/algorithms.jl     | 32 +++++++++++++++----
 lib/OrdinaryDiffEqBoris/src/boris_caches.jl   | 20 ++++++------
 .../src/boris_perform_step.jl                 | 30 +++++++++--------
 lib/OrdinaryDiffEqBoris/test/runtests.jl      | 11 +++----
 7 files changed, 71 insertions(+), 52 deletions(-)

diff --git a/lib/OrdinaryDiffEqBoris/Project.toml b/lib/OrdinaryDiffEqBoris/Project.toml
index a5521972a..486b7d568 100644
--- a/lib/OrdinaryDiffEqBoris/Project.toml
+++ b/lib/OrdinaryDiffEqBoris/Project.toml
@@ -3,19 +3,21 @@ uuid = "1e8b2ed2-ea1c-4b6e-b3f5-67c29d6bd197"
 version = "0.1.0"
 
 [deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
+MuladdMacro = "0.2"
 OrdinaryDiffEqCore = "3"
 RecursiveArrayTools = "3.24, 4"
+Reexport = "1.2.2"
 SciMLBase = "2, 3"
 StaticArrays = "1"
-MuladdMacro = "0.2"
 julia = "1.10"
 
 [extras]
diff --git a/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl b/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
index e2f0944f9..a3c661798 100644
--- a/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
+++ b/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
@@ -1,12 +1,10 @@
 module OrdinaryDiffEqBoris
 
-import OrdinaryDiffEqCore
+using Reexport
+@reexport using SciMLBase
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, OrdinaryDiffEqAdaptiveAlgorithm,
     OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
     alg_order, alg_cache, isfsal, initialize!, perform_step!
-
-
-using SciMLBase
 using RecursiveArrayTools
 using StaticArrays
 using MuladdMacro
@@ -18,9 +16,7 @@ include("boris_caches.jl")
 include("boris_perform_step.jl")
 
 export Boris
-export MultistepBoris
+export MultistepBoris2, MultistepBoris4, MultistepBoris6
 export AdaptiveBoris
-export boris_velocity_update # Exposing strictly for debugging if needed
-export update_velocity_multistep
 
 end
diff --git a/lib/OrdinaryDiffEqBoris/src/alg_utils.jl b/lib/OrdinaryDiffEqBoris/src/alg_utils.jl
index a8b8e9f0e..51ab170aa 100644
--- a/lib/OrdinaryDiffEqBoris/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqBoris/src/alg_utils.jl
@@ -1,8 +1,8 @@
-OrdinaryDiffEqCore.alg_order(alg::Boris) = 2
-OrdinaryDiffEqCore.isfsal(alg::Boris) = false
+alg_order(alg::Boris) = 2
+isfsal(alg::Boris) = false
 
-OrdinaryDiffEqCore.alg_order(alg::MultistepBoris) = 2
-OrdinaryDiffEqCore.isfsal(alg::MultistepBoris) = false
+alg_order(alg::MultistepBoris{N}) where {N} = 2
+isfsal(alg::MultistepBoris{N}) where {N} = false
 
-OrdinaryDiffEqCore.alg_order(alg::AdaptiveBoris) = 2
-OrdinaryDiffEqCore.isfsal(alg::AdaptiveBoris) = false
+alg_order(alg::AdaptiveBoris) = 2
+isfsal(alg::AdaptiveBoris) = false
diff --git a/lib/OrdinaryDiffEqBoris/src/algorithms.jl b/lib/OrdinaryDiffEqBoris/src/algorithms.jl
index 1ec0a4271..2fc362baf 100644
--- a/lib/OrdinaryDiffEqBoris/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqBoris/src/algorithms.jl
@@ -8,17 +8,37 @@ This solver expects a problem where `p` is structured as `(q2m, m, E, B, ...)`,
 struct Boris <: OrdinaryDiffEqAlgorithm end
 
 """
-    MultistepBoris(; n=1, N=2)
+    MultistepBoris{N}(; n=1)
 
-The Multistep/Hyper Boris method.
+The Multistep/Hyper Boris method of order `N`.
 `n` specifies the number of subcycles.
-`N` specifies the gyrophase correction order. `N=2` corresponds to the Multicycle solver, while `N=4` or `N=6` is the Hyper Boris solver.
+`N` specifies the gyrophase correction order. `N=2` corresponds to the Multicycle solver, while `N=4` or `N=6` are the Hyper Boris solvers.
 """
-struct MultistepBoris <: OrdinaryDiffEqAlgorithm
+struct MultistepBoris{N} <: OrdinaryDiffEqAlgorithm
     n::Int
-    N::Int
 end
-MultistepBoris(; n::Int = 1, N::Int = 2) = MultistepBoris(n, N)
+MultistepBoris{N}(; n::Int = 1) where {N} = MultistepBoris{N}(n)
+
+"""
+    MultistepBoris2(; n=1)
+
+The Multicycle Boris method (MultistepBoris with N=2).
+"""
+const MultistepBoris2 = MultistepBoris{2}
+
+"""
+    MultistepBoris4(; n=1)
+
+The 4th order Hyper Boris method (MultistepBoris with N=4).
+"""
+const MultistepBoris4 = MultistepBoris{4}
+
+"""
+    MultistepBoris6(; n=1)
+
+The 6th order Hyper Boris method (MultistepBoris with N=6).
+"""
+const MultistepBoris6 = MultistepBoris{6}
 
 """
     AdaptiveBoris(; safety=0.1)
diff --git a/lib/OrdinaryDiffEqBoris/src/boris_caches.jl b/lib/OrdinaryDiffEqBoris/src/boris_caches.jl
index 2e0ec7758..9578b3000 100644
--- a/lib/OrdinaryDiffEqBoris/src/boris_caches.jl
+++ b/lib/OrdinaryDiffEqBoris/src/boris_caches.jl
@@ -7,7 +7,7 @@ struct BorisCache{uType, rateType} <: OrdinaryDiffEqMutableCache
     k::rateType
 end
 
-function OrdinaryDiffEqCore.alg_cache(
+function alg_cache(
         alg::Boris, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
         ::Val{false}, args...; kwargs...
@@ -15,7 +15,7 @@ function OrdinaryDiffEqCore.alg_cache(
     return BorisConstantCache()
 end
 
-function OrdinaryDiffEqCore.alg_cache(
+function alg_cache(
         alg::Boris, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
         ::Val{true}, args...; kwargs...
@@ -32,19 +32,19 @@ struct MultistepBorisCache{uType, rateType} <: OrdinaryDiffEqMutableCache
     k::rateType
 end
 
-function OrdinaryDiffEqCore.alg_cache(
-        alg::MultistepBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
+function alg_cache(
+        alg::MultistepBoris{N}, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
         ::Val{false}, args...; kwargs...
-    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, N}
     return MultistepBorisConstantCache()
 end
 
-function OrdinaryDiffEqCore.alg_cache(
-        alg::MultistepBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
+function alg_cache(
+        alg::MultistepBoris{N}, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
         ::Val{true}, args...; kwargs...
-    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, N}
     return MultistepBorisCache(u, uprev, similar(u), similar(rate_prototype))
 end
 
@@ -57,7 +57,7 @@ struct AdaptiveBorisCache{uType, rateType} <: OrdinaryDiffEqMutableCache
     k::rateType
 end
 
-function OrdinaryDiffEqCore.alg_cache(
+function alg_cache(
         alg::AdaptiveBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
         ::Val{false}, args...; kwargs...
@@ -65,7 +65,7 @@ function OrdinaryDiffEqCore.alg_cache(
     return AdaptiveBorisConstantCache()
 end
 
-function OrdinaryDiffEqCore.alg_cache(
+function alg_cache(
         alg::AdaptiveBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
         ::Val{true}, args...; kwargs...
diff --git a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
index 0a0d9b499..200201b8a 100644
--- a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
+++ b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
@@ -14,14 +14,14 @@ const TN_MAG_THRESHOLD = 1.0e-4
     return v_new
 end
 
-function OrdinaryDiffEqCore.initialize!(integrator, cache::BorisConstantCache)
+function initialize!(integrator, cache::BorisConstantCache)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
     integrator.stats.nf += 1
     integrator.kshortsize = 0
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function OrdinaryDiffEqCore.perform_step!(integrator, cache::BorisConstantCache, repeat_step = false)
+function perform_step!(integrator, cache::BorisConstantCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
@@ -62,14 +62,14 @@ function OrdinaryDiffEqCore.perform_step!(integrator, cache::BorisConstantCache,
     return integrator.u = u_new
 end
 
-function OrdinaryDiffEqCore.initialize!(integrator, cache::BorisCache)
+function initialize!(integrator, cache::BorisCache)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
     integrator.stats.nf += 1
     integrator.kshortsize = 0
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function OrdinaryDiffEqCore.perform_step!(integrator, cache::BorisCache, repeat_step = false)
+function perform_step!(integrator, cache::BorisCache, repeat_step = false)
     # In-place version.
     # Since TestParticle.jl mostly focuses on SVector for small state vectors,
     # we can do a naive copy for the in-place cache.
@@ -183,14 +183,14 @@ end
     return v_new
 end
 
-function OrdinaryDiffEqCore.initialize!(integrator, cache::MultistepBorisConstantCache)
+function initialize!(integrator, cache::MultistepBorisConstantCache)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
     integrator.stats.nf += 1
     integrator.kshortsize = 0
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function OrdinaryDiffEqCore.perform_step!(integrator, cache::MultistepBorisConstantCache, repeat_step = false)
+function perform_step!(integrator, cache::MultistepBorisConstantCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
@@ -205,21 +205,22 @@ function OrdinaryDiffEqCore.perform_step!(integrator, cache::MultistepBorisConst
     t_half = t + dt / 2
 
     # Update velocity using multistep
-    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, alg.N, p)
+    N = typeof(alg).parameters[1]
+    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, N, p)
 
     r_new = r_half + v_new * (dt / 2)
 
     return integrator.u = vcat(r_new, v_new)
 end
 
-function OrdinaryDiffEqCore.initialize!(integrator, cache::MultistepBorisCache)
+function initialize!(integrator, cache::MultistepBorisCache)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
     integrator.stats.nf += 1
     integrator.kshortsize = 0
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function OrdinaryDiffEqCore.perform_step!(integrator, cache::MultistepBorisCache, repeat_step = false)
+function perform_step!(integrator, cache::MultistepBorisCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
@@ -232,7 +233,8 @@ function OrdinaryDiffEqCore.perform_step!(integrator, cache::MultistepBorisCache
     r_half = r + v * (dt / 2)
     t_half = t + dt / 2
 
-    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, alg.N, p)
+    N = typeof(alg).parameters[1]
+    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, N, p)
 
     r_new = r_half + v_new * (dt / 2)
 
@@ -244,14 +246,14 @@ function OrdinaryDiffEqCore.perform_step!(integrator, cache::MultistepBorisCache
     return integrator.u[6] = v_new[3]
 end
 
-function OrdinaryDiffEqCore.initialize!(integrator, cache::AdaptiveBorisConstantCache)
+function initialize!(integrator, cache::AdaptiveBorisConstantCache)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
     integrator.stats.nf += 1
     integrator.kshortsize = 0
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function OrdinaryDiffEqCore.perform_step!(integrator, cache::AdaptiveBorisConstantCache, repeat_step = false)
+function perform_step!(integrator, cache::AdaptiveBorisConstantCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
@@ -282,14 +284,14 @@ function OrdinaryDiffEqCore.perform_step!(integrator, cache::AdaptiveBorisConsta
     return set_proposed_dt!(integrator, dt_new)
 end
 
-function OrdinaryDiffEqCore.initialize!(integrator, cache::AdaptiveBorisCache)
+function initialize!(integrator, cache::AdaptiveBorisCache)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
     integrator.stats.nf += 1
     integrator.kshortsize = 0
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function OrdinaryDiffEqCore.perform_step!(integrator, cache::AdaptiveBorisCache, repeat_step = false)
+function perform_step!(integrator, cache::AdaptiveBorisCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
diff --git a/lib/OrdinaryDiffEqBoris/test/runtests.jl b/lib/OrdinaryDiffEqBoris/test/runtests.jl
index ec7e24175..3d5947389 100644
--- a/lib/OrdinaryDiffEqBoris/test/runtests.jl
+++ b/lib/OrdinaryDiffEqBoris/test/runtests.jl
@@ -1,6 +1,5 @@
 using Test
 using OrdinaryDiffEqBoris
-using SciMLBase
 using StaticArrays
 using LinearAlgebra: norm
 
@@ -34,12 +33,12 @@ using LinearAlgebra: norm
         prob = ODEProblem(dummy_f, u0, tspan, param)
 
         dt = 0.1
-        sol_multi_2 = solve(prob, MultistepBoris(n = 2, N = 2); dt, adaptive = false)
+        sol_multi_2 = solve(prob, MultistepBoris2(n = 2); dt, adaptive = false)
         @test sol_multi_2.u[end][1] ≈ 10.0 atol = 1.0e-6
         @test sol_multi_2.u[end][4] ≈ 1.0 atol = 1.0e-6
 
         # Test Hyper Boris N=4
-        sol_hyper_4 = solve(prob, MultistepBoris(n = 2, N = 4); dt, adaptive = false)
+        sol_hyper_4 = solve(prob, MultistepBoris4(n = 2); dt, adaptive = false)
         @test sol_hyper_4.u[end][1] ≈ 10.0 atol = 1.0e-6
     end
 
@@ -49,10 +48,10 @@ using LinearAlgebra: norm
         prob_gyro = ODEProblem(dummy_f, u0_gyro, (0.0, 2π), param_gyro)
 
         # After one period, should return to origin
-        sol_1step_gyro = solve(prob_gyro, MultistepBoris(n = 1, N = 2); dt = 0.1, adaptive = false)
+        sol_1step_gyro = solve(prob_gyro, MultistepBoris2(n = 1); dt = 0.1, adaptive = false)
         @test hypot(@views sol_1step_gyro.u[end][1:3]...) < 0.02
 
-        sol_hyper_6_gyro = solve(prob_gyro, MultistepBoris(n = 4, N = 6); dt = 0.1, adaptive = false)
+        sol_hyper_6_gyro = solve(prob_gyro, MultistepBoris6(n = 4); dt = 0.1, adaptive = false)
         @test hypot(@views sol_hyper_6_gyro.u[end][1:3]...) < 0.02
     end
 
@@ -136,7 +135,7 @@ using LinearAlgebra: norm
         @test sol.t[end] ≈ tspan[2]
 
         # Multistep Boris test flags
-        sol_ms = solve(prob, MultistepBoris(n = 2, N = 2); dt, save_everystep = false, save_start = true, save_on = false, adaptive = false)
+        sol_ms = solve(prob, MultistepBoris2(n = 2); dt, save_everystep = false, save_start = true, save_on = false, adaptive = false)
         @test length(sol_ms.t) == 2
         @test sol_ms.t[1] == tspan[1]
         @test sol_ms.t[end] ≈ tspan[2]

From 3671ca98cf1be3eefb30c3f4ff5a03b643b6b3ab Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 12:24:50 -0400
Subject: [PATCH 03/17] Unify adaptive and fixed timestep implementations

---
 .../src/OrdinaryDiffEqBoris.jl                |   1 -
 lib/OrdinaryDiffEqBoris/src/alg_utils.jl      |   3 -
 lib/OrdinaryDiffEqBoris/src/algorithms.jl     |  21 +--
 lib/OrdinaryDiffEqBoris/src/boris_caches.jl   |  25 ---
 .../src/boris_perform_step.jl                 | 155 +++++-------------
 lib/OrdinaryDiffEqBoris/test/runtests.jl      |  24 +--
 6 files changed, 56 insertions(+), 173 deletions(-)

diff --git a/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl b/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
index a3c661798..94bd1e1c8 100644
--- a/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
+++ b/lib/OrdinaryDiffEqBoris/src/OrdinaryDiffEqBoris.jl
@@ -17,6 +17,5 @@ include("boris_perform_step.jl")
 
 export Boris
 export MultistepBoris2, MultistepBoris4, MultistepBoris6
-export AdaptiveBoris
 
 end
diff --git a/lib/OrdinaryDiffEqBoris/src/alg_utils.jl b/lib/OrdinaryDiffEqBoris/src/alg_utils.jl
index 51ab170aa..932c26962 100644
--- a/lib/OrdinaryDiffEqBoris/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqBoris/src/alg_utils.jl
@@ -3,6 +3,3 @@ isfsal(alg::Boris) = false
 
 alg_order(alg::MultistepBoris{N}) where {N} = 2
 isfsal(alg::MultistepBoris{N}) where {N} = false
-
-alg_order(alg::AdaptiveBoris) = 2
-isfsal(alg::AdaptiveBoris) = false
diff --git a/lib/OrdinaryDiffEqBoris/src/algorithms.jl b/lib/OrdinaryDiffEqBoris/src/algorithms.jl
index 2fc362baf..c6c4c4138 100644
--- a/lib/OrdinaryDiffEqBoris/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqBoris/src/algorithms.jl
@@ -5,7 +5,10 @@ The standard Boris method for particle pushing in electric and magnetic fields.
 
 This solver expects a problem where `p` is structured as `(q2m, m, E, B, ...)`, which matches the signature used by TestParticle.jl's `TraceProblem`.
 """
-struct Boris <: OrdinaryDiffEqAlgorithm end
+struct Boris{T} <: OrdinaryDiffEqAlgorithm
+    safety::T
+end
+Boris(; safety = 0.0) = Boris(safety)
 
 """
     MultistepBoris{N}(; n=1)
@@ -14,10 +17,11 @@ The Multistep/Hyper Boris method of order `N`.
 `n` specifies the number of subcycles.
 `N` specifies the gyrophase correction order. `N=2` corresponds to the Multicycle solver, while `N=4` or `N=6` are the Hyper Boris solvers.
 """
-struct MultistepBoris{N} <: OrdinaryDiffEqAlgorithm
+struct MultistepBoris{N, T} <: OrdinaryDiffEqAlgorithm
     n::Int
+    safety::T
 end
-MultistepBoris{N}(; n::Int = 1) where {N} = MultistepBoris{N}(n)
+MultistepBoris{N}(; n::Int = 1, safety = 0.0) where {N} = MultistepBoris{N, typeof(safety)}(n, safety)
 
 """
     MultistepBoris2(; n=1)
@@ -39,14 +43,3 @@ const MultistepBoris4 = MultistepBoris{4}
 The 6th order Hyper Boris method (MultistepBoris with N=6).
 """
 const MultistepBoris6 = MultistepBoris{6}
-
-"""
-    AdaptiveBoris(; safety=0.1)
-
-Adaptive Boris method with adaptive time stepping based on the local gyroperiod.
-The time step is evaluated as `dt = safety * 2π / |qB/m|`.
-"""
-struct AdaptiveBoris{T} <: OrdinaryDiffEqAlgorithm
-    safety::T
-end
-AdaptiveBoris(; safety = 0.1) = AdaptiveBoris(safety)
diff --git a/lib/OrdinaryDiffEqBoris/src/boris_caches.jl b/lib/OrdinaryDiffEqBoris/src/boris_caches.jl
index 9578b3000..a28d40aee 100644
--- a/lib/OrdinaryDiffEqBoris/src/boris_caches.jl
+++ b/lib/OrdinaryDiffEqBoris/src/boris_caches.jl
@@ -47,28 +47,3 @@ function alg_cache(
     ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, N}
     return MultistepBorisCache(u, uprev, similar(u), similar(rate_prototype))
 end
-
-struct AdaptiveBorisConstantCache <: OrdinaryDiffEqConstantCache end
-
-struct AdaptiveBorisCache{uType, rateType} <: OrdinaryDiffEqMutableCache
-    u::uType
-    uprev::uType
-    tmp::uType
-    k::rateType
-end
-
-function alg_cache(
-        alg::AdaptiveBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
-        ::Val{false}, args...; kwargs...
-    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    return AdaptiveBorisConstantCache()
-end
-
-function alg_cache(
-        alg::AdaptiveBoris, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck,
-        ::Val{true}, args...; kwargs...
-    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    return AdaptiveBorisCache(u, uprev, similar(u), similar(rate_prototype))
-end
diff --git a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
index 200201b8a..d707a329c 100644
--- a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
+++ b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
@@ -21,34 +21,21 @@ function initialize!(integrator, cache::BorisConstantCache)
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function perform_step!(integrator, cache::BorisConstantCache, repeat_step = false)
+@muladd function perform_step!(integrator, cache::BorisConstantCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
     f = integrator.f
     p = integrator.p
 
-    # Extract particles from uprev
-    # TestParticle.jl standard: u[1:3] is r, u[4:6] is v
     r = uprev[SVector(1, 2, 3)]
     v = uprev[SVector(4, 5, 6)]
 
-    # Expected p structure in TestParticle.jl: (q2m, m, Efunc, Bfunc, Ffunc)
     q2m = p[1]
     Efunc = p[3]
     Bfunc = p[4]
 
-    # Evaluate fields at time t
-    # For a full time step, first we push velocity half step?
-    # Wait, the traditional Boris pushes v by half step, r by full step, v by half step.
-    # In TestParticle.jl, they do:
-    # v_new = velocity_updater(v, r, dt, t + dt, p) # actually it uses t + 0.5*dt somewhere
-
-    # We need to closely match what the Boris algorithm conceptually does.
-    # Actually, a standard SciML step goes from t to t+dt.
-    # If the user expects standard Boris, it is:
     r_half = r + v * (dt / 2)
-    # Evaluate fields at t + dt/2, r_half
     t_half = t + dt / 2
     E = Efunc(r_half, t_half)
     B = Bfunc(r_half, t_half)
@@ -58,8 +45,15 @@ function perform_step!(integrator, cache::BorisConstantCache, repeat_step = fals
 
     r_new = r_half + v_new * (dt / 2)
 
-    u_new = vcat(r_new, v_new)
-    return integrator.u = u_new
+    integrator.u = vcat(r_new, v_new)
+
+    alg = integrator.alg
+    if alg.safety > 0.0
+        Bmag = norm(Bfunc(r_new, t + dt))
+        dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
+        set_proposed_dt!(integrator, dt_new)
+    end
+    return integrator.u
 end
 
 function initialize!(integrator, cache::BorisCache)
@@ -69,10 +63,7 @@ function initialize!(integrator, cache::BorisCache)
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function perform_step!(integrator, cache::BorisCache, repeat_step = false)
-    # In-place version.
-    # Since TestParticle.jl mostly focuses on SVector for small state vectors,
-    # we can do a naive copy for the in-place cache.
+@muladd function perform_step!(integrator, cache::BorisCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
@@ -104,7 +95,9 @@ function perform_step!(integrator, cache::BorisCache, repeat_step = false)
     return integrator.u[6] = v_new[3]
 end
 
-@muladd function update_velocity_multistep(v, r, dt, t, n::Int, N::Int, param)
+_get_val_N(::MultistepBoris{N}) where {N} = Val{N}()
+
+@muladd function update_velocity_multistep(v, r, dt, t, n::Int, ::Val{N}, param) where {N}
     q2m = param[1]
     Efunc = param[3]
     Bfunc = param[4]
@@ -112,13 +105,11 @@ end
     E = Efunc(r, t)
     B = Bfunc(r, t)
 
-    # t_n and e_n vectors
     factor = q2m * dt / (2 * n)
 
-    t_n = factor * B # (q/m * dt/(2n)) * B
-    e_n = factor * E # (q/m * dt/(2n)) * E
+    t_n = factor * B
+    e_n = factor * E
 
-    # Hyper Boris N-th order gyrophase correction
     if N != 2
         t_mag2 = sum(abs2, t_n)
         if N == 4
@@ -137,10 +128,7 @@ end
     t_n_mag2 = sum(abs2, t_n)
     t_n_mag = sqrt(t_n_mag2)
 
-    # Calculate coefficients
-    # Check for small t_n to avoid division by zero or precision loss
     if t_n_mag < TN_MAG_THRESHOLD
-        # Taylor expansion limits as t_n -> 0
         c_n1 = 1 - 2 * n * n * t_n_mag2
 
         n_term1 = 2 * n
@@ -171,8 +159,6 @@ end
     v_cross_t = v × t_n
     e_cross_t = e_n × t_n
 
-    # Update velocity
-    # Equation 39:
     v_new = c_n1 * v +
         c_n2 * v_cross_t +
         c_n3 * v_dot_t * t_n +
@@ -190,7 +176,7 @@ function initialize!(integrator, cache::MultistepBorisConstantCache)
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function perform_step!(integrator, cache::MultistepBorisConstantCache, repeat_step = false)
+@muladd function perform_step!(integrator, cache::MultistepBorisConstantCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
@@ -200,98 +186,33 @@ function perform_step!(integrator, cache::MultistepBorisConstantCache, repeat_st
     r = uprev[SVector(1, 2, 3)]
     v = uprev[SVector(4, 5, 6)]
 
-    # Half step update
     r_half = r + v * (dt / 2)
     t_half = t + dt / 2
 
-    # Update velocity using multistep
-    N = typeof(alg).parameters[1]
-    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, N, p)
+    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, _get_val_N(alg), p)
 
     r_new = r_half + v_new * (dt / 2)
 
-    return integrator.u = vcat(r_new, v_new)
-end
-
-function initialize!(integrator, cache::MultistepBorisCache)
-    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
-    integrator.stats.nf += 1
-    integrator.kshortsize = 0
-    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
-end
-
-function perform_step!(integrator, cache::MultistepBorisCache, repeat_step = false)
-    t = integrator.t
-    dt = integrator.dt
-    uprev = integrator.uprev
-    p = integrator.p
-    alg = integrator.alg
-
-    r = SVector(uprev[1], uprev[2], uprev[3])
-    v = SVector(uprev[4], uprev[5], uprev[6])
-
-    r_half = r + v * (dt / 2)
-    t_half = t + dt / 2
-
-    N = typeof(alg).parameters[1]
-    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, N, p)
-
-    r_new = r_half + v_new * (dt / 2)
-
-    integrator.u[1] = r_new[1]
-    integrator.u[2] = r_new[2]
-    integrator.u[3] = r_new[3]
-    integrator.u[4] = v_new[1]
-    integrator.u[5] = v_new[2]
-    return integrator.u[6] = v_new[3]
-end
-
-function initialize!(integrator, cache::AdaptiveBorisConstantCache)
-    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
-    integrator.stats.nf += 1
-    integrator.kshortsize = 0
-    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
-end
-
-function perform_step!(integrator, cache::AdaptiveBorisConstantCache, repeat_step = false)
-    t = integrator.t
-    dt = integrator.dt
-    uprev = integrator.uprev
-    p = integrator.p
-    alg = integrator.alg
-
-    r = uprev[SVector(1, 2, 3)]
-    v = uprev[SVector(4, 5, 6)]
-
-    q2m = p[1]
-    Efunc = p[3]
-    Bfunc = p[4]
-
-    r_half = r + v * (dt / 2)
-    t_half = t + dt / 2
-    E = Efunc(r_half, t_half)
-    B = Bfunc(r_half, t_half)
-
-    qdt_2m = q2m * 0.5 * dt
-    v_new = boris_velocity_update(v, E, B, qdt_2m)
-
-    r_new = r_half + v_new * (dt / 2)
     integrator.u = vcat(r_new, v_new)
 
-    # Adaptive Step proposition based on local gyroperiod
-    Bmag = norm(Bfunc(r_new, t + dt))
-    dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
-    return set_proposed_dt!(integrator, dt_new)
+    if alg.safety > 0.0
+        q2m = p[1]
+        Bfunc = p[4]
+        Bmag = norm(Bfunc(r_new, t + dt))
+        dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
+        set_proposed_dt!(integrator, dt_new)
+    end
+    return integrator.u
 end
 
-function initialize!(integrator, cache::AdaptiveBorisCache)
+function initialize!(integrator, cache::MultistepBorisCache)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
     integrator.stats.nf += 1
     integrator.kshortsize = 0
     return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
 end
 
-function perform_step!(integrator, cache::AdaptiveBorisCache, repeat_step = false)
+@muladd function perform_step!(integrator, cache::MultistepBorisCache, repeat_step = false)
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
@@ -301,17 +222,10 @@ function perform_step!(integrator, cache::AdaptiveBorisCache, repeat_step = fals
     r = SVector(uprev[1], uprev[2], uprev[3])
     v = SVector(uprev[4], uprev[5], uprev[6])
 
-    q2m = p[1]
-    Efunc = p[3]
-    Bfunc = p[4]
-
     r_half = r + v * (dt / 2)
     t_half = t + dt / 2
-    E = Efunc(r_half, t_half)
-    B = Bfunc(r_half, t_half)
 
-    qdt_2m = q2m * 0.5 * dt
-    v_new = boris_velocity_update(v, E, B, qdt_2m)
+    v_new = update_velocity_multistep(v, r_half, dt, t_half, alg.n, _get_val_N(alg), p)
 
     r_new = r_half + v_new * (dt / 2)
 
@@ -322,7 +236,12 @@ function perform_step!(integrator, cache::AdaptiveBorisCache, repeat_step = fals
     integrator.u[5] = v_new[2]
     integrator.u[6] = v_new[3]
 
-    Bmag = norm(Bfunc(r_new, t + dt))
-    dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
-    return set_proposed_dt!(integrator, dt_new)
+    if alg.safety > 0.0
+        q2m = p[1]
+        Bfunc = p[4]
+        Bmag = norm(Bfunc(r_new, t + dt))
+        dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
+        set_proposed_dt!(integrator, dt_new)
+    end
+    return integrator.u[6]
 end
diff --git a/lib/OrdinaryDiffEqBoris/test/runtests.jl b/lib/OrdinaryDiffEqBoris/test/runtests.jl
index 3d5947389..5f439f664 100644
--- a/lib/OrdinaryDiffEqBoris/test/runtests.jl
+++ b/lib/OrdinaryDiffEqBoris/test/runtests.jl
@@ -21,7 +21,7 @@ using LinearAlgebra: norm
 
         dt = 0.1
         # Test standard Boris
-        sol_boris = solve(prob, Boris(); dt, adaptive = false)
+        sol_boris = solve(prob, Boris(); dt)
         @test sol_boris.u[end][1] ≈ 10.0 atol = 1.0e-6
         @test sol_boris.u[end][4] ≈ 1.0 atol = 1.0e-6
     end
@@ -33,12 +33,12 @@ using LinearAlgebra: norm
         prob = ODEProblem(dummy_f, u0, tspan, param)
 
         dt = 0.1
-        sol_multi_2 = solve(prob, MultistepBoris2(n = 2); dt, adaptive = false)
+        sol_multi_2 = solve(prob, MultistepBoris2(n = 2); dt)
         @test sol_multi_2.u[end][1] ≈ 10.0 atol = 1.0e-6
         @test sol_multi_2.u[end][4] ≈ 1.0 atol = 1.0e-6
 
         # Test Hyper Boris N=4
-        sol_hyper_4 = solve(prob, MultistepBoris4(n = 2); dt, adaptive = false)
+        sol_hyper_4 = solve(prob, MultistepBoris4(n = 2); dt)
         @test sol_hyper_4.u[end][1] ≈ 10.0 atol = 1.0e-6
     end
 
@@ -48,10 +48,10 @@ using LinearAlgebra: norm
         prob_gyro = ODEProblem(dummy_f, u0_gyro, (0.0, 2π), param_gyro)
 
         # After one period, should return to origin
-        sol_1step_gyro = solve(prob_gyro, MultistepBoris2(n = 1); dt = 0.1, adaptive = false)
+        sol_1step_gyro = solve(prob_gyro, MultistepBoris2(n = 1); dt = 0.1)
         @test hypot(@views sol_1step_gyro.u[end][1:3]...) < 0.02
 
-        sol_hyper_6_gyro = solve(prob_gyro, MultistepBoris6(n = 4); dt = 0.1, adaptive = false)
+        sol_hyper_6_gyro = solve(prob_gyro, MultistepBoris6(n = 4); dt = 0.1)
         @test hypot(@views sol_hyper_6_gyro.u[end][1:3]...) < 0.02
     end
 
@@ -68,7 +68,7 @@ using LinearAlgebra: norm
         init_Bmag = norm(gradient_B(u0[1:3], 0.0))
         dt_init = safety * 2π / (abs(param[1]) * init_Bmag)
 
-        sol = solve(prob, AdaptiveBoris(safety = safety); dt = dt_init, adaptive = false)
+        sol = solve(prob, Boris(safety = safety); dt = dt_init)
 
         # Check if step sizes are not uniform (hence adapted)
         dts = diff(sol.t)
@@ -97,12 +97,12 @@ using LinearAlgebra: norm
         prob = ODEProblem(dummy_f, u0, tspan, param)
 
         # Standard Boris
-        sol_boris = solve(prob, Boris(); dt = dt, adaptive = false)
+        sol_boris = solve(prob, Boris(); dt = dt)
         # If B was 0.01, vx should have deviated from initial 1.0e5
         @test abs(sol_boris.u[end][4] - 1.0e5) > 1.0e-4
 
         # Adaptive Boris
-        sol_adaptive = solve(prob, AdaptiveBoris(safety = 0.1); dt = dt, adaptive = false)
+        sol_adaptive = solve(prob, Boris(safety = 0.1); dt = dt)
         @test abs(sol_adaptive.u[end][4] - 1.0e5) > 1.0e-4
     end
 
@@ -120,22 +120,22 @@ using LinearAlgebra: norm
         prob = ODEProblem(dummy_f, stateinit, tspan, param)
 
         # Baseline: save_everystep=true (default)
-        sol = solve(prob, Boris(); dt, adaptive = false)
+        sol = solve(prob, Boris(); dt)
         @test length(sol.t) == 1001 # 3e-8 / 3e-11 = 1000 steps + start
 
         # Scenario 2: Only final state
-        sol = solve(prob, Boris(); dt, save_everystep = false, save_start = false, save_on = false, adaptive = false)
+        sol = solve(prob, Boris(); dt, save_everystep = false, save_start = false, save_on = false)
         @test length(sol.t) == 1
         @test sol.t[1] ≈ tspan[2]
 
         # Scenario 3: Start and End
-        sol = solve(prob, Boris(); dt, save_everystep = false, save_start = true, save_on = false, adaptive = false)
+        sol = solve(prob, Boris(); dt, save_everystep = false, save_start = true, save_on = false)
         @test length(sol.t) == 2
         @test sol.t[1] == tspan[1]
         @test sol.t[end] ≈ tspan[2]
 
         # Multistep Boris test flags
-        sol_ms = solve(prob, MultistepBoris2(n = 2); dt, save_everystep = false, save_start = true, save_on = false, adaptive = false)
+        sol_ms = solve(prob, MultistepBoris2(n = 2); dt, save_everystep = false, save_start = true, save_on = false)
         @test length(sol_ms.t) == 2
         @test sol_ms.t[1] == tspan[1]
         @test sol_ms.t[end] ≈ tspan[2]

From f49a140af3832d08c5dc87f0125575b03254b2a5 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 12:31:03 -0400
Subject: [PATCH 04/17] Format return

---
 lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
index d707a329c..8e50b1459 100644
--- a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
+++ b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
@@ -92,7 +92,8 @@ end
     integrator.u[3] = r_new[3]
     integrator.u[4] = v_new[1]
     integrator.u[5] = v_new[2]
-    return integrator.u[6] = v_new[3]
+    integrator.u[6] = v_new[3]
+    return
 end
 
 _get_val_N(::MultistepBoris{N}) where {N} = Val{N}()
@@ -209,7 +210,8 @@ function initialize!(integrator, cache::MultistepBorisCache)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
     integrator.stats.nf += 1
     integrator.kshortsize = 0
-    return integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+    return
 end
 
 @muladd function perform_step!(integrator, cache::MultistepBorisCache, repeat_step = false)
@@ -243,5 +245,5 @@ end
         dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
         set_proposed_dt!(integrator, dt_new)
     end
-    return integrator.u[6]
+    return
 end

From 5069434ffcb6e39719c37516d1262338313cbee0 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 12:34:14 -0400
Subject: [PATCH 05/17] Remove unused code

---
 lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
index 8e50b1459..dd5bf0740 100644
--- a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
+++ b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
@@ -25,7 +25,6 @@ end
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
-    f = integrator.f
     p = integrator.p
 
     r = uprev[SVector(1, 2, 3)]
@@ -67,7 +66,6 @@ end
     t = integrator.t
     dt = integrator.dt
     uprev = integrator.uprev
-    f = integrator.f
     p = integrator.p
 
     r = SVector(uprev[1], uprev[2], uprev[3])
@@ -99,9 +97,7 @@ end
 _get_val_N(::MultistepBoris{N}) where {N} = Val{N}()
 
 @muladd function update_velocity_multistep(v, r, dt, t, n::Int, ::Val{N}, param) where {N}
-    q2m = param[1]
-    Efunc = param[3]
-    Bfunc = param[4]
+    q2m, Efunc, Bfunc = param[1], param[3], param[4]
 
     E = Efunc(r, t)
     B = Bfunc(r, t)

From f103cdb81587d111bf1a720be6a1d1993b259bad Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 17:37:20 -0400
Subject: [PATCH 06/17] refactor: merge Boris and AdaptiveBoris solvers

---
 Project.toml          |   1 -
 src/TestParticle.jl   |   3 +-
 src/adaptive_boris.jl | 303 ----------------------------
 src/boris.jl          | 456 +++++++++++++++++++++---------------------
 src/precompile.jl     |   2 +-
 test/test_boris.jl    |  14 +-
 test/test_boundary.jl |   8 +-
 7 files changed, 246 insertions(+), 541 deletions(-)
 delete mode 100644 src/adaptive_boris.jl

diff --git a/Project.toml b/Project.toml
index 6c88be42b..ab595de9c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,7 +17,6 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
-OrdinaryDiffEqBoris = "1e8b2ed2-ea1c-4b6e-b3f5-67c29d6bd197"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
diff --git a/src/TestParticle.jl b/src/TestParticle.jl
index 0b034be71..06a696890 100644
--- a/src/TestParticle.jl
+++ b/src/TestParticle.jl
@@ -37,7 +37,7 @@ export trace!, trace_relativistic!, trace_normalized!, trace_relativistic_normal
     get_gc_velocity, full_to_gc, gc_to_full
 export Proton, Electron, Ion
 export Maxwellian, BiMaxwellian, Kappa, BiKappa
-export AdaptiveBoris, AdaptiveHybrid
+export AdaptiveHybrid, Boris, MultistepBoris, MultistepBoris2, MultistepBoris4, MultistepBoris6
 export get_gyrofrequency,
     get_gyroperiod, get_gyroradius, get_velocity, get_energy, get_mean_magnitude,
     energy2velocity, get_curvature_radius, get_adiabaticity,
@@ -63,7 +63,6 @@ include("gc_solver.jl")
 include("equations.jl")
 include("boris.jl")
 include("boris_kernel.jl")
-include("adaptive_boris.jl")
 include("hybrid.jl")
 include("fieldline.jl")
 
diff --git a/src/adaptive_boris.jl b/src/adaptive_boris.jl
deleted file mode 100644
index 8be6319c3..000000000
--- a/src/adaptive_boris.jl
+++ /dev/null
@@ -1,303 +0,0 @@
-# Adaptive Boris method
-
-struct AdaptiveBoris{T}
-    safety::T
-end
-
-"""
-    AdaptiveBoris(; safety=0.1)
-
-Adaptive Boris method with adaptive time stepping based on local gyroperiod.
-The time step is determined by `dt = safety * T_gyro = safety * 2π / |qB/m|`.
-"""
-function AdaptiveBoris(; safety = 0.1)
-    T = typeof(safety)
-    return AdaptiveBoris{T}(T(safety))
-end
-"""
-    solve(prob::TraceProblem, alg::AdaptiveBoris,
-        ensemblealg::BasicEnsembleAlgorithm=EnsembleSerial();
-        trajectories::Int=1, savestepinterval::Int=1,
-        isoutside::Function=ODE_DEFAULT_ISOUTOFDOMAIN,
-        save_start::Bool=true, save_end::Bool=true, save_everystep::Bool=true,
-        save_fields::Bool=false, save_work::Bool=false,
-        batch_size::Int = max(1, trajectories ÷ Threads.nthreads()))
-
-Trace particles using the Adaptive Boris method with specified `prob` and `alg`.
-
-# keywords
-
-  - `trajectories::Int`: number of trajectories to trace.
-  - `savestepinterval::Int`: saving output interval.
-  - `isoutside::Function`: a function with input of position and velocity vector `xv` that determines whether to stop tracing.
-  - `save_start::Bool`: save the initial condition. Default is `true`.
-  - `save_end::Bool`: save the final condition. Default is `true`.
-  - `save_everystep::Bool`: save the state at every `savestepinterval`. Default is `true`.
-  - `save_fields::Bool`: save the electric and magnetic fields. Default is `false`.
-  - `save_work::Bool`: save the work done by the electric field. Default is `false`.
-  - `batch_size::Int`: the number of trajectories to process per worker in `EnsembleDistributed` and `EnsembleSplitThreads`. Default is `max(1, trajectories ÷ nworkers())` for distributed and 1 for others.
-
-
-"""
-function solve(
-        prob::TraceProblem, alg::AdaptiveBoris,
-        ensemblealg::BasicEnsembleAlgorithm = EnsembleSerial();
-        trajectories::Int = 1, savestepinterval::Int = 1,
-        isoutside::F = ODE_DEFAULT_ISOUTOFDOMAIN,
-        save_start::Bool = true, save_end::Bool = true, save_everystep::Bool = true,
-        save_fields::Bool = false, save_work::Bool = false,
-        batch_size::Int = (ensemblealg isa EnsembleDistributed || ensemblealg isa EnsembleSplitThreads) ?
-            max(1, trajectories ÷ nworkers()) : 1
-    ) where {F}
-
-    return if save_fields
-        if save_work
-            return _solve(
-                ensemblealg, prob, trajectories, alg, savestepinterval, isoutside,
-                save_start, save_end, save_everystep, Val(true), Val(true), batch_size
-            )
-        else
-            return _solve(
-                ensemblealg, prob, trajectories, alg, savestepinterval, isoutside,
-                save_start, save_end, save_everystep, Val(true), Val(false), batch_size
-            )
-        end
-    else
-        if save_work
-            return _solve(
-                ensemblealg, prob, trajectories, alg, savestepinterval, isoutside,
-                save_start, save_end, save_everystep, Val(false), Val(true), batch_size
-            )
-        else
-            return _solve(
-                ensemblealg, prob, trajectories, alg, savestepinterval, isoutside,
-                save_start, save_end, save_everystep, Val(false), Val(false), batch_size
-            )
-        end
-    end
-end
-
-function _solve(
-        ::EnsembleSerial, prob, trajectories, alg::AdaptiveBoris, savestepinterval,
-        isoutside, save_start, save_end, save_everystep,
-        ::Val{SaveFields}, ::Val{SaveWork}, batch_size
-    ) where {SaveFields, SaveWork}
-    # We cannot precalculate nt for adaptive steps
-    sol_type = _get_sol_type(prob, zero(eltype(prob.tspan)), Val(SaveFields), Val(SaveWork))
-    sols = Vector{sol_type}(undef, trajectories)
-    irange = 1:trajectories
-
-    _adaptive_boris!(
-        sols, prob, irange, alg, savestepinterval, isoutside,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
-    )
-
-    return sols
-end
-
-function _solve(
-        ::EnsembleThreads, prob, trajectories, alg::AdaptiveBoris, savestepinterval,
-        isoutside, save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
-        batch_size
-    ) where {SaveFields, SaveWork}
-    sol_type = _get_sol_type(prob, zero(eltype(prob.tspan)), Val(SaveFields), Val(SaveWork))
-    sols = Vector{sol_type}(undef, trajectories)
-
-    nchunks = Threads.nthreads()
-    Threads.@threads for irange in index_chunks(1:trajectories; n = nchunks)
-        _adaptive_boris!(
-            sols, prob, irange, alg, savestepinterval, isoutside,
-            save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
-        )
-    end
-
-    return sols
-end
-
-"See `_solve_single_boris` for the rationale behind the `single_prob` construction."
-function _solve_single_adaptive_boris(
-        prob, i, alg::AdaptiveBoris, savestepinterval, isoutside,
-        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}
-    ) where {SaveFields, SaveWork}
-    new_prob = prob.prob_func(prob, i, false)
-    single_prob = TraceProblem(new_prob.u0, new_prob.tspan, new_prob.p)
-    sol_type = _get_sol_type(
-        single_prob, zero(eltype(single_prob.tspan)), Val(SaveFields), Val(SaveWork)
-    )
-    local_sols = Vector{sol_type}(undef, 1)
-    _adaptive_boris!(
-        local_sols, single_prob, 1:1, alg, savestepinterval, isoutside,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
-    )
-    return local_sols[1]
-end
-
-function _solve(
-        ::EnsembleDistributed, prob, trajectories, alg::AdaptiveBoris, savestepinterval,
-        isoutside, save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
-        batch_size
-    ) where {SaveFields, SaveWork}
-    return pmap(1:trajectories; batch_size = batch_size) do i
-        _solve_single_adaptive_boris(
-            prob, i, alg, savestepinterval, isoutside,
-            save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
-        )
-    end
-end
-
-function _solve(
-        ::EnsembleSplitThreads, prob, trajectories, alg::AdaptiveBoris, savestepinterval,
-        isoutside, save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
-        batch_size
-    ) where {SaveFields, SaveWork}
-    # _solve_single_adaptive_boris wraps each trajectory in a fresh TraceProblem(u0, tspan, p)
-    # with DEFAULT_PROB_FUNC. We get a sample problem from prob_func to ensure the
-    # u0 type (uType) matches.
-    sample_prob = prob.prob_func(prob, 1, false)
-    dummy_prob = TraceProblem(sample_prob.u0, sample_prob.tspan, sample_prob.p)
-    sol_type = _get_sol_type(
-        dummy_prob, zero(eltype(dummy_prob.tspan)), Val(SaveFields), Val(SaveWork)
-    )
-    ichunks = index_chunks(1:trajectories; size = batch_size)
-    results = pmap(ichunks) do irange
-        local_sols = Vector{sol_type}(undef, length(irange))
-        Threads.@threads for k in eachindex(irange)
-            i = irange[k]
-            local_sols[k] = _solve_single_adaptive_boris(
-                prob, i, alg, savestepinterval, isoutside,
-                save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
-            )
-        end
-        local_sols
-    end
-    return reduce(vcat, results)
-end
-
-
-@muladd function _adaptive_boris!(
-        sols, prob, irange, alg, savestepinterval, isoutside,
-        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}
-    ) where {SaveFields, SaveWork}
-    (; tspan, p, u0) = prob
-    q2m, _, _, Bfunc, _ = p
-    T = eltype(u0)
-    vars_dim = 6
-    if SaveFields
-        vars_dim += 6
-    end
-    if SaveWork
-        vars_dim += 4
-    end
-
-    # Pre-calculate common factors for time step calculation
-    # dt = safety * 2π / (abs(q2m) * Bmag)
-    C = (2π * alg.safety * sign(tspan[2] - tspan[1])) / abs(q2m)
-
-    @fastmath @inbounds for i in irange
-        # Initialize solution containers
-        initial_capacity = 1000
-        traj = Vector{SVector{vars_dim, T}}(undef, 0)
-        tsave = Vector{typeof(tspan[1])}(undef, 0)
-        sizehint!(traj, initial_capacity)
-        sizehint!(tsave, initial_capacity)
-
-        new_prob = prob.prob_func(prob, i, false)
-        # Load independent r and v SVector from u0
-        u0_i = SVector{6, T}(new_prob.u0)
-        r = u0_i[SVector(1, 2, 3)]
-        v = u0_i[SVector(4, 5, 6)]
-        t = tspan[1]
-        ttotal = tspan[2] - tspan[1]
-
-        if save_start
-            data = _prepare_saved_data(u0_i, p, t, Val(SaveFields), Val(SaveWork))
-            push!(traj, data)
-            push!(tsave, t)
-        end
-
-        # Initial dt calculation
-        Bmag = norm(Bfunc(r, t))
-        dt = C / Bmag
-
-        # Backstep velocity: v(0) -> v(-1/2) using dt
-        v = update_velocity(v, r, -0.5 * dt, t, p)
-        it = 1
-        should_save_final = save_end
-        retcode = ReturnCode.Success
-
-        while abs(t - tspan[1]) < abs(ttotal)
-            # Check if next step exceeds tspan[2]
-            if abs(t + dt - tspan[1]) > abs(ttotal)
-                dt_step = tspan[2] - t
-                # Resync v from `t - 0.5*dt` to `t - 0.5*dt_step`
-                v = update_velocity(v, r, 0.5 * dt, t, p)
-                v = update_velocity(v, r, -0.5 * dt_step, t, p)
-                dt = dt_step
-            end
-
-            # Saving logic (start of step)
-            if save_everystep && (it - 1) > 0 && (it - 1) % savestepinterval == 0
-                # Advance to t to get v_n
-                v_save = update_velocity(v, r, 0.5 * dt, t, p)
-
-                xv_s = vcat(r, v_save)
-                data = _prepare_saved_data(xv_s, p, t, Val(SaveFields), Val(SaveWork))
-                push!(traj, data)
-                push!(tsave, t)
-            end
-
-            # Update velocity to v_{n+1/2}
-            t_mid = t + 0.5 * dt
-            v_new = update_velocity(v, r, dt, t_mid, p)
-
-            # Update location x_{n} -> x_{n+1}
-            r_next = r + v_new * dt
-            t_next = t + dt
-
-            xv_new = vcat(r_next, v_new)
-            if isoutside(xv_new, p, t_next)
-                should_save_final = true
-                retcode = ReturnCode.Terminated
-                break
-            end
-
-            r = r_next
-            t = t_next
-            v = v_new
-
-            # New dt
-            Bmag = norm(Bfunc(r, t))
-            dt_new = C / Bmag
-
-            # Resync v_{n+1/2}(dt) to v_{n+1/2}(dt_new)
-            # v is at t_{new} - 0.5 * dt_old (relative to t_{new})
-            # i.e. it is v_{n+1/2} from step we just took.
-            v = update_velocity(v, r, 0.5 * dt, t, p)
-            v = update_velocity(v, r, -0.5 * dt_new, t, p)
-            dt = dt_new
-            it += 1
-            if save_everystep && (it - 1) % savestepinterval == 0
-                should_save_final = true
-            end
-        end
-
-        if should_save_final && (isempty(tsave) || tsave[end] != t)
-            # v is at t - 0.5*dt. To get v at t, advance by 0.5*dt
-            v_final = update_velocity(v, r, 0.5 * dt, t, p)
-
-            xv_f = vcat(r, v_final)
-            data = _prepare_saved_data(xv_f, p, t, Val(SaveFields), Val(SaveWork))
-            push!(traj, data)
-            push!(tsave, t)
-        end
-
-        # Construct solution
-        sol_alg = :adaptive_boris
-        interp = LinearInterpolation(tsave, traj)
-        stats = nothing
-
-        sols[i] = build_solution(prob, sol_alg, tsave, traj; interp, retcode, stats)
-    end
-
-    return
-end
diff --git a/src/boris.jl b/src/boris.jl
index c0c4f5079..e5f7d99b9 100644
--- a/src/boris.jl
+++ b/src/boris.jl
@@ -1,5 +1,40 @@
 # Native particle pusher
 
+"""
+    Boris(; safety=0.0)
+
+The standard Boris method for particle pushing in electric and magnetic fields.
+When `safety > 0.0`, it uses adaptive time stepping based on local gyroperiod.
+"""
+struct Boris{T}
+    safety::T
+end
+Boris(; safety = 0.0) = Boris(safety)
+
+"""
+    MultistepBoris{N}(; n=1, safety=0.0)
+
+The Multistep/Hyper Boris method of order `N`.
+`n` specifies the number of subcycles.
+`N` specifies the gyrophase correction order: 2 (standard), 4, or 6 (Hyper-Boris).
+When `safety > 0.0`, it uses adaptive time stepping.
+"""
+struct MultistepBoris{N, T}
+    n::Int
+    safety::T
+end
+@inline function MultistepBoris{N}(; n::Int = 1, safety = 0.0) where {N}
+    if N ∉ (2, 4, 6)
+        throw(ArgumentError("Multistep Boris order N must be 2, 4, or 6."))
+    end
+    return MultistepBoris{N, typeof(safety)}(n, safety)
+end
+
+const MultistepBoris2 = MultistepBoris{2}
+const MultistepBoris4 = MultistepBoris{4}
+const MultistepBoris6 = MultistepBoris{6}
+
+
 struct TraceProblem{uType, tType, isinplace, P, F <: AbstractODEFunction, PF} <:
     AbstractODEProblem{uType, tType, isinplace}
     f::F
@@ -108,9 +143,10 @@ In-place cross product.
 end
 
 """
-    solve(prob::TraceProblem; trajectories::Int=1, dt::AbstractFloat,
+    solve(prob::TraceProblem, alg::Union{Boris, MultistepBoris}=Boris();
+        trajectories::Int=1, dt=nothing,
         savestepinterval::Int=1, isoutside::Function=ODE_DEFAULT_ISOUTOFDOMAIN,
-        n::Int=1, save_start::Bool=true, save_end::Bool=true, save_everystep::Bool=true,
+        save_start::Bool=true, save_end::Bool=true, save_everystep::Bool=true,
         save_fields::Bool=false, save_work::Bool=false)
 
 Trace particles using the Boris method with specified `prob`.
@@ -118,92 +154,104 @@ Trace particles using the Boris method with specified `prob`.
 # keywords
 
   - `trajectories::Int`: number of trajectories to trace.
-  - `dt::AbstractFloat`: time step.
+  - `dt::AbstractFloat`: time step for fixed-step Boris. If not provided, `alg.safety` must be > 0.
   - `savestepinterval::Int`: saving output interval.
   - `isoutside::Function`: pinpointing impact or checking boundaries.
-  - `n::Int=1`: number of substeps for the Multistep Boris method. 1 is standard Boris.
-  - `N::Int=2`: order of the Hyper Boris gyrophase correction (2, 4, or 6). 2 is uncorrected.
   - `save_start::Bool=true`: save the initial condition.
   - `save_end::Bool=true`: save the final condition.
   - `save_everystep::Bool=true`: save the state at every `savestepinterval`.
   - `save_fields::Bool=false`: save the electric and magnetic fields.
   - `save_work::Bool=false`: save the work done by the electric field.
   - `batch_size::Int=max(1, trajectories ÷ nworkers())`: the number of trajectories to process per worker in `EnsembleDistributed` and `EnsembleSplitThreads`.
-
 """
 @inline function solve(
-        prob::TraceProblem, ensemblealg::EA = EnsembleSerial();
-        trajectories::Int = 1, savestepinterval::Int = 1, dt::AbstractFloat,
-        isoutside::F = ODE_DEFAULT_ISOUTOFDOMAIN, n::Int = 1, N::Int = 2,
+        prob::TraceProblem, alg::Union{Boris, MultistepBoris} = Boris(),
+        ensemblealg::EA = EnsembleSerial();
+        trajectories::Int = 1, savestepinterval::Int = 1, dt = nothing,
+        isoutside::F = ODE_DEFAULT_ISOUTOFDOMAIN,
         save_start::Bool = true, save_end::Bool = true, save_everystep::Bool = true,
         save_fields::Bool = false, save_work::Bool = false, maxiters::Int = 1_000_000,
         batch_size::Int = (ensemblealg isa EnsembleDistributed || ensemblealg isa EnsembleSplitThreads) ?
-            max(1, trajectories ÷ nworkers()) : 1
+            max(1, trajectories ÷ nworkers()) : 1,
+        n::Int = 1, N::Int = 2
     ) where {EA <: BasicEnsembleAlgorithm, F}
 
-    if N ∉ (2, 4, 6)
-        throw(ArgumentError("N must be 2, 4, or 6"))
+    # Backward compatibility: promote Boris to MultistepBoris if n > 1 or N != 2
+    if alg isa Boris && (n > 1 || N != 2)
+        if N ∉ (2, 4, 6)
+            throw(ArgumentError("Multistep Boris order N must be 2, 4, or 6."))
+        end
+        alg = MultistepBoris{N}(n = n, safety = alg.safety)
+    end
+
+    if !isnothing(dt)
+        if dt < eps(eltype(dt)) * 100
+            throw(ArgumentError("Time step dt is too small."))
+        end
+        ttotal = prob.tspan[2] - prob.tspan[1]
+        if abs(ttotal / dt) > maxiters
+            throw(ArgumentError("Total steps exceed maxiters. Increase maxiters or dt."))
+        end
     end
 
     return _solve(
-        ensemblealg, prob, trajectories, dt, savestepinterval, isoutside, n, N,
+        ensemblealg, prob, trajectories, alg, dt, savestepinterval, isoutside,
         save_start, save_end, save_everystep, Val(save_fields), Val(save_work), maxiters,
         batch_size
     )
 end
 
+function solve(prob::TraceProblem, ensemblealg::BasicEnsembleAlgorithm; kwargs...)
+    return solve(prob, Boris(), ensemblealg; kwargs...)
+end
+
 function _dispatch_boris!(
-        sols, prob::TraceProblem, irange, savestepinterval, dt, nt, nout, isoutside::F,
-        n, N, save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}
+        sols, prob::TraceProblem, irange, alg, dt, savestepinterval, isoutside::F,
+        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
+        maxiters
     ) where {SaveFields, SaveWork, F}
-    return if n == 1 && N == 2
-        _boris!(
-            sols, prob, irange, savestepinterval, dt, nt, nout, isoutside,
-            save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
-        )
-    else
-        _multistep_boris!(
-            sols, prob, irange, savestepinterval, dt, nt, nout, isoutside, n, N,
-            save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
-        )
-    end
+    return _generic_boris!(
+        sols, prob, irange, alg, dt, savestepinterval, isoutside,
+        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork),
+        maxiters
+    )
 end
 
 @inline function _solve(
-        ::EnsembleSerial, prob::TraceProblem, trajectories, dt, savestepinterval,
-        isoutside::F, n, N, save_start, save_end, save_everystep,
+        ::EnsembleSerial, prob::TraceProblem, trajectories, alg, dt, savestepinterval,
+        isoutside::F, save_start, save_end, save_everystep,
         ::Val{SaveFields}, ::Val{SaveWork}, maxiters, batch_size
     ) where {SaveFields, SaveWork, F}
-    sols, nt,
-        nout = _prepare(
-        prob, trajectories, dt, savestepinterval,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork), maxiters
+    sols = _prepare(
+        prob, trajectories, alg, dt, save_start, save_end, save_everystep,
+        Val(SaveFields), Val(SaveWork)
     )
     irange = 1:trajectories
     _dispatch_boris!(
-        sols, prob, irange, savestepinterval, dt, nt, nout, isoutside, n, N,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
+        sols, prob, irange, alg, dt, savestepinterval, isoutside,
+        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork),
+        maxiters
     )
 
     return sols
 end
 
 @inline function _solve(
-        ::EnsembleThreads, prob::TraceProblem, trajectories, dt, savestepinterval,
-        isoutside::F, n, N, save_start, save_end, save_everystep,
+        ::EnsembleThreads, prob::TraceProblem, trajectories, alg, dt, savestepinterval,
+        isoutside::F, save_start, save_end, save_everystep,
         ::Val{SaveFields}, ::Val{SaveWork}, maxiters, batch_size
     ) where {SaveFields, SaveWork, F}
-    sols, nt,
-        nout = _prepare(
-        prob, trajectories, dt, savestepinterval,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork), maxiters
+    sols = _prepare(
+        prob, trajectories, alg, dt, save_start, save_end, save_everystep,
+        Val(SaveFields), Val(SaveWork)
     )
 
     nchunks = Threads.nthreads()
     Threads.@threads for irange in index_chunks(1:trajectories; n = nchunks)
         _dispatch_boris!(
-            sols, prob, irange, savestepinterval, dt, nt, nout, isoutside, n, N,
-            save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
+            sols, prob, irange, alg, dt, savestepinterval, isoutside,
+            save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork),
+            maxiters
         )
     end
 
@@ -229,61 +277,59 @@ The `TraceProblem` construction is a negligible struct copy relative to the
 simulation cost and the serialization overhead inherent in `pmap`.
 """
 function _solve_single_boris(
-        prob::TraceProblem, i, savestepinterval, dt, nt, nout, isoutside::F, n, N,
-        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}
+        prob::TraceProblem, i, alg, dt, savestepinterval, isoutside::F,
+        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
+        maxiters
     ) where {SaveFields, SaveWork, F}
     new_prob = prob.prob_func(prob, i, false)
     single_prob = TraceProblem(new_prob.u0, new_prob.tspan, new_prob.p)
-    sol_type = _get_sol_type(single_prob, dt, Val(SaveFields), Val(SaveWork))
+    sol_type = _get_sol_type(
+        single_prob, isnothing(dt) ? zero(eltype(single_prob.tspan)) : dt,
+        Val(SaveFields), Val(SaveWork)
+    )
     local_sols = Vector{sol_type}(undef, 1)
     _dispatch_boris!(
-        local_sols, single_prob, 1:1, savestepinterval, dt, nt, nout,
-        isoutside, n, N, save_start, save_end, save_everystep,
-        Val(SaveFields), Val(SaveWork)
+        local_sols, single_prob, 1:1, alg, dt, savestepinterval,
+        isoutside, save_start, save_end, save_everystep,
+        Val(SaveFields), Val(SaveWork), maxiters
     )
     return local_sols[1]
 end
 
 @inline function _solve(
-        ::EnsembleDistributed, prob::TraceProblem, trajectories, dt, savestepinterval,
-        isoutside::F, n, N, save_start, save_end, save_everystep,
+        ::EnsembleDistributed, prob::TraceProblem, trajectories, alg, dt, savestepinterval,
+        isoutside::F, save_start, save_end, save_everystep,
         ::Val{SaveFields}, ::Val{SaveWork}, maxiters, batch_size
     ) where {SaveFields, SaveWork, F}
-    _, nt, nout = _prepare(
-        prob, trajectories, dt, savestepinterval,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork), maxiters
-    )
     return pmap(1:trajectories; batch_size = batch_size) do i
         _solve_single_boris(
-            prob, i, savestepinterval, dt, nt, nout, isoutside, n, N,
-            save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
+            prob, i, alg, dt, savestepinterval, isoutside,
+            save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork),
+            maxiters
         )
     end
 end
 
 @inline function _solve(
-        ::EnsembleSplitThreads, prob::TraceProblem, trajectories, dt, savestepinterval,
-        isoutside::F, n, N, save_start, save_end, save_everystep,
+        ::EnsembleSplitThreads, prob::TraceProblem, trajectories, alg, dt, savestepinterval,
+        isoutside::F, save_start, save_end, save_everystep,
         ::Val{SaveFields}, ::Val{SaveWork}, maxiters, batch_size
     ) where {SaveFields, SaveWork, F}
-    _, nt, nout = _prepare(
-        prob, trajectories, dt, savestepinterval,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork), maxiters
-    )
     # _solve_single_boris wraps each trajectory in a fresh TraceProblem(u0, tspan, p)
     # with DEFAULT_PROB_FUNC. We get a sample problem from prob_func to ensure the
     # u0 type (uType) matches.
     sample_prob = prob.prob_func(prob, 1, false)
     dummy_prob = TraceProblem(sample_prob.u0, sample_prob.tspan, sample_prob.p)
-    sol_type = _get_sol_type(dummy_prob, dt, Val(SaveFields), Val(SaveWork))
+    sol_type = _get_sol_type(dummy_prob, zero(eltype(dummy_prob.tspan)), Val(SaveFields), Val(SaveWork))
     ichunks = index_chunks(1:trajectories; size = batch_size)
     results = pmap(ichunks) do irange
         local_sols = Vector{sol_type}(undef, length(irange))
         Threads.@threads for k in eachindex(irange)
             i = irange[k]
             local_sols[k] = _solve_single_boris(
-                prob, i, savestepinterval, dt, nt, nout, isoutside, n, N,
-                save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork)
+                prob, i, alg, dt, savestepinterval, isoutside,
+                save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork),
+                maxiters
             )
         end
         local_sols
@@ -319,123 +365,129 @@ end
 Prepare for advancing.
 """
 function _prepare(
-        prob::TraceProblem, trajectories, dt, savestepinterval,
-        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}, maxiters
+        prob::TraceProblem, trajectories, alg, dt,
+        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}
     ) where {SaveFields, SaveWork}
-    if abs(dt) < 10 * eps(typeof(dt))
+    if !isnothing(dt) && abs(dt) < 10 * eps(typeof(dt))
         throw(ArgumentError("time step dt is too small, violating min_dt = 10 * eps(typeof(dt))"))
     end
-    ttotal = prob.tspan[2] - prob.tspan[1]
-    nt = round(Int, ttotal / dt) |> abs
-    if nt > maxiters
-        throw(ArgumentError("number of iterations nt ($nt) exceeds maxiters ($maxiters)"))
-    end
-
-    nout = 0
-    if save_start
-        nout += 1
-    end
-
-    if save_everystep
-        steps = nt ÷ savestepinterval
-        last_is_step = (nt > 0) && (nt % savestepinterval == 0)
-        nout += steps
-        if !save_end && last_is_step
-            nout -= 1
-        end
-        if save_end && !last_is_step
-            nout += 1
-        end
-    elseif save_end
-        nout += 1
-    end
-
-    sol_type = _get_sol_type(prob, dt, Val(SaveFields), Val(SaveWork))
+    
+    sol_type = _get_sol_type(
+        prob, isnothing(dt) ? zero(eltype(prob.tspan)) : dt,
+        Val(SaveFields), Val(SaveWork)
+    )
     sols = Vector{sol_type}(undef, trajectories)
 
-    return sols, nt, nout
+    return sols
 end
 
 @inline function _prepare_saved_data(xv, p, t, ::Val{SaveFields}, ::Val{SaveWork}) where {SaveFields, SaveWork}
     data = xv
-
-    # Pre-declare variables to share between blocks if both are true
-    local E_field, magnetic_props
-
-    if SaveFields
+    if SaveFields || SaveWork
         r = get_x(xv)
         T = eltype(xv)
+        Bfunc = get_BField(p)
         E_field = SVector{3, T}(get_EField(p)(r, t))
 
-        # We need magnetic properties for work, so if SaveWork is also true, compute them now
         if SaveWork
-            q2m, m, Efunc, Bfunc, _ = p
             # get_magnetic_properties returns (B, ∇B, κ, b̂, Bmag)
             magnetic_props = get_magnetic_properties(r, t, Bfunc)
-            B_vec = SVector{3, T}(magnetic_props[1])
-            data = vcat(data, E_field, B_vec)
-        else
-            B_vec = SVector{3, T}(get_BField(p)(r, t))
-            data = vcat(data, E_field, B_vec)
-        end
-    end
-
-    if SaveWork
-        # If SaveFields was true, we already computed E_field and magnetic_props
-        if SaveFields
+            if SaveFields
+                B_vec = SVector{3, T}(magnetic_props[1])
+                data = vcat(data, E_field, B_vec)
+            end
             work = get_work_rates(xv, p, t, magnetic_props, E_field)
-        else
-            work = get_work_rates(xv, p, t)
+            data = vcat(data, work)
+        elseif SaveFields
+            B_vec = SVector{3, T}(Bfunc(r, t))
+            data = vcat(data, E_field, B_vec)
         end
-        data = vcat(data, work)
     end
 
     return data
 end
 
-"""
-Apply Boris method for particles with index in `irange`.
-"""
 @inline @muladd function _boris_loop!(
-        traj, tsave, iout, r, v, p, dt, nt, tspan,
+        traj, tsave, r, v, p, dt, tspan,
         savestepinterval, save_everystep, isoutside::F1, velocity_updater::F2,
-        ::Val{SaveFields}, ::Val{SaveWork}
+        ::Val{SaveFields}, ::Val{SaveWork}, alg, maxiters
     ) where {F1, F2, SaveFields, SaveWork}
+    t = tspan[1]
+    ttotal = tspan[2] - tspan[1]
     it = 1
-    t = tspan[1] - 0.5 * dt
-    while it <= nt
-        v_prev = v
-        t += dt
-        v = velocity_updater(v, r, dt, t, p)
-
-        r_next = r + v * dt
-        t_next = t + 0.5 * dt
-        if isoutside(vcat(r_next, v), p, t_next)
-            return it - 1, iout, r, v_prev
+    retcode = ReturnCode.Success
+
+    C = zero(eltype(v))
+    if alg.safety > 0.0
+        q2m, _, _, Bfunc, _ = p
+        C = (2π * alg.safety * sign(ttotal)) / abs(q2m)
+        Bmag = norm(Bfunc(r, t))
+        dt = C / Bmag
+    end
+
+    # push velocity back in time by 1/2 dt
+    v = velocity_updater(v, r, -0.5 * dt, t, p)
+
+    while abs(t - tspan[1]) < abs(ttotal)
+        if abs(t + dt - tspan[1]) >= abs(ttotal) - 100 * eps(eltype(v)(abs(ttotal)))
+            dt_step = tspan[2] - t
+            if dt_step != dt
+                # Resync v from `t - 0.5*dt` to `t - 0.5*dt_step`
+                v = velocity_updater(v, r, 0.5 * dt, t, p)
+                v = velocity_updater(v, r, -0.5 * dt_step, t, p)
+                dt = dt_step
+            end
         end
 
         if save_everystep && (it - 1) > 0 && (it - 1) % savestepinterval == 0
-            iout += 1
-            if iout <= length(traj)
-                t_current = t - 0.5 * dt
-                v_save = velocity_updater(v_prev, r, 0.5 * dt, t_current, p)
-                data = vcat(r, v_save)
-                traj[iout] = _prepare_saved_data(data, p, t_current, Val(SaveFields), Val(SaveWork))
-                tsave[iout] = t_current
-            end
+            # Advance to t to get v_n
+            v_save = velocity_updater(v, r, 0.5 * dt, t, p)
+            data = vcat(r, v_save)
+            push!(traj, _prepare_saved_data(data, p, t, Val(SaveFields), Val(SaveWork)))
+            push!(tsave, t)
         end
 
-        r = r_next
+        v_new = velocity_updater(v, r, dt, t + 0.5 * dt, p)
+        r_next = r + v_new * dt
+        t_next = t + dt
+
+        # NaN check
+        if any(isnan, r_next) || any(isnan, v_new)
+            retcode = ReturnCode.Unstable
+            break
+        end
+
+        if isoutside(vcat(r_next, v_new), p, t_next)
+            retcode = ReturnCode.Terminated
+            break
+        end
+
+        r, v, t = r_next, v_new, t_next
+
+        if alg.safety > 0.0 && abs(t - tspan[1]) < abs(ttotal)
+            q2m, _, _, Bfunc, _ = p
+            Bmag = norm(Bfunc(r, t))
+            dt_new = C / Bmag
+            # Resync v_{n+1/2}(dt) to v_{n+1/2}(dt_new)
+            v = velocity_updater(v, r, 0.5 * dt, t, p)
+            v = velocity_updater(v, r, -0.5 * dt_new, t, p)
+            dt = dt_new
+        end
+        if it > maxiters
+            retcode = ReturnCode.MaxIters
+            break
+        end
         it += 1
     end
-    return it - 1, iout, r, v
+
+    return r, v, t, dt, retcode
 end
 
 @inline @muladd function _generic_boris!(
-        sols, prob::TraceProblem, irange, savestepinterval, dt, nt, nout, isoutside::F1,
+        sols, prob::TraceProblem, irange, alg, dt, savestepinterval, isoutside::F1,
         save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
-        velocity_updater::F2, alg_name
-    ) where {SaveFields, SaveWork, F1, F2}
+        maxiters
+    ) where {SaveFields, SaveWork, F1}
     (; tspan, p, u0) = prob
     T = eltype(u0)
 
@@ -447,101 +499,75 @@ end
         vars_dim += 4
     end
 
+    algorithm_name = alg isa Boris ? :boris : :multistep_boris
+
+    velocity_updater = if alg isa Boris
+        update_velocity
+    else
+        (v, r, dt, t, p) -> update_velocity_multistep(v, r, dt, t, alg.n, _get_val_N(alg), p)
+    end
+
     @inbounds for i in irange
-        traj = Vector{SVector{vars_dim, T}}(undef, nout)
-        tsave = Vector{typeof(tspan[1] + dt)}(undef, nout)
+        traj = Vector{SVector{vars_dim, T}}(undef, 0)
+        tsave = Vector{typeof(tspan[1] + (isnothing(dt) ? 0.0 : dt))}(undef, 0)
+        sizehint!(traj, 1000)
+        sizehint!(tsave, 1000)
 
         # set initial conditions for each trajectory i
-        iout = 0
         new_prob = prob.prob_func(prob, i, false)
         u0_i = SVector{6, T}(new_prob.u0)
         r = u0_i[SVector(1, 2, 3)]
         v = u0_i[SVector(4, 5, 6)]
 
+        # If dt is not provided, it must be adaptive (safety > 0)
+        _dt = isnothing(dt) ? zero(T) : T(dt)
+
         if save_start
-            iout += 1
-            traj[iout] = _prepare_saved_data(u0_i, p, tspan[1], Val(SaveFields), Val(SaveWork))
-            tsave[iout] = tspan[1]
+            push!(traj, _prepare_saved_data(u0_i, p, tspan[1], Val(SaveFields), Val(SaveWork)))
+            push!(tsave, tspan[1])
         end
 
-        # push velocity back in time by 1/2 dt
-        v = velocity_updater(v, r, -0.5 * dt, tspan[1], p)
-
-        it, iout, r, v = _boris_loop!(
-            traj, tsave, iout, r, v, p, dt, nt, tspan,
+        r, v, t, _dt_final, retcode = _boris_loop!(
+            traj, tsave, r, v, p, _dt, tspan,
             savestepinterval, save_everystep, isoutside, velocity_updater,
-            Val(SaveFields), Val(SaveWork)
+            Val(SaveFields), Val(SaveWork), alg, maxiters
         )
 
-        final_step = it
-        should_save_final = false
-        if save_end
-            should_save_final = true
-        elseif save_everystep && (final_step > 0) && (final_step % savestepinterval == 0)
-            should_save_final = true
-        end
-
-        if iout < nout && should_save_final
-            t_final = final_step == nt ? tspan[2] : tspan[1] + final_step * dt
-            if iout == 0 || tsave[iout] < t_final
-                iout += 1
-                dt_final = t_final - (tspan[1] + (final_step - 0.5) * dt)
-                v_final = velocity_updater(v, r, dt_final, t_final, p)
-
-                data = vcat(r, v_final)
-                traj[iout] = _prepare_saved_data(
-                    data, p, t_final, Val(SaveFields), Val(SaveWork)
-                )
-                tsave[iout] = t_final
-            end
-        end
-
-        if iout < nout
-            resize!(traj, iout)
-            resize!(tsave, iout)
-            retcode = ReturnCode.Terminated
-        else
-            retcode = ReturnCode.Success
+        should_save_final = save_end
+        if should_save_final && (isempty(tsave) || tsave[end] != t)
+            # v is at t - 0.5*_dt_final. To get v at t, advance by 0.5*_dt_final
+            v_save = velocity_updater(v, r, 0.5 * _dt_final, t, p)
+            data = vcat(r, v_save)
+            push!(traj, _prepare_saved_data(data, p, t, Val(SaveFields), Val(SaveWork)))
+            push!(tsave, t)
         end
-        traj_save = traj
-        t = tsave
 
-        alg = alg_name
-        interp = LinearInterpolation(t, traj_save)
+        alg_sol = algorithm_name
+        interp = LinearInterpolation(tsave, traj)
         stats = nothing
 
-        sols[i] = build_solution(prob, alg, t, traj_save; interp, retcode, stats)
+        sols[i] = build_solution(prob, alg_sol, tsave, traj; interp, retcode, stats)
     end
 
     return
 end
 
-"""
-Apply Boris method for particles with index in `irange`.
-"""
-@inline @muladd function _boris!(
-        sols, prob::TraceProblem, irange, savestepinterval, dt, nt, nout, isoutside::F,
-        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}
-    ) where {SaveFields, SaveWork, F}
-
-    _generic_boris!(
-        sols, prob, irange, savestepinterval, dt, nt, nout, isoutside,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork),
-        update_velocity, :boris
-    )
-
-    return
-end
+_get_val_N(::MultistepBoris{N}) where {N} = Val{N}()
 
 """
     update_velocity_multistep(v, r, dt, t, n, N, param)
+    update_velocity_multistep(v, r, dt, t, n, ::Val{N}, param)
 
 Update velocity using the Multistep/Hyper Boris method, returning the new velocity as an SVector.
 `n` specifies the number of subcycles.
 `N` specifies the gyrophase correction order. When N=2, it corresponds to the Multicycle solver. When N=4 or N=6, it is the Hyper Boris solver.
 Reference: [Zenitani & Kato 2025](https://arxiv.org/abs/2505.02270)
 """
-@muladd function update_velocity_multistep(v, r, dt, t, n::Int, N::Int, param)
+@inline function update_velocity_multistep(v, r, dt, t, n::Int, N::Int, param)
+    return update_velocity_multistep(v, r, dt, t, n, Val(N), param)
+end
+
+@muladd function update_velocity_multistep(v, r, dt, t, n::Int, ::Val{N}, param) where {N}
     q2m, _, Efunc, Bfunc, _ = param
     E = Efunc(r, t)
     B = Bfunc(r, t)
@@ -618,22 +644,6 @@ Reference: [Zenitani & Kato 2025](https://arxiv.org/abs/2505.02270)
     return v_new
 end
 
-@inline @muladd function _multistep_boris!(
-        sols, prob::TraceProblem, irange, savestepinterval, dt, nt, nout, isoutside::F,
-        n_steps::Int, N_order::Int, save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}
-    ) where {SaveFields, SaveWork, F}
-
-    velocity_updater = (v, r, dt, t, p) ->
-    update_velocity_multistep(v, r, dt, t, n_steps, N_order, p)
-
-    _generic_boris!(
-        sols, prob, irange, savestepinterval, dt, nt, nout, isoutside,
-        save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork),
-        velocity_updater, :multistep_boris
-    )
-
-    return
-end
 
 """
     get_fields(sol::AbstractODESolution)
diff --git a/src/precompile.jl b/src/precompile.jl
index 70dace791..af5ff1899 100644
--- a/src/precompile.jl
+++ b/src/precompile.jl
@@ -37,7 +37,7 @@
         sol_kernel = solve(prob, CPU(); dt, savestepinterval = 100)
 
         # Adaptive Boris
-        alg_adaptive = AdaptiveBoris(safety = 0.1)
+        alg_adaptive = Boris(safety = 0.1)
         sol_adaptive = solve(prob, alg_adaptive)[1]
 
         sol_adaptive_sf = solve(prob, alg_adaptive; save_fields = true, save_work = true)
diff --git a/test/test_boris.jl b/test/test_boris.jl
index 613018af8..7d270e9ae 100644
--- a/test/test_boris.jl
+++ b/test/test_boris.jl
@@ -188,7 +188,7 @@ using Distributed
 
     @testset "Adaptive Boris" begin
         # Check constructor default
-        alg1 = AdaptiveBoris(safety = 0.1)
+        alg1 = Boris(safety = 0.1)
         @test alg1.safety == 0.1
 
         x0 = [0.0, 0.0, 0.0]
@@ -198,7 +198,7 @@ using Distributed
         tperiod = abs(TP.get_gyroperiod(0.01; q = TP.qₑ, m = TP.mₑ))
         tspan = (0.0, 200 * tperiod)
 
-        alg_adaptive = AdaptiveBoris(safety = 0.1)
+        alg_adaptive = Boris(safety = 0.1)
 
         param = prepare(constant_E, gradient_B, species = Electron)
         prob = TraceProblem(stateinit, tspan, param)
@@ -243,7 +243,7 @@ using Distributed
         @test abs(sol_boris.u[end][4]) < 1.0e5 - 100
 
         # Adaptive Boris
-        alg_adaptive = AdaptiveBoris(safety = 0.1)
+        alg_adaptive = Boris(safety = 0.1)
         sol_adaptive = TP.solve(prob, alg_adaptive)[1]
         @test abs(sol_adaptive.u[end][4]) < 1.0e5 - 100
     end
@@ -408,7 +408,7 @@ using Distributed
 
         # Test Adaptive Boris with save_work
         # Use simple AdaptiveBoris
-        alg_adaptive = AdaptiveBoris(safety = 0.1)
+        alg_adaptive = Boris(safety = 0.1)
         sol_adaptive = TP.solve(prob, alg_adaptive; save_work = true, save_everystep = true)[1]
         @test length(sol_adaptive.u[1]) == 10
         work_adaptive = sol_adaptive.u[1][7:10]
@@ -523,10 +523,10 @@ using Distributed
                 end
             end
 
-            @testset "AdaptiveBoris" begin
+            @testset "BorisAdaptive" begin
                 tperiod = abs(TP.get_gyroperiod(0.01; q = TP.qₑ, m = TP.mₑ))
-                alg_adaptive = AdaptiveBoris(; safety = 0.1)
-
+                alg_adaptive = Boris(; safety = 0.1)
+                
                 sols_serial = TP.solve(
                     prob_dist, alg_adaptive, EnsembleSerial(); trajectories
                 )
diff --git a/test/test_boundary.jl b/test/test_boundary.jl
index b320c60af..508a205c8 100644
--- a/test/test_boundary.jl
+++ b/test/test_boundary.jl
@@ -45,7 +45,7 @@ using OrdinaryDiffEq: ReturnCode
         @testset "Adaptive Boris" begin
             prob = TraceProblem(u0, tspan, p)
             sol = TP.solve(
-                prob, AdaptiveBoris(safety = 0.05);
+                prob, Boris(safety = 0.05);
                 isoutside = callback.condition
             )[1]
 
@@ -82,13 +82,13 @@ using OrdinaryDiffEq: ReturnCode
         E_field_fixed(r, t = 0.0) = SA[0.0, 0.0, 0.0]
         param = prepare(E_field_fixed, B_field_fixed; species = Proton)
 
-        @testset "Spatial rejection (AdaptiveBoris)" begin
+        @testset "Spatial rejection (Boris)" begin
             isoutside(u, p, t) = u[1] > 0.5
             u0 = [0.0, 0.0, 0.0, 1.0e5, 0.0, 0.0]
             tspan = (0.0, 1.0e-5)
             prob = TraceProblem(u0, tspan, param)
-
-            alg = AdaptiveBoris(safety = 0.1)
+            
+            alg = Boris(safety = 0.1)
             sol = TP.solve(prob, alg; isoutside)[1]
             @test sol.u[end][1] <= 0.5
             @test sol.t[end] < tspan[2]

From fe240124ee754569aff730d1e72675bd3a8abea5 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 18:31:34 -0400
Subject: [PATCH 07/17] Update benchmarks, demos and docs

---
 README.md                                        |  2 +-
 benchmark/benchmarks.jl                          |  2 +-
 docs/examples/features/demo_adaptive_accuracy.jl | 12 ++++++------
 docs/examples/features/demo_boris.jl             |  4 ++--
 docs/src/tutorial/advanced_boris.md              |  6 +++---
 test/test_boris.jl                               |  3 +--
 test/test_boundary.jl                            |  2 +-
 7 files changed, 15 insertions(+), 16 deletions(-)

diff --git a/README.md b/README.md
index 36434fb25..274dada8d 100644
--- a/README.md
+++ b/README.md
@@ -56,7 +56,7 @@ sol = TestParticle.solve(prob; dt, savestepinterval)[1]
 
 Besides the standard Boris method, we also support various versions of Boris solvers including:
 - **Multistep Boris**: fixed time step with `n > 1` substeps.
-- **Adaptive Boris**: uses `AdaptiveBoris()` for automatic time step selection based on local gyroperiod.
+- **Adaptive Boris**: uses `Boris(safety = 0.1)` to automatically adjust time step based on local gyroperiod.
 
 For plotting with Makie,
 
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 106519682..1eceae722 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -168,7 +168,7 @@ SUITE["trace"]["numerical field"]["Hyper Boris (n=2, N=4)"] = @benchmarkable TP.
 SUITE["trace"]["numerical field"]["Hyper Boris (n=2, N=6)"] = @benchmarkable TP.solve(
     $prob_boris; dt = 1 / 7, savestepinterval = 10, n = 2, N = 6
 )
-alg_adaptive = AdaptiveBoris(safety = 0.1)
+alg_adaptive = Boris(safety = 0.1)
 SUITE["trace"]["numerical field"]["Adaptive Boris"] = @benchmarkable TP.solve(
     $prob_boris, $alg_adaptive
 )
diff --git a/docs/examples/features/demo_adaptive_accuracy.jl b/docs/examples/features/demo_adaptive_accuracy.jl
index cea0b4418..58418cbe9 100644
--- a/docs/examples/features/demo_adaptive_accuracy.jl
+++ b/docs/examples/features/demo_adaptive_accuracy.jl
@@ -92,15 +92,15 @@ prob_ode1 = ODEProblem(trace_normalized!, u0_1, tspan1, param1)
 adaptive_solvers_1 = [
     (
         "AdaptiveBoris (η=1/16)",
-        AdaptiveBoris(; safety = 1 / 16),
+        Boris(; safety = 1 / 16),
     ),
     (
         "AdaptiveBoris (η=1/8)",
-        AdaptiveBoris(; safety = 1 / 8),
+        Boris(; safety = 1 / 8),
     ),
     (
         "AdaptiveBoris (η=1/4)",
-        AdaptiveBoris(; safety = 1 / 4),
+        Boris(; safety = 1 / 4),
     ),
 ]
 
@@ -201,7 +201,7 @@ plot_table(results1) #hide
 #
 # With E = [0, 0.5, 0.1] and B = [0, 0, 1], an exact velocity
 # solution exists (Section 6, Zenitani & Kato 2025). We sweep the
-# `safety` parameter η (representing dt / T_gyro) for `AdaptiveBoris`
+# `safety` parameter η (representing dt / T_gyro) for `Boris(safety = η)`
 # and the tolerance for the ODE solvers, then plot maximum velocity
 # error vs. number of time steps (computational cost).
 
@@ -336,11 +336,11 @@ prob_ode3 = ODEProblem(trace_normalized!, u0_2, tspan3, param3)
 adaptive_solvers_3 = [
     (
         "AdaptiveBoris (η=0.1)",
-        AdaptiveBoris(; safety = 0.1),
+        Boris(; safety = 0.1),
     ),
     (
         "AdaptiveBoris (η=0.2)",
-        AdaptiveBoris(; safety = 0.2),
+        Boris(; safety = 0.2),
     ),
 ]
 
diff --git a/docs/examples/features/demo_boris.jl b/docs/examples/features/demo_boris.jl
index 6afc13454..3ff5b806e 100644
--- a/docs/examples/features/demo_boris.jl
+++ b/docs/examples/features/demo_boris.jl
@@ -99,7 +99,7 @@ sol_boris_2 = TP.solve(prob; dt, n = 2)[1];
 sol_boris_4 = TP.solve(prob; dt, n = 4)[1];
 sol_boris_hyper = TP.solve(prob; dt, n = 2, N = 4)[1];
 
-alg_adaptive = AdaptiveBoris(safety = 0.1)
+alg_adaptive = Boris(safety = 0.1)
 sol_boris_adaptive = TP.solve(prob, alg_adaptive)[1];
 
 # Let's compare against the default ODE solver `Tsit5` from DifferentialEquations.jl, in both fixed time step mode and adaptive mode:
@@ -148,7 +148,7 @@ sol_boris_4 = TP.solve(prob_boris; dt, savestepinterval = 36, n = 4)[1];
 sol_boris_hyper = TP.solve(prob_boris; dt, savestepinterval = 36, n = 2, N = 4)[1];
 sol_boris_adaptive = TP.solve(
     prob_boris,
-    AdaptiveBoris(safety = 0.1)
+    Boris(safety = 0.1)
 )[1]
 sol1 = solve(prob, Tsit5(); adaptive = false, dt, dense = false, saveat = dt);
 sol2 = solve(prob, Tsit5());
diff --git a/docs/src/tutorial/advanced_boris.md b/docs/src/tutorial/advanced_boris.md
index 7fffc373d..48adf1084 100644
--- a/docs/src/tutorial/advanced_boris.md
+++ b/docs/src/tutorial/advanced_boris.md
@@ -74,7 +74,7 @@ Combining both $n > 1$ and $N > 2$ ensures ultra-high stability tracking over dr
 
 ## 3. Adaptive Boris Method
 
-The `AdaptiveBoris` solver adjusts the time step $\Delta t$ dynamically based on the local cyclotron frequency $\Omega_c = |q B / m|$. This is particularly useful in systems with strong magnetic field gradients, such as magnetic mirrors or planetary magnetospheres, where the required resolution varies significantly along the particle's trajectory.
+The `Boris(safety = η)` solver adjusts the time step $\Delta t$ dynamically based on the local cyclotron frequency $\Omega_c = |q B / m|$. This is particularly useful in systems with strong magnetic field gradients, such as magnetic mirrors or planetary magnetospheres, where the required resolution varies significantly along the particle's trajectory.
 
 The time step is determined by:
 ```math
@@ -113,10 +113,10 @@ sol = TestParticle.solve(prob; dt, n=2, N=4)
 
 ### Adaptive Boris
 
-You can use the adaptive solver by passing an `AdaptiveBoris` object as the second argument to `solve`.
+You can use the adaptive solver by passing `Boris(safety = η)` as the second argument to `solve`.
 
 ```julia
 # Adaptive Boris with safety factor 0.05 (20 steps per period)
-alg = AdaptiveBoris(safety=0.05)
+alg = Boris(safety=0.05)
 sol = TestParticle.solve(prob, alg)[1]
 ```
diff --git a/test/test_boris.jl b/test/test_boris.jl
index 7d270e9ae..301f6e0ac 100644
--- a/test/test_boris.jl
+++ b/test/test_boris.jl
@@ -407,7 +407,6 @@ using Distributed
         @test abs(work_ms[3]) > 0.0
 
         # Test Adaptive Boris with save_work
-        # Use simple AdaptiveBoris
         alg_adaptive = Boris(safety = 0.1)
         sol_adaptive = TP.solve(prob, alg_adaptive; save_work = true, save_everystep = true)[1]
         @test length(sol_adaptive.u[1]) == 10
@@ -526,7 +525,7 @@ using Distributed
             @testset "BorisAdaptive" begin
                 tperiod = abs(TP.get_gyroperiod(0.01; q = TP.qₑ, m = TP.mₑ))
                 alg_adaptive = Boris(; safety = 0.1)
-                
+
                 sols_serial = TP.solve(
                     prob_dist, alg_adaptive, EnsembleSerial(); trajectories
                 )
diff --git a/test/test_boundary.jl b/test/test_boundary.jl
index 508a205c8..9255a14f8 100644
--- a/test/test_boundary.jl
+++ b/test/test_boundary.jl
@@ -87,7 +87,7 @@ using OrdinaryDiffEq: ReturnCode
             u0 = [0.0, 0.0, 0.0, 1.0e5, 0.0, 0.0]
             tspan = (0.0, 1.0e-5)
             prob = TraceProblem(u0, tspan, param)
-            
+
             alg = Boris(safety = 0.1)
             sol = TP.solve(prob, alg; isoutside)[1]
             @test sol.u[end][1] <= 0.5

From 7dc9a98f08b710c0a1f3b2632cbe7b3897a1aaa7 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 18:36:14 -0400
Subject: [PATCH 08/17] Formatting

---
 src/boris.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/boris.jl b/src/boris.jl
index e5f7d99b9..23fb3d6c3 100644
--- a/src/boris.jl
+++ b/src/boris.jl
@@ -371,7 +371,7 @@ function _prepare(
     if !isnothing(dt) && abs(dt) < 10 * eps(typeof(dt))
         throw(ArgumentError("time step dt is too small, violating min_dt = 10 * eps(typeof(dt))"))
     end
-    
+
     sol_type = _get_sol_type(
         prob, isnothing(dt) ? zero(eltype(prob.tspan)) : dt,
         Val(SaveFields), Val(SaveWork)

From a7fe7ccf7334ac44e97f07ffaeb98b70a3936e67 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 19:06:02 -0400
Subject: [PATCH 09/17] Update text in demo

---
 docs/examples/features/demo_adaptive_accuracy.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/examples/features/demo_adaptive_accuracy.jl b/docs/examples/features/demo_adaptive_accuracy.jl
index 58418cbe9..05ec650b9 100644
--- a/docs/examples/features/demo_adaptive_accuracy.jl
+++ b/docs/examples/features/demo_adaptive_accuracy.jl
@@ -1,7 +1,7 @@
 # # Adaptive Solver Accuracy
 #
 # This example compares the accuracy and energy conservation of solvers that
-# use **changing time steps**: the native `AdaptiveBoris` solver and the
+# use **changing time steps**: the native `Boris(safety = η)` solver and the
 # adaptive ODE solvers `Tsit5` and `Vern6`.
 #
 # Unlike the fixed-step tests in [Solver Accuracy Analysis](@ref) and

From a923db6edbd028371f7e469457d75292e9c7cb92 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 19:18:02 -0400
Subject: [PATCH 10/17] doc: update demo

---
 docs/examples/features/demo_adaptive_accuracy.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/docs/examples/features/demo_adaptive_accuracy.jl b/docs/examples/features/demo_adaptive_accuracy.jl
index 05ec650b9..7b21f5523 100644
--- a/docs/examples/features/demo_adaptive_accuracy.jl
+++ b/docs/examples/features/demo_adaptive_accuracy.jl
@@ -91,15 +91,15 @@ prob_ode1 = ODEProblem(trace_normalized!, u0_1, tspan1, param1)
 
 adaptive_solvers_1 = [
     (
-        "AdaptiveBoris (η=1/16)",
+        "Boris (η=1/16)",
         Boris(; safety = 1 / 16),
     ),
     (
-        "AdaptiveBoris (η=1/8)",
+        "Boris (η=1/8)",
         Boris(; safety = 1 / 8),
     ),
     (
-        "AdaptiveBoris (η=1/4)",
+        "Boris (η=1/4)",
         Boris(; safety = 1 / 4),
     ),
 ]
@@ -134,7 +134,7 @@ ax1b = Axis(
     f1[1, 2];
     xlabel = "Time [Gyroperiod]",
     ylabel = L"\Delta t / T_\mathrm{gyro}",
-    title = "AdaptiveBoris: Time Step History",
+    title = "Adaptive Boris: Time Step History",
 )
 
 for (i, (name, alg)) in enumerate(adaptive_solvers_1)
@@ -245,7 +245,7 @@ boris_steps = Int[]
 boris_errors = Float64[]
 
 for s in safety_values
-    alg = AdaptiveBoris(; safety = s)
+    alg = Boris(; safety = s)
     sol = TP.solve(prob_tp2, alg)[1]
     push!(boris_steps, length(sol.t))
     push!(boris_errors, max_velocity_error(sol))
@@ -288,7 +288,7 @@ ax2 = Axis(
 
 scatterlines!(
     ax2, boris_steps, boris_errors;
-    label = "AdaptiveBoris",
+    label = "Adaptive Boris",
     marker = :circle, linewidth = 2,
 )
 
@@ -335,11 +335,11 @@ prob_ode3 = ODEProblem(trace_normalized!, u0_2, tspan3, param3)
 
 adaptive_solvers_3 = [
     (
-        "AdaptiveBoris (η=0.1)",
+        "Boris (η=0.1)",
         Boris(; safety = 0.1),
     ),
     (
-        "AdaptiveBoris (η=0.2)",
+        "Boris (η=0.2)",
         Boris(; safety = 0.2),
     ),
 ]

From fac42da26043e73d7d75ea8c069ac7661fabccba Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 20:38:32 -0400
Subject: [PATCH 11/17] Fix memory allocation regression for the Boris fixed
 time step solver

---
 src/boris.jl | 24 ++++++++++++++++++++++--
 1 file changed, 22 insertions(+), 2 deletions(-)

diff --git a/src/boris.jl b/src/boris.jl
index 23fb3d6c3..8e5911175 100644
--- a/src/boris.jl
+++ b/src/boris.jl
@@ -507,11 +507,31 @@ end
         (v, r, dt, t, p) -> update_velocity_multistep(v, r, dt, t, alg.n, _get_val_N(alg), p)
     end
 
+    # Calculate exact nout for fixed-step to avoid memory regression
+    ttotal = tspan[2] - tspan[1]
+    nout_fixed = if !isnothing(dt)
+        nt = round(Int, abs(ttotal / dt))
+        nsteps = 0
+        if save_start; nsteps += 1; end
+        if save_everystep
+            steps = nt ÷ savestepinterval
+            last_is_step = (nt > 0) && (nt % savestepinterval == 0)
+            nsteps += steps
+            if !save_end && last_is_step; nsteps -= 1; end
+            if save_end && !last_is_step; nsteps += 1; end
+        elseif save_end
+            nsteps += 1
+        end
+        nsteps
+    else
+        1000 # Default capacity for adaptive Boris
+    end
+
     @inbounds for i in irange
         traj = Vector{SVector{vars_dim, T}}(undef, 0)
         tsave = Vector{typeof(tspan[1] + (isnothing(dt) ? 0.0 : dt))}(undef, 0)
-        sizehint!(traj, 1000)
-        sizehint!(tsave, 1000)
+        sizehint!(traj, nout_fixed)
+        sizehint!(tsave, nout_fixed)
 
         # set initial conditions for each trajectory i
         new_prob = prob.prob_func(prob, i, false)

From 05b909889f19f7087ab198bf25e08b56f0e0fa6a Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 20:47:45 -0400
Subject: [PATCH 12/17] formatting

---
 src/boris.jl | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/boris.jl b/src/boris.jl
index 8e5911175..c98fc0a7b 100644
--- a/src/boris.jl
+++ b/src/boris.jl
@@ -512,13 +512,19 @@ end
     nout_fixed = if !isnothing(dt)
         nt = round(Int, abs(ttotal / dt))
         nsteps = 0
-        if save_start; nsteps += 1; end
+        if save_start
+            nsteps += 1
+        end
         if save_everystep
             steps = nt ÷ savestepinterval
             last_is_step = (nt > 0) && (nt % savestepinterval == 0)
             nsteps += steps
-            if !save_end && last_is_step; nsteps -= 1; end
-            if save_end && !last_is_step; nsteps += 1; end
+            if !save_end && last_is_step
+                nsteps -= 1
+            end
+            if save_end && !last_is_step
+                nsteps += 1
+            end
         elseif save_end
             nsteps += 1
         end

From 51b5c28fa348fba272ddb9985187c7a32794b4f0 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 21:38:16 -0400
Subject: [PATCH 13/17] Reintroduce AdaptiveBoris for backward compatibility;
 Support backward tracing in OrdinaryDiffEqBoris

---
 .../src/boris_perform_step.jl                 | 13 +++++++---
 src/TestParticle.jl                           |  2 +-
 src/boris.jl                                  | 25 ++++++++++++++++++-
 3 files changed, 35 insertions(+), 5 deletions(-)

diff --git a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
index dd5bf0740..c0ad1a7e5 100644
--- a/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
+++ b/lib/OrdinaryDiffEqBoris/src/boris_perform_step.jl
@@ -49,7 +49,7 @@ end
     alg = integrator.alg
     if alg.safety > 0.0
         Bmag = norm(Bfunc(r_new, t + dt))
-        dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
+        dt_new = integrator.tdir * (2π * alg.safety) / (abs(q2m) * Bmag)
         set_proposed_dt!(integrator, dt_new)
     end
     return integrator.u
@@ -91,6 +91,13 @@ end
     integrator.u[4] = v_new[1]
     integrator.u[5] = v_new[2]
     integrator.u[6] = v_new[3]
+
+    alg = integrator.alg
+    if alg.safety > 0.0
+        Bmag = norm(Bfunc(r_new, t + dt))
+        dt_new = integrator.tdir * (2π * alg.safety) / (abs(q2m) * Bmag)
+        set_proposed_dt!(integrator, dt_new)
+    end
     return
 end
 
@@ -196,7 +203,7 @@ end
         q2m = p[1]
         Bfunc = p[4]
         Bmag = norm(Bfunc(r_new, t + dt))
-        dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
+        dt_new = integrator.tdir * (2π * alg.safety) / (abs(q2m) * Bmag)
         set_proposed_dt!(integrator, dt_new)
     end
     return integrator.u
@@ -238,7 +245,7 @@ end
         q2m = p[1]
         Bfunc = p[4]
         Bmag = norm(Bfunc(r_new, t + dt))
-        dt_new = (2π * alg.safety) / (abs(q2m) * Bmag)
+        dt_new = integrator.tdir * (2π * alg.safety) / (abs(q2m) * Bmag)
         set_proposed_dt!(integrator, dt_new)
     end
     return
diff --git a/src/TestParticle.jl b/src/TestParticle.jl
index 06a696890..8b2f57d67 100644
--- a/src/TestParticle.jl
+++ b/src/TestParticle.jl
@@ -37,7 +37,7 @@ export trace!, trace_relativistic!, trace_normalized!, trace_relativistic_normal
     get_gc_velocity, full_to_gc, gc_to_full
 export Proton, Electron, Ion
 export Maxwellian, BiMaxwellian, Kappa, BiKappa
-export AdaptiveHybrid, Boris, MultistepBoris, MultistepBoris2, MultistepBoris4, MultistepBoris6
+export AdaptiveHybrid, AdaptiveBoris, Boris, MultistepBoris, MultistepBoris2, MultistepBoris4, MultistepBoris6
 export get_gyrofrequency,
     get_gyroperiod, get_gyroradius, get_velocity, get_energy, get_mean_magnitude,
     energy2velocity, get_curvature_radius, get_adiabaticity,
diff --git a/src/boris.jl b/src/boris.jl
index c98fc0a7b..5a718a0b7 100644
--- a/src/boris.jl
+++ b/src/boris.jl
@@ -11,6 +11,20 @@ struct Boris{T}
 end
 Boris(; safety = 0.0) = Boris(safety)
 
+"""
+    AdaptiveBoris(; safety=0.1)
+
+The adaptive Boris method with adaptive time stepping based on local gyroperiod.
+Returns a `Boris` instance with the specified `safety` factor.
+"""
+struct AdaptiveBoris{T}
+    safety::T
+end
+function AdaptiveBoris(; safety = 0.1)
+    @warn "AdaptiveBoris is deprecated. Use Boris(safety=$safety) instead." maxlog=1
+    return AdaptiveBoris(safety)
+end
+
 """
     MultistepBoris{N}(; n=1, safety=0.0)
 
@@ -165,7 +179,7 @@ Trace particles using the Boris method with specified `prob`.
   - `batch_size::Int=max(1, trajectories ÷ nworkers())`: the number of trajectories to process per worker in `EnsembleDistributed` and `EnsembleSplitThreads`.
 """
 @inline function solve(
-        prob::TraceProblem, alg::Union{Boris, MultistepBoris} = Boris(),
+        prob::TraceProblem, alg::Union{Boris, MultistepBoris, AdaptiveBoris} = Boris(),
         ensemblealg::EA = EnsembleSerial();
         trajectories::Int = 1, savestepinterval::Int = 1, dt = nothing,
         isoutside::F = ODE_DEFAULT_ISOUTOFDOMAIN,
@@ -176,6 +190,11 @@ Trace particles using the Boris method with specified `prob`.
         n::Int = 1, N::Int = 2
     ) where {EA <: BasicEnsembleAlgorithm, F}
 
+    # Backward compatibility: promote AdaptiveBoris to Boris
+    if alg isa AdaptiveBoris
+        alg = Boris(alg.safety)
+    end
+
     # Backward compatibility: promote Boris to MultistepBoris if n > 1 or N != 2
     if alg isa Boris && (n > 1 || N != 2)
         if N ∉ (2, 4, 6)
@@ -184,6 +203,10 @@ Trace particles using the Boris method with specified `prob`.
         alg = MultistepBoris{N}(n = n, safety = alg.safety)
     end
 
+    if isnothing(dt) && alg.safety <= 0.0
+        throw(ArgumentError("Time step dt must be provided for fixed-step Boris solver."))
+    end
+
     if !isnothing(dt)
         if dt < eps(eltype(dt)) * 100
             throw(ArgumentError("Time step dt is too small."))

From 7844b7263ced30fd7d3c8e4dc527d635d899b62a Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 21:51:26 -0400
Subject: [PATCH 14/17] Formatting

---
 src/boris.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/boris.jl b/src/boris.jl
index 5a718a0b7..ce4bc3121 100644
--- a/src/boris.jl
+++ b/src/boris.jl
@@ -21,7 +21,7 @@ struct AdaptiveBoris{T}
     safety::T
 end
 function AdaptiveBoris(; safety = 0.1)
-    @warn "AdaptiveBoris is deprecated. Use Boris(safety=$safety) instead." maxlog=1
+    @warn "AdaptiveBoris is deprecated. Use Boris(safety=$safety) instead." maxlog = 1
     return AdaptiveBoris(safety)
 end
 

From a3b1bc73f68647571fa360a7595daff037016769 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 23:12:18 -0400
Subject: [PATCH 15/17] Fix backward tracing bug

---
 src/boris.jl       |  2 +-
 test/test_boris.jl | 20 ++++++++++++++++++++
 2 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/src/boris.jl b/src/boris.jl
index ce4bc3121..4286ab984 100644
--- a/src/boris.jl
+++ b/src/boris.jl
@@ -208,7 +208,7 @@ Trace particles using the Boris method with specified `prob`.
     end
 
     if !isnothing(dt)
-        if dt < eps(eltype(dt)) * 100
+        if abs(dt) < eps(eltype(dt)) * 100
             throw(ArgumentError("Time step dt is too small."))
         end
         ttotal = prob.tspan[2] - prob.tspan[1]
diff --git a/test/test_boris.jl b/test/test_boris.jl
index 301f6e0ac..5d7761791 100644
--- a/test/test_boris.jl
+++ b/test/test_boris.jl
@@ -248,6 +248,26 @@ using Distributed
         @test abs(sol_adaptive.u[end][4]) < 1.0e5 - 100
     end
 
+    @testset "Backward tracing" begin
+        # B = [0, 0, 0.01], E = 0
+        param = prepare(ZeroField(), uniform_B2, species = Proton)
+        u0 = SA[0.0, 0.0, 0.0, 1.0e5, 0.0, 0.0]
+        # Forward tspan
+        tspan_fw = (0.0, 1.0e-5)
+        dt = 1.0e-7
+        sol_fw = TP.solve(TraceProblem(u0, tspan_fw, param); dt)[1]
+        
+        # Backward tspan
+        u_end = sol_fw.u[end]
+        tspan_bw = (1.0e-5, 0.0)
+        prob_bw = TraceProblem(u_end, tspan_bw, param)
+        sol_bw = TP.solve(prob_bw; dt = -dt)[1]
+        
+        # Should return close to original IC
+        @test sol_bw.u[end][1:3] ≈ u0[1:3] atol = 1.0e-4
+        @test sol_bw.u[end][4:6] ≈ u0[4:6] atol = 1.0e-4
+    end
+
     @testset "Output saving flags" begin
         # Setup
         x0 = [0.0, 0.0, 0.0]

From 4b2ba280a98977f01c7b9b964721c0d85ca5439c Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Thu, 16 Apr 2026 23:16:39 -0400
Subject: [PATCH 16/17] Formatting

---
 test/test_boris.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_boris.jl b/test/test_boris.jl
index 5d7761791..7a700f392 100644
--- a/test/test_boris.jl
+++ b/test/test_boris.jl
@@ -256,13 +256,13 @@ using Distributed
         tspan_fw = (0.0, 1.0e-5)
         dt = 1.0e-7
         sol_fw = TP.solve(TraceProblem(u0, tspan_fw, param); dt)[1]
-        
+
         # Backward tspan
         u_end = sol_fw.u[end]
         tspan_bw = (1.0e-5, 0.0)
         prob_bw = TraceProblem(u_end, tspan_bw, param)
         sol_bw = TP.solve(prob_bw; dt = -dt)[1]
-        
+
         # Should return close to original IC
         @test sol_bw.u[end][1:3] ≈ u0[1:3] atol = 1.0e-4
         @test sol_bw.u[end][4:6] ≈ u0[4:6] atol = 1.0e-4

From be80053eb6d1454200bc8eaae8bbfd47b8678b46 Mon Sep 17 00:00:00 2001
From: Hongyang Zhou <hyzhou@umich.edu>
Date: Fri, 17 Apr 2026 11:52:41 -0400
Subject: [PATCH 17/17] refactor the Boris solvers

---
 src/boris.jl | 601 +++++++++++++++++++++++----------------------------
 1 file changed, 266 insertions(+), 335 deletions(-)

diff --git a/src/boris.jl b/src/boris.jl
index 4286ab984..f6b978104 100644
--- a/src/boris.jl
+++ b/src/boris.jl
@@ -1,12 +1,14 @@
 # Native particle pusher
 
+abstract type AbstractBoris end
+
 """
     Boris(; safety=0.0)
 
 The standard Boris method for particle pushing in electric and magnetic fields.
 When `safety > 0.0`, it uses adaptive time stepping based on local gyroperiod.
 """
-struct Boris{T}
+struct Boris{T} <: AbstractBoris
     safety::T
 end
 Boris(; safety = 0.0) = Boris(safety)
@@ -17,7 +19,7 @@ Boris(; safety = 0.0) = Boris(safety)
 The adaptive Boris method with adaptive time stepping based on local gyroperiod.
 Returns a `Boris` instance with the specified `safety` factor.
 """
-struct AdaptiveBoris{T}
+struct AdaptiveBoris{T} <: AbstractBoris
     safety::T
 end
 function AdaptiveBoris(; safety = 0.1)
@@ -33,7 +35,7 @@ The Multistep/Hyper Boris method of order `N`.
 `N` specifies the gyrophase correction order: 2 (standard), 4, or 6 (Hyper-Boris).
 When `safety > 0.0`, it uses adaptive time stepping.
 """
-struct MultistepBoris{N, T}
+struct MultistepBoris{N, T} <: AbstractBoris
     n::Int
     safety::T
 end
@@ -145,6 +147,252 @@ Update location in one timestep `dt`.
     return
 end
 
+@inline _get_val_N(::MultistepBoris{N}) where {N} = Val(N)
+
+"""
+    update_velocity_multistep(v, r, dt, t, n, ::Val{N}, param) where {N}
+
+Update velocity using the Multistep/Hyper Boris method.
+"""
+@muladd function update_velocity_multistep(v, r, dt, t, n::Int, ::Val{N}, param) where {N}
+    q2m, _, Efunc, Bfunc, _ = param
+    E = Efunc(r, t)
+    B = Bfunc(r, t)
+
+    # t_n and e_n vectors
+    factor = q2m * dt / (2 * n)
+
+    t_n = factor * B # (q/m * dt/(2n)) * B
+    e_n = factor * E # (q/m * dt/(2n)) * E
+
+    # Hyper Boris N-th order gyrophase correction
+    if N != 2
+        t_mag2 = sum(abs2, t_n)
+        if N == 4
+            f_N = 1 + t_mag2 / 3
+            e_corr_factor = -1 / 3
+        else # N == 6
+            f_N = 1 + t_mag2 / 3 + 2 * t_mag2 * t_mag2 / 15
+            e_corr_factor = -1 / 3 - 2 * t_mag2 / 15
+        end
+
+        e_dot_t = e_n ⋅ t_n
+        e_n = f_N * e_n + (e_corr_factor * e_dot_t) * t_n
+        t_n = f_N * t_n
+    end
+
+    t_n_mag2 = sum(abs2, t_n)
+    t_n_mag = sqrt(t_n_mag2)
+
+    # Calculate coefficients
+    if t_n_mag < TN_MAG_THRESHOLD
+        c_n1 = 1 - 2 * n * n * t_n_mag2
+        n_term1 = 2 * n
+        n_term3 = 4 * n * n * n
+        c_n2 = n_term1 - (n_term1 + n_term3) / 3 * t_n_mag2
+        c_n3 = 2 * n * n - (4 * n * n + 2 * n * n * n * n) / 3 * t_n_mag2
+        c_n6 = (n_term1 + n_term3) / 3
+    else
+        alpha_n = atan(t_n_mag)
+        n_alpha_n = n * alpha_n
+        sin_n_alpha, cos_n_alpha = sincos(n_alpha_n)
+        sin_2n_alpha = 2 * sin_n_alpha * cos_n_alpha
+        cos_2n_alpha = cos_n_alpha * cos_n_alpha - sin_n_alpha * sin_n_alpha
+
+        c_n1 = cos_2n_alpha
+        c_n2 = sin_2n_alpha / t_n_mag
+        c_n3 = 2 * sin_n_alpha * sin_n_alpha / t_n_mag2
+        c_n6 = (2 * n - c_n2) / t_n_mag2
+    end
+
+    c_n4 = c_n2
+    c_n5 = c_n3
+
+    v_dot_t = v ⋅ t_n
+    e_dot_t = e_n ⋅ t_n
+    v_cross_t = v × t_n
+    e_cross_t = e_n × t_n
+
+    v_new = c_n1 * v +
+        c_n2 * v_cross_t +
+        c_n3 * v_dot_t * t_n +
+        c_n4 * e_n +
+        c_n5 * e_cross_t +
+        c_n6 * e_dot_t * t_n
+
+    return v_new
+end
+
+@inline _update_velocity(v, r, dt, t, p, ::Boris) = update_velocity(v, r, dt, t, p)
+@inline function _update_velocity(v, r, dt, t, p, alg::MultistepBoris)
+    return update_velocity_multistep(v, r, dt, t, alg.n, _get_val_N(alg), p)
+end
+
+@inline @muladd function _boris_loop!(
+        traj, tsave, r, v, p, dt::T, tspan,
+        savestepinterval::Int, save_everystep::Bool, isoutside::F1,
+        alg, ::Val{SaveFields}, ::Val{SaveWork}, maxiters::Int
+    ) where {T, F1, SaveFields, SaveWork}
+    t = tspan[1]
+    ttotal = tspan[2] - tspan[1]
+    it = 1
+    retcode = ReturnCode.Success
+
+    C = zero(eltype(v))
+    if alg.safety > 0.0
+        q2m, _, _, Bfunc, _ = p
+        C = (2π * alg.safety * sign(ttotal)) / abs(q2m)
+        Bmag = norm(Bfunc(r, t))
+        dt = C / Bmag
+    end
+
+    # push velocity back in time by 1/2 dt
+    v = _update_velocity(v, r, -0.5 * dt, t, p, alg)
+
+    while abs(t - tspan[1]) < abs(ttotal)
+        if abs(t + dt - tspan[1]) >= abs(ttotal) - 100 * eps(eltype(v)(abs(ttotal)))
+            dt_step = tspan[2] - t
+            if dt_step != dt
+                # Resync v from `t - 0.5*dt` to `t - 0.5*dt_step`
+                v = _update_velocity(v, r, 0.5 * dt, t, p, alg)
+                v = _update_velocity(v, r, -0.5 * dt_step, t, p, alg)
+                dt = dt_step
+            end
+        end
+
+        if save_everystep && (it - 1) > 0 && (it - 1) % savestepinterval == 0
+            # Advance to t to get v_n
+            v_save = _update_velocity(v, r, 0.5 * dt, t, p, alg)
+            data = vcat(r, v_save)
+            push!(traj, _prepare_saved_data(data, p, t, Val(SaveFields), Val(SaveWork)))
+            push!(tsave, t)
+        end
+
+        v_new = _update_velocity(v, r, dt, t + 0.5 * dt, p, alg)
+        r_next = r + v_new * dt
+        t_next = t + dt
+
+        # NaN check
+        if any(isnan, r_next) || any(isnan, v_new)
+            retcode = ReturnCode.Unstable
+            break
+        end
+
+        if isoutside(vcat(r_next, v_new), p, t_next)
+            retcode = ReturnCode.Terminated
+            break
+        end
+
+        r, v, t = r_next, v_new, t_next
+
+        if alg.safety > 0.0 && abs(t - tspan[1]) < abs(ttotal)
+            q2m, _, _, Bfunc, _ = p
+            Bmag = norm(Bfunc(r, t))
+            dt_new = C / Bmag
+            # Resync v_{n+1/2}(dt) to v_{n+1/2}(dt_new)
+            v = _update_velocity(v, r, 0.5 * dt, t, p, alg)
+            v = _update_velocity(v, r, -0.5 * dt_new, t, p, alg)
+            dt = dt_new
+        end
+        if it > maxiters
+            retcode = ReturnCode.MaxIters
+            break
+        end
+        it += 1
+    end
+
+    return r, v, t, dt, retcode
+end
+
+@inline @muladd function boris_generic_loop!(
+        sols, prob::TraceProblem, irange, alg, dt,
+        savestepinterval::Int, isoutside::F1, save_start, save_end,
+        save_everystep, ::Val{SaveFields}, ::Val{SaveWork}, maxiters::Int
+    ) where {SaveFields, SaveWork, F1}
+    (; tspan, p, u0) = prob
+    T = eltype(u0)
+
+    vars_dim = 6
+    if SaveFields
+        vars_dim += 6
+    end
+    if SaveWork
+        vars_dim += 4
+    end
+
+    algorithm_name = alg isa Boris ? :boris : :multistep_boris
+
+    # Calculate exact nout for fixed-step to avoid memory regression
+    ttotal = tspan[2] - tspan[1]
+    nout_fixed = if !iszero(dt)
+        nt = round(Int, abs(ttotal / dt))
+        nsteps = 0
+        if save_start
+            nsteps += 1
+        end
+        if save_everystep
+            steps = nt ÷ savestepinterval
+            last_is_step = (nt > 0) && (nt % savestepinterval == 0)
+            nsteps += steps
+            if !save_end && last_is_step
+                nsteps -= 1
+            end
+            if save_end && !last_is_step
+                nsteps += 1
+            end
+        elseif save_end
+            nsteps += 1
+        end
+        nsteps
+    else
+        1000 # Default capacity for adaptive Boris
+    end
+
+    @inbounds for i in irange
+        traj = Vector{SVector{vars_dim, T}}(undef, 0)
+        tsave = Vector{typeof(tspan[1] + (iszero(dt) ? 0.0 : dt))}(undef, 0)
+        sizehint!(traj, nout_fixed)
+        sizehint!(tsave, nout_fixed)
+
+        # set initial conditions for each trajectory i
+        new_prob = prob.prob_func(prob, i, false)
+        u0_i = SVector{6, T}(new_prob.u0)
+        r = u0_i[SVector(1, 2, 3)]
+        v = u0_i[SVector(4, 5, 6)]
+
+        # If dt is not provided, it must be adaptive (safety > 0)
+        _dt = iszero(dt) ? zero(T) : T(dt)
+
+        if save_start
+            push!(traj, _prepare_saved_data(u0_i, p, tspan[1], Val(SaveFields), Val(SaveWork)))
+            push!(tsave, tspan[1])
+        end
+
+        r, v, t, _dt_final, retcode = _boris_loop!(
+            traj, tsave, r, v, p, _dt, tspan,
+            savestepinterval, save_everystep, isoutside, alg,
+            Val(SaveFields), Val(SaveWork), maxiters
+        )
+
+        should_save_final = save_end
+        if should_save_final && (isempty(tsave) || tsave[end] != t)
+            # v is at t - 0.5*_dt_final. To get v at t, advance by 0.5*_dt_final
+            v_save = _update_velocity(v, r, 0.5 * _dt_final, t, p, alg)
+            data = vcat(r, v_save)
+            push!(traj, _prepare_saved_data(data, p, t, Val(SaveFields), Val(SaveWork)))
+            push!(tsave, t)
+        end
+
+        alg_sol = algorithm_name
+        interp = LinearInterpolation(tsave, traj)
+        stats = nothing
+
+        sols[i] = build_solution(prob, alg_sol, tsave, traj; interp, retcode, stats)
+    end
+
+    return
+end
+
 """
 In-place cross product.
 """
@@ -164,24 +412,11 @@ end
         save_fields::Bool=false, save_work::Bool=false)
 
 Trace particles using the Boris method with specified `prob`.
-
-# keywords
-
-  - `trajectories::Int`: number of trajectories to trace.
-  - `dt::AbstractFloat`: time step for fixed-step Boris. If not provided, `alg.safety` must be > 0.
-  - `savestepinterval::Int`: saving output interval.
-  - `isoutside::Function`: pinpointing impact or checking boundaries.
-  - `save_start::Bool=true`: save the initial condition.
-  - `save_end::Bool=true`: save the final condition.
-  - `save_everystep::Bool=true`: save the state at every `savestepinterval`.
-  - `save_fields::Bool=false`: save the electric and magnetic fields.
-  - `save_work::Bool=false`: save the work done by the electric field.
-  - `batch_size::Int=max(1, trajectories ÷ nworkers())`: the number of trajectories to process per worker in `EnsembleDistributed` and `EnsembleSplitThreads`.
 """
 @inline function solve(
         prob::TraceProblem, alg::Union{Boris, MultistepBoris, AdaptiveBoris} = Boris(),
         ensemblealg::EA = EnsembleSerial();
-        trajectories::Int = 1, savestepinterval::Int = 1, dt = nothing,
+        trajectories::Int = 1, savestepinterval::Int = 1, dt = 0.0,
         isoutside::F = ODE_DEFAULT_ISOUTOFDOMAIN,
         save_start::Bool = true, save_end::Bool = true, save_everystep::Bool = true,
         save_fields::Bool = false, save_work::Bool = false, maxiters::Int = 1_000_000,
@@ -203,11 +438,11 @@ Trace particles using the Boris method with specified `prob`.
         alg = MultistepBoris{N}(n = n, safety = alg.safety)
     end
 
-    if isnothing(dt) && alg.safety <= 0.0
+    if iszero(dt) && alg.safety <= 0.0
         throw(ArgumentError("Time step dt must be provided for fixed-step Boris solver."))
     end
 
-    if !isnothing(dt)
+    if !iszero(dt)
         if abs(dt) < eps(eltype(dt)) * 100
             throw(ArgumentError("Time step dt is too small."))
         end
@@ -229,11 +464,11 @@ function solve(prob::TraceProblem, ensemblealg::BasicEnsembleAlgorithm; kwargs..
 end
 
 function _dispatch_boris!(
-        sols, prob::TraceProblem, irange, alg, dt, savestepinterval, isoutside::F,
-        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
-        maxiters
+        sols, prob::TraceProblem, irange, alg, dt, savestepinterval::Int,
+        isoutside::F, save_start, save_end, save_everystep, ::Val{SaveFields},
+        ::Val{SaveWork}, maxiters::Int
     ) where {SaveFields, SaveWork, F}
-    return _generic_boris!(
+    return boris_generic_loop!(
         sols, prob, irange, alg, dt, savestepinterval, isoutside,
         save_start, save_end, save_everystep, Val(SaveFields), Val(SaveWork),
         maxiters
@@ -241,9 +476,9 @@ function _dispatch_boris!(
 end
 
 @inline function _solve(
-        ::EnsembleSerial, prob::TraceProblem, trajectories, alg, dt, savestepinterval,
-        isoutside::F, save_start, save_end, save_everystep,
-        ::Val{SaveFields}, ::Val{SaveWork}, maxiters, batch_size
+        ::EnsembleSerial, prob::TraceProblem, trajectories, alg, dt,
+        savestepinterval::Int, isoutside::F, save_start, save_end, save_everystep,
+        ::Val{SaveFields}, ::Val{SaveWork}, maxiters::Int, batch_size
     ) where {SaveFields, SaveWork, F}
     sols = _prepare(
         prob, trajectories, alg, dt, save_start, save_end, save_everystep,
@@ -260,9 +495,9 @@ end
 end
 
 @inline function _solve(
-        ::EnsembleThreads, prob::TraceProblem, trajectories, alg, dt, savestepinterval,
-        isoutside::F, save_start, save_end, save_everystep,
-        ::Val{SaveFields}, ::Val{SaveWork}, maxiters, batch_size
+        ::EnsembleThreads, prob::TraceProblem, trajectories, alg, dt,
+        savestepinterval::Int, isoutside::F, save_start, save_end, save_everystep,
+        ::Val{SaveFields}, ::Val{SaveWork}, maxiters::Int, batch_size
     ) where {SaveFields, SaveWork, F}
     sols = _prepare(
         prob, trajectories, alg, dt, save_start, save_end, save_everystep,
@@ -281,24 +516,6 @@ end
     return sols
 end
 
-"""
-    _solve_single_boris(prob, i, ...)
-
-Solve a single trajectory `i` of `prob` for use with `EnsembleDistributed`.
-
-`_generic_boris!` uses the loop index `i` for two purposes simultaneously:
-applying `prob_func(prob, i, false)` to select per-particle initial conditions,
-and storing the result at `sols[i]`. For a distributed worker handling only one
-trajectory at a time, a 1-element `local_sols` would be out of bounds if `i > 1`
-were passed directly. Decoupling the two uses would require threading a storage
-offset through the entire `_dispatch_boris!` → `_generic_boris!` call chain.
-
-Instead, we pre-apply `prob_func` to get the correct IC for trajectory `i` and
-wrap the result in a fresh `TraceProblem` with the default (identity) `prob_func`,
-so `_generic_boris!` can safely iterate `1:1` without applying `prob_func` again.
-The `TraceProblem` construction is a negligible struct copy relative to the
-simulation cost and the serialization overhead inherent in `pmap`.
-"""
 function _solve_single_boris(
         prob::TraceProblem, i, alg, dt, savestepinterval, isoutside::F,
         save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
@@ -307,7 +524,7 @@ function _solve_single_boris(
     new_prob = prob.prob_func(prob, i, false)
     single_prob = TraceProblem(new_prob.u0, new_prob.tspan, new_prob.p)
     sol_type = _get_sol_type(
-        single_prob, isnothing(dt) ? zero(eltype(single_prob.tspan)) : dt,
+        single_prob, iszero(dt) ? zero(eltype(single_prob.tspan)) : dt,
         Val(SaveFields), Val(SaveWork)
     )
     local_sols = Vector{sol_type}(undef, 1)
@@ -338,9 +555,6 @@ end
         isoutside::F, save_start, save_end, save_everystep,
         ::Val{SaveFields}, ::Val{SaveWork}, maxiters, batch_size
     ) where {SaveFields, SaveWork, F}
-    # _solve_single_boris wraps each trajectory in a fresh TraceProblem(u0, tspan, p)
-    # with DEFAULT_PROB_FUNC. We get a sample problem from prob_func to ensure the
-    # u0 type (uType) matches.
     sample_prob = prob.prob_func(prob, 1, false)
     dummy_prob = TraceProblem(sample_prob.u0, sample_prob.tspan, sample_prob.p)
     sol_type = _get_sol_type(dummy_prob, zero(eltype(dummy_prob.tspan)), Val(SaveFields), Val(SaveWork))
@@ -365,7 +579,6 @@ function _get_sol_type(prob, dt, ::Val{SaveFields}, ::Val{SaveWork}) where {Save
     tspan = prob.tspan
     T_t = typeof(tspan[1] + dt)
     t = Vector{T_t}(undef, 0)
-    # Force u to be Vector{SVector{6, T}} as used in _boris!
     T = eltype(u0)
 
     n_vars = 6
@@ -384,19 +597,12 @@ function _get_sol_type(prob, dt, ::Val{SaveFields}, ::Val{SaveWork}) where {Save
     return typeof(sol)
 end
 
-"""
-Prepare for advancing.
-"""
 function _prepare(
         prob::TraceProblem, trajectories, alg, dt,
         save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork}
     ) where {SaveFields, SaveWork}
-    if !isnothing(dt) && abs(dt) < 10 * eps(typeof(dt))
-        throw(ArgumentError("time step dt is too small, violating min_dt = 10 * eps(typeof(dt))"))
-    end
-
     sol_type = _get_sol_type(
-        prob, isnothing(dt) ? zero(eltype(prob.tspan)) : dt,
+        prob, iszero(dt) ? zero(eltype(prob.tspan)) : dt,
         Val(SaveFields), Val(SaveWork)
     )
     sols = Vector{sol_type}(undef, trajectories)
@@ -430,275 +636,6 @@ end
     return data
 end
 
-@inline @muladd function _boris_loop!(
-        traj, tsave, r, v, p, dt, tspan,
-        savestepinterval, save_everystep, isoutside::F1, velocity_updater::F2,
-        ::Val{SaveFields}, ::Val{SaveWork}, alg, maxiters
-    ) where {F1, F2, SaveFields, SaveWork}
-    t = tspan[1]
-    ttotal = tspan[2] - tspan[1]
-    it = 1
-    retcode = ReturnCode.Success
-
-    C = zero(eltype(v))
-    if alg.safety > 0.0
-        q2m, _, _, Bfunc, _ = p
-        C = (2π * alg.safety * sign(ttotal)) / abs(q2m)
-        Bmag = norm(Bfunc(r, t))
-        dt = C / Bmag
-    end
-
-    # push velocity back in time by 1/2 dt
-    v = velocity_updater(v, r, -0.5 * dt, t, p)
-
-    while abs(t - tspan[1]) < abs(ttotal)
-        if abs(t + dt - tspan[1]) >= abs(ttotal) - 100 * eps(eltype(v)(abs(ttotal)))
-            dt_step = tspan[2] - t
-            if dt_step != dt
-                # Resync v from `t - 0.5*dt` to `t - 0.5*dt_step`
-                v = velocity_updater(v, r, 0.5 * dt, t, p)
-                v = velocity_updater(v, r, -0.5 * dt_step, t, p)
-                dt = dt_step
-            end
-        end
-
-        if save_everystep && (it - 1) > 0 && (it - 1) % savestepinterval == 0
-            # Advance to t to get v_n
-            v_save = velocity_updater(v, r, 0.5 * dt, t, p)
-            data = vcat(r, v_save)
-            push!(traj, _prepare_saved_data(data, p, t, Val(SaveFields), Val(SaveWork)))
-            push!(tsave, t)
-        end
-
-        v_new = velocity_updater(v, r, dt, t + 0.5 * dt, p)
-        r_next = r + v_new * dt
-        t_next = t + dt
-
-        # NaN check
-        if any(isnan, r_next) || any(isnan, v_new)
-            retcode = ReturnCode.Unstable
-            break
-        end
-
-        if isoutside(vcat(r_next, v_new), p, t_next)
-            retcode = ReturnCode.Terminated
-            break
-        end
-
-        r, v, t = r_next, v_new, t_next
-
-        if alg.safety > 0.0 && abs(t - tspan[1]) < abs(ttotal)
-            q2m, _, _, Bfunc, _ = p
-            Bmag = norm(Bfunc(r, t))
-            dt_new = C / Bmag
-            # Resync v_{n+1/2}(dt) to v_{n+1/2}(dt_new)
-            v = velocity_updater(v, r, 0.5 * dt, t, p)
-            v = velocity_updater(v, r, -0.5 * dt_new, t, p)
-            dt = dt_new
-        end
-        if it > maxiters
-            retcode = ReturnCode.MaxIters
-            break
-        end
-        it += 1
-    end
-
-    return r, v, t, dt, retcode
-end
-
-@inline @muladd function _generic_boris!(
-        sols, prob::TraceProblem, irange, alg, dt, savestepinterval, isoutside::F1,
-        save_start, save_end, save_everystep, ::Val{SaveFields}, ::Val{SaveWork},
-        maxiters
-    ) where {SaveFields, SaveWork, F1}
-    (; tspan, p, u0) = prob
-    T = eltype(u0)
-
-    vars_dim = 6
-    if SaveFields
-        vars_dim += 6
-    end
-    if SaveWork
-        vars_dim += 4
-    end
-
-    algorithm_name = alg isa Boris ? :boris : :multistep_boris
-
-    velocity_updater = if alg isa Boris
-        update_velocity
-    else
-        (v, r, dt, t, p) -> update_velocity_multistep(v, r, dt, t, alg.n, _get_val_N(alg), p)
-    end
-
-    # Calculate exact nout for fixed-step to avoid memory regression
-    ttotal = tspan[2] - tspan[1]
-    nout_fixed = if !isnothing(dt)
-        nt = round(Int, abs(ttotal / dt))
-        nsteps = 0
-        if save_start
-            nsteps += 1
-        end
-        if save_everystep
-            steps = nt ÷ savestepinterval
-            last_is_step = (nt > 0) && (nt % savestepinterval == 0)
-            nsteps += steps
-            if !save_end && last_is_step
-                nsteps -= 1
-            end
-            if save_end && !last_is_step
-                nsteps += 1
-            end
-        elseif save_end
-            nsteps += 1
-        end
-        nsteps
-    else
-        1000 # Default capacity for adaptive Boris
-    end
-
-    @inbounds for i in irange
-        traj = Vector{SVector{vars_dim, T}}(undef, 0)
-        tsave = Vector{typeof(tspan[1] + (isnothing(dt) ? 0.0 : dt))}(undef, 0)
-        sizehint!(traj, nout_fixed)
-        sizehint!(tsave, nout_fixed)
-
-        # set initial conditions for each trajectory i
-        new_prob = prob.prob_func(prob, i, false)
-        u0_i = SVector{6, T}(new_prob.u0)
-        r = u0_i[SVector(1, 2, 3)]
-        v = u0_i[SVector(4, 5, 6)]
-
-        # If dt is not provided, it must be adaptive (safety > 0)
-        _dt = isnothing(dt) ? zero(T) : T(dt)
-
-        if save_start
-            push!(traj, _prepare_saved_data(u0_i, p, tspan[1], Val(SaveFields), Val(SaveWork)))
-            push!(tsave, tspan[1])
-        end
-
-        r, v, t, _dt_final, retcode = _boris_loop!(
-            traj, tsave, r, v, p, _dt, tspan,
-            savestepinterval, save_everystep, isoutside, velocity_updater,
-            Val(SaveFields), Val(SaveWork), alg, maxiters
-        )
-
-        should_save_final = save_end
-        if should_save_final && (isempty(tsave) || tsave[end] != t)
-            # v is at t - 0.5*_dt_final. To get v at t, advance by 0.5*_dt_final
-            v_save = velocity_updater(v, r, 0.5 * _dt_final, t, p)
-            data = vcat(r, v_save)
-            push!(traj, _prepare_saved_data(data, p, t, Val(SaveFields), Val(SaveWork)))
-            push!(tsave, t)
-        end
-
-        alg_sol = algorithm_name
-        interp = LinearInterpolation(tsave, traj)
-        stats = nothing
-
-        sols[i] = build_solution(prob, alg_sol, tsave, traj; interp, retcode, stats)
-    end
-
-    return
-end
-
-_get_val_N(::MultistepBoris{N}) where {N} = Val{N}()
-
-"""
-    update_velocity_multistep(v, r, dt, t, n, N, param)
-    update_velocity_multistep(v, r, dt, t, n, ::Val{N}, param)
-
-Update velocity using the Multistep/Hyper Boris method, returning the new velocity as an SVector.
-`n` specifies the number of subcycles.
-`N` specifies the gyrophase correction order. When N=2, it corresponds to the Multicycle solver. When N=4 or N=6, it is the Hyper Boris solver.
-Reference: [Zenitani & Kato 2025](https://arxiv.org/abs/2505.02270)
-"""
-@inline function update_velocity_multistep(v, r, dt, t, n::Int, N::Int, param)
-    return update_velocity_multistep(v, r, dt, t, n, Val(N), param)
-end
-
-@muladd function update_velocity_multistep(v, r, dt, t, n::Int, ::Val{N}, param) where {N}
-    q2m, _, Efunc, Bfunc, _ = param
-    E = Efunc(r, t)
-    B = Bfunc(r, t)
-
-    # t_n and e_n vectors
-    factor = q2m * dt / (2 * n)
-
-    t_n = factor * B # (q/m * dt/(2n)) * B
-    e_n = factor * E # (q/m * dt/(2n)) * E
-
-    # Hyper Boris N-th order gyrophase correction
-    if N != 2
-        t_mag2 = sum(abs2, t_n)
-        if N == 4
-            f_N = 1 + t_mag2 / 3
-            e_corr_factor = -1 / 3
-        else # N == 6
-            f_N = 1 + t_mag2 / 3 + 2 * t_mag2 * t_mag2 / 15
-            e_corr_factor = -1 / 3 - 2 * t_mag2 / 15
-        end
-
-        e_dot_t = e_n ⋅ t_n
-        e_n = f_N * e_n + (e_corr_factor * e_dot_t) * t_n
-        t_n = f_N * t_n
-    end
-
-    t_n_mag2 = sum(abs2, t_n)
-    t_n_mag = sqrt(t_n_mag2)
-
-    # Calculate coefficients
-    # Check for small t_n to avoid division by zero or precision loss
-    if t_n_mag < TN_MAG_THRESHOLD
-        # Taylor expansion limits as t_n -> 0
-        c_n1 = 1 - 2 * n * n * t_n_mag2
-
-        n_term1 = 2 * n
-        n_term3 = 4 * n * n * n
-
-        c_n2 = n_term1 - (n_term1 + n_term3) / 3 * t_n_mag2
-        c_n3 = 2 * n * n - (4 * n * n + 2 * n * n * n * n) / 3 * t_n_mag2
-        c_n6 = (n_term1 + n_term3) / 3
-    else
-        alpha_n = atan(t_n_mag)
-        n_alpha_n = n * alpha_n
-        sin_n_alpha, cos_n_alpha = sincos(n_alpha_n)
-        sin_2n_alpha = 2 * sin_n_alpha * cos_n_alpha
-        cos_2n_alpha = cos_n_alpha * cos_n_alpha - sin_n_alpha * sin_n_alpha
-
-        c_n1 = cos_2n_alpha
-        c_n2 = sin_2n_alpha / t_n_mag
-        c_n3 = 2 * sin_n_alpha * sin_n_alpha / t_n_mag2
-        c_n6 = (2 * n - c_n2) / t_n_mag2
-    end
-
-    c_n4 = c_n2
-    c_n5 = c_n3
-
-    v_dot_t = v ⋅ t_n
-    e_dot_t = e_n ⋅ t_n
-
-    v_cross_t = v × t_n
-    e_cross_t = e_n × t_n
-
-    # Update velocity
-    # Equation 39:
-    # v_new = c_n1*v + c_n2*(v x t_n) + c_n3*(v . t_n)t_n + c_n4*e_n + c_n5*(e_n x t_n) + c_n6*(e_n . t_n)t_n
-    v_new = c_n1 * v +
-        c_n2 * v_cross_t +
-        c_n3 * v_dot_t * t_n +
-        c_n4 * e_n +
-        c_n5 * e_cross_t +
-        c_n6 * e_dot_t * t_n
-
-    return v_new
-end
-
-
-"""
-    get_fields(sol::AbstractODESolution)
-
-Return the electric and magnetic fields from the solution `sol`.
-"""
 function get_fields(sol::AbstractODESolution)
     Efunc, Bfunc = _get_field_funcs(sol.prob)
 
@@ -709,7 +646,6 @@ function get_fields(sol::AbstractODESolution)
 end
 
 function _get_field_funcs(prob::TraceGCProblem)
-    # p = (q, q2m, μ, Efunc, Bfunc)
     p = prob.p
     return p[4], p[5]
 end
@@ -719,11 +655,6 @@ function _get_field_funcs(prob)
     return get_EField(p), get_BField(p)
 end
 
-"""
-    get_work(sol::AbstractODESolution)
-
-Return the work done by the electric field from the solution `sol`.
-"""
 function get_work(sol::AbstractODESolution)
     return _get_work(sol, sol.prob)
 end