Relativistic momentum #200
Merged
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## master #200 +/- ##
==========================================
+ Coverage 83.48% 84.39% +0.90%
==========================================
Files 9 9
Lines 660 692 +32
==========================================
+ Hits 551 584 +33
+ Misses 109 108 -1 ☔ View full report in Codecov by Sentry. |
Benchmark resultJudge resultBenchmark Report for /home/runner/work/TestParticle.jl/TestParticle.jlJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/TestParticle.jl/TestParticle.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/TestParticle.jl/TestParticle.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
|
|
I would also prefer solving momentum (which is also more common in PIC codes). PS: If you are concerned with concerned about velocity accuracy, PPS: this may be faster. |
|
There should be no
|
|
I feel like in the normalized relativistic case, the only reasonable velocity normalization factor is c; otherwise, in the Lorentz factor |
|
BTW why use v̂ and I feel like v̂ is an intermediate variable, creating it any way would slow the code if not optimized by the compiler. |
Modern compilers are very good at branch prediction. In our case, small γ²v² value is a rare case, so this branch will almost never be executed. Even if you write vx, vy, vz = vmag * normalize(γv)without v̂, it will still allocate because of Let's compare them side-by-side: function trace_relativistic_normalized!(dy, y, p::TestParticle.TPNormalizedTuple, t)
_, E, B = p
Ex, Ey, Ez = E(y, t)
Bx, By, Bz = B(y, t)
γv = @view y[4:6]
γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2
if γ²v² > eps(eltype(dy))
v̂ = SVector{3, eltype(dy)}(normalize(γv))
else # no velocity
v̂ = SVector{3, eltype(dy)}(0, 0, 0)
end
vmag = √(γ²v² / (1 + γ²v²))
vx, vy, vz = vmag * v̂[1], vmag * v̂[2], vmag * v̂[3]
dy[1], dy[2], dy[3] = vx, vy, vz
dy[4] = vy*Bz - vz*By + Ex
dy[5] = vz*Bx - vx*Bz + Ey
dy[6] = vx*By - vy*Bx + Ez
return
end
function trace_relativistic_normalized2!(dy, y, p::TestParticle.TPNormalizedTuple, t)
_, E, B = p
Ex, Ey, Ez = E(y, t)
Bx, By, Bz = B(y, t)
γv = @view y[4:6]
γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2
T = eltype(dy)
if γ²v² > eps(T)
vmag = √(γ²v² / (1 + γ²v²))
vx, vy, vz = vmag * normalize(γv)
else # no velocity
vx, vy, vz = zero(T), zero(T), zero(T)
end
dy[1], dy[2], dy[3] = vx, vy, vz
dy[4] = vy*Bz - vz*By + Ex
dy[5] = vz*Bx - vx*Bz + Ey
dy[6] = vx*By - vy*Bx + Ez
return
end
using BenchmarkTools
# Tracing relativistic particle in dimensionless units
param = prepare(xu -> SA[0.0, 0.0, 0.0], xu -> SA[0.0, 0.0, 1.0]; species=User)
tspan = (0.0, 1.0) # 1/2π period
stateinit = [0.0, 0.0, 0.0, 0.5, 0.0, 0.0]
prob = ODEProblem(trace_relativistic_normalized!, stateinit, tspan, param)julia> @benchmark TestParticle.trace_relativistic_normalized!(stateinit, stateinit, prob.p, 0.0)
BenchmarkTools.Trial: 10000 samples with 927 evaluations.
Range (min … max): 109.385 ns … 30.531 μs ┊ GC (min … max): 0.00% … 99.22%
Time (median): 132.039 ns ┊ GC (median): 0.00%
Time (mean ± σ): 157.137 ns ± 382.613 ns ┊ GC (mean ± σ): 7.44% ± 3.89%
▃█▆▃
▃▄████▇▆▅▅▄▄▃▃▄▅▄▅▅▄▄▃▃▂▂▂▂▂▂▂▂▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
109 ns Histogram: frequency by time 289 ns <
Memory estimate: 96 bytes, allocs estimate: 3.
julia> @benchmark TestParticle.trace_relativistic_normalized2!(stateinit, stateinit, prob.p, 0.0)
BenchmarkTools.Trial: 10000 samples with 875 evaluations.
Range (min … max): 129.257 ns … 32.808 μs ┊ GC (min … max): 0.00% … 99.35%
Time (median): 143.086 ns ┊ GC (median): 0.00%
Time (mean ± σ): 177.729 ns ± 448.856 ns ┊ GC (mean ± σ): 10.46% ± 5.04%
▂▅█▁
▅████▆▄▃▂▂▂▂▂▂▂▂▂▁▁▂▂▂▂▂▂▁▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
129 ns Histogram: frequency by time 344 ns <
Memory estimate: 176 bytes, allocs estimate: 5.Regarding 1e-20, I agree it's an arbitrarily chosen bad value. I will replace it with |
|
Your benchmark is surprising. I did a simple test. Results |
I followed your advice in the new commit. This can now remove all the allocations by using SVector before passing to julia> @benchmark TestParticle.trace_relativistic_normalized!($stateinit, $stateinit, $prob.p, 0.0)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
Range (min … max): 17.900 ns … 222.200 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 18.000 ns ┊ GC (median): 0.00%
Time (mean ± σ): 18.883 ns ± 4.304 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▂▂▁ ▁
█▇▆▆▇▇▇▇█████████▇▇▆▇▆▆▇▇▅▄▄▅▆▄▅▅▅▄▄▅▄▄▄▄▄▄▄▄▅▅▅▄▄▁▅▃▄▄▃▁▃▁▄ █
17.9 ns Histogram: log(frequency) by time 36.6 ns <
Memory estimate: 0 bytes, allocs estimate: 0.P.S. Previously we did not exclude the allocation from passing the input arrays. |
Benchmark resultJudge resultBenchmark Report for /home/runner/work/TestParticle.jl/TestParticle.jlJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/TestParticle.jl/TestParticle.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/TestParticle.jl/TestParticle.jlJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
|
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Handle #198 by switching to solve "momentum"$\vec{p}/m = \gamma \vec{v}$ instead of $\vec{v}$ in the relativistic case. In this way we do not need to check whether the computed velocity is larger than the speed of light, since the derivation guarantees smaller than c speed. A small caveat is that when velocity is 0, the direction is undetermined, which requires an additional branch.
As quoted from the original discussion note, this form is a bit unintuitive and requires a conversion from relativistic momentum to velocity in certain cases. However, I feel like it is more natural with relativity. Further discussions are welcomed! @Beforerr
The normalization case may need further testing.