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.

using TestParticle
using StaticArrays
using Chairmarks

function setup_spherical_field(ns = 16)
   r = logrange(0.1, 10.0, length = ns)
   r_uniform = range(0.1, 10.0, length = ns)
   θ = range(0, π, length = ns)
   ϕ = range(0, 2π, length = ns)

   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 = build_interpolator(StructuredGrid, B, r, θ, ϕ)
   A_field_nu = build_interpolator(StructuredGrid, A, r, θ, ϕ)
   B_field = build_interpolator(StructuredGrid, B, r_uniform, θ, ϕ)
   A_field = build_interpolator(StructuredGrid, A, r_uniform, θ, ϕ)

   return B_field_nu, A_field_nu, B_field, A_field
end

function setup_cartesian_field(ns = 16)
   x = range(-10, 10, length = ns)
   y = range(-10, 10, length = ns)
   z = range(-10, 10, length = ns)
   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 = build_interpolator(B, x, y, z)
   A_field = 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 = build_interpolator(RectilinearGrid, B, x, y, z)
   A_field = build_interpolator(RectilinearGrid, A, x, y, z)

   return B_field, A_field
end

function setup_time_dependent_field(ns = 16)
   x = range(-10, 10, length = ns)
   y = range(-10, 10, length = ns)
   z = range(-10, 10, length = ns)

   # 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 build_interpolator(CartesianGrid, B0, x, y, z)
       elseif i == 2
           return build_interpolator(CartesianGrid, B1, x, y, z)
       else
           error("Index out of bounds")
       end
   end

   # B_field_t(x, t)
   B_field_t = LazyTimeInterpolator(times, loader)

   return B_field_t
end

function setup_mixed_precision_field(ns = 11, order = 1, bc = FillExtrap(NaN); coeffs = OnTheFly())
   x = range(0.0f0, 10.0f0, length = ns)
   y = range(0.0f0, 10.0f0, length = ns)
   z = range(0.0f0, 10.0f0, length = ns)
   B = fill(0.0f0, 3, ns, ns, ns)
   B[3, :, :, :] .= 1.0f-8

   itp = build_interpolator(B, x, y, z, order, bc; coeffs)
   return itp
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();
itp_f32 = setup_mixed_precision_field();

loc = SA[1.0, 1.0, 1.0];
loc_f32 = SA[1.0f0, 1.0f0, 1.0f0];
loc_f64 = SA[1.0, 1.0, 1.0];

3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 1.0
 1.0
 1.0

Gridded spherical interpolation

Input Location

For spherical data, the input location is still in Cartesian coordinates!

julia> @be B_sph_nu($loc)
Benchmark: 3121 samples with 225 evaluations
 min    128.729 ns (2 allocs: 64 bytes)
 median 130.689 ns (2 allocs: 64 bytes)
 mean   132.770 ns (2 allocs: 64 bytes)
 max    215.031 ns (2 allocs: 64 bytes)

julia> @be A_sph_nu($loc)
Benchmark: 3162 samples with 287 evaluations
 min    100.362 ns (2 allocs: 48 bytes)
 median 101.934 ns (2 allocs: 48 bytes)
 mean   103.322 ns (2 allocs: 48 bytes)
 max    186.484 ns (2 allocs: 48 bytes)

Uniform spherical interpolation

julia> @be B_sph($loc)
Benchmark: 3142 samples with 226 evaluations
 min    128.381 ns (2 allocs: 64 bytes)
 median 130.159 ns (2 allocs: 64 bytes)
 mean   131.854 ns (2 allocs: 64 bytes)
 max    240.854 ns (2 allocs: 64 bytes)

julia> @be A_sph($loc)
Benchmark: 3144 samples with 282 evaluations
 min    102.461 ns (2 allocs: 48 bytes)
 median 103.954 ns (2 allocs: 48 bytes)
 mean   105.206 ns (2 allocs: 48 bytes)
 max    185.918 ns (2 allocs: 48 bytes)

Uniform Cartesian interpolation

julia> @be B_car($loc)
Benchmark: 1937 samples with 481 evaluations
 min    58.800 ns (2 allocs: 64 bytes)
 median 79.505 ns (2 allocs: 64 bytes)
 mean   224.065 ns (2 allocs: 64 bytes, 0.05% gc time)
 max    289.638 μs (2 allocs: 64 bytes, 99.86% gc time)

