TestParticle.jl

Appearance

Field Interpolation ​

A robust field interpolation is the prerequisite for pushing particles. This example demonstrates the construction of scalar/vector field interpolators for Cartesian/Spherical grids. If the field is analytic, you can directly pass the generated function to prepare.

julia
import TestParticle as TP
using StaticArrays
using Chairmarks

function setup_spherical_field()
   r = logrange(0.1, 10.0, length = 16)
   r_uniform = range(0.1, 10.0, length = 16)
   θ = range(0, π, length = 16)
   ϕ = range(0, 2π, length = 16)

   B₀ = 1e-8 # [nT]
   B = zeros(3, length(r), length(θ), length(ϕ)) # vector
   A = zeros(length(r), length(θ), length(ϕ)) # scalar

   for (iθ, θ_val) in enumerate(θ)
      sinθ, cosθ = sincos(θ_val)
      B[1, :, iθ, :] .= B₀ * cosθ
      B[2, :, iθ, :] .= -B₀ * sinθ
      A[:, iθ, :] .= B₀ * sinθ
   end

   B_field_nu = TP.getinterp(TP.StructuredGrid, B, r, θ, ϕ)
   A_field_nu = TP.getinterp_scalar(TP.StructuredGrid, A, r, θ, ϕ)
   B_field = TP.getinterp(TP.StructuredGrid, B, r_uniform, θ, ϕ)
   A_field = TP.getinterp_scalar(TP.StructuredGrid, A, r_uniform, θ, ϕ)

   return B_field_nu, A_field_nu, B_field, A_field
end

function setup_cartesian_field()
   x = range(-10, 10, length = 16)
   y = range(-10, 10, length = 16)
   z = range(-10, 10, length = 16)
   B = zeros(3, length(x), length(y), length(z)) # vector
   B[3, :, :, :] .= 10e-9
   A = zeros(length(x), length(y), length(z)) # scalar
   A[:, :, :] .= 10e-9

   B_field = TP.getinterp(B, x, y, z)
   A_field = TP.getinterp_scalar(A, x, y, z)

   return B_field, A_field
end

function setup_cartesian_nonuniform_field()
   x = logrange(0.1, 10.0, length = 16)
   y = range(-10, 10, length = 16)
   z = range(-10, 10, length = 16)
   B = zeros(3, length(x), length(y), length(z)) # vector
   B[3, :, :, :] .= 10e-9
   A = zeros(length(x), length(y), length(z)) # scalar
   A[:, :, :] .= 10e-9

   B_field = TP.getinterp(TP.RectilinearGrid, B, x, y, z)
   A_field = TP.getinterp_scalar(TP.RectilinearGrid, A, x, y, z)

   return B_field, A_field
end

B_sph_nu, A_sph_nu, B_sph, A_sph = setup_spherical_field();
B_car, A_car = setup_cartesian_field();
B_car_nu, A_car_nu = setup_cartesian_nonuniform_field();
loc = SA[1.0, 1.0, 1.0];
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 1.0
 1.0
 1.0

Gridded spherical interpolation ​

julia
julia> @be B_sph_nu($loc)
Benchmark: 3128 samples with 78 evaluations
 min    374.423 ns (2 allocs: 64 bytes)
 median 377.756 ns (2 allocs: 64 bytes)
 mean   381.865 ns (2 allocs: 64 bytes)
 max    644.154 ns (2 allocs: 64 bytes)

julia> @be A_sph_nu($loc)
Benchmark: 3124 samples with 88 evaluations
 min    332.557 ns (2 allocs: 48 bytes)
 median 334.841 ns (2 allocs: 48 bytes)
 mean   339.021 ns (2 allocs: 48 bytes)
 max    911.830 ns (2 allocs: 48 bytes)

Uniform spherical interpolation ​

julia
julia> @be B_sph($loc)
Benchmark: 3154 samples with 213 evaluations
 min    135.606 ns (2 allocs: 64 bytes)
 median 137.300 ns (2 allocs: 64 bytes)
 mean   139.019 ns (2 allocs: 64 bytes)
 max    296.329 ns (2 allocs: 64 bytes)

julia> @be A_sph($loc)
Benchmark: 3180 samples with 285 evaluations
 min    100.470 ns (2 allocs: 48 bytes)
 median 102.789 ns (2 allocs: 48 bytes)
 mean   103.867 ns (2 allocs: 48 bytes)
 max    192.782 ns (2 allocs: 48 bytes)

Uniform Cartesian interpolation ​

