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: 3124 samples with 73 evaluations
min 393.877 ns (2 allocs: 64 bytes)
median 402.247 ns (2 allocs: 64 bytes)
mean 407.136 ns (2 allocs: 64 bytes)
max 799.151 ns (2 allocs: 64 bytes)
julia> @be A_sph_nu($loc)
Benchmark: 3135 samples with 81 evaluations
min 358.679 ns (2 allocs: 48 bytes)
median 361.160 ns (2 allocs: 48 bytes)
mean 366.060 ns (2 allocs: 48 bytes)
max 567.025 ns (2 allocs: 48 bytes)Uniform spherical interpolation
julia> @be B_sph($loc)
Benchmark: 3139 samples with 194 evaluations
min 148.881 ns (2 allocs: 64 bytes)
median 151.773 ns (2 allocs: 64 bytes)
mean 153.992 ns (2 allocs: 64 bytes)
max 261.371 ns (2 allocs: 64 bytes)
julia> @be A_sph($loc)
Benchmark: 3145 samples with 257 evaluations
min 111.097 ns (2 allocs: 48 bytes)
median 113.790 ns (2 allocs: 48 bytes)
mean 115.174 ns (2 allocs: 48 bytes)
max 186.856 ns (2 allocs: 48 bytes)Uniform Cartesian interpolation
julia> @be B_car($loc)
Benchmark: 3152 samples with 361 evaluations
min 78.983 ns (2 allocs: 64 bytes)
median 80.759 ns (2 allocs: 64 bytes)
mean 82.138 ns (2 allocs: 64 bytes)
max 128.753 ns (2 allocs: 64 bytes)
julia> @be A_car($loc)
Benchmark: 3131 samples with 435 evaluations
min 64.970 ns (2 allocs: 48 bytes)
median 67.503 ns (2 allocs: 48 bytes)
mean 68.302 ns (2 allocs: 48 bytes)
max 122.460 ns (2 allocs: 48 bytes)Non-uniform Cartesian interpolation
julia> @be B_car_nu($loc)
Benchmark: 3107 samples with 91 evaluations
min 321.143 ns (2 allocs: 64 bytes)
median 324.560 ns (2 allocs: 64 bytes)
mean 328.594 ns (2 allocs: 64 bytes)
max 626.209 ns (2 allocs: 64 bytes)
julia> @be A_car_nu($loc)
Benchmark: 3134 samples with 92 evaluations
min 315.913 ns (2 allocs: 48 bytes)
median 317.870 ns (2 allocs: 48 bytes)
mean 322.674 ns (2 allocs: 48 bytes)
max 628.120 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: 3131 samples with 120 evaluations
min 235.933 ns (2 allocs: 64 bytes)
median 245.708 ns (2 allocs: 64 bytes)
mean 248.730 ns (2 allocs: 64 bytes)
max 718.492 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.