From 393ca596d4e818c1e37ba1b6745333906ad5ce47 Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Sun, 27 Oct 2024 17:58:26 -0400 Subject: [PATCH 01/13] solve momentum in the relativistic case #198 --- src/equations.jl | 96 +++++++++++++++++++++--------------------------- test/runtests.jl | 37 ++++++++----------- 2 files changed, 57 insertions(+), 76 deletions(-) diff --git a/src/equations.jl b/src/equations.jl index 2065369a..758a2433 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -80,43 +80,30 @@ function trace(y, p::FullTPTuple, t) SVector{6}(dx, dy, dz, dux, duy, duz) end -const FTLError = """ -Particle faster than the speed of light! - -If the initial velocity is slower than light and -adaptive timestepping of the solver is turned on, it -is better to set a small initial stepsize (dt) or -maximum dt for adaptive timestepping (dtmax). - -More details about the keywords of initial stepsize -can be found in this documentation page: -https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Stepsize-Control -""" - """ trace_relativistic!(dy, y, p::TPTuple, t) -ODE equations for relativistic charged particle moving in static EM field with in-place -form. +ODE equations for relativistic charged particle (x, γv) moving in static EM field with in-place form. """ function trace_relativistic!(dy, y, p::TPTuple, t) q2m, E, B = p - - u2 = y[4]^2 + y[5]^2 + y[6]^2 - c2 = c^2 - if u2 ≥ c2 - throw(DomainError(u2, FTLError)) - end - - γInv3 = √(1.0 - u2/c2)^3 - vx, vy, vz = @view y[4:6] 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² < 1e-20 # no velocity + v̂ = SVector{3, Float64}(0, 0, 0) + else + v̂ = SVector{3, Float64}(normalize(γv)) + end + vmag = √(γ²v² / (1 + γ²v²/c^2)) + vx, vy, vz = vmag * v̂[1], vmag * v̂[2], vmag * v̂[3] + dy[1], dy[2], dy[3] = vx, vy, vz - dy[4] = q2m*γInv3*(vy*Bz - vz*By + Ex) - dy[5] = q2m*γInv3*(vz*Bx - vx*Bz + Ey) - dy[6] = q2m*γInv3*(vx*By - vy*Bx + Ez) + dy[4] = q2m * (vy*Bz - vz*By + Ex) + dy[5] = q2m * (vz*Bx - vx*Bz + Ey) + dy[6] = q2m * (vx*By - vy*Bx + Ez) return end @@ -124,27 +111,27 @@ end """ trace_relativistic(y, p::TPTuple, t) -> SVector{6, Float64} -ODE equations for relativistic charged particle moving in static EM field with out-of-place -form. +ODE equations for relativistic charged particle (x, γv) moving in static EM field with out-of-place form. """ function trace_relativistic(y, p::TPTuple, t) q2m, E, B = p - - u2 = y[4]^2 + y[5]^2 + y[6]^2 - c2 = c^2 - if u2 ≥ c2 - throw(DomainError(u2, FTLError)) - end - - γInv3 = √(1.0 - u2/c2)^3 - vx, vy, vz = @view y[4:6] 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² < 1e-20 # no velocity + v̂ = SVector{3, Float64}(0, 0, 0) + else + v̂ = SVector{3, Float64}(normalize(γv)) + end + vmag = √(γ²v² / (1 + γ²v²/c^2)) + vx, vy, vz = vmag * v̂[1], vmag * v̂[2], vmag * v̂[3] + dx, dy, dz = vx, vy, vz - dux = q2m*γInv3*(vy*Bz - vz*By + Ex) - duy = q2m*γInv3*(vz*Bx - vx*Bz + Ey) - duz = q2m*γInv3*(vx*By - vy*Bx + Ez) + dux = q2m * (vy*Bz - vz*By + Ex) + duy = q2m * (vz*Bx - vx*Bz + Ey) + duz = q2m * (vx*By - vy*Bx + Ez) SVector{6}(dx, dy, dz, dux, duy, duz) end @@ -175,26 +162,27 @@ end """ trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) -Normalized ODE equations for relativistic charged particle moving in static EM field with -in-place form. +Normalized ODE equations for relativistic charged particle (x, γv) moving in static EM field with in-place form. """ function trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) Ω, E, B = p - - u2 = y[4]^2 + y[5]^2 + y[6]^2 - if u2 ≥ 1 - throw(DomainError(u2, FTLError)) - end - - γInv3 = √(1.0 - u2)^3 - vx, vy, vz = @view y[4:6] 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² < 1e-20 # no velocity + v̂ = SVector{3, Float64}(0, 0, 0) + else + v̂ = SVector{3, Float64}(normalize(γv)) + end + vmag = √(γ²v² / (1 + γ²v²/c^2)) + vx, vy, vz = vmag * v̂[1], vmag * v̂[2], vmag * v̂[3] + dy[1], dy[2], dy[3] = vx, vy, vz - dy[4] = Ω*γInv3*(vy*Bz - vz*By + Ex) - dy[5] = Ω*γInv3*(vz*Bx - vx*Bz + Ey) - dy[6] = Ω*γInv3*(vx*By - vy*Bx + Ez) + dy[4] = Ω * (vy*Bz - vz*By + Ex) + dy[5] = Ω * (vz*Bx - vx*Bz + Ey) + dy[6] = Ω * (vx*By - vy*Bx + Ez) return end diff --git a/test/runtests.jl b/test/runtests.jl index 9b9560fb..bad2b51d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -275,18 +275,6 @@ end F_field(r, v, t) = SA[r, v, t] @test typeof(Field(F_field)).parameters[1] == false - - x0 = [10.0, 10.0, 0.0] # initial position, [m] - u0 = [1e10, 0.0, 0.0] # initial velocity, [m/s] - stateinit = [x0..., u0...] - tspan = (0.0, 2e-7) - - param = prepare(E_field, E_field; species=Electron) - prob = ODEProblem(trace_relativistic!, stateinit, tspan, param) - @test_throws DomainError solve(prob, Tsit5()) - - prob = ODEProblem(trace_relativistic, SA[stateinit...], tspan, param) - @test_throws DomainError solve(prob, Tsit5()) end @testset "relativistic particle" begin @@ -299,37 +287,42 @@ end # calculate the energy [eV] of a electron function calc_energy(sol) - v = hypot(sol.u[end][4:6]...) - γ = 1/sqrt(1-(v/c)^2) + γv = @view sol.u[end][4:6] + γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 + v² = γ²v² / (1 + γ²v²/c^2) + γ = 1 / √(1 - v²/c^2) + return -(γ-1)*mₑ*c^2/qₑ end x0 = [10.0, 10.0, 0.0] # initial position, [m] - u0 = [0.0, 0.0, 0.0] # initial velocity, [m/s] + u0 = [0.0, 0.0, 0.0] # initial velocity γv, [m/s] tspan = (0.0, 2e-7) stateinit = [x0..., u0...] param = prepare(E_field, B_field; species=Electron) prob = ODEProblem(trace_relativistic!, stateinit, tspan, param) - sol = solve(prob, Vern6(); dtmax=1e-10) + sol = solve(prob, Vern6(); dt=1e-10, adaptive=false) - x = sol.u[end][1:3] # Test whether the kinetic energy [eV] of the electron # is equal to the electric potential energy gained. - @test calc_energy(sol)/(x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 + x = sol.u[end][1:3] + @test calc_energy(sol) / (x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 + # Convert to velocity + @test get_velocity(sol)[1,end] ≈ -1.6588694948554998e7 prob = ODEProblem(trace_relativistic, SA[stateinit...], tspan, param) - sol = solve(prob, Vern6(); dtmax=1e-10) + sol = solve(prob, Vern6(); dt=1e-10, adaptive=false) + x = sol.u[end][1:3] - - @test calc_energy(sol)/(x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 + @test calc_energy(sol) / (x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 # 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) sol = solve(prob, Vern6()) - @test sol.u[end][1] ≈ 0.46557792820784516 + @test sol.u[end][1] ≈ 0.42073549161010637 end @testset "normalized fields" begin From e900cd80e655f0eca51ef8b4f9609c05262576e6 Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Sun, 27 Oct 2024 17:58:45 -0400 Subject: [PATCH 02/13] Add get_velocity --- src/TestParticle.jl | 2 +- src/utility/utility.jl | 16 ++++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/TestParticle.jl b/src/TestParticle.jl index 0fbef57f..45599983 100644 --- a/src/TestParticle.jl +++ b/src/TestParticle.jl @@ -19,7 +19,7 @@ export trace!, trace_relativistic!, trace_normalized!, trace, trace_relativistic trace_relativistic_normalized!, trace_gc!, trace_gc_1st!, trace_gc_drifts! export Proton, Electron, Ion, User export Maxwellian, BiMaxwellian -export get_gyrofrequency, get_gyroperiod, get_gyroradius +export get_gyrofrequency, get_gyroperiod, get_gyroradius, get_velocity export orbit, monitor export TraceProblem diff --git a/src/utility/utility.jl b/src/utility/utility.jl index 5614e769..1c590aad 100644 --- a/src/utility/utility.jl +++ b/src/utility/utility.jl @@ -104,4 +104,20 @@ end function get_gyroperiod(B::AbstractFloat=5e-9; q::AbstractFloat=qᵢ, m::AbstractFloat=mᵢ) ω = get_gyrofrequency(B; q, m) 2π / ω +end + +"Return velocity from relativistic γv in `sol`." +function get_velocity(sol) + v = Array{eltype(sol.u[1]), 2}(undef, 3, length(sol)) + for is in axes(v, 2) + γv = @view sol[4:6, is] + γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 + v² = γ²v² / (1 + γ²v²/c^2) + γ = 1 / √(1 - v²/c^2) + for i in axes(v, 1) + v[i,is] = γv[i] / γ + end + end + + v end \ No newline at end of file From e2b3f61b748eb1621c2ed3d7e650378b3d2494b0 Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Sun, 27 Oct 2024 18:58:53 -0400 Subject: [PATCH 03/13] Add relativistic benchmark --- benchmark/benchmarks.jl | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 83961ec0..7e984f62 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -43,25 +43,36 @@ stateinit = [x0..., u0...] param_analytic = prepare(E_analytic, B_analytic) prob_ip = ODEProblem(trace!, stateinit, tspan, param_analytic) # in place +prob_rel_ip = ODEProblem(trace_relativistic!, stateinit, tspan, param_analytic) # in place prob_oop = ODEProblem(trace, SA[stateinit...], tspan, param_analytic) # out of place -SUITE["trace"]["analytic field"]["in place"] = @benchmarkable solve($prob_ip, Tsit5(); save_idxs=[1,2,3]) -SUITE["trace"]["analytic field"]["out of place"] = @benchmarkable solve($prob_oop, Tsit5(); save_idxs=[1,2,3]) +SUITE["trace"]["analytic field"]["in place"] = + @benchmarkable solve($prob_ip, Tsit5(); save_idxs=[1,2,3]) +SUITE["trace"]["analytic field"]["out of place"] = + @benchmarkable solve($prob_oop, Tsit5(); save_idxs=[1,2,3]) +SUITE["trace"]["analytic field"]["in place relativistic"] = + @benchmarkable solve($prob_rel_ip, Tsit5(); save_idxs=[1,2,3]) param_numeric = prepare(mesh, E_numeric, B_numeric) prob_ip = ODEProblem(trace!, stateinit, tspan, param_numeric) # in place prob_oop = ODEProblem(trace, SA[stateinit...], tspan, param_numeric) # out of place prob_boris = TraceProblem(stateinit, tspan, param_numeric) -SUITE["trace"]["numerical field"]["in place"] = @benchmarkable solve($prob_ip, Tsit5(); save_idxs=[1,2,3]) -SUITE["trace"]["numerical field"]["out of place"] = @benchmarkable solve($prob_oop, Tsit5(); save_idxs=[1,2,3]) -SUITE["trace"]["numerical field"]["Boris"] = @benchmarkable TestParticle.solve($prob_boris; dt=1/7, savestepinterval=10) -SUITE["trace"]["numerical field"]["Boris ensemble"] = @benchmarkable TestParticle.solve($prob_boris; dt=1/7, savestepinterval=10, trajectories=2) +SUITE["trace"]["numerical field"]["in place"] = + @benchmarkable solve($prob_ip, Tsit5(); save_idxs=[1,2,3]) +SUITE["trace"]["numerical field"]["out of place"] = + @benchmarkable solve($prob_oop, Tsit5(); save_idxs=[1,2,3]) +SUITE["trace"]["numerical field"]["Boris"] = + @benchmarkable TestParticle.solve($prob_boris; dt=1/7, savestepinterval=10) +SUITE["trace"]["numerical field"]["Boris ensemble"] = + @benchmarkable TestParticle.solve($prob_boris; dt=1/7, savestepinterval=10, trajectories=2) param_td = prepare(E_td, B_td, F_td) prob_ip = ODEProblem(trace!, stateinit, tspan, param_td) # in place prob_oop = ODEProblem(trace, SA[stateinit...], tspan, param_td) # out of place -SUITE["trace"]["time-dependent field"]["in place"] = @benchmarkable solve($prob_ip, Tsit5(); save_idxs=[1,2,3]) -SUITE["trace"]["time-dependent field"]["out of place"] = @benchmarkable solve($prob_oop, Tsit5(); save_idxs=[1,2,3]) +SUITE["trace"]["time-dependent field"]["in place"] = + @benchmarkable solve($prob_ip, Tsit5(); save_idxs=[1,2,3]) +SUITE["trace"]["time-dependent field"]["out of place"] = + @benchmarkable solve($prob_oop, Tsit5(); save_idxs=[1,2,3]) stateinit_gc, param_gc = TestParticle.prepare_gc(stateinit, E_analytic, B_analytic, species=Proton, removeExB=true) From 6641b6c8ca993098ff4df930afc35cf5427fcddf Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Sun, 27 Oct 2024 20:36:45 -0400 Subject: [PATCH 04/13] Update sol extraction; add utility methods --- docs/examples/advanced/demo_proton_dipole.jl | 6 +- src/TestParticle.jl | 3 +- src/utility/utility.jl | 28 +++++- test/runtests.jl | 94 +++++--------------- 4 files changed, 56 insertions(+), 75 deletions(-) diff --git a/docs/examples/advanced/demo_proton_dipole.jl b/docs/examples/advanced/demo_proton_dipole.jl index fbb05684..6ca28e72 100644 --- a/docs/examples/advanced/demo_proton_dipole.jl +++ b/docs/examples/advanced/demo_proton_dipole.jl @@ -62,9 +62,9 @@ f = DisplayAs.PNG(f) #hide # In the above we used Verner's “Most Efficient” 9/8 Runge-Kutta method. Let's check other algorithms. function get_energy_ratio(sol) - vx = getindex.(sol.u, 4) - vy = getindex.(sol.u, 5) - vz = getindex.(sol.u, 6) + vx = @view sol[4,:] + vy = @view sol[5,:] + vz = @view sol[6,:] Einit = vx[1]^2 + vy[1]^2 + vz[1]^2 Eend = vx[end]^2 + vy[end]^2 + vz[end]^2 diff --git a/src/TestParticle.jl b/src/TestParticle.jl index 45599983..d5015d28 100644 --- a/src/TestParticle.jl +++ b/src/TestParticle.jl @@ -19,7 +19,8 @@ export trace!, trace_relativistic!, trace_normalized!, trace, trace_relativistic trace_relativistic_normalized!, trace_gc!, trace_gc_1st!, trace_gc_drifts! export Proton, Electron, Ion, User export Maxwellian, BiMaxwellian -export get_gyrofrequency, get_gyroperiod, get_gyroradius, get_velocity +export get_gyrofrequency, get_gyroperiod, get_gyroradius, get_velocity, get_energy, + energy2velocity export orbit, monitor export TraceProblem diff --git a/src/utility/utility.jl b/src/utility/utility.jl index 1c590aad..7e0c7113 100644 --- a/src/utility/utility.jl +++ b/src/utility/utility.jl @@ -120,4 +120,30 @@ function get_velocity(sol) end v -end \ No newline at end of file +end + +"Return the energy [eV] from relativistic `sol`." +function get_energy(sol::AbstractODESolution; m=mᵢ, q=qᵢ) + e = Vector{eltype(sol.u[1])}(undef, length(sol)) + for i in eachindex(e) + γv = @view sol[4:6, i] + γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 + v² = γ²v² / (1 + γ²v²/c^2) + γ = 1 / √(1 - v²/c^2) + e[i] = (γ-1)*m*c^2/abs(q) + end + + e +end + +"Calculate the energy [eV] of a relativistic particle from γv." +function get_energy(γv; m=mᵢ, q=qᵢ) + γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 + v² = γ²v² / (1 + γ²v²/c^2) + γ = 1 / √(1 - v²/c^2) + + (γ-1)*m*c^2/abs(q) +end + +"Return velocity magnitude from energy in [eV]." +energy2velocity(Ek; m=mᵢ, q=qᵢ) = c*sqrt(1 - 1/(1+Ek*abs(q)/(m*c^2))^2) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index bad2b51d..59b6595b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -92,17 +92,17 @@ end prob = ODEProblem(trace!, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1], isoutofdomain, verbose=false) # There are numerical differences on x86 and ARM platforms! - @test getindex.(sol.u, 1)[end] ≈ 0.7388945226814018 + @test sol[1,end] ≈ 0.7388945226814018 # Because the field is uniform, the order of interpolation does not matter. param = prepare(x, y, z, E, B; order=2) prob = remake(prob; p=param) sol = solve(prob, Tsit5(); save_idxs=[1], isoutofdomain, verbose=false) - @test getindex.(sol.u, 1)[end] ≈ 0.7388945226814018 + @test sol[1,end] ≈ 0.7388945226814018 param = prepare(x, y, z, E, B; order=3) prob = remake(prob; p=param) sol = solve(prob, Tsit5(); save_idxs=[1], isoutofdomain, verbose=false) - @test getindex.(sol.u, 1)[end] ≈ 0.7388945226814018 + @test sol[1,end] ≈ 0.7388945226814018 # GC prepare stateinit_gc, param_gc = prepare_gc(stateinit, x, y, z, E, B; @@ -115,19 +115,14 @@ end param = prepare(grid, E, B) prob = ODEProblem(trace!, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1]) - - x = getindex.(sol.u, 1) - - @test length(x) == 8 && x[end] ≈ 0.8540967226885379 + @test length(sol) == 8 && sol[1, end] ≈ 0.8540967226885379 trajectories = 10 prob = ODEProblem(trace!, stateinit, tspan, param) ensemble_prob = EnsembleProblem(prob, prob_func=prob_func) sol = solve(ensemble_prob, Tsit5(), EnsembleThreads(); trajectories=trajectories, save_idxs=[1]) - x = getindex.(sol.u[10].u, 1) - @test x[7] ≈ 0.08230289216655486 rtol=1e-6 stateinit = SA[x0..., u0...] @@ -136,10 +131,7 @@ end param = prepare(grid, E, B) prob = ODEProblem(trace, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1]) - - x = getindex.(sol.u, 1) - - @test length(x) == 8 && x[end] ≈ 0.8540967226885379 + @test length(sol) == 8 && sol[1,end] ≈ 0.8540967226885379 end @testset "analytical field" begin @@ -152,7 +144,7 @@ end Rₑ = TestParticle.Rₑ # initial velocity, [m/s] - v₀ = TestParticle.sph2cart(c*sqrt(1-1/(1+Ek*q/(m*c^2))^2), 0.0, π/4) + v₀ = TestParticle.sph2cart(energy2velocity(Ek; q, m), 0.0, π/4) # initial position, [m] r₀ = TestParticle.sph2cart(2.5*Rₑ, 0.0, π/2) stateinit = [r₀..., v₀...] @@ -162,21 +154,16 @@ end prob = ODEProblem(trace!, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1]) - x = getindex.(sol.u, 1) - @test guiding_center([stateinit..., 0.0], param)[1] == 1.59275e7 @test get_gc(param) isa Function - @test x[300] ≈ 1.2563192407332942e7 rtol=1e-6 + @test sol[1, 300] ≈ 1.2563192407332942e7 rtol=1e-6 # static array version (results not identical with above: maybe some bugs?) stateinit = SA[r₀..., v₀...] prob = ODEProblem(trace, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1]) - - x = getindex.(sol.u, 1) - - @test x[306] ≈ 1.2440619301099773e7 rtol=1e-5 + @test sol[1,306] ≈ 1.2440619301099773e7 rtol=1e-5 end @testset "mixed type fields" begin @@ -206,11 +193,7 @@ end param = prepare(grid, E_field, B, F; species=Electron) prob = ODEProblem(trace!, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1,2]) - - x = getindex.(sol.u, 1) - y = getindex.(sol.u, 2) - - @test x[end] ≈ 1.5324506965560782 && y[end] ≈ -2.8156470047903706 + @test sol[1,end] ≈ 1.5324506965560782 && sol[2,end] ≈ -2.8156470047903706 end @testset "time-independent fields" begin @@ -239,41 +222,30 @@ end param = prepare(grid, E_field, B_field, F; species=Electron) prob = ODEProblem(trace!, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1,2,3]) - - x = getindex.(sol.u, 1) - y = getindex.(sol.u, 2) - z = getindex.(sol.u, 3) - - @test x[end] ≈ -1.2828663442681638 && y[end] ≈ 1.5780464321537067 && z[end] ≈ 1.0 + @test sol[1,end] ≈ -1.2828663442681638 && sol[2,end] ≈ 1.5780464321537067 && + sol[3,end] ≈ 1.0 @test guiding_center([stateinit..., 0.0], param)[1] == -0.5685630064930044 F_field(r) = SA[0, 9.10938356e-42, 0] # [N] param = prepare(E_field, B_field, F_field; species=Electron) _, _, _, _, F = param - @test F(x0)[2] ≈ 9.10938356e-42 stateinit = SA[x0..., u0...] prob = ODEProblem(trace, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1,2,3]) - - x = getindex.(sol.u, 1) - y = getindex.(sol.u, 2) - z = getindex.(sol.u, 3) - - @test x[end] ≈ -1.2828663442681638 && y[end] ≈ 1.5780464321537067 && z[end] ≈ 1.0 + @test sol[1,end] ≈ -1.2828663442681638 && sol[2,end] ≈ 1.5780464321537067 && + sol[3,end] ≈ 1.0 end @testset "Exceptions" begin E_field(r, t) = SA[5e-11*sin(2π*t), 0, 0] E = Field(E_field) - @test_throws ArgumentError E([0, 0, 0]) # Test unsupported function form F_field(r, v, t) = SA[r, v, t] - @test typeof(Field(F_field)).parameters[1] == false end @@ -285,16 +257,6 @@ end return SA[Ex, Ey, 0.0] end - # calculate the energy [eV] of a electron - function calc_energy(sol) - γv = @view sol.u[end][4:6] - γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 - v² = γ²v² / (1 + γ²v²/c^2) - γ = 1 / √(1 - v²/c^2) - - return -(γ-1)*mₑ*c^2/qₑ - end - x0 = [10.0, 10.0, 0.0] # initial position, [m] u0 = [0.0, 0.0, 0.0] # initial velocity γv, [m/s] tspan = (0.0, 2e-7) @@ -307,15 +269,16 @@ end # Test whether the kinetic energy [eV] of the electron # is equal to the electric potential energy gained. x = sol.u[end][1:3] - @test calc_energy(sol) / (x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 + @test get_energy(sol[4:6,end]; m=mₑ, q=qₑ) / (x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 + @test get_energy(sol; m=mₑ, q=qₑ)[end] / (x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 # Convert to velocity @test get_velocity(sol)[1,end] ≈ -1.6588694948554998e7 prob = ODEProblem(trace_relativistic, SA[stateinit...], tspan, param) sol = solve(prob, Vern6(); dt=1e-10, adaptive=false) - x = sol.u[end][1:3] - @test calc_energy(sol) / (x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 + x = sol[1:3, end] + @test get_energy(sol[4:6,end]; m=mₑ, q=qₑ) / (x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 # 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 @@ -353,10 +316,7 @@ end prob = ODEProblem(trace_normalized!, stateinit, tspan, param) sol = solve(prob, Vern9(); save_idxs=[1]) - - x = getindex.(sol.u, 1) - - @test length(x) == 7 && x[end] ≈ 1.0 + @test length(sol) == 7 && sol[1,end] ≈ 1.0 # 2D x = range(-10, 10, length=15) @@ -383,23 +343,19 @@ end param = prepare(x, y, E, B; species=Proton, bc=2) prob = ODEProblem(trace_normalized!, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1]) - - xs = getindex.(sol.u, 1) - @test length(xs) == 9 && xs[end] ≈ 0.9999998697180689 + @test length(sol) == 9 && sol[1,end] ≈ 0.9999998697180689 # Because the field is uniform, the order of interpolation does not matter. param = prepare(grid, E, B; order=2) prob = remake(prob; p=param) sol = solve(prob, Tsit5(); save_idxs=[1]) - xs = getindex.(sol.u, 1) - @test length(xs) == 9 && xs[end] ≈ 0.9999998697180689 + @test length(sol) == 9 && sol[1,end] ≈ 0.9999998697180689 # Because the field is uniform, the order of interpolation does not matter. param = prepare(grid, E, B; order=3) prob = remake(prob; p=param) sol = solve(prob, Tsit5(); save_idxs=[1]) - xs = getindex.(sol.u, 1) - @test length(xs) == 9 && xs[end] ≈ 0.9999998697180689 + @test length(sol) == 9 && sol[1,end] ≈ 0.9999998697180689 # 1D x = range(-10, 10, length=15) @@ -418,9 +374,7 @@ end param = prepare(x, E, B; species=Proton, bc=3) prob = ODEProblem(trace_normalized!, stateinit, tspan, param) sol = solve(prob, Tsit5(); save_idxs=[1]) - - xs = getindex.(sol.u, 1) - @test length(xs) == 9 && xs[end] ≈ 0.9999998697180689 + @test length(sol) == 9 && sol[1,end] ≈ 0.9999998697180689 end @testset "Boris pusher" begin @@ -483,8 +437,8 @@ end gc_x0 = gc(stateinit) prob_gc_analytic = ODEProblem(trace_gc_drifts!, gc_x0, tspan, (param..., sol)) sol_gc_analytic = solve(prob_gc_analytic, Vern9(); save_idxs=[1,2,3]) - @test sol_gc.u[end][1] ≈ 0.9896155284173717 - @test sol_gc_analytic.u[end][1] ≈ 0.9906923500002904 rtol=1e-5 + @test sol_gc[1,end] ≈ 0.9896155284173717 + @test sol_gc_analytic[1,end] ≈ 0.9906923500002904 rtol=1e-5 stateinit_gc, param_gc = TestParticle.prepare_gc(stateinit, uniform_E, curved_B, species=Proton, removeExB=false) From fd18a87d827f6b1a803ec0cf1dac6184c7afd8a2 Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Sun, 27 Oct 2024 20:41:32 -0400 Subject: [PATCH 05/13] Add more coverage --- test/runtests.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 59b6595b..edf1c366 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -286,6 +286,10 @@ end prob = ODEProblem(trace_relativistic_normalized!, stateinit, tspan, param) sol = solve(prob, Vern6()) @test sol.u[end][1] ≈ 0.42073549161010637 + stateinit = zeros(6) + prob = ODEProblem(trace_relativistic_normalized!, stateinit, tspan, param) + sol = solve(prob, Vern6()) + @test sol[1,end] == 0.0 end @testset "normalized fields" begin From 6c5dedca2bb57cef86ab23d13c858cb3c80e65bb Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 28 Oct 2024 10:48:54 -0400 Subject: [PATCH 06/13] Fix normalized relativistic eq --- .../basics/demo_dimensionless_dimensional.jl | 8 +++-- src/equations.jl | 34 +++++++++---------- 2 files changed, 22 insertions(+), 20 deletions(-) diff --git a/docs/examples/basics/demo_dimensionless_dimensional.jl b/docs/examples/basics/demo_dimensionless_dimensional.jl index 2f839cfb..4f59c79f 100644 --- a/docs/examples/basics/demo_dimensionless_dimensional.jl +++ b/docs/examples/basics/demo_dimensionless_dimensional.jl @@ -83,8 +83,8 @@ sol2 = solve(prob2, Vern9(); reltol=1e-4, abstol=1e-6) ### Visualization f = Figure(fontsize=18) ax = Axis(f[1, 1], - xlabel = "x", - ylabel = "y", + xlabel = "x [km]", + ylabel = "y [km]", aspect = DataAspect(), ) @@ -93,7 +93,9 @@ lines!(ax, sol1, idxs=(1, 2)) xp, yp = let trange = range(tspan2..., length=40) sol2.(trange, idxs=1) .* l₀, sol2.(trange, idxs=2) .* l₀ end -lines!(ax, xp, yp, linestyle=:dashdot, linewidth=5) +lines!(ax, xp, yp, linestyle=:dashdot, linewidth=5, color=Makie.wong_colors()[2]) +invL = inv(1e3) +scale!(ax.scene, invrL, invrL) f = DisplayAs.PNG(f) #hide diff --git a/src/equations.jl b/src/equations.jl index 758a2433..fdd22cd3 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -92,10 +92,10 @@ function trace_relativistic!(dy, y, p::TPTuple, t) γv = @view y[4:6] γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 - if γ²v² < 1e-20 # no velocity - v̂ = SVector{3, Float64}(0, 0, 0) - else - v̂ = SVector{3, Float64}(normalize(γv)) + if γ²v² > 1e-20 + v̂ = SVector{3, eltype(dy)}(normalize(γv)) + else # no velocity + v̂ = SVector{3, eltype(dy)}(0, 0, 0) end vmag = √(γ²v² / (1 + γ²v²/c^2)) vx, vy, vz = vmag * v̂[1], vmag * v̂[2], vmag * v̂[3] @@ -109,7 +109,7 @@ function trace_relativistic!(dy, y, p::TPTuple, t) end """ - trace_relativistic(y, p::TPTuple, t) -> SVector{6, Float64} + trace_relativistic(y, p::TPTuple, t) -> SVector{6} ODE equations for relativistic charged particle (x, γv) moving in static EM field with out-of-place form. """ @@ -120,10 +120,10 @@ function trace_relativistic(y, p::TPTuple, t) γv = @view y[4:6] γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 - if γ²v² < 1e-20 # no velocity - v̂ = SVector{3, Float64}(0, 0, 0) - else - v̂ = SVector{3, Float64}(normalize(γv)) + if γ²v² > 1e-20 + v̂ = SVector{3, eltype(y)}(normalize(γv)) + else # no velocity + v̂ = SVector{3, eltype(y)}(0, 0, 0) end vmag = √(γ²v² / (1 + γ²v²/c^2)) vx, vy, vz = vmag * v̂[1], vmag * v̂[2], vmag * v̂[3] @@ -165,24 +165,24 @@ end Normalized ODE equations for relativistic charged particle (x, γv) moving in static EM field with in-place form. """ function trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) - Ω, E, B = p + _, 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² < 1e-20 # no velocity - v̂ = SVector{3, Float64}(0, 0, 0) - else - v̂ = SVector{3, Float64}(normalize(γv)) + if γ²v² > 1e-20 + v̂ = SVector{3, eltype(dy)}(normalize(γv)) + else # no velocity + v̂ = SVector{3, eltype(dy)}(0, 0, 0) end vmag = √(γ²v² / (1 + γ²v²/c^2)) 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) + dy[4] = vy*Bz - vz*By + Ex + dy[5] = vz*Bx - vx*Bz + Ey + dy[6] = vx*By - vy*Bx + Ez return end From 3720118dedbfa0060565dc6b1a2c72f5c428e017 Mon Sep 17 00:00:00 2001 From: Beforerr <58776897+Beforerr@users.noreply.github.com> Date: Mon, 28 Oct 2024 07:56:55 -0700 Subject: [PATCH 07/13] correct relativistic normalization equation (#202) --- src/equations.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/equations.jl b/src/equations.jl index fdd22cd3..c271284f 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -163,6 +163,9 @@ end trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) Normalized ODE equations for relativistic charged particle (x, γv) moving in static EM field with in-place form. + +The velocity is normalized by the characteristic velocity v0 (c); the magnetic field is normalized by the characteristic magnetic field B0; + the electric field is normalized by E0=v0 B0; the position is normalized by the L = v0 / Ω0, where Ω0 = q B0 / m and the time is normalized by Ω0^-1. """ function trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) _, E, B = p @@ -176,7 +179,7 @@ function trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) else # no velocity v̂ = SVector{3, eltype(dy)}(0, 0, 0) end - vmag = √(γ²v² / (1 + γ²v²/c^2)) + 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 From 9922059d0a51c525de41df2b805539a30beda62e Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 28 Oct 2024 11:44:03 -0400 Subject: [PATCH 08/13] Add normalized relativistic demo --- docs/examples/basics/demo_dimensionless.jl | 58 ++++++++++++++++++++-- src/equations.jl | 3 -- 2 files changed, 54 insertions(+), 7 deletions(-) diff --git a/docs/examples/basics/demo_dimensionless.jl b/docs/examples/basics/demo_dimensionless.jl index a5c16348..c08df3e4 100644 --- a/docs/examples/basics/demo_dimensionless.jl +++ b/docs/examples/basics/demo_dimensionless.jl @@ -1,9 +1,9 @@ # --- # title: Dimensionless Units # id: demo_dimensionless -# date: 2023-12-14 +# date: 2024-10-28 # author: "[Hongyang Zhou](https://github.com/henry2004y)" -# julia: 1.9.4 +# julia: 1.11.1 # description: Tracing charged particle in dimensionless units # --- @@ -15,7 +15,7 @@ # Tracing in dimensionless units is beneficial for many scenarios. For example, MHD simulations do not have intrinsic scales. Therefore, we can do dimensionless particle tracing in MHD fields, and then convert to any scale we would like. -# Now let's demonstrate this with `trace_normalized!`. +# Now let's demonstrate this with `trace_normalized!` and `trace_relativistic_normalized!`. import DisplayAs #hide using TestParticle @@ -67,6 +67,56 @@ ax = Axis(f[1, 1], aspect = DataAspect() ) -lines!(ax, sol, vars=(1,2)) +lines!(ax, sol, idxs=(1,2)) + +f = DisplayAs.PNG(f) #hide + +# In the relativistic case, +# * Velocity is normalized by speed of light c, $V = V^\prime c$; +# * Magnetic field is normalized by the characteristic magnetic field, $B = B^\prime B\_0$; +# * Electric field is normalized by $E_0 = c\,B_0$, $E = E^\prime E_0$; +# * Location is normalized by the $L = c / \Omega_0$, where $\Omega_0 = q\,B_0 / m$, and +# * Time is normalized by $\Omega_0^{-1}$, $t = t^\prime * \Omega_0^{-1}$. + +# In the small velocity scenario, it should behave similar to the non-relativistic case: + +param = prepare(xu -> SA[0.0, 0.0, 0.0], xu -> SA[0.0, 0.0, 1.0]; species=User) +tspan = (0.0, π) # half period +stateinit = [0.0, 0.0, 0.0, 0.01, 0.0, 0.0] +prob = ODEProblem(trace_relativistic_normalized!, stateinit, tspan, param) +sol = solve(prob, Vern9()) + +### Visualization +f = Figure(fontsize = 18) +ax = Axis(f[1, 1], + title = "Relativistic particle trajectory", + xlabel = "X", + ylabel = "Y", + #limits = (-0.6, 0.6, -1.1, 0.1), + aspect = DataAspect() +) + +lines!(ax, sol, idxs=(1,2)) + +f = DisplayAs.PNG(f) #hide + +# In the large velocity scenario, relativistic effect takes place: + +param = prepare(xu -> SA[0.0, 0.0, 0.0], xu -> SA[0.0, 0.0, 1.0]; species=User) +tspan = (0.0, π) # half period +stateinit = [0.0, 0.0, 0.0, 0.9, 0.0, 0.0] +prob = ODEProblem(trace_relativistic_normalized!, stateinit, tspan, param) +sol = solve(prob, Vern9()) + +### Visualization +f = Figure(fontsize = 18) +ax = Axis(f[1, 1], + title = "Relativistic particle trajectory", + xlabel = "X", + ylabel = "Y", + aspect = DataAspect() +) + +lines!(ax, sol, idxs=(1,2)) f = DisplayAs.PNG(f) #hide \ No newline at end of file diff --git a/src/equations.jl b/src/equations.jl index c271284f..7ddb88bb 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -163,9 +163,6 @@ end trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) Normalized ODE equations for relativistic charged particle (x, γv) moving in static EM field with in-place form. - -The velocity is normalized by the characteristic velocity v0 (c); the magnetic field is normalized by the characteristic magnetic field B0; - the electric field is normalized by E0=v0 B0; the position is normalized by the L = v0 / Ω0, where Ω0 = q B0 / m and the time is normalized by Ω0^-1. """ function trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) _, E, B = p From 08e2aa696b01d5d6c4dd2ffcda58e246b5fa909b Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 28 Oct 2024 12:07:15 -0400 Subject: [PATCH 09/13] Remove magic numbers --- src/equations.jl | 6 +++--- test/runtests.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations.jl b/src/equations.jl index 7ddb88bb..2e5178a3 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -92,7 +92,7 @@ function trace_relativistic!(dy, y, p::TPTuple, t) γv = @view y[4:6] γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 - if γ²v² > 1e-20 + if γ²v² > eps(eltype(dy)) v̂ = SVector{3, eltype(dy)}(normalize(γv)) else # no velocity v̂ = SVector{3, eltype(dy)}(0, 0, 0) @@ -120,7 +120,7 @@ function trace_relativistic(y, p::TPTuple, t) γv = @view y[4:6] γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 - if γ²v² > 1e-20 + if γ²v² > eps(eltype(dy)) v̂ = SVector{3, eltype(y)}(normalize(γv)) else # no velocity v̂ = SVector{3, eltype(y)}(0, 0, 0) @@ -171,7 +171,7 @@ function trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) γv = @view y[4:6] γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 - if γ²v² > 1e-20 + if γ²v² > eps(eltype(dy)) v̂ = SVector{3, eltype(dy)}(normalize(γv)) else # no velocity v̂ = SVector{3, eltype(dy)}(0, 0, 0) diff --git a/test/runtests.jl b/test/runtests.jl index edf1c366..637964d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -285,7 +285,7 @@ end stateinit = [0.0, 0.0, 0.0, 0.5, 0.0, 0.0] prob = ODEProblem(trace_relativistic_normalized!, stateinit, tspan, param) sol = solve(prob, Vern6()) - @test sol.u[end][1] ≈ 0.42073549161010637 + @test sol.u[end][1] ≈ 0.38992532495827226 stateinit = zeros(6) prob = ODEProblem(trace_relativistic_normalized!, stateinit, tspan, param) sol = solve(prob, Vern6()) From 556ea6c46c190a0094416b5edc786b21bc65622b Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 28 Oct 2024 12:16:27 -0400 Subject: [PATCH 10/13] Fix typo --- src/equations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations.jl b/src/equations.jl index 2e5178a3..6ffe8ade 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -120,7 +120,7 @@ function trace_relativistic(y, p::TPTuple, t) γv = @view y[4:6] γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 - if γ²v² > eps(eltype(dy)) + if γ²v² > eps(eltype(y)) v̂ = SVector{3, eltype(y)}(normalize(γv)) else # no velocity v̂ = SVector{3, eltype(y)}(0, 0, 0) From 31ffce17ad1d73b498d0190fd7e3037938e256ed Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 28 Oct 2024 12:38:47 -0400 Subject: [PATCH 11/13] Import StaticArrays --- docs/examples/basics/demo_dimensionless.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/docs/examples/basics/demo_dimensionless.jl b/docs/examples/basics/demo_dimensionless.jl index c08df3e4..36395283 100644 --- a/docs/examples/basics/demo_dimensionless.jl +++ b/docs/examples/basics/demo_dimensionless.jl @@ -20,8 +20,7 @@ import DisplayAs #hide using TestParticle using TestParticle: qᵢ, mᵢ -using OrdinaryDiffEq - +using OrdinaryDiffEq, StaticArrays using CairoMakie CairoMakie.activate!(type = "png") #hide From 298da3a98a6e16d7146db7127585e6c579e645dd Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 28 Oct 2024 12:59:41 -0400 Subject: [PATCH 12/13] Update dimensionless demos --- docs/examples/basics/demo_dimensionless.jl | 2 +- docs/examples/basics/demo_dimensionless_dimensional.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/examples/basics/demo_dimensionless.jl b/docs/examples/basics/demo_dimensionless.jl index 36395283..8559919e 100644 --- a/docs/examples/basics/demo_dimensionless.jl +++ b/docs/examples/basics/demo_dimensionless.jl @@ -11,7 +11,7 @@ # After normalization, ``q=1, B=1, m=1`` so that the gyroradius ``r_L = mv_\perp/qB = v_\perp``. All the quantities are given in dimensionless units. # If the magnetic field is homogeneous and the initial perpendicular velocity is 4, then the gyroradius is 4. # To convert them to the original units, ``v_\perp = 4*U_0`` and ``r_L = 4*l_0``. -# Check [Demo: single tracing with additional diagnostics](@ref demo_savingcallback) for explaining the unit conversion. +# Check [Demo: single tracing with additional diagnostics](@ref demo_savingcallback) and [Demo: Dimensionless and Dimensional Tracing](@ref demo_dimensionless_dimensional) for explaining the unit conversion. # Tracing in dimensionless units is beneficial for many scenarios. For example, MHD simulations do not have intrinsic scales. Therefore, we can do dimensionless particle tracing in MHD fields, and then convert to any scale we would like. diff --git a/docs/examples/basics/demo_dimensionless_dimensional.jl b/docs/examples/basics/demo_dimensionless_dimensional.jl index 4f59c79f..a9687680 100644 --- a/docs/examples/basics/demo_dimensionless_dimensional.jl +++ b/docs/examples/basics/demo_dimensionless_dimensional.jl @@ -3,7 +3,7 @@ # id: demo_dimensionless_dimensional # date: 2023-12-19 # author: "[Hongyang Zhou](https://github.com/henry2004y)" -# julia: 1.9.4 +# julia: 1.11.1 # description: Tracing charged particle in both dimensional and dimensionless units. # --- @@ -95,7 +95,7 @@ xp, yp = let trange = range(tspan2..., length=40) end lines!(ax, xp, yp, linestyle=:dashdot, linewidth=5, color=Makie.wong_colors()[2]) invL = inv(1e3) -scale!(ax.scene, invrL, invrL) +scale!(ax.scene, invL, invL) f = DisplayAs.PNG(f) #hide From aeaae8d37f9fb7933c94efd2b7a8e4d0ed3daeb2 Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 28 Oct 2024 13:27:15 -0400 Subject: [PATCH 13/13] Optimization --- src/equations.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/equations.jl b/src/equations.jl index 6ffe8ade..7a15a160 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -90,15 +90,15 @@ function trace_relativistic!(dy, y, p::TPTuple, t) Ex, Ey, Ez = E(y, t) Bx, By, Bz = B(y, t) - γv = @view y[4:6] + γv = @views SVector{3, eltype(dy)}(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)) + v̂ = normalize(γv) else # no velocity v̂ = SVector{3, eltype(dy)}(0, 0, 0) end vmag = √(γ²v² / (1 + γ²v²/c^2)) - vx, vy, vz = vmag * v̂[1], vmag * v̂[2], vmag * v̂[3] + vx, vy, vz = vmag * v̂ dy[1], dy[2], dy[3] = vx, vy, vz dy[4] = q2m * (vy*Bz - vz*By + Ex) @@ -118,15 +118,15 @@ function trace_relativistic(y, p::TPTuple, t) Ex, Ey, Ez = E(y, t) Bx, By, Bz = B(y, t) - γv = @view y[4:6] + γv = @views SVector{3, eltype(y)}(y[4:6]) γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 if γ²v² > eps(eltype(y)) - v̂ = SVector{3, eltype(y)}(normalize(γv)) + v̂ = normalize(γv) else # no velocity v̂ = SVector{3, eltype(y)}(0, 0, 0) end vmag = √(γ²v² / (1 + γ²v²/c^2)) - vx, vy, vz = vmag * v̂[1], vmag * v̂[2], vmag * v̂[3] + vx, vy, vz = vmag * v̂ dx, dy, dz = vx, vy, vz dux = q2m * (vy*Bz - vz*By + Ex) @@ -169,15 +169,15 @@ function trace_relativistic_normalized!(dy, y, p::TPNormalizedTuple, t) Ex, Ey, Ez = E(y, t) Bx, By, Bz = B(y, t) - γv = @view y[4:6] + γv = @views SVector{3, eltype(dy)}(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)) + v̂ = 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] + vx, vy, vz = vmag * v̂ dy[1], dy[2], dy[3] = vx, vy, vz dy[4] = vy*Bz - vz*By + Ex @@ -197,7 +197,7 @@ requires the full particle trajectory `p.sol`. function trace_gc_drifts!(dx, x, p, t) q2m, E, B, sol = p xu = sol(t) - v = @view xu[4:6] + v = @views SVector{3, eltype(dx)}(xu[4:6]) abs_B(x) = norm(B(x)) gradient_B = ForwardDiff.gradient(abs_B, x) Bv = B(x)