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: 3093 samples with 226 evaluations
min 129.442 ns (2 allocs: 64 bytes)
median 131.305 ns (2 allocs: 64 bytes)
mean 133.212 ns (2 allocs: 64 bytes)
max 299.186 ns (2 allocs: 64 bytes)
julia> @be A_sph_nu($loc)
Benchmark: 3134 samples with 288 evaluations
min 100.778 ns (2 allocs: 48 bytes)
median 102.448 ns (2 allocs: 48 bytes)
mean 103.813 ns (2 allocs: 48 bytes)
max 219.056 ns (2 allocs: 48 bytes)Uniform spherical interpolation
julia> @be B_sph($loc)
Benchmark: 3113 samples with 219 evaluations
min 132.160 ns (2 allocs: 64 bytes)
median 134.863 ns (2 allocs: 64 bytes)
mean 136.892 ns (2 allocs: 64 bytes)
max 372.566 ns (2 allocs: 64 bytes)
julia> @be A_sph($loc)
Benchmark: 3151 samples with 278 evaluations
min 103.428 ns (2 allocs: 48 bytes)
median 105.665 ns (2 allocs: 48 bytes)
mean 107.137 ns (2 allocs: 48 bytes)
max 205.633 ns (2 allocs: 48 bytes)Uniform Cartesian interpolation
julia> @be B_car($loc)
Benchmark: 3179 samples with 415 evaluations
min 68.730 ns (2 allocs: 64 bytes)
median 69.937 ns (2 allocs: 64 bytes)
mean 70.924 ns (2 allocs: 64 bytes)
max 170.366 ns (2 allocs: 64 bytes)
julia> @be A_car($loc)
Benchmark: 3191 samples with 439 evaluations
min 65.269 ns (2 allocs: 48 bytes)
median 66.754 ns (2 allocs: 48 bytes)
mean 67.652 ns (2 allocs: 48 bytes)
max 130.722 ns (2 allocs: 48 bytes)Non-uniform Cartesian interpolation
julia> @be B_car_nu($loc)
Benchmark: 2619 samples with 438 evaluations
min 65.237 ns (2 allocs: 64 bytes)
median 68.347 ns (2 allocs: 64 bytes)
mean 81.642 ns (2 allocs: 64 bytes)
max 609.790 ns (2 allocs: 64 bytes)
julia> @be A_car_nu($loc)
Benchmark: 3224 samples with 463 evaluations
min 60.782 ns (2 allocs: 48 bytes)
median 62.416 ns (2 allocs: 48 bytes)
mean 63.313 ns (2 allocs: 48 bytes)
max 176.933 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: 3122 samples with 141 evaluations
min 205.915 ns (2 allocs: 64 bytes)
median 209.184 ns (2 allocs: 64 bytes)
mean 211.473 ns (2 allocs: 64 bytes)
max 461.362 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.