julia
julia> @be B_car($loc)
Benchmark: 3207 samples with 401 evaluations
 min    71.130 ns (2 allocs: 64 bytes)
 median 72.828 ns (2 allocs: 64 bytes)
 mean   73.627 ns (2 allocs: 64 bytes)
 max    131.968 ns (2 allocs: 64 bytes)

julia> @be A_car($loc)
Benchmark: 2581 samples with 463 evaluations
 min    61.886 ns (2 allocs: 48 bytes)
 median 80.430 ns (2 allocs: 48 bytes)
 mean   78.848 ns (2 allocs: 48 bytes)
 max    539.674 ns (2 allocs: 48 bytes)

Non-uniform Cartesian interpolation ​

julia
julia> @be B_car_nu($loc)
Benchmark: 3129 samples with 98 evaluations
 min    298.827 ns (2 allocs: 64 bytes)
 median 300.969 ns (2 allocs: 64 bytes)
 mean   304.155 ns (2 allocs: 64 bytes)
 max    461.469 ns (2 allocs: 64 bytes)

julia> @be A_car_nu($loc)
Benchmark: 3154 samples with 100 evaluations
 min    291.950 ns (2 allocs: 48 bytes)
 median 293.260 ns (2 allocs: 48 bytes)
 mean   296.872 ns (2 allocs: 48 bytes)
 max    509.050 ns (2 allocs: 48 bytes)

Based on the benchmarks, for the same grid size, gridded interpolation (StructuredGrid with non-uniform ranges, RectilinearGrid) is 2x slower than uniform mesh interpolation (StructuredGrid with uniform ranges, CartesianGrid).

TestParticle.getinterp Function
julia
 getinterp(::Type{<:CartesianGrid}, A, gridx, gridy, gridz, order::Int=1, bc::Int=1)

Return a function for interpolating field array A on the grid given by gridx, gridy, and gridz.

Arguments

  • order::Int=1: order of interpolation in [1,2,3].

  • bc::Int=1: type of boundary conditions, 1 -> NaN, 2 -> periodic, 3 -> Flat.

  • dir::Int: 1/2/3, representing x/y/z direction.

Notes

The input array A may be modified in-place for memory optimization.

source
TestParticle.getinterp_scalar Function
julia
 getinterp_scalar(::Type{<:CartesianGrid}, A, gridx, gridy, gridz, order::Int=1, bc::Int=1)

Return a function for interpolating scalar array A on the grid given by gridx, gridy, and gridz. Currently only 3D arrays are supported.

Arguments

  • order::Int=1: order of interpolation in [1,2,3].

  • bc::Int=1: type of boundary conditions, 1 -> NaN, 2 -> periodic, 3 -> Flat.

  • dir::Int: 1/2/3, representing x/y/z direction.

source
TestParticle.prepare Function
julia
prepare(args...; kwargs...) -> (q2m, m, E, B, F)
prepare(E, B, F = ZeroField(); kwargs...)
prepare(grid::CartesianGrid, E, B, F = ZeroField(); kwargs...)
prepare(x, E, B, F = ZeroField(); dir = 1, kwargs...)
prepare(x, y, E, B, F = ZeroField(); kwargs...)
prepare(x, y, z, E, B, F = ZeroField(); kwargs...)
prepare(B; E = ZeroField(), F = ZeroField(), kwargs...)

Return a tuple consists of particle charge-mass ratio for a prescribed species of charge q and mass m, mass m for a prescribed species, analytic/interpolated EM field functions, and external force F.

Prescribed species are Electron and Proton; other species can be manually specified with m and q keywords or species = Ion(m̄, q̄), where m̄ and q̄ are the mass and charge numbers respectively.

Direct range input for uniform grid in 1/2/3D is supported. For 1D grid, an additional keyword dir is used for specifying the spatial direction, 1 -> x, 2 -> y, 3 -> z. For 3D grid, the default grid type is CartesianGrid. To use StructuredGrid (spherical) grid, an additional keyword gridtype is needed. For StructuredGrid (spherical) grid, dimensions of field arrays should be (Br, Bθ, Bϕ).

Keywords

  • order::Int=1: order of interpolation in [1,2,3].

  • bc::Int=1: type of boundary conditions, 1 -> NaN, 2 -> periodic.

  • species=Proton: particle species.

  • q=nothing: particle charge.

  • m=nothing: particle mass.

  • gridtype: CartesianGrid, RectilinearGrid, StructuredGrid.

source