diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml
index 39e9d56..23c0134 100644
--- a/.github/workflows/Tests.yml
+++ b/.github/workflows/Tests.yml
@@ -5,14 +5,14 @@ on:
     branches:
       - main
     paths-ignore:
-      - 'docs/**'
+      - "docs/**"
   push:
     branches:
       - main
     paths-ignore:
-      - 'docs/**'
+      - "docs/**"
   schedule:
-    - cron: '59 12 * * 2'
+    - cron: "59 12 * * 2"
 
 concurrency:
   group: ${{ github.workflow }}-${{ github.ref }}
@@ -26,6 +26,7 @@ jobs:
       matrix:
         group:
           - Core
+          - Extensions
           - QA
         version:
           - "1"
@@ -40,6 +41,7 @@ jobs:
             version: "pre"
     uses: "SciML/.github/.github/workflows/tests.yml@v1"
     with:
+      group: "${{ matrix.group }}"
       julia-version: "${{ matrix.version }}"
       os: "${{ matrix.os }}"
     secrets: "inherit"
diff --git a/Project.toml b/Project.toml
index 9090b47..eab00fe 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,13 +22,14 @@ DataInterpolations = "8"
 EllipsisNotation = "1.8.0"
 ForwardDiff = "0"
 KernelAbstractions = "0.9.34"
+Pkg = "1"
 PrecompileTools = "1.0"
 Random = "1"
 RecipesBase = "1.3.4"
-SafeTestsets = "0.1"
 SymbolicUtils = "4"
 Symbolics = "7"
 Test = "1"
+TestItemRunner = "1"
 julia = "1"
 
 [extras]
@@ -36,10 +37,10 @@ DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["DataInterpolations", "ForwardDiff", "Pkg", "Random", "SafeTestsets", "Symbolics", "Test", "SymbolicUtils"]
+test = ["DataInterpolations", "ForwardDiff", "Pkg", "Random", "Symbolics", "Test", "TestItemRunner", "SymbolicUtils"]
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index b0c2db9..b425ac8 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -51,24 +51,18 @@ function validate_size_u(
     return @assert expected_size == size(u)[1:N_in] "Expected the size of the first N_in dimensions of u to be $expected_size based on the BSplineInterpolation properties."
 end
 
-function validate_cache(
-        ::EmptyCache, ::NTuple{N_in, ID}, ::AbstractArray
-    ) where {N_in, ID}
-    return nothing
-end
-
-function validate_cache(
-        nurbs_weights::NURBSWeights,
-        ::NTuple{N_in, BSplineInterpolationDimension},
-        u::AbstractArray
-    ) where {N_in}
+function validate_cache(nurbs_weights::NURBSWeights, interp_dims, u::AbstractArray)
+    N_in = length(interp_dims)
     size_expected = size(u)[1:N_in]
     return @assert size(nurbs_weights.weights) == size_expected "The size of the weights array must match the length of the first N_in dimensions of u ($size_expected)."
 end
 
+validate_cache(cache::AbstractInterpolationCache, interp_dims, u::AbstractArray) = validate_cache(cache, first(interp_dims))
+validate_cache(::EmptyCache, ::AbstractInterpolationDimension) = nothing
+
 function validate_cache(
-        ::gType, ::NTuple{N_in, ID}, ::AbstractArray
-    ) where {gType, N_in, ID}
+        ::gType, ::ID
+    ) where {gType <: AbstractInterpolationCache, ID <: AbstractInterpolationDimension}
     return @error("Interpolation dimension type $ID is not compatible with global cache type $gType.")
 end
 
