Fix dimensionless tracing (#122)

* Bump minor version

* Add Boris method to benchmark

* Update prepare args

* Clarify spatial coordinates in dimensionless units

Project.toml

@@ -1,7 +1,7 @@
 name = "TestParticle"
 uuid = "953b605b-f162-4481-8f7f-a191c2bb40e3"
 authors = ["Hongyang Zhou <[email protected]>, and Tiancheng Liu <[email protected]>"]
-version = "0.5.1"
+version = "0.6.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

benchmark/benchmarks.jl

@@ -50,8 +50,12 @@ SUITE["trace"]["analytic field"]["out of place"] = @benchmarkable solve($prob_oo
 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 = let dt = 1/7, paramBoris = BorisMethod(param_numeric)
+    TraceProblem(stateinit, tspan, dt, paramBoris)
+end
 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 trace_trajectory($prob_boris; savestepinterval=10)
 
 param_td = prepare(E_td, B_td, F_td)
 prob_ip = ODEProblem(trace!, stateinit, tspan, param_td) # in place

docs/examples/basics/demo_dimensionless.jl

@@ -8,60 +8,47 @@
 # ---
 
 # This example shows how to trace charged particles in dimensionless units.
-# Let ``B_0`` be the reference magnetic field strength and ``U_0`` the reference velocity. The Lorentz equation without the electric field can be written as
-#
-# ```math
-# \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = \Omega \mathbf{u}\times\mathbf{b}
-# ```
-#
-# where ``\Omega = qB_0/m`` is the reference gyrofrequency and  ``\mathbf{b} = \mathbf{B}/B_0`` is the normalized magnetic field.
-#
-# Then we can normalize the time scale to the gyroperiod ``T=2\pi/\Omega``. We use the gyroradius for the spatial scale, which is proportional to the velocity scale
-#
-# ```math
-# r_L = \frac{U_0}{\Omega}
-# ```
-#
-# For instance, if the magnetic field is homogeneous and the initial perpendicular velocity is 4 (in the normalized units), then the gyroradius is 4.
+# 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.
+
 # 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!`.
 
 import DisplayAs #hide
 using TestParticle
+using TestParticle: qᵢ, mᵢ
 using OrdinaryDiffEq
-using Meshes
 
 using CairoMakie
 CairoMakie.activate!(type = "png")
 
-x = range(-10, 10, length=15)
-y = range(-10, 10, length=20)
-z = range(-10, 10, length=25)
-B = fill(0.0, 3, length(x), length(y), length(z)) # [B₀]
-E = fill(0.0, 3, length(x), length(y), length(z)) # [E₀]
-
-B₀ = 10e-9
+## Number of cells for the field along each dimension
+nx, ny, nz = 4, 6, 8
+## Unit conversion factors between SI and dimensionless units
+B₀ = 10e-9            # [T]
+Ω = abs(qᵢ) * B₀ / mᵢ # [1/s]
+t₀ = 1 / Ω            # [s]
+U₀ = 1.0              # [m/s]
+l₀ = U₀ * t₀          # [m]
+E₀ = U₀*B₀            # [V/m]
+## All quantities are in dimensionless units
+x = range(-10, 10, length=nx) # [l₀]
+y = range(-10, 10, length=ny) # [l₀]
+z = range(-10, 10, length=nz) # [l₀]
 
+B = fill(0.0, 3, nx, ny, nz) # [B₀]
 B[3,:,:,:] .= 1.0
+E = fill(0.0, 3, nx, ny, nz) # [E₀]
 
-Δx = x[2] - x[1]
-Δy = y[2] - y[1]
-Δz = z[2] - z[1]
-
-mesh = CartesianGrid((length(x)-1, length(y)-1, length(z)-1),
-   (x[1], y[1], z[1]),
-   (Δx, Δy, Δz))
-
-param = prepare(mesh, E, B, B₀; species=Proton)
-
-Ω = param[1]
-U₀ = 1.0
+param = prepare(x, y, z, E, B; species=User)
 
 x0 = [0.0, 0.0, 0.0] # initial position [l₀]
-u0 = [4*Ω*U₀, 0.0, 0.0] # initial velocity [v₀]
+u0 = [4.0, 0.0, 0.0] # initial velocity [v₀]
 stateinit = [x0..., u0...]
-tspan = (0.0, 2π/Ω) # one gyroperiod
+tspan = (0.0, π) # half gyroperiod
 
 prob = ODEProblem(trace_normalized!, stateinit, tspan, param)
 sol = solve(prob, Vern9())

docs/examples/basics/demo_dimensionless_dimensional.jl

@@ -43,37 +43,37 @@ using StaticArrays
 using CairoMakie
 CairoMakie.activate!(type = "png")
 
-## Unit conversion factors
-const B₀ = 1e-8 # [T]
-const U₀ = c    # [m/s]
-const E₀ = U₀*B₀ # [V/m]
-const Ω = abs(qᵢ) * B₀ / mᵢ
-const t₀ = 1 / Ω  # [s]
-const l₀ = U₀ * t₀ # [m]
-
-const Emag = 1e-8 # [V/m]
+## Unit conversion factors between SI and dimensionless units
+const B₀ = 1e-8             # [T]
+const U₀ = c                # [m/s]
+const E₀ = U₀*B₀            # [V/m]
+const Ω = abs(qᵢ) * B₀ / mᵢ # [1/s]
+const t₀ = 1 / Ω            # [s]
+const l₀ = U₀ * t₀          # [m]
+## Electric field magnitude in SI units
+const Emag = 1e-8           # [V/m]
 ### Solving in SI units
 B(x) = SA[0, 0, B₀]
 E(x) = SA[Emag, 0.0, 0.0]
 
 x0 = [0.0, 0.0, 0.0] # [m]
 v0 = [0.0, 0.01c, 0.0] # [m/s]
 stateinit1 = [x0..., v0...]
-tspan1 = (0, 2π/Ω) # [s]
+tspan1 = (0, 2π*t₀) # [s]
 
 param1 = prepare(E, B, species=Proton)
 prob1 = ODEProblem(trace!, stateinit1, tspan1, param1)
 sol1 = solve(prob1, Vern9(); reltol=1e-4, abstol=1e-6)
 
 ### Solving in dimensionless units
-B_normalize(x) = SA[0, 0, 1.0]
+B_normalize(x) = SA[0, 0, B₀/B₀]
 E_normalize(x) = SA[Emag/E₀, 0.0, 0.0]
-## Default User type has q=1, m=1. Here we also set B=0. 
-param2 = prepare(E_normalize, B_normalize, 1.0; species=User)
+## For full EM problems, the normalization of E and B should be done separately.
+param2 = prepare(E_normalize, B_normalize; species=User)
 ## Scale initial quantities by the conversion factors
 x0 ./= l₀
 v0 ./= U₀
-tspan2 = (0, 2π/Ω/t₀)
+tspan2 = (0, 2π)
 stateinit2 = [x0..., v0...]
 
 prob2 = ODEProblem(trace_normalized!, stateinit2, tspan2, param2)

docs/examples/basics/demo_dimensionless_periodic.jl

@@ -7,47 +7,47 @@
 # description: Tracing charged particle in dimensionless units and periodic boundary
 # ---
 
-# This example shows how to trace charged particles in dimensionless units and EM fields with periodic boundaries.
-# For details about dimensionless units, please check another [example](@ref demo_dimensionless).
+# This example shows how to trace charged particles in dimensionless units and EM fields with periodic boundaries in a 2D spatial domain.
+# For details about dimensionless units, please check [Demo: dimensionless tracing](@ref demo_dimensionless).
 
 # Now let's demonstrate this with `trace_normalized!`.
 
 import DisplayAs #hide
 using TestParticle
+using TestParticle: qᵢ, mᵢ
 using OrdinaryDiffEq
-using Meshes
 using StaticArrays
 using CairoMakie
 CairoMakie.activate!(type = "png")
 
-x = range(-10, 10, length=15)
-y = range(-10, 10, length=20)
-B = fill(0.0, 3, length(x), length(y)) # [B₀]
+## Number of cells for the field along each dimension
+nx, ny = 4, 6
+## Unit conversion factors between SI and dimensionless units
+B₀ = 10e-9            # [T]
+Ω = abs(qᵢ) * B₀ / mᵢ # [1/s]
+t₀ = 1 / Ω            # [s]
+U₀ = 1.0              # [m/s]
+l₀ = U₀ * t₀          # [m]
+E₀ = U₀*B₀            # [V/m]
 
-E(x) = SA[0.0, 0.0, 0.0]
-
-B₀ = 10e-9
+x = range(-10, 10, length=nx) # [l₀]
+y = range(-10, 10, length=ny) # [l₀]
 
+B = fill(0.0, 3, nx, ny) # [B₀]
 B[3,:,:] .= 1.0
 
-Δx = x[2] - x[1]
-Δy = y[2] - y[1]
-
-grid = CartesianGrid((length(x)-1, length(y)-1), (x[1], y[1]), (Δx, Δy))
+E(x) = SA[0.0, 0.0, 0.0] # [E₀]
 
 ## If bc == 1, we set a NaN value outside the domain (default);
 ## If bc == 2, we set periodic boundary conditions.
-param = prepare(grid, E, B, B₀; species=Proton, bc=2)
-
-Ω = param[1]
-U₀ = 1.0
+param = prepare(x, y, E, B; species=User, bc=2)
 
 # Note that we set a radius of 10, so the trajectory extent from -20 to 0 in y, which is beyond the original y range.
 
 x0 = [0.0, 0.0, 0.0] # initial position [l₀]
-u0 = [10*Ω*U₀, 0.0, 0.0] # initial velocity [v₀]
+u0 = [10.0, 0.0, 0.0] # initial velocity [v₀]
 stateinit = [x0..., u0...]
-tspan = (0.0, 1.5π/Ω) # 3/4 gyroperiod
+tspan = (0.0, 1.5π) # 3/4 gyroperiod
 
 prob = ODEProblem(trace_normalized!, stateinit, tspan, param)
 sol = solve(prob, Vern9())

docs/examples/basics/demo_output_func.jl

@@ -8,21 +8,26 @@
 # ---
 
 # This example demonstrates tracing multiple protons in an analytic E field and numerical B field.
-# It also combines one type of normalization using a reference velocity `U₀`, a reference magnetic field `B₀`, and an initial reference gyroradius `rL`.
-# One way to think of this is that a proton with initial perpendicular velocity `U₀` will gyrate with a radius `rL`.
-# Check [demo_ensemble](@ref demo_ensemble) for basic usages of the ensemble problem.
+# See [Demo: single tracing with additional diagnostics](@ref demo_savingcallback) for explaining the unit conversion.
+# Also check [demo_ensemble](@ref demo_ensemble) for basic usages of the ensemble problem.
+
 # The `output_func` argument can be used to change saving outputs. It works as a reduction function, but here we demonstrate how to add additional outputs.
 # Besides the regular outputs, we also save the magnetic field along the trajectory, together with the parallel velocity.
 
 import DisplayAs #hide
 using TestParticle
+using TestParticle: qᵢ, mᵢ
 using OrdinaryDiffEq
 using StaticArrays
 using Statistics
 using LinearAlgebra
+using Random
 using CairoMakie
 CairoMakie.activate!(type = "png")
 
+seed = 1 # seed for random number
+Random.seed!(seed)
+
 "Set initial state for EnsembleProblem."
 function prob_func(prob, i, repeat)
    B0 = prob.p[3](prob.u0)
@@ -42,39 +47,45 @@ function prob_func(prob, i, repeat)
    prob
 end
 
-## Analytical electric field
-E(x) = SA[0.0, 0.0, 0.0]
 ## Number of cells for the field along each dimension
 nx, ny, nz = 4, 6, 8
-## Spatial extent along each dimension
-x = range(-10.0, 10.0, length=nx)
-y = range(-10.0, 10.0, length=ny)
-z = range(-10.0, 10.0, length=nz)
-
-## Numerical magnetic field
+## Spatial coordinates given in customized units
+x = range(0, 1, length=nx)
+y = range(0, 1, length=ny)
+z = range(0, 1, length=nz)
+## Numerical magnetic field given in customized units
 B = Array{Float32, 4}(undef, 3, nx, ny, nz)
 
 B[1,:,:,:] .= 0.0
 B[2,:,:,:] .= 0.0
-B[3,:,:,:] .= 1.0
-
-Bmag = @views hypot.(B[1,:,:,:], B[2,:,:,:], B[3,:,:,:])
-
-B₀ = sqrt(mean(vec(Bmag) .^ 2))
+B[3,:,:,:] .= 2.0
 
+## Reference values for unit conversions between the customized and dimensionless units
+const B₀ = let Bmag = @views hypot.(B[1,:,:,:], B[2,:,:,:], B[3,:,:,:])
+   sqrt(mean(vec(Bmag) .^ 2))
+end
 const U₀ = 1.0
-const rL = 4.0
-
-Ω = q2mc = U₀ / (rL*B₀)
+const l₀ = 2*nx
+const t₀ = l₀ / U₀
+const E₀ = U₀ * B₀
+
+### Convert from customized to default dimensionless units
+## Dimensionless spatial extents [l₀]
+x /= l₀
+y /= l₀
+z /= l₀
+## For full EM problems, the normalization of E and B should be done separately.
+B ./= B₀
+E(x) = SA[0.0/E₀, 0.0/E₀, 0.0/E₀]
 
 ## By default User type assumes q=1, m=1
 ## bc=2 uses periodic boundary conditions
-param = prepare(x, y, z, E, B, Ω; species=User, bc=2)
+param = prepare(x, y, z, E, B; species=User, bc=2)
 
 x0 = [0.0, 0.0, 0.0] # initial position [l₀]
-u0 = [1U₀, 0.0, 0.0] # initial velocity [v₀]
+u0 = [1.0, 0.0, 0.0] # initial velocity [v₀]
 stateinit = [x0..., u0...]
-tspan = (0.0, 2π/Ω) # one averaged gyroperiod based on B₀
+tspan = (0.0, 2π) # one averaged gyroperiod based on B₀
 
 saveat = tspan[2] / 40 # save interval