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.build_interpolator(TP.StructuredGrid, B, r, θ, ϕ)
A_field_nu = TP.build_interpolator(TP.StructuredGrid, A, r, θ, ϕ)
B_field = TP.build_interpolator(TP.StructuredGrid, B, r_uniform, θ, ϕ)
A_field = TP.build_interpolator(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.build_interpolator(B, x, y, z)
A_field = TP.build_interpolator(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.build_interpolator(TP.RectilinearGrid, B, x, y, z)
A_field = TP.build_interpolator(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.build_interpolator(TP.CartesianGrid, B0, x, y, z)
elseif i == 2
return TP.build_interpolator(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
Input Location
For spherical data, the input location is still in Cartesian coordinates!
julia> @be B_sph_nu($loc)
Benchmark: 3026 samples with 224 evaluations
min 130.424 ns (2 allocs: 64 bytes)
median 132.254 ns (2 allocs: 64 bytes)
mean 138.165 ns (2 allocs: 64 bytes)
max 295.152 ns (2 allocs: 64 bytes)
julia> @be A_sph_nu($loc)
Benchmark: 3182 samples with 286 evaluations
min 101.273 ns (2 allocs: 48 bytes)
median 103.024 ns (2 allocs: 48 bytes)
mean 104.342 ns (2 allocs: 48 bytes)
max 204.122 ns (2 allocs: 48 bytes)Uniform spherical interpolation
julia> @be B_sph($loc)
Benchmark: 3100 samples with 216 evaluations
min 134.231 ns (2 allocs: 64 bytes)
median 137.106 ns (2 allocs: 64 bytes)
mean 139.810 ns (2 allocs: 64 bytes)
max 283.958 ns (2 allocs: 64 bytes)
julia> @be A_sph($loc)
Benchmark: 3162 samples with 273 evaluations
min 105.103 ns (2 allocs: 48 bytes)
median 106.978 ns (2 allocs: 48 bytes)
mean 108.531 ns (2 allocs: 48 bytes)
max 186.176 ns (2 allocs: 48 bytes)Uniform Cartesian interpolation
julia> @be B_car($loc)
Benchmark: 3296 samples with 392 evaluations
min 69.212 ns (2 allocs: 64 bytes)
median 70.490 ns (2 allocs: 64 bytes)
mean 71.964 ns (2 allocs: 64 bytes)
max 177.296 ns (2 allocs: 64 bytes)
julia> @be A_car($loc)
Benchmark: 3171 samples with 436 evaluations
min 65.720 ns (2 allocs: 48 bytes)
median 67.165 ns (2 allocs: 48 bytes)
mean 68.103 ns (2 allocs: 48 bytes)
max 219.908 ns (2 allocs: 48 bytes)Non-uniform Cartesian interpolation
julia> @be B_car_nu($loc)
Benchmark: 3178 samples with 448 evaluations
min 64.004 ns (2 allocs: 64 bytes)
median 65.143 ns (2 allocs: 64 bytes)
mean 66.086 ns (2 allocs: 64 bytes)
max 148.717 ns (2 allocs: 64 bytes)
julia> @be A_car_nu($loc)
Benchmark: 2065 samples with 345 evaluations
min 62.670 ns (2 allocs: 48 bytes)
median 83.957 ns (2 allocs: 48 bytes)
mean 206.563 ns (2 allocs: 48 bytes, 0.05% gc time)
max 246.359 μs (2 allocs: 48 bytes, 99.87% gc time)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: 3091 samples with 129 evaluations
min 216.302 ns (2 allocs: 64 bytes)
median 227.636 ns (2 allocs: 64 bytes)
mean 233.946 ns (2 allocs: 64 bytes)
max 798.078 ns (2 allocs: 64 bytes)Related API
TestParticle.build_interpolator Function
build_interpolator(gridtype, A, grids..., order::Int=1, bc::Int=1)
build_interpolator(A, grids..., order::Int=1, bc::Int=1)Return a function for interpolating field array A on the given grids.
Arguments
gridtype:CartesianGrid,RectilinearGridorStructuredGrid. Usually determined by the number of grids.A: field array. For vector field, the first dimension should be 3 if it's not an SVector wrapper.order::Int=1: order of interpolation in [1,2,3].bc::Int=1: type of boundary conditions, 1 -> NaN, 2 -> periodic, 3 -> Clamp (flat extrapolation).
Notes
The input array A may be modified in-place for memory optimization.
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. The grid vectors must be sorted. 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.