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: 3102 samples with 78 evaluations
min 376.603 ns (2 allocs: 64 bytes)
median 379.423 ns (2 allocs: 64 bytes)
mean 383.978 ns (2 allocs: 64 bytes)
max 711.077 ns (2 allocs: 64 bytes)
julia> @be A_sph_nu($loc)
Benchmark: 3140 samples with 87 evaluations
min 335.575 ns (2 allocs: 48 bytes)
median 337.874 ns (2 allocs: 48 bytes)
mean 341.672 ns (2 allocs: 48 bytes)
max 700.621 ns (2 allocs: 48 bytes)Uniform spherical interpolation
julia> @be B_sph($loc)
Benchmark: 3153 samples with 212 evaluations
min 136.429 ns (2 allocs: 64 bytes)
median 138.467 ns (2 allocs: 64 bytes)
mean 140.028 ns (2 allocs: 64 bytes)
max 256.755 ns (2 allocs: 64 bytes)
julia> @be A_sph($loc)
Benchmark: 3173 samples with 287 evaluations
min 100.362 ns (2 allocs: 48 bytes)
median 102.456 ns (2 allocs: 48 bytes)
mean 103.562 ns (2 allocs: 48 bytes)
max 203.446 ns (2 allocs: 48 bytes)Uniform Cartesian interpolation
julia> @be B_car($loc)
Benchmark: 3171 samples with 381 evaluations
min 74.837 ns (2 allocs: 64 bytes)
median 76.756 ns (2 allocs: 64 bytes)
mean 77.713 ns (2 allocs: 64 bytes)
max 139.341 ns (2 allocs: 64 bytes)
julia> @be A_car($loc)
Benchmark: 3181 samples with 468 evaluations
min 60.863 ns (2 allocs: 48 bytes)
median 62.553 ns (2 allocs: 48 bytes)
mean 63.267 ns (2 allocs: 48 bytes)
max 108.024 ns (2 allocs: 48 bytes)Non-uniform Cartesian interpolation
julia> @be B_car_nu($loc)
Benchmark: 3119 samples with 98 evaluations
min 299.847 ns (2 allocs: 64 bytes)
median 301.582 ns (2 allocs: 64 bytes)
mean 304.910 ns (2 allocs: 64 bytes)
max 530.582 ns (2 allocs: 64 bytes)
julia> @be A_car_nu($loc)
Benchmark: 3171 samples with 99 evaluations
min 291.455 ns (2 allocs: 48 bytes)
median 293.273 ns (2 allocs: 48 bytes)
mean 297.520 ns (2 allocs: 48 bytes)
max 712.646 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: 3165 samples with 143 evaluations
min 201.986 ns (2 allocs: 64 bytes)
median 204.720 ns (2 allocs: 64 bytes)
mean 206.949 ns (2 allocs: 64 bytes)
max 381.622 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.