refactor Boris kernel

[henry2004y]

• Remove temporary field arrays
• Reuse CPU buffers for saving states
• Optimize KernelAbstraction usage for CPU: benchmark in the suite shows 10x improvement
• Use FieldInterpolator types instead of function capture for the interpolation functions, which can be adapted to both CPU and GPU
• Add multithreading support for the kernel Boris solver
• Refactor isoutofdomain and velocity_updater function calls in the native Boris solver to make them type stable. This solves the large allocations in the native solver.
• Clean up package imports

CPU Backend vs Original

Comparing Script

using TestParticle
import TestParticle: solve
import SciMLBase: remake
using KernelAbstractions
using StaticArrays
using BenchmarkTools
using Printf
using Random
using LinearAlgebra

# Ensure we are using the same BLAS threads settings if that mattered, but here we are primarily scalar.
# The user will run this with julia -t 1

println("="^70)
println("Benchmark: Original (Serial) vs KA (CPU Backend)")
println("Threads: ", Threads.nthreads())
println("="^70)

# Setup
uniform_B(x) = SA[0.0, 0.0, 1.0e-8]
uniform_E(x) = SA[0.0, 0.0, 0.0]
param = prepare(uniform_E, uniform_B; species = Proton)

x0 = [0.0, 0.0, 0.0]
v0 = [1.0e5, 0.0, 0.0]
stateinit = [x0..., v0...]

println("\n## Verification")
# Run a short simulation to verify equivalence
tspan_verify = (0.0, 1.0e-5)
dt = 1.0e-9
prob_verify = TraceProblem(stateinit, tspan_verify, param)

println("Running Original Solver (Verification)...")
sol_orig = solve(prob_verify; dt = dt)

println("Running KA Solver (Verification)...")
backend_cpu = KernelAbstractions.CPU()
sol_ka = solve(prob_verify, backend_cpu; dt = dt)

# Compare final states
u_orig = sol_orig[1].u[end]
u_ka = sol_ka[1].u[end]
diff = norm(u_orig - u_ka)
println("Difference in final state: $diff")

if diff > 1.0e-12
    error("Verification FAILED: Results differ by $diff!")
else
    println("Verification PASSED: Results are equivalent.")
end

# Benchmark 1: Single Particle, Long Trace
# -----------------------------------------------------------------------------
println("\n## Benchmark 1: Single Particle, Long Trace (Low overhead check)")
println("Time span: 1ms, dt=1ns => 1,000,000 steps")
tspan_long = (0.0, 1.0e-3)
dt = 1.0e-9
prob_long = TraceProblem(stateinit, tspan_long, param)

println("\nRunning Original Solver (1 particle)...")
# Bench
b_orig_long = @benchmark solve($prob_long; dt = $dt, savestepinterval = 100000)
display(b_orig_long)

println("\nRunning KA Solver (CPU) (1 particle)...")
backend = KernelAbstractions.CPU()
# Bench
b_ka_long = @benchmark solve($prob_long, $backend; dt = $dt, trajectories = 1, savestepinterval = 100000)
display(b_ka_long)

# Benchmark 2: Ensemble (Throughput)
# -----------------------------------------------------------------------------
println("\n" * "-"^70)
println("## Benchmark 2: Ensemble 1000 Particles")
println("Time span: 10us, dt=1ns => 10,000 steps")
tspan_multi = (0.0, 1.0e-5)
n_particles = 1000

prob_func(prob, i, repeat) = remake(prob; u0 = [prob.u0[1:3]..., (i / 1000.0) * 1.0e5, 0.0, 0.0])
prob_multi = TraceProblem(stateinit, tspan_multi, param; prob_func = prob_func)

println("\nRunning Original Solver ($n_particles particles)...")
# Bench
b_orig_multi = @benchmark solve($prob_multi; dt = $dt, trajectories = $n_particles, savestepinterval = 10000)
display(b_orig_multi)

println("\nRunning KA Solver (CPU) ($n_particles particles)...")
# Bench
b_ka_multi = @benchmark solve($prob_multi, $backend; dt = $dt, trajectories = $n_particles, savestepinterval = 10000)
display(b_ka_multi)

# Comparison Summary
# -----------------------------------------------------------------------------
println("\n" * "="^70)
println("Summary Results (Median Times):")

