feat: add out-of-place form relativistic equations

[Beforerr]

This pull adds out-of-place form relativistic equations.

PS: I think if γ²v² underflow, γ²v² would be zero, and the result would be zero. So, I think it is safe to remove the if statement.

PPS: Additionally, I’d like to share a simple test setup. In this setup, test resembles the current code, test2 uses the simplest form, and test3 is a midpoint between the two. Interestingly, test3 appears to run approximately 10% faster than test1.”

using LinearAlgebra
using BenchmarkTools
using StaticArrays

u = rand(6)
b = rand(3)
E = rand(3)

function test(u, b, E)
    v = @view u[4:6]
    vx, vy, vz = v
    Bx, By, Bz = b
    Ex, Ey, Ez = E

    dx, dy, dz = vx, vy, vz
    dux = vy * Bz - vz * By + Ex
    duy = vz * Bx - vx * Bz + Ey
    duz = vx * By - vy * Bx + Ez
    return SVector{6}(dx, dy, dz, dux, duy, duz)
end

function test2(u, b, E)
    v = @view u[4:6]
    dv = v × b + E
    return SVector{6}(v..., dv...)
end


function test3(u, b, E)
    v = @view u[4:6]
    v = SVector{3}(v)
    b = SVector{3}(b)
    E = SVector{3}(E)
    dv = v × b + E
    return SVector{6}(v..., dv...)
end

test(u, b, E) == test2(u, b, E) == test3(u, b, E)

@benchmark test($u, $b, $E)
@benchmark test2($u, $b, $E)
@benchmark test3($u, $b, $E

Results

julia> @benchmark test($u, $b, $E)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  3.000 ns … 33.084 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.125 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.116 ns ±  0.329 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

           ▂         █        █         ▂         ▁        ▁ ▂
  ▄▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁█ █
  3 ns         Histogram: log(frequency) by time     3.25 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark test2($u, $b, $E)
BenchmarkTools.Trial: 10000 samples with 60 evaluations.
 Range (min … max):  859.017 ns …  37.826 μs  ┊ GC (min … max): 0.00% … 96.42%
 Time  (median):     876.383 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   910.540 ns ± 715.482 ns  ┊ GC (mean ± σ):  2.63% ±  3.39%

     ▁█▆▇▆▇▁                                                     
  ▂▃▅███████▆▄▃▃▂▂▂▂▂▂▃▃▃▃▂▃▃▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  859 ns           Histogram: frequency by time          996 ns <

 Memory estimate: 688 bytes, allocs estimate: 20.

julia> @benchmark test3($u, $b, $E)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  2.708 ns … 27.500 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.792 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.837 ns ±  0.644 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

          ▂       █       ▅▄       ▃       ▂       ▁         ▁
  ▄▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█▁▁▁▁▁▁▁██▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█ █
  2.71 ns      Histogram: log(frequency) by time        3 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 84.43%. Comparing base (f1f6197) to head (5082c88).
Report is 6 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #203      +/-   ##
==========================================
+ Coverage   84.39%   84.43%   +0.04%     
==========================================
  Files           9        9              
  Lines         692      694       +2     
==========================================
+ Hits          584      586       +2     
  Misses        108      108              

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

[henry2004y]

Thanks! Could you please add a new test for trace_relativistic_normalized?

[Beforerr]

Sure, would add the test later this week.

[Beforerr]

Also slightly changetrace_normalized!. In my application, this improve the speed by about 25%.

[henry2004y]

Also slightly changetrace_normalized!. In my application, this improve the speed by about 25%.

My guess is that in your application, functions E and B were not returning static vectors? If wrapping the return values from these function as SVectors do not trigger any regression, I am considering modifying other similar equations as well.

[henry2004y]

I will merge this first. The benchmark failure seems not related to the PR itself. The rest of the discussed optimization will be added in another PR.

[Beforerr]

Also slightly changetrace_normalized!. In my application, this improve the speed by about 25%.

My guess is that in your application, functions E and B were not returning static vectors? If wrapping the return values from these function as SVectors do not trigger any regression, I am considering modifying other similar equations as well.

Actually not, E and B is returning static vectors. I think this speedup mainly comes from ensuring v, B, E are SVector{3}, and thus utiliziing a faster × for SVector. From my experience the regression is negligible.

[Beforerr]

Also for all cases I tries, out-of-place with SVector is about 20% faster than the correspondong in-place form even with save_everystep = false. And significantly faster if we save everystep.