From f1f61979cdc3bf12a6251ff745a6f375e5f9abe2 Mon Sep 17 00:00:00 2001 From: Hongyang Zhou Date: Mon, 28 Oct 2024 16:37:34 -0400 Subject: [PATCH] Relativistic momentum (#200) * Solve momentum in the relativistic case #198 * Add relativistic benchmark * Add utility methods * Add normalized relativistic demo * Update dimensionless demos --------- Co-authored-by: Beforerr <58776897+Beforerr@users.noreply.github.com> --- benchmark/benchmarks.jl | 27 ++-- docs/examples/advanced/demo_proton_dipole.jl | 6 +- docs/examples/basics/demo_dimensionless.jl | 63 ++++++++- .../basics/demo_dimensionless_dimensional.jl | 10 +- src/TestParticle.jl | 3 +- src/equations.jl | 102 +++++++-------- src/utility/utility.jl | 44 ++++++- test/runtests.jl | 121 ++++++------------ 8 files changed, 210 insertions(+), 166 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) 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/docs/examples/basics/demo_dimensionless.jl b/docs/examples/basics/demo_dimensionless.jl index a5c16348..8559919e 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 # --- @@ -11,17 +11,16 @@ # 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. -# 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 using TestParticle: qᵢ, mᵢ -using OrdinaryDiffEq - +using OrdinaryDiffEq, StaticArrays using CairoMakie CairoMakie.activate!(type = "png") #hide @@ -67,6 +66,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/docs/examples/basics/demo_dimensionless_dimensional.jl b/docs/examples/basics/demo_dimensionless_dimensional.jl index 2f839cfb..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. # --- @@ -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, invL, invL) f = DisplayAs.PNG(f) #hide diff --git a/src/TestParticle.jl b/src/TestParticle.jl index 0fbef57f..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 +export get_gyrofrequency, get_gyroperiod, get_gyroradius, get_velocity, get_energy, + energy2velocity export orbit, monitor export TraceProblem diff --git a/src/equations.jl b/src/equations.jl index 2065369a..7a15a160 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -80,71 +80,58 @@ 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 = @views SVector{3, eltype(dy)}(y[4:6]) + γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 + if γ²v² > eps(eltype(dy)) + 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̂ + 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 """ - trace_relativistic(y, p::TPTuple, t) -> SVector{6, Float64} + trace_relativistic(y, p::TPTuple, t) -> SVector{6} -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 = @views SVector{3, eltype(y)}(y[4:6]) + γ²v² = γv[1]^2 + γv[2]^2 + γv[3]^2 + if γ²v² > eps(eltype(y)) + 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̂ + 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] + _, E, B = p Ex, Ey, Ez = E(y, t) Bx, By, Bz = B(y, t) + γ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̂ = normalize(γv) + else # no velocity + v̂ = SVector{3, eltype(dy)}(0, 0, 0) + end + vmag = √(γ²v² / (1 + γ²v²)) + vx, vy, vz = vmag * v̂ + 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 @@ -209,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) diff --git a/src/utility/utility.jl b/src/utility/utility.jl index 5614e769..7e0c7113 100644 --- a/src/utility/utility.jl +++ b/src/utility/utility.jl @@ -104,4 +104,46 @@ end function get_gyroperiod(B::AbstractFloat=5e-9; q::AbstractFloat=qᵢ, m::AbstractFloat=mᵢ) ω = get_gyrofrequency(B; q, m) 2π / ω -end \ No newline at end of file +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 + +"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 9b9560fb..637964d6 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,54 +222,31 @@ 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 - - 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 @@ -297,39 +257,39 @@ end return SA[Ex, Ey, 0.0] 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) - 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 - - prob = ODEProblem(trace_relativistic, SA[stateinit...], tspan, param) - sol = solve(prob, Vern6(); dtmax=1e-10) x = sol.u[end][1:3] + @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 - @test calc_energy(sol)/(x[1]-x0[1]+x[2]-x0[2]) ≈ 1e5 + prob = ODEProblem(trace_relativistic, SA[stateinit...], tspan, param) + sol = solve(prob, Vern6(); dt=1e-10, adaptive=false) + + 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 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.38992532495827226 + 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 @@ -360,10 +320,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) @@ -390,23 +347,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) @@ -425,9 +378,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 @@ -490,8 +441,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)