t_orig_long = median(b_orig_long).time / 1.0e6 # ms
t_ka_long = median(b_ka_long).time / 1.0e6 # ms
println(@sprintf("Single Particle Long Trace: Original=%.2f ms, KA=%.2f ms, Factor=%.2fx", t_orig_long, t_ka_long, t_ka_long / t_orig_long))

t_orig_multi = median(b_orig_multi).time / 1.0e6 # ms
t_ka_multi = median(b_ka_multi).time / 1.0e6 # ms
println(@sprintf("Ensemble %d Short Trace:    Original=%.2f ms, KA=%.2f ms, Factor=%.2fx", n_particles, t_orig_multi, t_ka_multi, t_ka_multi / t_orig_multi))

if minimum(b_ka_multi).allocs > minimum(b_orig_multi).allocs * 1.5
    println("WARNING: KA version has significantly more allocations!")
    println("Original Allocs: ", minimum(b_orig_multi).allocs)
    println("KA Allocs:       ", minimum(b_ka_multi).allocs)
end

We verified that the results are identical from the two methods. The AbstractKernel based Boris solvers somehow shows crazy performance improvements:

Benchmark Case	Description	Original Solver (Serial)	KA Solver (CPU)	Factor
Single Particle Long Trace	1 particle, 1ms duration (1,000,000 steps).	~6.33 ms Mem: 928 bytes Allocs: 6	~6.40 ms Mem: 1.88 KiB Allocs: 38	~1.01x (KA is roughly equal)
Ensemble Short Trace	1000 particles, 10µs duration (10,000 steps each). Total 10,000,000 steps.	~49.06 ms Mem: 1.32 MiB Allocs: 34,003	~50.23 ms Mem: 1.01 MiB Allocs: 20,052	~1.02x (KA is roughly equal)

Note on Original Solver Performance: The original Boris solver previously shows significantly higher memory allocations (~80 million) in the ensemble case. It is related to the type instability in the function call. Now it's fixed.

CUDA Backend vs Original

CUDA script

# GPU vs CPU Boris Solver Benchmark (CUDA Version)

using TestParticle
import TestParticle: solve
import SciMLBase: remake
using KernelAbstractions
import KernelAbstractions as KA
using CUDA
using StaticArrays
using Printf

println("="^70)
println("GPU Boris Solver Benchmarks (CUDA Backend)")
println("="^70)

if !CUDA.functional()
    println("Error: CUDA is not functional.")
    exit(1)
end

# Setup
uniform_B(x) = SA[0.0, 0.0, 1.0e-8]
uniform_E(x) = SA[0.0, 0.0, 0.0]
param = prepare(uniform_E, uniform_B; species = Proton)

using LinearAlgebra

backend_cuda = CUDA.CUDABackend()

# Benchmark 1: Single particle
println("\n## Benchmark 1: Single Particle, Long Tracing")
println("-"^70)

x0 = [0.0, 0.0, 0.0]
v0 = [1.0e5, 0.0, 0.0]
stateinit = [x0..., v0...]

# Verification
println("\n## Verification")
tspan_verify = (0.0, 1.0e-5)
dt_verify = 1.0e-9
prob_verify = TraceProblem(stateinit, tspan_verify, param)

println("Running Original Solver (Verification)...")
sol_orig = solve(prob_verify; dt = dt_verify)

println("Running KA Solver (CUDA) (Verification)...")
sol_ka = solve(prob_verify, backend_cuda; dt = dt_verify)

u_orig = sol_orig[1].u[end]
u_ka = sol_ka[1].u[end]
diff = norm(u_orig - u_ka)
println("Difference in final state: $diff")

if diff > 1.0e-12
    error("Verification FAILED: Results differ by $diff!")
else
    println("Verification PASSED: Results are equivalent.")
end

tspan_long = (0.0, 1.0e-3)
dt = 1.0e-9

prob_long = TraceProblem(stateinit, tspan_long, param)

println("Running CPU solver...")
cpu_time = @elapsed begin
    sol_cpu = solve(prob_long; dt, savestepinterval = 100000)
end
println(@sprintf("CPU Time: %.4f seconds", cpu_time))

