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: 3124 samples with 78 evaluations
 min    374.795 ns (2 allocs: 64 bytes)
 median 378.000 ns (2 allocs: 64 bytes)
 mean   382.635 ns (2 allocs: 64 bytes)
 max    815.231 ns (2 allocs: 64 bytes)

julia> @be A_sph_nu($loc)
Benchmark: 3095 samples with 84 evaluations
 min    336.345 ns (2 allocs: 48 bytes)
 median 356.143 ns (2 allocs: 48 bytes)
 mean   356.738 ns (2 allocs: 48 bytes)
 max    1.011 μs (2 allocs: 48 bytes)

Uniform spherical interpolation ​

julia
julia> @be B_sph($loc)
Benchmark: 1216 samples with 214 evaluations
 min    134.500 ns (2 allocs: 64 bytes)
 median 136.093 ns (2 allocs: 64 bytes)
 mean   867.748 ns (2.00 allocs: 64.006 bytes, 0.08% gc time)
 max    888.047 μs (2.15 allocs: 71.850 bytes, 99.81% gc time)

julia> @be A_sph($loc)
Benchmark: 3210 samples with 276 evaluations
 min    102.362 ns (2 allocs: 48 bytes)
 median 104.652 ns (2 allocs: 48 bytes)
 mean   105.935 ns (2 allocs: 48 bytes)
 max    257.906 ns (2 allocs: 48 bytes)

Uniform Cartesian interpolation ​

julia
julia> @be B_car($loc)
Benchmark: 3196 samples with 392 evaluations
 min    71.867 ns (2 allocs: 64 bytes)
 median 73.913 ns (2 allocs: 64 bytes)
 mean   74.856 ns (2 allocs: 64 bytes)
 max    145.704 ns (2 allocs: 64 bytes)

julia> @be A_car($loc)
Benchmark: 3177 samples with 441 evaluations
 min    64.587 ns (2 allocs: 48 bytes)
 median 66.088 ns (2 allocs: 48 bytes)
 mean   67.153 ns (2 allocs: 48 bytes)
 max    138.102 ns (2 allocs: 48 bytes)

Non-uniform Cartesian interpolation ​

julia
julia> @be B_car_nu($loc)
Benchmark: 3216 samples with 98 evaluations
 min    297.592 ns (2 allocs: 64 bytes)
 median 299.531 ns (2 allocs: 64 bytes)
 mean   303.470 ns (2 allocs: 64 bytes)
 max    540.490 ns (2 allocs: 64 bytes)

julia> @be A_car_nu($loc)
Benchmark: 3142 samples with 97 evaluations
 min    299.938 ns (2 allocs: 48 bytes)
 median 302.320 ns (2 allocs: 48 bytes)
 mean   305.921 ns (2 allocs: 48 bytes)
 max    569.825 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