diff --git a/test/qa/aqua.jl b/test/qa/aqua.jl
deleted file mode 100644
index 8728c99..0000000
--- a/test/qa/aqua.jl
+++ /dev/null
@@ -1,6 +0,0 @@
-using Aqua
-using DataInterpolationsND
-
-@testset "Aqua" begin
-    Aqua.test_all(DataInterpolationsND)
-end
diff --git a/test/qa/runtests.jl b/test/qa/runtests.jl
index a0846f9..d95fee9 100644
--- a/test/qa/runtests.jl
+++ b/test/qa/runtests.jl
@@ -1,5 +1,5 @@
-using SafeTestsets
+using TestItemRunner
 
-@safetestset "Aqua" include("aqua.jl")
-@safetestset "ExplicitImports" include("explicit_imports.jl")
-@safetestset "JET" include("jet.jl")
+include("test_qa_aqua.jl")
+include("test_qa_explicit_imports.jl")
+include("test_qa_jet.jl")
diff --git a/test/qa/test_qa_aqua.jl b/test/qa/test_qa_aqua.jl
new file mode 100644
index 0000000..fe01f71
--- /dev/null
+++ b/test/qa/test_qa_aqua.jl
@@ -0,0 +1,5 @@
+@testitem "Aqua" begin
+    using Aqua
+    using DataInterpolationsND
+    Aqua.test_all(DataInterpolationsND)
+end
diff --git a/test/qa/explicit_imports.jl b/test/qa/test_qa_explicit_imports.jl
similarity index 61%
rename from test/qa/explicit_imports.jl
rename to test/qa/test_qa_explicit_imports.jl
index f5de23f..55fbec9 100644
--- a/test/qa/explicit_imports.jl
+++ b/test/qa/test_qa_explicit_imports.jl
@@ -1,8 +1,6 @@
-using ExplicitImports
-using DataInterpolationsND
-using Test
-
-@testset "ExplicitImports" begin
+@testitem "ExplicitImports" begin
+    using ExplicitImports
+    using DataInterpolationsND
     @test check_no_implicit_imports(DataInterpolationsND) === nothing
     @test check_no_stale_explicit_imports(DataInterpolationsND) === nothing
 end
diff --git a/test/qa/jet.jl b/test/qa/test_qa_jet.jl
similarity index 96%
rename from test/qa/jet.jl
rename to test/qa/test_qa_jet.jl
index 3e11630..5a8fd16 100644
--- a/test/qa/jet.jl
+++ b/test/qa/test_qa_jet.jl
@@ -1,8 +1,6 @@
-using JET
-using DataInterpolationsND
-using Test
-
-@testset "JET static analysis" begin
+@testitem "JET static analysis" begin
+    using JET
+    using DataInterpolationsND
     # Linear interpolation setup
     t1 = [-3.14, 1.0, 3.0, 7.6, 12.8]
     t2 = [-2.71, 1.41, 12.76, 50.2, 120.0]
@@ -118,7 +116,9 @@ using Test
     end
 end
 
-@testset "JET type optimization" begin
+@testitem "JET type optimization" begin
+    using JET
+    using DataInterpolationsND
     # Linear interpolation setup
     t1 = [-3.14, 1.0, 3.0, 7.6, 12.8]
     t2 = [-2.71, 1.41, 12.76, 50.2, 120.0]
diff --git a/test/runtests.jl b/test/runtests.jl
index d58239b..9ceca8b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-using SafeTestsets, Pkg
+using TestItemRunner, Pkg
 const GROUP = get(ENV, "GROUP", "All")
 
 function activate_gpu_env()
@@ -14,15 +14,15 @@ function activate_qa_env()
 end
 
 if GROUP == "All" || GROUP == "Core"
-    @safetestset "Interpolations" include("test_interpolations.jl")
-    @safetestset "Derivatives" include("test_derivatives.jl")
-    @safetestset "DataInterpolations" include("test_datainterpolations_comparison.jl")
-    @safetestset "Interface" include("test_interface.jl")
+    core_files = ("test_interpolations.jl", "test_derivatives.jl", "test_datainterpolations_comparison.jl", "test_interface.jl")
+    @run_package_tests filter = ti -> any(endswith(ti.filename, file) for file in core_files)
+
 elseif GROUP == "Extensions"