println("\nRunning GPU solver with CUDA backend...")
gpu_cuda_time = @elapsed begin
    sol_gpu_cuda = solve(prob_long, backend_cuda; dt, trajectories = 1, savestepinterval = 100000)
end
println(@sprintf("GPU (CUDA backend) Time: %.4f seconds", gpu_cuda_time))
println(@sprintf("Speedup: %.2fx", cpu_time / gpu_cuda_time))

# Benchmark 2: Multi-particle
println("\n## Benchmark 2: Multi-Particle Ensemble")
println("-"^70)

tspan_multi = (0.0, 1.0e-5)
n_particles_list = [1000, 10000, 100000]

for n_particles in n_particles_list
    println(@sprintf("\n### %d Particles", n_particles))

    # Define prob_func closure capturing n_particles
    prob_func(prob, i, repeat) = remake(prob; u0 = [prob.u0[1:3]..., (i / n_particles) * 1.0e5, 0.0, 0.0])
    prob_multi = TraceProblem(stateinit, tspan_multi, param; prob_func = prob_func)

    print("  CPU: ")
    cpu_multi_time = @elapsed begin
        sols_cpu = solve(prob_multi; dt, trajectories = n_particles, savestepinterval = 10000)
    end
    println(@sprintf("%.4f s", cpu_multi_time))

    print("  GPU (CUDA backend): ")
    gpu_multi_time = @elapsed begin
        sols_gpu = solve(prob_multi, backend_cuda; dt, trajectories = n_particles, savestepinterval = 10000)
    end
    println(@sprintf("%.4f s", gpu_multi_time))
    println(@sprintf("  Speedup: %.2fx", cpu_multi_time / gpu_multi_time))
end

# Benchmark 4: Numerical
println("\n## Benchmark 4: Numerical Field (Interpolated)")
println("-"^70)

x = range(-0.1, 0.1, length = 10)
y = range(-0.1, 0.1, length = 10)
z = range(-1.0, 1.0, length = 10)
B_data = [SA[0.0, 0.0, 1.0e-8] for xi in x, yi in y, zi in z]
E_data = [SA[0.0, 0.0, 0.0] for xi in x, yi in y, zi in z]
param_num = prepare(x, y, z, E_data, B_data; species = Proton)

println("\n### 10000 Particles (Numerical)")
# Use slightly more particles than CPU bench to respect GPU parallelism
n_num = 10000
prob_func_num(prob, i, repeat) = remake(prob; u0 = [prob.u0[1:3]..., (i / n_num) * 1.0e5, 0.0, 0.0])
prob_num_multi = TraceProblem(stateinit, tspan_multi, param_num; prob_func = prob_func_num)

print("  CPU: ")
cpu_num_multi_time = @elapsed begin
    sol_cpu_num = solve(prob_num_multi; dt, trajectories = n_num, savestepinterval = 100)
end
println(@sprintf("%.4f s", cpu_num_multi_time))

print("  GPU (CUDA backend): ")
gpu_num_multi_time = @elapsed begin
    sol_gpu_num = solve(prob_num_multi, backend_cuda; dt, trajectories = n_num, savestepinterval = 100)
end
println(@sprintf("%.4f s", gpu_num_multi_time))
println(@sprintf("  Speedup: %.2fx", cpu_num_multi_time / gpu_num_multi_time))

Note: GPU times include some compilation/overhead in the first run, especially for the single particle case. Speedups are most significant for large ensembles, but apparently still not ideal.

Benchmark Case	Description	Original Solver (CPU)	KA Solver (CUDA)	Speedup
Single Particle	1 particle, 1ms duration.	~0.011 s	~20.6 s (Compilation)	0.00x
Ensemble 1k	1000 particles, 10µs duration.	~0.52 s	~0.53 s	0.97x
Ensemble 10k	10,000 particles, 10µs duration.	~0.64 s	~0.22 s	2.9x
Ensemble 100k	100,000 particles, 10µs duration.	~5.61 s	~0.54 s	10.4x
Numerical Field	10,000 particles, 10µs (Numerical).	~4.54 s	~2.17 s	2.1x

---

Updated performance comparison:

Particles	Threads	Native (Serial)	Native (Threaded)	KA (Serial)	KA (Threaded)
1,000	1	49.05 ms	49.03 ms	44.28 ms	48.22 ms
1,000	2	49.23 ms	25.23 ms	44.33 ms	22.02 ms
1,000	4	49.16 ms	12.68 ms	44.36 ms	10.67 ms
10,000	1	491.30 ms	491.75 ms	504.56 ms	498.62 ms
10,000	2	491.40 ms	250.45 ms	473.18 ms	265.85 ms
10,000	4	497.61 ms	125.24 ms	472.31 ms	133.42 ms

Updated memory allocations:

Solver	Particles	Memory (Serial)	Memory (Threaded)
Native	1,000	1.32 MiB	1.32 MiB
KA (CPU)	1,000	1.01 MiB	1.13 MiB
Native	10,000	13.20 MiB	13.20 MiB
KA (CPU)	10,000	10.08 MiB	11.23 MiB

Updated scaling analysis

Particles	Threads	Native Speedup	KA Speedup
1,000	2	1.95x	2.01x
1,000	4	3.88x	4.16x
10,000	2	1.96x	1.78x
10,000	4	3.97x	3.54x

[gemini-code-assist]

Summary of Changes

Hello @henry2004y, I'm Gemini Code Assist¹! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed!

This pull request refactors the Boris kernel implementation by shifting from a pre-calculation and storage approach for electric and magnetic fields to an on-the-fly evaluation directly within the GPU kernels. This change streamlines field handling, significantly reduces memory overhead by eliminating temporary field arrays, and optimizes CPU buffer usage for state saving. The overall impact is a cleaner, more efficient, and potentially faster simulation process, particularly beneficial for complex or analytic field definitions.

Highlights

• Direct Field Evaluation in Kernels: Electric and magnetic fields are now evaluated directly within the boris_push_kernel! and velocity_back_kernel! functions using the provided field functions (Efunc, Bfunc), eliminating the need for pre-calculated field arrays.
• Removal of Temporary Field Arrays: The pull request removes the Ex_arr, Ey_arr, Ez_arr, Bx_arr, By_arr, Bz_arr arrays and associated field evaluation functions (evaluate_interp_fields_kernel!, evaluate_fields_on_particles!, gpu_field_evaluation!). This significantly reduces memory overhead.
• Optimized CPU Buffer Usage: A dedicated xv_cpu_buffer is now consistently used for saving particle states to the CPU, avoiding redundant array allocations and improving memory management.
• FieldInterpolator Wrapper: A FieldInterpolator wrapper is introduced during GPU adaptation of interpolation objects to ensure consistent calling conventions for field functions.

🧠 New Feature in Public Preview: You can now enable Memory to help Gemini Code Assist learn from your team's feedback. This makes future code reviews more consistent and personalized to your project's style. Click here to enable Memory in your admin console.

Changelog

• src/boris_kernel.jl
  • Removed evaluate_interp_fields_kernel!, evaluate_fields_on_particles!, and gpu_field_evaluation! functions.
  • Modified boris_push_kernel! and velocity_back_kernel! to accept field functions (Efunc, Bfunc) and time (t) directly, instead of pre-computed field arrays.
  • Updated adapt_field_to_gpu to wrap adapted interpolation functions in FieldInterpolator to maintain calling conventions.
  • Eliminated the creation and usage of temporary Ex_arr, Ey_arr, Ez_arr, Bx_arr, By_arr, Bz_arr arrays.
  • Refactored CPU buffer management for saving particle states, reusing a single xv_cpu_buffer.

Activity

• No specific activity (comments, reviews, or progress updates) has been recorded for this pull request yet.

Using Gemini Code Assist

The full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips.

Invoking Gemini

You can request assistance from Gemini at any point by creating a comment using either /gemini <command> or @gemini-code-assist <command>. Below is a summary of the supported commands on the current page.

Feature	Command	Description
Code Review	/gemini review	Performs a code review for the current pull request in its current state.
Pull Request Summary	/gemini summary	Provides a summary of the current pull request in its current state.
Comment	@gemini-code-assist	Responds in comments when explicitly tagged, both in pull request comments and review comments.
Help	/gemini help	Displays a list of available commands.

Customization

To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a .gemini/ folder in the base of the repository. Detailed instructions can be found here.

