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.
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
function setup_time_dependent_field()
x = range(-10, 10, length = 16)
y = range(-10, 10, length = 16)
z = range(-10, 10, length = 16)
# Create two time snapshots
B0 = zeros(3, length(x), length(y), length(z))
B0[3, :, :, :] .= 1.0 # Bz = 1 at t=0
B1 = zeros(3, length(x), length(y), length(z))
B1[3, :, :, :] .= 2.0 # Bz = 2 at t=1
times = [0.0, 1.0]
function loader(i)
if i == 1
# For demonstration, we assume we load from disk here
return TP.getinterp(TP.CartesianGrid, B0, x, y, z)
elseif i == 2
return TP.getinterp(TP.CartesianGrid, B1, x, y, z)
else
error("Index out of bounds")
end
end
# B_field_t(x, t)
B_field_t = TP.LazyTimeInterpolator(times, loader)
return B_field_t
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();
B_td = setup_time_dependent_field();
loc = SA[1.0, 1.0, 1.0];3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
1.0
1.0
1.0Gridded spherical interpolation
julia> @be B_sph_nu($loc)
Benchmark: 3134 samples with 78 evaluations
min 373.526 ns (2 allocs: 64 bytes)
median 376.474 ns (2 allocs: 64 bytes)
mean 380.388 ns (2 allocs: 64 bytes)
max 778.115 ns (2 allocs: 64 bytes)
julia> @be A_sph_nu($loc)
Benchmark: 3077 samples with 86 evaluations
min 335.628 ns (2 allocs: 48 bytes)
median 338.186 ns (2 allocs: 48 bytes)
mean 350.849 ns (2 allocs: 48 bytes)
max 2.031 μs (2 allocs: 48 bytes)Uniform spherical interpolation
julia> @be B_sph($loc)
Benchmark: 3148 samples with 215 evaluations
min 134.670 ns (2 allocs: 64 bytes)
median 136.488 ns (2 allocs: 64 bytes)
mean 137.896 ns (2 allocs: 64 bytes)
max 309.414 ns (2 allocs: 64 bytes)
julia> @be A_sph($loc)
Benchmark: 3117 samples with 262 evaluations
min 109.481 ns (2 allocs: 48 bytes)
median 111.775 ns (2 allocs: 48 bytes)
mean 113.953 ns (2 allocs: 48 bytes)
max 985.580 ns (2 allocs: 48 bytes)Uniform Cartesian interpolation
julia> @be B_car($loc)
Benchmark: 2370 samples with 305 evaluations
min 93.354 ns (2 allocs: 64 bytes)
median 100.089 ns (2 allocs: 64 bytes)
mean 328.534 ns (2 allocs: 64 bytes, 0.04% gc time)
max 534.172 μs (2 allocs: 64 bytes, 99.87% gc time)
julia> @be A_car($loc)
Benchmark: 3208 samples with 460 evaluations
min 61.330 ns (2 allocs: 48 bytes)
median 62.638 ns (2 allocs: 48 bytes)
mean 63.294 ns (2 allocs: 48 bytes)
max 96.593 ns (2 allocs: 48 bytes)Non-uniform Cartesian interpolation
julia> @be B_car_nu($loc)
Benchmark: 3122 samples with 98 evaluations
min 298.918 ns (2 allocs: 64 bytes)
median 301.378 ns (2 allocs: 64 bytes)
mean 304.486 ns (2 allocs: 64 bytes)
max 445.520 ns (2 allocs: 64 bytes)
julia> @be A_car_nu($loc)
Benchmark: 3142 samples with 100 evaluations
min 291.640 ns (2 allocs: 48 bytes)
median 293.250 ns (2 allocs: 48 bytes)
mean 296.444 ns (2 allocs: 48 bytes)
max 424.190 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).
Time-dependent field interpolation
For time-dependent fields, we can use LazyTimeInterpolator. It takes a list of time points and a loader function that returns a spatial interpolator for a given time index. The interpolator will linearly interpolate between the two nearest time points.
julia> @be B_td($loc, 0.5)
Benchmark: 3153 samples with 143 evaluations
min 202.685 ns (2 allocs: 64 bytes)
median 205.420 ns (2 allocs: 64 bytes)
mean 207.518 ns (2 allocs: 64 bytes)
max 308.406 ns (2 allocs: 64 bytes)Related API
TestParticle.getinterp Function
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.
TestParticle.getinterp_scalar Function
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.
TestParticle.prepare Function
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.