-    @safetestset "Symbolics Extension" include("test_symbolics_ext.jl")
+    extension_files = ("test_symbolics_ext.jl",)
+    @run_package_tests filter = ti -> any(endswith(ti.filename, file) for file in extension_files)
 elseif GROUP == "QA"
     activate_qa_env()
-    include("qa/runtests.jl")
+    @run_package_tests filter = ti -> contains(ti.filename, r"test_qa_")
 elseif GROUP == "GPU"
     activate_gpu_env()
     # TODO: Add GPU tests
diff --git a/test/test_datainterpolations_comparison.jl b/test/test_datainterpolations_comparison.jl
index 5897736..ba9bda3 100644
--- a/test/test_datainterpolations_comparison.jl
+++ b/test/test_datainterpolations_comparison.jl
@@ -1,75 +1,76 @@
-using DataInterpolationsND
-using DataInterpolations
+@testitem "DataInterpolations comparison" begin
+    using DataInterpolations
 
-isapprox_or_nan(v1, v2) = (v1 ≈ v2) || (all(isnan, v1) && all(isnan, v2))
+    isapprox_or_nan(v1, v2) = (v1 ≈ v2) || (all(isnan, v1) && all(isnan, v2))
 
-function interpolation_comparison(
-        itp_dim_type::Type{<:DataInterpolationsND.AbstractInterpolationDimension},
-        itp_type::Type{<:DataInterpolations.AbstractInterpolation}
-    )
-    t = [0.0, 0.2, 0.5, 0.7, 0.9, 1.1]
-    u = [1.5, -5.3, 0.5, 3.3, 3.6, 1.0]
+    function interpolation_comparison(
+            itp_dim_type::Type{<:DataInterpolationsND.AbstractInterpolationDimension},
+            itp_type::Type{<:DataInterpolations.AbstractInterpolation}
+        )
+        t = [0.0, 0.2, 0.5, 0.7, 0.9, 1.1]
+        u = [1.5, -5.3, 0.5, 3.3, 3.6, 1.0]
 
-    n = length(t)
+        n = length(t)
 
-    itp_DI = itp_type(u, t)
-    itp_NDI = NDInterpolation(u, itp_dim_type(t))
+        itp_DI = itp_type(u, t)
+        itp_NDI = NDInterpolation(u, itp_dim_type(t))
 
-    # Interior non data points
-    for i in 1:(n - 1)
-        t_eval = range(t[i], t[i + 1]; length = 25)[2:24]
+        # Interior non data points
+        for i in 1:(n - 1)
+            t_eval = range(t[i], t[i + 1]; length = 25)[2:24]
+            @test itp_DI(t_eval) ≈ itp_NDI.(t_eval)
+        end
+
+        # Interior data points
+        t_eval = view(t, 2:(n - 1))
         @test itp_DI(t_eval) ≈ itp_NDI.(t_eval)
+
+        # Boundary points
+        t_eval = [first(t), last(t)]
+        return @test itp_DI(t_eval) ≈ itp_NDI.(t_eval)
     end
 
-    # Interior data points
-    t_eval = view(t, 2:(n - 1))
-    @test itp_DI(t_eval) ≈ itp_NDI.(t_eval)
+    function interpolation_derivative_comparison(
+            itp_dim_type::Type{<:DataInterpolationsND.AbstractInterpolationDimension},
+            itp_type::Type{<:DataInterpolations.AbstractInterpolation}
+        )
+        t = [0.0, 0.2, 0.5, 0.7, 0.9, 1.1]
+        u = [1.5, -5.3, 0.5, 3.3, 3.6, 1.0]
 
-    # Boundary points
-    t_eval = [first(t), last(t)]
-    return @test itp_DI(t_eval) ≈ itp_NDI.(t_eval)
-end
+        n = length(t)
 
