TestParticle.jl

Appearance

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 OrdinaryDiffEq
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, Vern7())
sol[end]
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, Vern7())
sol_p[end]
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, Vern7())
sol_SI[end]
0.0

Compare performance ​

julia
julia> using Chairmarks

julia> @b solve(prob_SI, Vern7()) # "SI units (unitless)"
13.040 μs (503 allocs: 30.164 KiB)

julia> @b solve(prob_p, Vern7()) # "Partitioned (unitful)"
117.486 μs (2361 allocs: 96.922 KiB)

julia> @b solve(prob, Vern7()) # "Basic (unitful)"
733.121 μs (11454 allocs: 231.297 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