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.4088344909683774e-10
 1.925859561433506e-9
 1.3194349022483877e-8
 8.812966040768264e-8
 5.494574873428513e-7
 3.4629406848752854e-6
 2.148429027387439e-5
 0.00012918278550991505
 0.000730736126206529
 0.004081541295463183
 0.022937388071097683
 0.12877413196982104
 0.6213053233871128
 1.5037537833398693
 2.6140719028720536
 3.793135407570003
 5.061853931919584
 6.437803884405918
 6.5594474868589705
u: 20-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 0.0, 2.99792458e6, 0.0]
 [8.331388431966723e-14, 0.0007221504129625887, 0.0, 0.0006917360626646803, 2.99792458e6, 0.0]
 [5.3254055810762926e-12, 0.005773581716849528, 0.0, 0.005530419442539567, 2.99792458e6, 0.0]
 [2.499651115455538e-10, 0.03955566325160339, 0.0, 0.03788972250462677, 2.9979245799999996e6, 0.0]
 [1.1151868872734071e-8, 0.26420607516324424, 0.0, 0.25307867569547343, 2.9979245799999894e6, 0.0]
 [4.3348258297628074e-7, 1.647232106970097, 0.0, 1.5778566784941657, 2.9979245799995847e6, 0.0]
 [1.7218430678788385e-5, 10.381634998250615, 0.0, 9.94439827050789, 2.997924579983507e6, 0.0]
 [0.0006627435105106637, 64.40828189135667, 0.0, 61.69563918927235, 2.997924579365169e6, 0.0]
 [0.023961430944091422, 387.2802470046967, 0.0, 370.96941080709803, 2.997924557047745e6, 0.0]
 [0.7666984765891037, 2190.691615362459, 0.0, 2098.4276479281075, 2.9979238455923214e6, 0.0]
 [23.919455956688612, 12236.121801781845, 0.0, 11720.780831147582, 2.9979016679522214e6, 0.0]
 [755.3946673017898, 68759.02706562518, 0.0, 65863.14679400499, 2.997201000051503e6, 0.0]
 [23779.898179761392, 385076.88449089293, 0.0, 368858.8462254769, 2.9751462122870744e6, 0.0]
 [538092.6069289364, 1.7546036335232435e6, 0.0, 1.6807061411711655e6, 2.4824946495780963e6, 0.0]
 [2.722852914234935e6, 3.103176622232147e6, 0.0, 2.9724825521993074e6, 389749.2026368658, 0.0]
 [5.644538065119795e6, 1.8630757004553287e6, 0.0, 1.7846116331313036e6, -2.4088842568155425e6, 0.0]
 [5.888573180714667e6, -1.4778767264858738e6, 0.0, -1.4156298254249124e6, -2.6426414224630124e6, 0.0]
 [2.7045489242454525e6, -3.1007275429247064e6, 0.0, -2.9701303309060545e6, 407282.28925659746, 0.0]
 [21223.86367614255, -363859.6938069101, 0.0, -348528.9632499821, 2.9775945950921746e6, 0.0]
 [1.7093031248411283, -6.190710236927165, 0.0, 0.35320726048385764, 2.997922942687056e6, 0.0]

Compare performance

julia

julia> using Chairmarks

julia> @b solve(prob_SI, Vern9()) # "SI units (unitless)"
15.609 μs (528 allocs: 31.445 KiB)

julia> @b solve(prob_p, Vern9()) # "Partitioned (unitful)"
167.222 μs (3437 allocs: 136.211 KiB)

julia> @b solve(prob, Vern9()) # "Basic (unitful)"
1.319 ms (20053 allocs: 384.453 KiB)

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

Unit support ​

Heterogeneous arrays with ArrayPartition ​

Tracing in standard SI units ​

Compare performance ​

Related API ​

Unit support

Heterogeneous arrays with `ArrayPartition`

Tracing in standard SI units

Compare performance

Related API