julia> @be A_car($loc)
Benchmark: 3162 samples with 495 evaluations
 min    57.420 ns (2 allocs: 48 bytes)
 median 59.081 ns (2 allocs: 48 bytes)
 mean   59.903 ns (2 allocs: 48 bytes)
 max    107.376 ns (2 allocs: 48 bytes)

Non-uniform Cartesian interpolation

julia> @be B_car_nu($loc)
Benchmark: 3147 samples with 507 evaluations
 min    56.241 ns (2 allocs: 64 bytes)
 median 57.899 ns (2 allocs: 64 bytes)
 mean   58.811 ns (2 allocs: 64 bytes)
 max    139.335 ns (2 allocs: 64 bytes)

julia> @be A_car_nu($loc)
Benchmark: 3198 samples with 523 evaluations
 min    54.040 ns (2 allocs: 48 bytes)
 median 55.572 ns (2 allocs: 48 bytes)
 mean   56.308 ns (2 allocs: 48 bytes)
 max    93.694 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).

Mixed precision interpolation

Numerical field data from files is often stored in Float32. TestParticle supports constructing interpolators from Float32 data and ranges, which can then be queried with both Float32 and Float64 location vectors.

julia> @be itp_f32($loc_f32)
Benchmark: 3158 samples with 480 evaluations
 min    59.006 ns (2 allocs: 64 bytes)
 median 60.508 ns (2 allocs: 64 bytes)
 mean   61.438 ns (2 allocs: 64 bytes)
 max    121.810 ns (2 allocs: 64 bytes)

julia> @be itp_f32($loc_f64)
Benchmark: 2214 samples with 463 evaluations
 min    61.281 ns (2 allocs: 64 bytes)
 median 83.892 ns (2 allocs: 64 bytes)
 mean   177.166 ns (2 allocs: 64 bytes, 0.05% gc time)
 max    210.332 μs (2 allocs: 64 bytes, 99.89% gc time)

Memory usage analysis

Large numerical field arrays can consume significant amounts of memory. It's important that interpolators are memory-efficient during both construction and storage.

For linear interpolation (order=1), construction is near zero-allocation because it creates a wrapper around a reinterpreted view of your existing array.

For higher-order interpolation (order = 3), FastInterpolations.jl provides two modes: OnTheFly() and PreCompute(). By default, OnTheFly() is used, which calculates interpolation coefficients during each query. This approach is memory-efficient as it does not require additional storage beyond the input data. Alternatively, PreCompute() precalculates and stores coefficients in an additional array of the same size, enabling faster $O(1)$ lookup performance at the cost of higher memory usage.

We can measure this difference using @be.

julia> # Order 1: Minimal allocations (Uses a view)
       @be setup_mixed_precision_field(11, 1)
Benchmark: 1268 samples with 7 evaluations
 min    2.970 μs (8 allocs: 31.781 KiB)
 median 3.562 μs (8 allocs: 31.781 KiB)
 mean   10.735 μs (8 allocs: 31.781 KiB, 0.23% gc time)
 max    4.562 ms (8 allocs: 31.781 KiB, 98.60% gc time)

julia> # Order 3: Minimal allocations (Uses a view)
       @be setup_mixed_precision_field(11, 3)
Benchmark: 1528 samples with 6 evaluations
 min    3.338 μs (14 allocs: 32.109 KiB)
 median 3.971 μs (14 allocs: 32.109 KiB)
 mean   10.639 μs (14 allocs: 32.109 KiB, 0.19% gc time)
 max    4.961 ms (14 allocs: 32.109 KiB, 98.76% gc time)

Comparing the ratios relative to the original array size illustrates the overhead:

julia> B = fill(0.0f0, 3, 11, 11, 11);

julia> size_B = Base.summarysize(B) # Original 4D field array size
16036

julia> size_itp1 = Base.summarysize(itp_f32) # Total size for order=1 (Essentially the input array size)
16196

