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.

julia

import TestParticle as TP
using StaticArrays

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.SphericalNonUniformR(), B, r, θ, ϕ)
   A_field_nu = TP.getinterp_scalar(TP.SphericalNonUniformR(), A, r, θ, ϕ)
   B_field = TP.getinterp(TP.Spherical(), B, r_uniform, θ, ϕ)
   A_field = TP.getinterp_scalar(TP.Spherical(), 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

B_sph_nu, A_sph_nu, B_sph, A_sph = setup_spherical_field()
B_car, A_car = setup_cartesian_field();

Gridded spherical interpolation:

julia

loc = SA[1.0, 1.0, 1.0];

@time B_sph_nu(loc); # precompilation
@time B_sph_nu(loc);
@time A_sph_nu(loc); # precompilation
@time A_sph_nu(loc);

  0.045479 seconds (44.77 k allocations: 1.840 MiB, 99.95% compilation time)
  0.000011 seconds (1 allocation: 32 bytes)
  0.108227 seconds (381.23 k allocations: 17.740 MiB, 99.98% compilation time)
  0.000008 seconds (1 allocation: 16 bytes)

Uniform spherical interpolation:

julia

@time B_sph(loc); # precompilation
@time B_sph(loc);
@time A_sph(loc); # precompilation
@time A_sph(loc);

  0.097449 seconds (395.09 k allocations: 18.077 MiB, 99.98% compilation time)
  0.000009 seconds (3 allocations: 80 bytes)
  0.055758 seconds (265.50 k allocations: 12.266 MiB, 99.96% compilation time)
  0.000007 seconds (1 allocation: 16 bytes)

Uniform Cartesian interpolation:

julia

@time B_car(loc); # precompilation
@time B_car(loc);
@time A_car(loc); # precompilation
@time A_car(loc);

  0.000014 seconds (5 allocations: 224 bytes)
  0.000004 seconds (1 allocation: 32 bytes)
  0.047614 seconds (186.60 k allocations: 8.481 MiB, 99.96% compilation time)
  0.000004 seconds (1 allocation: 16 bytes)

Based on the benchmarks, for the same grid size, gridded interpolation (SphericalNonuniformR()) is 2x slower than uniform mesh interpolation (Spherical(), Cartesian()).

Field Interpolation ​

Field Interpolation