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
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.0

Gridded spherical interpolation

julia

julia> @be B_sph_nu($loc)
Benchmark: 3096 samples with 79 evaluations
 min    372.595 ns (2 allocs: 64 bytes)
 median 376.392 ns (2 allocs: 64 bytes)
 mean   380.401 ns (2 allocs: 64 bytes)
 max    974.873 ns (2 allocs: 64 bytes)

julia> @be A_sph_nu($loc)
Benchmark: 3124 samples with 87 evaluations
 min    337.184 ns (2 allocs: 48 bytes)
 median 338.920 ns (2 allocs: 48 bytes)
 mean   342.528 ns (2 allocs: 48 bytes)
 max    531.575 ns (2 allocs: 48 bytes)

Uniform spherical interpolation

julia

julia> @be B_sph($loc)
Benchmark: 3141 samples with 213 evaluations
 min    135.559 ns (2 allocs: 64 bytes)
 median 137.681 ns (2 allocs: 64 bytes)
 mean   139.101 ns (2 allocs: 64 bytes)
 max    232.268 ns (2 allocs: 64 bytes)

julia> @be A_sph($loc)
Benchmark: 3195 samples with 270 evaluations
 min    103.081 ns (2 allocs: 48 bytes)
 median 105.644 ns (2 allocs: 48 bytes)
 mean   108.544 ns (2 allocs: 48 bytes)
 max    490.956 ns (2 allocs: 48 bytes)

Uniform Cartesian interpolation

julia

julia> @be B_car($loc)
Benchmark: 3007 samples with 399 evaluations
 min    71.987 ns (2 allocs: 64 bytes)
 median 73.925 ns (2 allocs: 64 bytes)
 mean   78.240 ns (2 allocs: 64 bytes)
 max    265.283 ns (2 allocs: 64 bytes)

julia> @be A_car($loc)
Benchmark: 3247 samples with 402 evaluations
 min    62.930 ns (2 allocs: 48 bytes)
 median 64.423 ns (2 allocs: 48 bytes)
 mean   71.241 ns (2 allocs: 48 bytes)
 max    451.716 ns (2 allocs: 48 bytes)

Non-uniform Cartesian interpolation

julia

julia> @be B_car_nu($loc)
Benchmark: 3071 samples with 90 evaluations
 min    317.033 ns (2 allocs: 64 bytes)
 median 323.167 ns (2 allocs: 64 bytes)
 mean   333.109 ns (2 allocs: 64 bytes)
 max    1.979 μs (2 allocs: 64 bytes)

julia> @be A_car_nu($loc)
Benchmark: 3136 samples with 100 evaluations
 min    292.750 ns (2 allocs: 48 bytes)
 median 294.350 ns (2 allocs: 48 bytes)
 mean   297.190 ns (2 allocs: 48 bytes)
 max    474.690 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

julia> @be B_td($loc, 0.5)
Benchmark: 1 sample with 3 evaluations
        437.667 ns (2 allocs: 64 bytes)

TestParticle.getinterp Function

julia

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.

source

TestParticle.getinterp_scalar Function

julia

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.

source

TestParticle.prepare Function

julia

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.

source

Field Interpolation ​

Gridded spherical interpolation ​

Uniform spherical interpolation ​

Uniform Cartesian interpolation ​

Non-uniform Cartesian interpolation ​

Time-dependent field interpolation ​

Related API ​

Field Interpolation

Gridded spherical interpolation

Uniform spherical interpolation

Uniform Cartesian interpolation

Non-uniform Cartesian interpolation

Time-dependent field interpolation

Related API