Limitations & Feedback

Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here.

You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension.

Footnotes

1. Review the Privacy Notices, Generative AI Prohibited Use Policy, Terms of Service, and learn how to configure Gemini Code Assist in GitHub here. Gemini can make mistakes, so double check it and use code with caution. ↩

[gemini-code-assist]

Code Review

This pull request is a great refactoring of the Boris kernel. It significantly simplifies the code by removing temporary field arrays and evaluating fields directly inside the GPU kernels. Reusing CPU buffers for saving particle states is also a good performance optimization. Overall, these changes should improve performance and reduce memory usage. I've found one potential issue regarding time-dependent interpolated fields which I've commented on.

[codecov]

Codecov Report

❌ Patch coverage is 77.50000% with 54 lines in your changes missing coverage. Please review.
✅ Project coverage is 80.93%. Comparing base (9dc631e) to head (e01136a).
⚠️ Report is 2 commits behind head on master.

Files with missing lines	Patch %	Lines
src/boris_kernel.jl	75.00%	41 Missing ⚠️
src/utility/interpolation.jl	60.00%	12 Missing ⚠️
src/prepare.jl	80.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #452      +/-   ##
==========================================
- Coverage   81.40%   80.93%   -0.48%     
==========================================
  Files          20       20              
  Lines        2044     2025      -19     
==========================================
- Hits         1664     1639      -25     
- Misses        380      386       +6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[github-actions]

Benchmark Results (Julia v1)

Time benchmarks

	master	e01136a...	master / e01136a...
interpolation/cartesian	0.12 ± 0.001 μs	0.12 ± 0.01 μs	1 ± 0.084
interpolation/spherical	0.411 ± 0.001 μs	0.411 ± 0.01 μs	1 ± 0.024
interpolation/time-dependent	4.22 ± 3.1 μs	4.19 ± 3.1 μs	1.01 ± 1.1
trace/GC/DiffEq Vern6	5.65 ± 1.4 μs	5.6 ± 1.4 μs	1.01 ± 0.35
trace/GC/Native RK4	0.137 ± 0.00019 s	0.135 ± 0.0013 s	1.01 ± 0.01
trace/GC/Native RK45	3.35 ± 0.08 μs	3.37 ± 0.13 μs	0.994 ± 0.045
trace/Hybrid/Sheared	3.74 ± 0.081 μs	3.42 ± 0.08 μs	1.09 ± 0.035
trace/analytic field/in place	0.0561 ± 0.03 ms	0.0585 ± 0.032 ms	0.958 ± 0.74
trace/analytic field/in place relativistic	0.0852 ± 0.032 ms	0.0803 ± 0.034 ms	1.06 ± 0.6
trace/analytic field/out of place	0.0454 ± 0.026 ms	0.0454 ± 0.032 ms	0.999 ± 0.92
trace/normalized/out of place	15.8 ± 8 μs	15.6 ± 7.8 μs	1.01 ± 0.72
trace/numerical field/Adaptive Boris	4.27 ± 0.056 ms	4.22 ± 0.082 ms	1.01 ± 0.024
trace/numerical field/Boris	7.63 ± 0.07 μs	7.35 ± 0.071 μs	1.04 ± 0.014
trace/numerical field/Boris ensemble	15.2 ± 0.14 μs	14.6 ± 0.95 μs	1.04 ± 0.068
trace/numerical field/Boris kernel	0.135 ± 0.044 ms	9.44 ± 0.12 μs	14.3 ± 4.7
trace/numerical field/Boris with fields	8.28 ± 0.08 μs	7.99 ± 0.08 μs	1.04 ± 0.014
trace/numerical field/Multistep Boris	11.4 ± 0.081 μs	11 ± 0.08 μs	1.04 ± 0.011
trace/numerical field/in place	27 ± 4.9 μs	27.2 ± 4.8 μs	0.993 ± 0.25
trace/numerical field/out of place	19.1 ± 4.7 μs	19.3 ± 4.9 μs	0.99 ± 0.35
trace/time-dependent field/in place	0.135 ± 0.0031 ms	0.135 ± 0.0034 ms	1 ± 0.034
trace/time-dependent field/out of place	0.107 ± 0.003 ms	0.106 ± 0.0032 ms	1.01 ± 0.041
time_to_load	2.16 ± 0.012 s	2.15 ± 0.0036 s	1 ± 0.0057