julia> itp_f32_q = setup_mixed_precision_field(11, 3) # Total size for order=3
(::TestParticle.FieldInterpolator{FastInterpolations.HeteroInterpolantND{Float32, StaticArraysCore.SVector{3, Float32}, 3, Tuple{FastInterpolations._CachedRange{Float32}, FastInterpolations._CachedRange{Float32}, FastInterpolations._CachedRange{Float32}}, Tuple{FastInterpolations.ScalarSpacing{Float32}, FastInterpolations.ScalarSpacing{Float32}, FastInterpolations.ScalarSpacing{Float32}}, Tuple{FastInterpolations.CardinalInterp{Float64}, FastInterpolations.CardinalInterp{Float64}, FastInterpolations.CardinalInterp{Float64}}, Tuple{FillExtrap{StaticArraysCore.SVector{3, Float32}}, FillExtrap{StaticArraysCore.SVector{3, Float32}}, FillExtrap{StaticArraysCore.SVector{3, Float32}}}, Tuple{FastInterpolations.AutoSearch, FastInterpolations.AutoSearch, FastInterpolations.AutoSearch}, Array{StaticArraysCore.SVector{3, Float32}, 3}}}) (generic function with 2 methods)

julia> size_itp3 = Base.summarysize(itp_f32_q)
16220

julia> # Ratios relative to raw data
       size_itp1 / size_B
1.0099775505113495

julia> size_itp3 / size_B
1.0114741830880518

As a rule of thumb, linear interpolation and cubic interpolation with OnTheFly() coefficients have nearly zero memory overhead (ratio ≈ 1.0), as they both operate directly on the input data. When supported, cubic interpolation with PreCompute() coefficients increases the memory footprint by $2^{\mathrm{DIM}}$ to store the extra coefficients, where $\mathrm{DIM}$ is the dimension of the field.

On-the-fly vs Precomputed coefficients

Cubic interpolation (order = 3) requires high-order coefficients. By default, TestParticle uses OnTheFly() coefficients, which are calculated at query time. This saves memory but increases evaluation time. For maximum performance, you can use PreCompute(), which stores the coefficients in an additional array.

Status in v0.4.8

Support for PreCompute() with local Hermite cubic interpolation is currently under development and not yet available in FastInterpolations.jl v0.4.8.

julia> # Benchmark evaluation time
       itp_fly = setup_mixed_precision_field(11, 3; coeffs = OnTheFly());

julia> # itp_pre = setup_mixed_precision_field(11, 3; coeffs = PreCompute()); # Not yet supported in v0.4.8
       
       @be itp_fly($loc_f64)
Benchmark: 3177 samples with 43 evaluations
 min    664.512 ns (2 allocs: 64 bytes)
 median 674.512 ns (2 allocs: 64 bytes)
 mean   682.713 ns (2 allocs: 64 bytes)
 max    1.434 μs (2 allocs: 64 bytes)

julia> # Compare total object size
       Base.summarysize(itp_fly)
16220

As shown, OnTheFly() preserves memory efficiency while providing higher-order accuracy. Once available, PreCompute() will offer a faster alternative for memory-abundant systems.

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: 3151 samples with 163 evaluations
 min    176.773 ns (2 allocs: 64 bytes)
 median 179.902 ns (2 allocs: 64 bytes)
 mean   182.034 ns (2 allocs: 64 bytes)
 max    304.681 ns (2 allocs: 64 bytes)

Related API

TestParticle.build_interpolator Function

build_interpolator(gridtype, A, grids..., order::Int=1, bc=FillExtrap(NaN))
build_interpolator(A, grids..., order::Int=1, bc=FillExtrap(NaN))

Return a function for interpolating field array A on the given grids.

Arguments

• gridtype: CartesianGrid, RectilinearGrid or StructuredGrid. 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 [0,1,3].

• bc=FillExtrap(NaN): boundary condition type from FastInterpolations.jl.

  • FillExtrap(NaN): Fill with NaN (default).

  • ClampExtrap(): Clamp (flat extrapolation).

  • WrapExtrap(): Exclusive periodic wrapping ($L=N\mathrm{\Delta }x$).

• coeffs=OnTheFly(): coefficient strategy for cubic interpolation (order=3). Default is OnTheFly().

Notes

• The input array A may be modified in-place for memory optimization.

source

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 [0,1,3].

• bc=FillExtrap(NaN): boundary condition type from FastInterpolations.jl.

• species=Proton: particle species.

• q=nothing: particle charge.

• m=nothing: particle mass.

• gridtype: CartesianGrid, RectilinearGrid, StructuredGrid.

source