-function interpolation_derivative_comparison(
-        itp_dim_type::Type{<:DataInterpolationsND.AbstractInterpolationDimension},
-        itp_type::Type{<:DataInterpolations.AbstractInterpolation}
-    )
-    t = [0.0, 0.2, 0.5, 0.7, 0.9, 1.1]
-    u = [1.5, -5.3, 0.5, 3.3, 3.6, 1.0]
+        itp_DI = itp_type(u, t)
+        itp_NDI = NDInterpolation(u, itp_dim_type(t))
 
-    n = length(t)
+        # Interior non data points
+        for i in 1:(n - 1)
+            t_eval = range(t[i], t[i + 1]; length = 25)[2:24]
+            @test DataInterpolations.derivative.(Ref(itp_DI), t_eval) ≈
+                itp_NDI.(t_eval; derivative_orders = (1,))
+        end
 
-    itp_DI = itp_type(u, t)
-    itp_NDI = NDInterpolation(u, itp_dim_type(t))
+        # Interior data points
+        t_eval = view(t, 2:(n - 1))
+        @test isapprox_or_nan(
+            DataInterpolations.derivative.(Ref(itp_DI), t_eval),
+            itp_NDI.(t_eval; derivative_orders = (1,))
+        )
 
-    # Interior non data points
-    for i in 1:(n - 1)
-        t_eval = range(t[i], t[i + 1]; length = 25)[2:24]
-        @test DataInterpolations.derivative.(Ref(itp_DI), t_eval) ≈
+        # Boundary points
+        t_eval = [first(t), last(t)]
+        return @test isapprox_or_nan(
+            DataInterpolations.derivative.(Ref(itp_DI), t_eval),
             itp_NDI.(t_eval; derivative_orders = (1,))
+        )
     end
 
-    # Interior data points
-    t_eval = view(t, 2:(n - 1))
-    @test isapprox_or_nan(
-        DataInterpolations.derivative.(Ref(itp_DI), t_eval),
-        itp_NDI.(t_eval; derivative_orders = (1,))
-    )
-
-    # Boundary points
-    t_eval = [first(t), last(t)]
-    return @test isapprox_or_nan(
-        DataInterpolations.derivative.(Ref(itp_DI), t_eval),
-        itp_NDI.(t_eval; derivative_orders = (1,))
-    )
-end
-
-for (itp_dim_type, itp_type) in (
-        (LinearInterpolationDimension, LinearInterpolation),
-        (ConstantInterpolationDimension, ConstantInterpolation),
-    )
-    @testset "$itp_type" begin
-        interpolation_comparison(itp_dim_type, itp_type)
-        interpolation_derivative_comparison(itp_dim_type, itp_type)
+    for (itp_dim_type, itp_type) in (
+            (LinearInterpolationDimension, LinearInterpolation),
+            (ConstantInterpolationDimension, ConstantInterpolation),
+        )
+        @testset "$itp_type" begin
+            interpolation_comparison(itp_dim_type, itp_type)
+            interpolation_derivative_comparison(itp_dim_type, itp_type)
+        end
     end
 end
diff --git a/test/test_derivatives.jl b/test/test_derivatives.jl
index 667ba16..e12a4c1 100644
--- a/test/test_derivatives.jl
+++ b/test/test_derivatives.jl
@@ -1,79 +1,80 @@
-using DataInterpolationsND
-using ForwardDiff
-using Random
+@testitem "derivatives" begin
+    using Random
+    using ForwardDiff
 
-function test_deriv_10(itp, t)
-    return @test itp(t; derivative_orders = (1, 0)) ≈
-        ForwardDiff.derivative(t₁ -> itp(t₁, t[2]), t[1])
-end
+    function test_deriv_10(itp, t)
+        return @test itp(t; derivative_orders = (1, 0)) ≈
+            ForwardDiff.derivative(t₁ -> itp(t₁, t[2]), t[1])
+    end
 