Memory benchmarks

	master	e01136a...	master / e01136a...
interpolation/cartesian	2 allocs: 0.234 kB	2 allocs: 0.234 kB	1
interpolation/spherical	2 allocs: 0.203 kB	2 allocs: 0.203 kB	1
interpolation/time-dependent	0.044 k allocs: 9.62 kB	0.044 k allocs: 9.62 kB	1
trace/GC/DiffEq Vern6	0.177 k allocs: 12.2 kB	0.177 k allocs: 12.2 kB	1
trace/GC/Native RK4	13 allocs: 0.382 MB	13 allocs: 0.382 MB	1
trace/GC/Native RK45	16 allocs: 2.44 kB	16 allocs: 2.44 kB	1
trace/Hybrid/Sheared	8 allocs: 2.98 kB	8 allocs: 2.98 kB	1
trace/analytic field/in place	2.07 k allocs: 0.0883 MB	2.07 k allocs: 0.0883 MB	1
trace/analytic field/in place relativistic	2.07 k allocs: 0.0883 MB	2.07 k allocs: 0.0883 MB	1
trace/analytic field/out of place	2.03 k allocs: 0.0866 MB	2.03 k allocs: 0.0866 MB	1
trace/normalized/out of place	0.753 k allocs: 0.0321 MB	0.753 k allocs: 0.0321 MB	1
trace/numerical field/Adaptive Boris	27 allocs: 1.5 MB	27 allocs: 1.5 MB	1
trace/numerical field/Boris	6 allocs: 1.16 kB	6 allocs: 1.16 kB	1
trace/numerical field/Boris ensemble	10 allocs: 2.28 kB	10 allocs: 2.28 kB	1
trace/numerical field/Boris kernel	2.03 k allocs: 0.113 MB	0.038 k allocs: 2.12 kB	54.3
trace/numerical field/Boris with fields	6 allocs: 1.52 kB	6 allocs: 1.52 kB	1
trace/numerical field/Multistep Boris	6 allocs: 1.16 kB	6 allocs: 1.16 kB	1
trace/numerical field/in place	0.432 k allocs: 20.8 kB	0.432 k allocs: 20.8 kB	1
trace/numerical field/out of place	0.393 k allocs: 18.6 kB	0.393 k allocs: 18.6 kB	1
trace/time-dependent field/in place	3.97 k allocs: 0.177 MB	3.97 k allocs: 0.177 MB	1
trace/time-dependent field/out of place	3.93 k allocs: 0.175 MB	3.93 k allocs: 0.175 MB	1
time_to_load	0.149 k allocs: 11.1 kB	0.145 k allocs: 11 kB	1.02

[henry2004y]

/gemini review

[gemini-code-assist]

Code Review

This pull request is a significant and well-executed refactoring of the Boris kernel. The main changes, which involve removing temporary field arrays and evaluating fields directly within the GPU kernels, are excellent for performance and code simplification. Reusing CPU buffers for saving states is another great optimization that reduces memory allocations. The changes to adapt_field_to_gpu and FieldInterpolator to unify the calling conventions are also well thought out. Overall, this is a high-quality contribution that improves both performance and maintainability. I have one minor suggestion for consistency.

[henry2004y]

/gemini review

[gemini-code-assist]

Code Review

This is an excellent refactoring of the Boris kernel solver. The changes significantly improve performance, especially on the CPU backend, by removing temporary field arrays, reusing buffers, and using specialized dispatches for CPU loops vs. GPU kernels. The code is now much cleaner, more modular, and easier to follow.

I've found a critical issue related to the state-saving logic on the CPU path, which leads to incorrect results being saved. The optimization to avoid copying data on the CPU is clever, but the aliased buffer becomes stale after the ping-pong buffer swap in the main loop. I've provided a detailed comment and suggestion to fix this while preserving the performance benefits. Once this is addressed, this will be a fantastic improvement to the library.

[henry2004y]

Hmm this is strange. On my local PC, the ensemble serial version is slower than the 1-thread version, based on the scripts Gemini wrote. However, it is not the case in the CI benchmark. I must have missed something.