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 Meshes
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: 3126 samples with 239 evaluations
 min    120.184 ns (2 allocs: 64 bytes)
 median 122.782 ns (2 allocs: 64 bytes)
 mean   124.816 ns (2 allocs: 64 bytes)
 max    276.251 ns (2 allocs: 64 bytes)

julia> @be A_sph_nu($loc)
Benchmark: 3157 samples with 311 evaluations
 min    92.428 ns (2 allocs: 48 bytes)
 median 94.325 ns (2 allocs: 48 bytes)
 mean   95.699 ns (2 allocs: 48 bytes)
 max    172.997 ns (2 allocs: 48 bytes)

Uniform spherical interpolation

julia> @be B_sph($loc)
Benchmark: 3128 samples with 246 evaluations
 min    118.350 ns (2 allocs: 64 bytes)
 median 119.654 ns (2 allocs: 64 bytes)
 mean   121.276 ns (2 allocs: 64 bytes)
 max    216.996 ns (2 allocs: 64 bytes)

julia> @be A_sph($loc)
Benchmark: 2819 samples with 298 evaluations
 min    94.406 ns (2 allocs: 48 bytes)
 median 106.879 ns (2 allocs: 48 bytes)
 mean   111.081 ns (2 allocs: 48 bytes)
 max    902.302 ns (2 allocs: 48 bytes)

Uniform Cartesian interpolation

julia> @be B_car($loc)
Benchmark: 3201 samples with 366 evaluations
 min    49.683 ns (2 allocs: 64 bytes)
 median 75.937 ns (2 allocs: 64 bytes)
 mean   78.832 ns (2 allocs: 64 bytes)
 max    746.678 ns (2 allocs: 64 bytes)

julia> @be A_car($loc)
Benchmark: 3402 samples with 576 evaluations
 min    45.467 ns (2 allocs: 48 bytes)
 median 46.929 ns (2 allocs: 48 bytes)
 mean   47.687 ns (2 allocs: 48 bytes)
 max    95.179 ns (2 allocs: 48 bytes)

Non-uniform Cartesian interpolation

julia> @be B_car_nu($loc)
Benchmark: 4053 samples with 445 evaluations
 min    49.578 ns (2 allocs: 64 bytes)
 median 51.333 ns (2 allocs: 64 bytes)
 mean   52.075 ns (2 allocs: 64 bytes)
 max    102.575 ns (2 allocs: 64 bytes)

julia> @be A_car_nu($loc)
Benchmark: 3659 samples with 514 evaluations
 min    47.113 ns (2 allocs: 48 bytes)
 median 48.630 ns (2 allocs: 48 bytes)
 mean   49.761 ns (2 allocs: 48 bytes)
 max    164.103 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: 2733 samples with 565 evaluations
 min    46.742 ns (2 allocs: 64 bytes)
 median 48.729 ns (2 allocs: 64 bytes)
 mean   60.830 ns (2 allocs: 64 bytes)
 max    139.821 ns (2 allocs: 64 bytes)

julia> @be itp_f32($loc_f64)
Benchmark: 447 samples with 345 evaluations
 min    72.951 ns (2 allocs: 64 bytes)
 median 78.551 ns (2 allocs: 64 bytes)
 mean   957.506 ns (2 allocs: 64 bytes, 0.22% gc time)
 max    392.252 μs (2 allocs: 64 bytes, 99.91% 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: 1504 samples with 6 evaluations
 min    3.275 μs (8 allocs: 31.750 KiB)
 median 4.082 μs (8 allocs: 31.750 KiB)
 mean   11.880 μs (8 allocs: 31.750 KiB, 0.13% gc time)
 max    6.276 ms (8 allocs: 31.750 KiB, 99.71% gc time)

julia> # Order 3: Minimal allocations (Uses a view)
       @be setup_mixed_precision_field(11, 3)
Benchmark: 1381 samples with 5 evaluations
 min    3.667 μs (14 allocs: 32.047 KiB)
 median 4.575 μs (14 allocs: 32.047 KiB)
 mean   14.217 μs (14 allocs: 32.047 KiB, 0.14% gc time)
 max    7.131 ms (14 allocs: 32.047 KiB, 99.72% 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)
16172

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, Float32, FastInterpolations._Generic}, FastInterpolations._CachedRange{Float32, Float32, FastInterpolations._Generic}, FastInterpolations._CachedRange{Float32, Float32, FastInterpolations._Generic}}, Tuple{FastInterpolations.CardinalInterp{Float64, FastInterpolations.NoBC}, FastInterpolations.CardinalInterp{Float64, FastInterpolations.NoBC}, FastInterpolations.CardinalInterp{Float64, FastInterpolations.NoBC}}, 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)
16196

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

julia> size_itp3 / size_B
1.0099775505113495

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 FastInterpolations v0.4.15

PreCompute() is now available for the global natural cubic spline (CubicInterp) in ND as of FastInterpolations.jl v0.4.15. However, it is not yet supported for the local Hermite cubic spline (CardinalInterp, i.e. order = 3 in TestParticle) in ND — constructing such an interpolator with coeffs = PreCompute() raises an ArgumentError. OnTheFly() therefore remains the only option for cubic interpolation in TestParticle.

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()); # unsupported for ND cardinal cubic in v0.4.15
       
       @be itp_fly($loc_f64)
Benchmark: 3310 samples with 80 evaluations
 min    341.763 ns (2 allocs: 64 bytes)
 median 347.913 ns (2 allocs: 64 bytes)
 mean   352.383 ns (2 allocs: 64 bytes)
 max    607.763 ns (2 allocs: 64 bytes)

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

As shown, OnTheFly() preserves memory efficiency while providing higher-order accuracy. Once PreCompute() support for ND local Hermite cubic interpolation lands, it 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: 3219 samples with 187 evaluations
 min    151.193 ns (2 allocs: 64 bytes)
 median 153.925 ns (2 allocs: 64 bytes)
 mean   155.900 ns (2 allocs: 64 bytes)
 max    270.722 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