-function test_deriv_01(itp, t)
-    return @test itp(t; derivative_orders = (0, 1)) ≈
-        ForwardDiff.derivative(t₂ -> itp(t[1], t₂), t[2])
-end
+    function test_deriv_01(itp, t)
+        return @test itp(t; derivative_orders = (0, 1)) ≈
+            ForwardDiff.derivative(t₂ -> itp(t[1], t₂), t[2])
+    end
 
-function test_deriv_11(itp, t)
-    return @test itp(t; derivative_orders = (1, 1)) ≈
-        ForwardDiff.derivative(
-        t₂ -> ForwardDiff.derivative(
-            t₁ -> itp(t₁, t₂), t[1]
-        ),
-        t[2]
-    )
-end
+    function test_deriv_11(itp, t)
+        return @test itp(t; derivative_orders = (1, 1)) ≈
+            ForwardDiff.derivative(
+            t₂ -> ForwardDiff.derivative(
+                t₁ -> itp(t₁, t₂), t[1]
+            ),
+            t[2]
+        )
+    end
 
-function test_derivatives(itp::NDInterpolation{2})
-    t1 = itp.interp_dims[1].t
-    t2 = itp.interp_dims[2].t
+    function test_derivatives(itp::NDInterpolation{2})
+        t1 = itp.interp_dims[1].t
+        t2 = itp.interp_dims[2].t
 
-    # Derivatives in data points
-    for t in Iterators.product(t1, t2)
-        test_deriv_10(itp, t)
-        test_deriv_01(itp, t)
-        test_deriv_11(itp, t)
+        # Derivatives in data points
+        for t in Iterators.product(t1, t2)
+            test_deriv_10(itp, t)
+            test_deriv_01(itp, t)
+            test_deriv_11(itp, t)
+        end
+
+        # Derivatives between data points
+        for t in Iterators.product(
+                t1[1:(end - 1)] + diff(t1) / 2, t2[1:(end - 1)] + diff(t2) / 2
+            )
+            test_deriv_10(itp, t)
+            test_deriv_01(itp, t)
+            test_deriv_11(itp, t)
+        end
+        return
     end
 
