Fast interpolations

Project.toml

@@ -11,8 +11,8 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
@@ -38,8 +38,9 @@ Adapt = "4.4"
 ChunkSplitters = "3"
 DiffResults = "1"
 Distributed = "1"
+FastInterpolations = "0.4"
 ForwardDiff = "1"
-Interpolations = "0.15, 0.16"
+Interpolations = "0.16"
 KernelAbstractions = "0.9"
 LinearAlgebra = "1"
 Meshes = "0.55, 0.56"

src/TestParticle.jl

@@ -1,9 +1,9 @@
 module TestParticle
 
 using LinearAlgebra: norm, ×, ⋅, diag, normalize
-using Interpolations: interpolate, interpolate!, extrapolate, scale, BSpline, Linear,
-    Quadratic, Cubic,
-    Line, OnCell, Periodic, Flat, Gridded
+using FastInterpolations: linear_interp, quadratic_interp, cubic_interp, constant_interp,
+    Extrap, NoExtrap, PeriodicBC, ZeroCurvBC
+import Interpolations
 using SciMLBase: AbstractODEProblem, AbstractODEFunction, AbstractODESolution, ReturnCode,
     BasicEnsembleAlgorithm, EnsembleThreads, EnsembleSerial, EnsembleDistributed,
     DEFAULT_SPECIALIZATION, ODEFunction, ODEProblem, remake,
@@ -41,11 +41,14 @@ export get_gyrofrequency,
 export orbit, monitor
 export get_fields, get_work
 export LazyTimeInterpolator
+export AbstractInterpolationBackend, FastInterpolationsBackend, InterpolationsBackend
 export TraceProblem, TraceGCProblem, TraceHybridProblem, CartesianGrid, RectilinearGrid, StructuredGrid
 export EnsembleSerial, EnsembleThreads, EnsembleDistributed, remake
 
 include("types.jl")
 include("utility/utility.jl")
+include("utility/interp_backends.jl")
+include("utility/fastinterpolation.jl")
 include("utility/interpolation.jl")
 include("sampler.jl")
 include("prepare.jl")

src/prepare.jl

@@ -64,8 +64,12 @@ get_EField(param) = param[3]
 prepare_field(f, args...; kwargs...) = Field(f)
 prepare_field(f::ZeroField, args...; kwargs...) = f
 
-function prepare_field(f::AbstractArray, x...; gridtype, order, bc, kw...)
-    return Field(build_interpolator(gridtype, f, x..., order, bc; kw...))
+function prepare_field(
+        f::AbstractArray, x...;
+        gridtype, order, bc,
+        backend::AbstractInterpolationBackend = FastInterpolationsBackend(), kw...
+    )
+    return Field(build_interpolator(backend, gridtype, f, x..., order, bc; kw...))
 end
 
 function _prepare(