Unit support
This example shows how to trace charged particles with Unitful units.
TL;DR: Tracing with units is convenient although not as performant as tracing in dimensionless units.
julia
using TestParticle, StaticArrays
using Unitful
using OrdinaryDiffEqVerner
import TestParticle as TP
const Bmag = 1e-8u"T"
const Ω = abs(Unitful.q) * Bmag / Unitful.mp
const t_max = 1 / Ω |> u"s"
const Emag = 1e-8u"V/m"
B(x) = SA[0.0u"T", 0.0u"T", Bmag]
E(x) = SA[Emag, 0.0u"V/m", 0.0u"V/m"]
x0 = [0.0, 0.0, 0.0] * u"m" # [m]
v0 = [0.0, 0.01, 0.0] * Unitful.c0 # [m/s]
u0 = [x0..., v0...]
tspan = (0.0u"s", 2π * t_max) # [s]
param = prepare(E, B, q = Unitful.q, m = Unitful.mp)
prob = ODEProblem(trace!, u0, tspan, param)
sol = solve(prob, Vern9())retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 20-element Vector{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}:
0.0 s
2.4088344909683774e-10 s
1.925859561433506e-9 s
1.3194349022483877e-8 s
8.812966040768264e-8 s
5.494574873428513e-7 s
3.4629406848752854e-6 s
2.148429027387439e-5 s
0.00012918278550991505 s
0.000730736126206529 s
0.004081541295463183 s
0.022937388071097683 s
0.12877413196982104 s
0.6213053233871128 s
1.5037537833398693 s
2.6140719028720536 s
3.793135407570003 s
5.061853931919584 s
6.437803884405918 s
6.5594474868589705 s
u: 20-element Vector{Vector{Unitful.Quantity{Float64}}}:
[0.0 m, 0.0 m, 0.0 m, 0.0 m s^-1, 2.99792458e6 m s^-1, 0.0 m s^-1]
[8.331388431966723e-14 m, 0.0007221504129625887 m, 0.0 m, 0.0006917360626646803 m s^-1, 2.99792458e6 m s^-1, 0.0 m s^-1]
[5.3254055810762926e-12 m, 0.005773581716849528 m, 0.0 m, 0.005530419442539567 m s^-1, 2.99792458e6 m s^-1, 0.0 m s^-1]
[2.499651115455538e-10 m, 0.03955566325160339 m, 0.0 m, 0.03788972250462677 m s^-1, 2.9979245799999996e6 m s^-1, 0.0 m s^-1]
[1.1151868872734071e-8 m, 0.26420607516324424 m, 0.0 m, 0.25307867569547343 m s^-1, 2.9979245799999894e6 m s^-1, 0.0 m s^-1]
[4.3348258297628074e-7 m, 1.647232106970097 m, 0.0 m, 1.5778566784941657 m s^-1, 2.9979245799995847e6 m s^-1, 0.0 m s^-1]
[1.7218430678788385e-5 m, 10.381634998250615 m, 0.0 m, 9.94439827050789 m s^-1, 2.997924579983507e6 m s^-1, 0.0 m s^-1]
[0.0006627435105106637 m, 64.40828189135667 m, 0.0 m, 61.69563918927235 m s^-1, 2.997924579365169e6 m s^-1, 0.0 m s^-1]
[0.023961430944091422 m, 387.2802470046967 m, 0.0 m, 370.96941080709803 m s^-1, 2.997924557047745e6 m s^-1, 0.0 m s^-1]
[0.7666984765891037 m, 2190.691615362459 m, 0.0 m, 2098.4276479281075 m s^-1, 2.9979238455923214e6 m s^-1, 0.0 m s^-1]
[23.919455956688612 m, 12236.121801781845 m, 0.0 m, 11720.780831147582 m s^-1, 2.9979016679522214e6 m s^-1, 0.0 m s^-1]
[755.3946673017898 m, 68759.02706562518 m, 0.0 m, 65863.14679400499 m s^-1, 2.997201000051503e6 m s^-1, 0.0 m s^-1]
[23779.898179761392 m, 385076.88449089293 m, 0.0 m, 368858.8462254769 m s^-1, 2.9751462122870744e6 m s^-1, 0.0 m s^-1]
[538092.6069289364 m, 1.7546036335232435e6 m, 0.0 m, 1.6807061411711655e6 m s^-1, 2.4824946495780963e6 m s^-1, 0.0 m s^-1]
[2.722852914234935e6 m, 3.103176622232147e6 m, 0.0 m, 2.9724825521993074e6 m s^-1, 389749.2026368658 m s^-1, 0.0 m s^-1]
[5.644538065119795e6 m, 1.8630757004553287e6 m, 0.0 m, 1.7846116331313036e6 m s^-1, -2.4088842568155425e6 m s^-1, 0.0 m s^-1]
[5.888573180714667e6 m, -1.4778767264858738e6 m, 0.0 m, -1.4156298254249124e6 m s^-1, -2.6426414224630124e6 m s^-1, 0.0 m s^-1]
[2.7045489242454525e6 m, -3.1007275429247064e6 m, 0.0 m, -2.9701303309060545e6 m s^-1, 407282.28925659746 m s^-1, 0.0 m s^-1]
[21223.86367614255 m, -363859.6938069101 m, 0.0 m, -348528.9632499821 m s^-1, 2.9775945950921746e6 m s^-1, 0.0 m s^-1]
[1.7093031248411283 m, -6.190710236927165 m, 0.0 m, 0.35320726048385764 m s^-1, 2.997922942687056e6 m s^-1, 0.0 m s^-1]Heterogeneous arrays with ArrayPartition
ArrayPartitions in DiffEq can be used for heterogeneous arrays. See https://docs.sciml.ai/DiffEqDocs/stable/features/diffeq_arrays for more details.
julia
using RecursiveArrayTools
u0_p = ArrayPartition(x0, v0)
prob_p = ODEProblem(trace!, u0_p, tspan, param)
sol_p = solve(prob_p, Vern9())retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 20-element Vector{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}:
0.0 s
2.4088344909683774e-10 s
1.925859561433506e-9 s
1.3194349022483877e-8 s
8.812966040768264e-8 s
5.494574873428513e-7 s
3.4629406848752854e-6 s
2.148429027387439e-5 s
0.00012918278550991505 s
0.000730736126206529 s
0.004081541295463183 s
0.022937388071097683 s
0.12877413196982104 s
0.6213053233871128 s
1.5037537833398693 s
2.6140719028720536 s
3.793135407570003 s
5.061853931919584 s
6.437803884405918 s
6.5594474868589705 s
u: 20-element Vector{RecursiveArrayTools.ArrayPartition{Unitful.Quantity{Float64}, Tuple{Vector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, Vector{Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}}}}}:
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[0.0 m, 0.0 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[0.0 m s^-1, 2.99792458e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[8.331388431966723e-14 m, 0.0007221504129625887 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[0.0006917360626646803 m s^-1, 2.99792458e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[5.3254055810762926e-12 m, 0.005773581716849528 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[0.005530419442539567 m s^-1, 2.99792458e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[2.499651115455538e-10 m, 0.03955566325160339 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[0.03788972250462677 m s^-1, 2.9979245799999996e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[1.1151868872734071e-8 m, 0.26420607516324424 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[0.25307867569547343 m s^-1, 2.9979245799999894e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[4.3348258297628074e-7 m, 1.647232106970097 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[1.5778566784941657 m s^-1, 2.9979245799995847e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[1.7218430678788385e-5 m, 10.381634998250615 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[9.94439827050789 m s^-1, 2.997924579983507e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[0.0006627435105106637 m, 64.40828189135667 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[61.69563918927235 m s^-1, 2.997924579365169e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[0.023961430944091422 m, 387.2802470046967 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[370.96941080709803 m s^-1, 2.997924557047745e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[0.7666984765891037 m, 2190.691615362459 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[2098.4276479281075 m s^-1, 2.9979238455923214e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[23.919455956688612 m, 12236.121801781845 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[11720.780831147582 m s^-1, 2.9979016679522214e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[755.3946673017898 m, 68759.02706562518 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[65863.14679400499 m s^-1, 2.997201000051503e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[23779.898179761392 m, 385076.88449089293 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[368858.8462254769 m s^-1, 2.9751462122870744e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[538092.6069289364 m, 1.7546036335232435e6 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[1.6807061411711655e6 m s^-1, 2.4824946495780963e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[2.722852914234935e6 m, 3.103176622232147e6 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[2.9724825521993074e6 m s^-1, 389749.2026368658 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[5.644538065119795e6 m, 1.8630757004553287e6 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[1.7846116331313036e6 m s^-1, -2.4088842568155425e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[5.888573180714667e6 m, -1.4778767264858738e6 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[-1.4156298254249124e6 m s^-1, -2.6426414224630124e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[2.7045489242454525e6 m, -3.1007275429247064e6 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[-2.9701303309060545e6 m s^-1, 407282.28925659746 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[21223.86367614255 m, -363859.6938069101 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[-348528.9632499821 m s^-1, 2.9775945950921746e6 m s^-1, 0.0 m s^-1])
(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[1.7093031248411283 m, -6.190710236927165 m, 0.0 m], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}[0.35320726048385764 m s^-1, 2.997922942687056e6 m s^-1, 0.0 m s^-1])Tracing in standard SI units
julia
const Bmag_SI = Bmag / u"T"
const Emag_SI = Emag / u"V/m"
### Solving in SI units
B_SI(x) = SA[0, 0, Bmag_SI]
E_SI(x) = SA[Emag_SI, 0.0, 0.0]
## Initial conditions
x0_SI = [0.0, 0.0, 0.0]
v0_SI = [0.0, 0.01TP.c, 0.0]
u0_SI = [x0_SI..., v0_SI...]
tspan_SI = tspan ./ u"s" # [s]
param_SI = prepare(E_SI, B_SI, species = Proton)
prob_SI = ODEProblem(trace!, u0_SI, tspan_SI, param_SI)
sol_SI = solve(prob_SI, Vern9())retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 20-element Vector{Float64}:
0.0
2.4094789678540764e-10
1.8595945654928207e-9
1.2255312606307927e-8
7.842088982695753e-8
4.5778877484385586e-7
2.934447935719967e-6
1.823935528733788e-5
0.00011466221019483176
0.000678969151658004
0.003863557038823506
0.02245766265908049
0.13392290942149107
0.63407355466219
1.5224182814287581
2.6353557630902236
3.815720891592856
5.085210308104797
6.461798280084288
6.5594474868589705
u: 20-element Vector{Vector{Float64}}:
[0.0, 0.0, 0.0, 0.0, 2.99792458e6, 0.0]
[8.33118677981114e-14, 0.0007223436222722765, 0.0, 0.0006915343018935781, 2.99792458e6, 0.0]
[4.962461247962274e-12, 0.005574924256725347, 0.0, 0.0053371432031983235, 2.99792458e6, 0.0]
[2.1553081889055898e-10, 0.0367405028980344, 0.0, 0.03517345102720973, 2.99792458e6, 0.0]
[8.825191203083897e-9, 0.23509991319770776, 0.0, 0.22507245767186368, 2.9979245799999917e6, 0.0]
[3.007397750335094e-7, 1.3724162205524373, 0.0, 1.3138800755264561, 2.9979245799997123e6, 0.0]
[1.2357010299517912e-5, 8.79725359521358, 0.0, 8.422034106717247, 2.9979245799881704e6, 0.0]
[0.00047739685971894246, 54.680211536484535, 0.0, 52.34799719464109, 2.9979245795429656e6, 0.0]
[0.018866928733539685, 343.7486576498617, 0.0, 329.0871278782458, 2.997924561937788e6, 0.0]
[0.6615469360336905, 2035.4981654799858, 0.0, 1948.6803225907272, 2.9979239466695026e6, 0.0]
[21.420753968438788, 11582.626202777306, 0.0, 11088.605309536182, 2.9979040728885042e6, 0.0]
[723.7239926329265, 67321.19217623271, 0.0, 64449.81612626845, 2.997231724434262e6, 0.0]
[25702.498947924723, 400391.7337537003, 0.0, 383314.2701085212, 2.9733183494055616e6, 0.0]
[559451.955712179, 1.786294570069213e6, 0.0, 1.7101057680006837e6, 2.4623344859529277e6, 0.0]
[2.7774122436741334e6, 3.1114067892124127e6, 0.0, 2.978699874457961e6, 338975.10344794666, 0.0]
[5.682614222362476e6, 1.8160371979311747e6, 0.0, 1.7385817230099782e6, -2.4423135936496574e6, 0.0]
[5.862370322100434e6, -1.5324943456566713e6, 0.0, -1.4671264237689467e6, -2.61440269867679e6, 0.0]
[2.6451863465969297e6, -3.0935042509279093e6, 0.0, -2.96155459317695e6, 465561.2739456176, 0.0]
[14721.308087227339, -303275.0403806198, 0.0, -290333.5038263785, 2.9838311682200776e6, 0.0]
[21.031791049507152, -11006.288245703789, 0.0, -10530.566068053173, 2.997904445261294e6, 0.0]Compare performance
julia
julia> using Chairmarks
julia> @b solve(prob_SI, Vern9()) # "SI units (unitless)"
15.819 μs (534 allocs: 31.766 KiB)
julia> @b solve(prob_p, Vern9()) # "Partitioned (unitful)"
118.381 μs (2511 allocs: 97.727 KiB)
julia> @b solve(prob, Vern9()) # "Basic (unitful)"
708.423 μs (21072 allocs: 376.484 KiB)Related API
TestParticle.trace! Function
julia
trace!(dy, y, p, t)ODE equations for charged particle moving in EM field and external force field with in-place form.
source