-    # Derivatives between data points
-    for t in Iterators.product(
-            t1[1:(end - 1)] + diff(t1) / 2, t2[1:(end - 1)] + diff(t2) / 2
+    @testset "Linear Interpolation" begin
+        Random.seed!(1)
+        t1 = [1.0, 2.0, 3.5, 4.0, 10.0]
+        t2 = [-1.5, 0.0, 2.5, 7.8, 12.9]
+        u = rand(5, 5)
+        itp_dims = (
+            LinearInterpolationDimension(t1),
+            LinearInterpolationDimension(t2),
         )
-        test_deriv_10(itp, t)
-        test_deriv_01(itp, t)
-        test_deriv_11(itp, t)
+        itp = NDInterpolation(u, itp_dims)
+        test_derivatives(itp)
     end
-    return
-end
-
-@testset "Linear Interpolation" begin
-    Random.seed!(1)
-    t1 = [1.0, 2.0, 3.5, 4.0, 10.0]
-    t2 = [-1.5, 0.0, 2.5, 7.8, 12.9]
-    u = rand(5, 5)
-    itp_dims = (
-        LinearInterpolationDimension(t1),
-        LinearInterpolationDimension(t2),
-    )
-    itp = NDInterpolation(u, itp_dims)
-    test_derivatives(itp)
-end
 
-@testset "BSpline interpolation" begin
-    Random.seed!(1)
-    t1 = [1.0, 2.0, 3.5, 4.0, 10.0]
-    t2 = [-1.5, 0.0, 2.5, 7.8, 12.9]
-    u = rand(6, 6)
-    itp_dims = (
-        BSplineInterpolationDimension(t1, 2),
-        BSplineInterpolationDimension(t2, 2),
-    )
-    itp = NDInterpolation(u, itp_dims)
-    test_derivatives(itp)
+    @testset "BSpline interpolation" begin
+        Random.seed!(1)
+        t1 = [1.0, 2.0, 3.5, 4.0, 10.0]
+        t2 = [-1.5, 0.0, 2.5, 7.8, 12.9]
+        u = rand(6, 6)
+        itp_dims = (
+            BSplineInterpolationDimension(t1, 2),
+            BSplineInterpolationDimension(t2, 2),
+        )
+        itp = NDInterpolation(u, itp_dims)
+        test_derivatives(itp)
 
-    u = rand(9, 9)
-    itp_dims = (
-        BSplineInterpolationDimension(t1, 5),
-        BSplineInterpolationDimension(t2, 5),
-    )
-    itp = NDInterpolation(u, itp_dims)
-    test_derivatives(itp)
+        u = rand(9, 9)
+        itp_dims = (
+            BSplineInterpolationDimension(t1, 5),
+            BSplineInterpolationDimension(t2, 5),
+        )
+        itp = NDInterpolation(u, itp_dims)
+        test_derivatives(itp)
+    end
 end
diff --git a/test/test_interface.jl b/test/test_interface.jl
index 3f8426d..55dd3e2 100644
--- a/test/test_interface.jl
+++ b/test/test_interface.jl
@@ -4,7 +4,7 @@ using Test
 # Interface tests to verify the package properly adheres to Julia's standard interfaces
 # and SciML's array/number interfaces. Uses BigFloat to test numeric type genericity.
 
-@testset "BigFloat Support" begin
+@testitem "BigFloat Support" begin
     @testset "LinearInterpolationDimension" begin
         t = BigFloat[1.0, 2.0, 3.0, 4.0]
         t_eval = BigFloat[1.5, 2.5]
@@ -141,7 +141,7 @@ using Test
     end
 end
 
-@testset "Float32 Support" begin
+@testitem "Float32 Support" begin
     @testset "NDInterpolation with Linear" begin
         t1 = Float32[1.0, 2.0, 3.0, 4.0]
         t2 = Float32[1.0, 2.0, 3.0, 4.0]
@@ -180,7 +180,7 @@ end
     end
 end
 
-@testset "Type Promotion" begin
+@testitem "Type Promotion" begin
     @testset "Mixed precision t and u" begin
         t1 = Float32[1.0, 2.0, 3.0, 4.0]
         t2 = Float32[1.0, 2.0, 3.0, 4.0]
diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl
index 1a3075e..66e1aaf 100644
--- a/test/test_interpolations.jl
+++ b/test/test_interpolations.jl
@@ -1,70 +1,6 @@
-using DataInterpolationsND: AbstractInterpolationDimension, EmptyCache
-using DataInterpolationsND
-using Random
-
-function test_globally_constant(
-        ID::Type{<:AbstractInterpolationDimension}; args1 = (), args2 = (), kwargs1 = (),
-        kwargs2 = (), cache = EmptyCache(), test_derivatives = true
-    )
-    t1 = [-3.14, 1.0, 3.0, 7.6, 12.8]
-    t2 = [-2.71, 1.41, 12.76, 50.2, 120.0]
-
-    u = if ID == BSplineInterpolationDimension
-        fill(2.0, 4 + args1[1], 4 + args2[1])
-    else
-        fill(2.0, 5, 5)
-    end
-
-    # Evaluation in data points
-    itp_dims = (
-        ID(t1, args1...; t_eval = t1, kwargs1...),
-        ID(t2, args2...; t_eval = t2, kwargs2...),
-    )
-
-    itp = NDInterpolation(u, itp_dims; cache)
-    @test all(x -> isapprox(x, 2.0; atol = 1.0e-10), eval_grid(itp))
-    if test_derivatives
-        @test all(
-            x -> isapprox(x, 0.0; atol = 1.0e-10), eval_grid(itp, derivative_orders = (1, 0))
-        )
-        @test all(
-            x -> isapprox(x, 0.0; atol = 1.0e-10), eval_grid(itp, derivative_orders = (0, 1))
-        )
-    end
-
-    # Evaluation between data points
-    itp_dims = (
-        ID(t1, args1...; t_eval = t1[1:(end - 1)] + diff(t1) / 2, kwargs1...),
-        ID(t2, args2...; t_eval = t2[1:(end - 1)] + diff(t2) / 2, kwargs2...),
-    )
-    itp = NDInterpolation(u, itp_dims)
-    @test all(x -> isapprox(x, 2.0; atol = 1.0e-10), eval_grid(itp))
-    return if test_derivatives
-        @test all(
-            x -> isapprox(x, 0.0; atol = 1.0e-10), eval_grid(itp, derivative_orders = (1, 0))
-        )
-        @test all(
-            x -> isapprox(x, 0.0; atol = 1.0e-10), eval_grid(itp, derivative_orders = (0, 1))
-        )
-    end
-end
-
-function test_analytic(itp::NDInterpolation{N_in}, f) where {N_in}
-    # Evaluation in data points
-    ts = ntuple(dim_in -> itp.interp_dims[dim_in].t, N_in)
-    for t in Iterators.product(ts...)
-        @test itp(t) ≈ f(t...)
-    end
-
-    # Evaluation between data points
-    ts_ = ntuple(dim_in -> ts[dim_in][1:(end - 1)] + diff(ts[dim_in]) / 2, N_in)
-    for t in Iterators.product(ts_...)
-        @test itp(t) ≈ f(t...)
-    end
-    return
-end
-
-@testset "Linear Interpolation" begin
+@testitem "Linear Interpolation" begin
+    using Random
+    include("utils.jl")
     test_globally_constant(LinearInterpolationDimension)
 
     f(t1, t2) = 3.0 + 2.3t1 - 4.7t2
@@ -82,7 +18,8 @@ end
     test_analytic(itp, f)
 end
 
-@testset "BSpline Interpolation" begin
+@testitem "BSpline Interpolation" begin
+    include("utils.jl")
     test_globally_constant(
         BSplineInterpolationDimension, args1 = (2,), args2 = (3,),
         kwargs1 = (:max_derivative_order_eval => 1,),
@@ -102,8 +39,8 @@ end
     test_analytic(itp, f)
 end
 
-@testset "NURBS Interpolation" begin
-    Random.seed!(10)
+@testitem "NURBS Interpolation" begin
+    include("utils.jl")
 
     test_globally_constant(
         BSplineInterpolationDimension; args1 = (3,), args2 = (1,),
diff --git a/test/test_symbolics_ext.jl b/test/test_symbolics_ext.jl
index 1b63629..35ada07 100644
--- a/test/test_symbolics_ext.jl
+++ b/test/test_symbolics_ext.jl
@@ -1,23 +1,21 @@
-using DataInterpolationsND
-using Symbolics
-import SymbolicUtils as SU
-using Symbolics: unwrap
-using Test
+@testitem "Basics" begin
+    using Symbolics
+    using Symbolics: unwrap
+    import SymbolicUtils as SU
 
-t1 = cumsum(rand(5))
-t2 = cumsum(rand(7))
+    t1 = cumsum(rand(5))
+    t2 = cumsum(rand(7))
 
-interpolation_dimensions = (
-    LinearInterpolationDimension(t1),
-    LinearInterpolationDimension(t2),
-)
+    interpolation_dimensions = (
+        LinearInterpolationDimension(t1),
+        LinearInterpolationDimension(t2),
+    )
 
-u = rand(5, 7, 2)
+    u = rand(5, 7, 2)
+    interp = NDInterpolation(u, interpolation_dimensions)
 
-interp = NDInterpolation(u, interpolation_dimensions)
-@variables x y
+    @variables x y
 
-@testset "Basics" begin
     ex = interp(x, y)
     @test ex isa Symbolics.Arr
     @test size(ex) == (2,)
@@ -36,7 +34,23 @@ interp = NDInterpolation(u, interpolation_dimensions)
     @test ex isa SU.BasicSymbolic{Vector{Real}}
 end
 
-@testset "Differentiation" begin
+@testitem "Differentiation" begin
+    using Symbolics
+    import SymbolicUtils as SU
+
+    t1 = cumsum(rand(5))
+    t2 = cumsum(rand(7))
+
+    interpolation_dimensions = (
+        LinearInterpolationDimension(t1),
+        LinearInterpolationDimension(t2),
+    )
+
+    u = rand(5, 7, 2)
+    interp = NDInterpolation(u, interpolation_dimensions)
+
+    @variables x y
+
     ex = interp(x, y)
     der = Symbolics.derivative(ex[1], x)
     @test size(der) == ()
diff --git a/test/utils.jl b/test/utils.jl
new file mode 100644
index 0000000..37a18a3
--- /dev/null
+++ b/test/utils.jl
@@ -0,0 +1,67 @@
+using DataInterpolationsND: AbstractInterpolationDimension, EmptyCache
+
+###
+#### Interpolations
+###
+
+function test_globally_constant(
+        ID::Type{<:AbstractInterpolationDimension}; args1 = (), args2 = (), kwargs1 = (),
+        kwargs2 = (), cache = EmptyCache(), test_derivatives = true
+    )
+    t1 = [-3.14, 1.0, 3.0, 7.6, 12.8]
+    t2 = [-2.71, 1.41, 12.76, 50.2, 120.0]
+
+    u = if ID == BSplineInterpolationDimension
+        fill(2.0, 4 + args1[1], 4 + args2[1])
+    else
+        fill(2.0, 5, 5)
+    end
+
+    # Evaluation in data points
+    itp_dims = (
+        ID(t1, args1...; t_eval = t1, kwargs1...),
+        ID(t2, args2...; t_eval = t2, kwargs2...),
+    )
+
+    itp = NDInterpolation(u, itp_dims; cache)
+    @test all(x -> isapprox(x, 2.0; atol = 1.0e-10), eval_grid(itp))
+    if test_derivatives
+        @test all(
+            x -> isapprox(x, 0.0; atol = 1.0e-10), eval_grid(itp, derivative_orders = (1, 0))
+        )
+        @test all(
+            x -> isapprox(x, 0.0; atol = 1.0e-10), eval_grid(itp, derivative_orders = (0, 1))
+        )
+    end
+
+    # Evaluation between data points
+    itp_dims = (
+        ID(t1, args1...; t_eval = t1[1:(end - 1)] + diff(t1) / 2, kwargs1...),
+        ID(t2, args2...; t_eval = t2[1:(end - 1)] + diff(t2) / 2, kwargs2...),
+    )
+    itp = NDInterpolation(u, itp_dims)
+    @test all(x -> isapprox(x, 2.0; atol = 1.0e-10), eval_grid(itp))
+    return if test_derivatives
+        @test all(
+            x -> isapprox(x, 0.0; atol = 1.0e-10), eval_grid(itp, derivative_orders = (1, 0))
+        )
+        @test all(
+            x -> isapprox(x, 0.0; atol = 1.0e-10), eval_grid(itp, derivative_orders = (0, 1))
+        )
+    end
+end
+
+function test_analytic(itp::NDInterpolation{N_in}, f) where {N_in}
+    # Evaluation in data points
+    ts = ntuple(dim_in -> itp.interp_dims[dim_in].t, N_in)
+    for t in Iterators.product(ts...)
+        @test itp(t) ≈ f(t...)
+    end
+
+    # Evaluation between data points
+    ts_ = ntuple(dim_in -> ts[dim_in][1:(end - 1)] + diff(ts[dim_in]) / 2, N_in)
+    for t in Iterators.product(ts_...)
+        @test itp(t) ≈ f(t...)
+    end
+    return
+end