diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 2dffa0582..b9a145400 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -18,7 +18,7 @@ steps:
     if: build.message !~ /\[skip tests\]/
 
 env:
-  GROUP: GPU
+  GROUP: CUDA
   JULIA_PKG_SERVER: "" # it often struggles with our large artifacts
   RETESTITEMS_NWORKERS: 4
   RETESTITEMS_NWORKER_THREADS: 2
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 2dfaa4f65..21263f5f3 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -20,6 +20,8 @@ jobs:
         group:
           - Core
           - Downstream
+          - Misc
+          - Wrappers
         version:
           - "1"
         os:
diff --git a/Project.toml b/Project.toml
index cb34c6fb5..f8e7162b8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -76,9 +76,10 @@ LazyArrays = "1.8.2"
 LeastSquaresOptim = "0.8.5"
 LineSearches = "7.2"
 LinearAlgebra = "1.10"
-LinearSolve = "2.29"
+LinearSolve = "2.30"
 MINPACK = "1.2"
 MaybeInplace = "0.1.1"
+ModelingToolkit = "9.13.0"
 NLSolvers = "0.5"
 NLsolve = "4.5"
 NaNMath = "1"
@@ -89,11 +90,10 @@ PrecompileTools = "1.2"
 Preferences = "1.4"
 Printf = "1.10"
 Random = "1.91"
-ReTestItems = "1"
+ReTestItems = "1.24"
 RecursiveArrayTools = "3.8"
 Reexport = "1.2"
 SIAMFANLEquations = "1.0.1"
-SafeTestsets = "0.1"
 SciMLBase = "2.34.0"
 SimpleNonlinearSolve = "1.8"
 SparseArrays = "1.10"
@@ -104,7 +104,7 @@ StaticArrays = "1.9"
 StaticArraysCore = "1.4"
 Sundials = "4.23.1"
 SymbolicIndexingInterface = "0.3.15"
-Symbolics = "5.13"
+Symbolics = "5.26"
 Test = "1.10"
 TimerOutputs = "0.5.23"
 Zygote = "0.6.69"
@@ -130,7 +130,6 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 SIAMFANLEquations = "084e46ad-d928-497d-ad5e-07fa361a48c4"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -140,4 +139,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "BandedMatrices", "BenchmarkTools", "CUDA", "Enzyme", "FastLevenbergMarquardt", "FixedPointAcceleration", "LeastSquaresOptim", "MINPACK", "ModelingToolkit", "NLSolvers", "NLsolve", "NaNMath", "NonlinearProblemLibrary", "OrdinaryDiffEq", "Pkg", "Random", "ReTestItems", "SIAMFANLEquations", "SafeTestsets", "SparseDiffTools", "SpeedMapping", "StableRNGs", "StaticArrays", "Sundials", "Symbolics", "Test", "Zygote"]
+test = ["Aqua", "BandedMatrices", "BenchmarkTools", "CUDA", "Enzyme", "FastLevenbergMarquardt", "FixedPointAcceleration", "LeastSquaresOptim", "MINPACK", "ModelingToolkit", "NLSolvers", "NLsolve", "NaNMath", "NonlinearProblemLibrary", "OrdinaryDiffEq", "Pkg", "Random", "ReTestItems", "SIAMFANLEquations", "SpeedMapping", "StableRNGs", "StaticArrays", "Sundials", "Symbolics", "Test", "Zygote"]
diff --git a/test/core/23_test_problems_tests.jl b/test/core/23_test_problems_tests.jl
index b0f17113a..7305af96f 100644
--- a/test/core/23_test_problems_tests.jl
+++ b/test/core/23_test_problems_tests.jl
@@ -40,7 +40,7 @@ end
 export test_on_library, problems, dicts
 end
 
-@testitem "PolyAlgorithms" setup=[RobustnessTesting] begin
+@testitem "PolyAlgorithms" setup=[RobustnessTesting] tags=[:core] begin
     alg_ops = (RobustMultiNewton(), FastShortcutNonlinearPolyalg())
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
@@ -50,7 +50,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testitem "NewtonRaphson" setup=[RobustnessTesting] begin
+@testitem "NewtonRaphson" setup=[RobustnessTesting] tags=[:core] begin
     alg_ops = (NewtonRaphson(),)
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
@@ -59,7 +59,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testitem "TrustRegion" setup=[RobustnessTesting] begin
+@testitem "TrustRegion" setup=[RobustnessTesting] tags=[:core] begin
     alg_ops = (TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Simple),
         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan),
         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei),
@@ -78,7 +78,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testitem "LevenbergMarquardt" setup=[RobustnessTesting] begin
+@testitem "LevenbergMarquardt" setup=[RobustnessTesting] tags=[:core] begin
     using LinearSolve
 
     alg_ops = (LevenbergMarquardt(), LevenbergMarquardt(; α_geodesic = 0.1),
@@ -92,7 +92,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testitem "DFSane" setup=[RobustnessTesting] begin
+@testitem "DFSane" setup=[RobustnessTesting] tags=[:core] begin
     alg_ops = (DFSane(),)
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
@@ -101,21 +101,28 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testitem "Broyden" setup=[RobustnessTesting] begin
+@testitem "Broyden" setup=[RobustnessTesting] tags=[:core] begin
     alg_ops = (Broyden(), Broyden(; init_jacobian = Val(:true_jacobian)),
         Broyden(; update_rule = Val(:bad_broyden)),
         Broyden(; init_jacobian = Val(:true_jacobian), update_rule = Val(:bad_broyden)))
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
-    broken_tests[alg_ops[1]] = [1, 5, 11, 15]
-    broken_tests[alg_ops[2]] = [1, 5, 8, 11, 18]
-    broken_tests[alg_ops[3]] = [1, 5, 9, 11]
-    broken_tests[alg_ops[4]] = [5, 6, 8, 11]
+    if Sys.isapple()
+        broken_tests[alg_ops[1]] = [1, 5, 11]
+        broken_tests[alg_ops[2]] = [1, 5, 8, 11, 18]
+        broken_tests[alg_ops[3]] = [1, 5, 6, 9, 11]
+        broken_tests[alg_ops[4]] = [5, 6, 8, 11]
+    else
+        broken_tests[alg_ops[1]] = [1, 5, 11, 15]
+        broken_tests[alg_ops[2]] = [1, 5, 8, 11, 18]
+        broken_tests[alg_ops[3]] = [1, 5, 9, 11]
+        broken_tests[alg_ops[4]] = [5, 6, 8, 11]
+    end
 
-    test_on_library(problems, dicts, alg_ops, broken_tests)
+    test_on_library(problems, dicts, alg_ops, broken_tests, Sys.isapple() ? 1e-3 : 1e-4)
 end
 
-@testitem "Klement" retries=5 setup=[RobustnessTesting] begin
+@testitem "Klement" retries=5 setup=[RobustnessTesting] tags=[:core] begin
     alg_ops = (Klement(), Klement(; init_jacobian = Val(:true_jacobian_diagonal)))
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
@@ -125,7 +132,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testitem "PseudoTransient" setup=[RobustnessTesting] begin
+@testitem "PseudoTransient" setup=[RobustnessTesting] tags=[:core] begin
     # PT relies on the root being a stable equilibrium for convergence, so it won't work on
     # most problems
     alg_ops = (PseudoTransient(),)
diff --git a/test/core/forward_ad_tests.jl b/test/core/forward_ad_tests.jl
index 46034d844..24711b9d9 100644
--- a/test/core/forward_ad_tests.jl
+++ b/test/core/forward_ad_tests.jl
@@ -62,12 +62,14 @@ __compatible(::KINSOL, ::Val{:oop_cache}) = false
 export test_f!, test_f, jacobian_f, solve_with, __compatible
 end
 
-@testitem "ForwardDiff.jl Integration" setup=[ForwardADTesting] begin
+@testitem "ForwardDiff.jl Integration" setup=[ForwardADTesting] tags=[:core] begin
     for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(),
         PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(), DFSane(),
         nothing, NLsolveJL(), CMINPACK(), KINSOL(; globalization_strategy = :LineSearch))
         us = (2.0, @SVector[1.0, 1.0], [1.0, 1.0], ones(2, 2), @SArray ones(2, 2))
 
+        alg isa CMINPACK && Sys.isapple() && continue
+
         @testset "Scalar AD" begin
             for p in 1.0:0.1:100.0, u0 in us, mode in (:iip, :oop, :iip_cache, :oop_cache)
                 __compatible(u0, alg) || continue
diff --git a/test/core/nlls_tests.jl b/test/core/nlls_tests.jl
index a1506d356..483107f69 100644
--- a/test/core/nlls_tests.jl
+++ b/test/core/nlls_tests.jl
@@ -48,7 +48,7 @@ end
 export solvers, θ_init, x, y_target, true_function, θ_true, loss_function
 end
 
-@testitem "General NLLS Solvers" setup=[CoreNLLSTesting] begin
+@testitem "General NLLS Solvers" setup=[CoreNLLSTesting] tags=[:core] begin
     prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x)
     prob_iip = NonlinearLeastSquaresProblem(
         NonlinearFunction(loss_function; resid_prototype = zero(y_target)), θ_init, x)
@@ -62,7 +62,7 @@ end
     end
 end
 
-@testitem "Custom VJP" setup=[CoreNLLSTesting] begin
+@testitem "Custom VJP" setup=[CoreNLLSTesting] tags=[:core] begin
     # This is just for testing that we can use vjp provided by the user
     function vjp(v, θ, p)
         resid = zeros(length(p))
@@ -96,7 +96,7 @@ end
     end
 end
 
-@testitem "NLLS Analytic Jacobian" begin
+@testitem "NLLS Analytic Jacobian" tags=[:core] begin
     dataIn = 1:10
     f(x, p) = x[1] * dataIn .^ 2 .+ x[2] * dataIn .+ x[3]
     dataOut = f([1, 2, 3], nothing) + 0.1 * randn(10, 1)
diff --git a/test/core/rootfind_tests.jl b/test/core/rootfind_tests.jl
index c68fc8280..f72c9cc6c 100644
--- a/test/core/rootfind_tests.jl
+++ b/test/core/rootfind_tests.jl
@@ -20,7 +20,6 @@ function newton_fails(u, p)
 end
 
 const TERMINATION_CONDITIONS = [
-    SteadyStateDiffEqTerminationMode(), SimpleNonlinearSolveTerminationMode(),
     NormTerminationMode(), RelTerminationMode(), RelNormTerminationMode(),
     AbsTerminationMode(), AbsNormTerminationMode(), RelSafeTerminationMode(),
     AbsSafeTerminationMode(), RelSafeBestTerminationMode(), AbsSafeBestTerminationMode()]
@@ -53,7 +52,7 @@ end
 
 # --- NewtonRaphson tests ---
 
-@testitem "NewtonRaphson" setup=[CoreRootfindTesting] begin
+@testitem "NewtonRaphson" setup=[CoreRootfindTesting] tags=[:core] begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (
             Static(), StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()),
         ad in (AutoFiniteDiff(), AutoZygote())
@@ -119,7 +118,7 @@ end
 
 # --- TrustRegion tests ---
 
-@testitem "TrustRegion" setup=[CoreRootfindTesting] begin
+@testitem "TrustRegion" setup=[CoreRootfindTesting] tags=[:core] timeout=3600 begin
     radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.NocedalWright,
         RadiusUpdateSchemes.NLsolve, RadiusUpdateSchemes.Hei,
         RadiusUpdateSchemes.Yuan, RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin]
@@ -237,7 +236,7 @@ end
 
 # --- LevenbergMarquardt tests ---
 
-@testitem "LevenbergMarquardt" setup=[CoreRootfindTesting] begin
+@testitem "LevenbergMarquardt" setup=[CoreRootfindTesting] tags=[:core] begin
     u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
     @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
         sol = benchmark_nlsolve_oop(quadratic_f, u0; solver = LevenbergMarquardt())
@@ -323,7 +322,7 @@ end
 
 # --- DFSane tests ---
 
-@testitem "DFSane" setup=[CoreRootfindTesting] begin
+@testitem "DFSane" setup=[CoreRootfindTesting] tags=[:core] begin
     u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
 
     @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
@@ -394,7 +393,7 @@ end
 
 # --- PseudoTransient tests ---
 
-@testitem "PseudoTransient" setup=[CoreRootfindTesting] begin
+@testitem "PseudoTransient" setup=[CoreRootfindTesting] tags=[:core] begin
     # These are tests for NewtonRaphson so we should set alpha_initial to be high so that we
     # converge quickly
     @testset "PT: alpha_initial = 10.0 PT AD: $(ad)" for ad in (
@@ -463,7 +462,7 @@ end
 
 # --- Broyden tests ---
 
-@testitem "Broyden" setup=[CoreRootfindTesting] begin
+@testitem "Broyden" setup=[CoreRootfindTesting] tags=[:core] timeout=3600 begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian) Update Rule: $(update_rule)" for lsmethod in (
             Static(), StrongWolfe(), BackTracking(),
             HagerZhang(), MoreThuente(), LiFukushimaLineSearch()),
@@ -513,7 +512,7 @@ end
 
 # --- Klement tests ---
 
-@testitem "Klement" setup=[CoreRootfindTesting] begin
+@testitem "Klement" setup=[CoreRootfindTesting] tags=[:core] begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian)" for lsmethod in (
             Static(), StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()),
         ad in (AutoFiniteDiff(), AutoZygote()),
@@ -562,7 +561,7 @@ end
 
 # --- LimitedMemoryBroyden tests ---
 
-@testitem "LimitedMemoryBroyden" setup=[CoreRootfindTesting] begin
+@testitem "LimitedMemoryBroyden" setup=[CoreRootfindTesting] tags=[:core] begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (
             Static(), StrongWolfe(), BackTracking(),
             HagerZhang(), MoreThuente(), LiFukushimaLineSearch()),
@@ -612,7 +611,7 @@ end
 end
 
 # Miscellaneous Tests
-@testitem "Custom JVP" setup=[CoreRootfindTesting] begin
+@testitem "Custom JVP" setup=[CoreRootfindTesting] tags=[:core] begin
     function F(u::Vector{Float64}, p::Vector{Float64})
         Δ = Tridiagonal(-ones(99), 2 * ones(100), -ones(99))
         return u + 0.1 * u .* Δ * u - p
diff --git a/test/downstream/mtk_cache_indexing_tests.jl b/test/downstream/mtk_cache_indexing_tests.jl
index a7770f32e..64a073612 100644
--- a/test/downstream/mtk_cache_indexing_tests.jl
+++ b/test/downstream/mtk_cache_indexing_tests.jl
@@ -1,46 +1,48 @@
-using ModelingToolkit, NonlinearSolve
-using ModelingToolkit: t_nounits as t
-
-@parameters p d
-@variables X(t)
-eqs = [0 ~ sin(X + p) - d * sqrt(X + 1)]
-@mtkbuild nlsys = NonlinearSystem(eqs, [X], [p, d])
-
-# Creates an integrator.
-nlprob = NonlinearProblem(nlsys, [X => 1.0], [p => 2.0, d => 3.0])
-
-@testset "GeneralizedFirstOrderAlgorithmCache" begin
-    nint = init(nlprob, NewtonRaphson())
-    @test nint isa NonlinearSolve.GeneralizedFirstOrderAlgorithmCache
-
-    @test nint[X] == 1.0
-    @test nint[nlsys.X] == 1.0
-    @test nint[:X] == 1.0
-    @test nint.ps[p] == 2.0
-    @test nint.ps[nlsys.p] == 2.0
-    @test nint.ps[:p] == 2.0
-end
-
-@testset "NonlinearSolvePolyAlgorithmCache" begin
-    nint = init(nlprob, FastShortcutNonlinearPolyalg())
-    @test nint isa NonlinearSolve.NonlinearSolvePolyAlgorithmCache
-
-    @test nint[X] == 1.0
-    @test nint[nlsys.X] == 1.0
-    @test nint[:X] == 1.0
-    @test nint.ps[p] == 2.0
-    @test nint.ps[nlsys.p] == 2.0
-    @test nint.ps[:p] == 2.0
-end
-
-@testset "NonlinearSolveNoInitCache" begin
-    nint = init(nlprob, SimpleNewtonRaphson())
-    @test nint isa NonlinearSolve.NonlinearSolveNoInitCache
-
-    @test nint[X] == 1.0
-    @test nint[nlsys.X] == 1.0
-    @test nint[:X] == 1.0
-    @test nint.ps[p] == 2.0
-    @test nint.ps[nlsys.p] == 2.0
-    @test nint.ps[:p] == 2.0
+@testitem "Modeling Toolkit Cache Indexing" tags=[:downstream] begin
+    using ModelingToolkit
+    using ModelingToolkit: t_nounits as t
+
+    @parameters p d
+    @variables X(t)
+    eqs = [0 ~ sin(X + p) - d * sqrt(X + 1)]
+    @mtkbuild nlsys = NonlinearSystem(eqs, [X], [p, d])
+
+    # Creates an integrator.
+    nlprob = NonlinearProblem(nlsys, [X => 1.0], [p => 2.0, d => 3.0])
+
+    @testset "GeneralizedFirstOrderAlgorithmCache" begin
+        nint = init(nlprob, NewtonRaphson())
+        @test nint isa NonlinearSolve.GeneralizedFirstOrderAlgorithmCache
+
+        @test nint[X] == 1.0
+        @test nint[nlsys.X] == 1.0
+        @test nint[:X] == 1.0
+        @test nint.ps[p] == 2.0
+        @test nint.ps[nlsys.p] == 2.0
+        @test nint.ps[:p] == 2.0
+    end
+
+    @testset "NonlinearSolvePolyAlgorithmCache" begin
+        nint = init(nlprob, FastShortcutNonlinearPolyalg())
+        @test nint isa NonlinearSolve.NonlinearSolvePolyAlgorithmCache
+
+        @test nint[X] == 1.0
+        @test nint[nlsys.X] == 1.0
+        @test nint[:X] == 1.0
+        @test nint.ps[p] == 2.0
+        @test nint.ps[nlsys.p] == 2.0
+        @test nint.ps[:p] == 2.0
+    end
+
+    @testset "NonlinearSolveNoInitCache" begin
+        nint = init(nlprob, SimpleNewtonRaphson())
+        @test nint isa NonlinearSolve.NonlinearSolveNoInitCache
+
+        @test nint[X] == 1.0
+        @test nint[nlsys.X] == 1.0
+        @test nint[:X] == 1.0
+        @test nint.ps[p] == 2.0
+        @test nint.ps[nlsys.p] == 2.0
+        @test nint.ps[:p] == 2.0
+    end
 end
diff --git a/test/gpu/core_tests.jl b/test/gpu/core_tests.jl
index 6e9806498..79a20cc97 100644
--- a/test/gpu/core_tests.jl
+++ b/test/gpu/core_tests.jl
@@ -1,4 +1,4 @@
-@testitem "CUDA Tests" begin
+@testitem "CUDA Tests" tags=[:cuda] skip=:(using CUDA; !CUDA.functional()) begin
     using CUDA, NonlinearSolve, LinearSolve, StableRNGs
 
     CUDA.allowscalar(false)
diff --git a/test/misc/aliasing_tests.jl b/test/misc/aliasing_tests.jl
index 857233029..653490dfa 100644
--- a/test/misc/aliasing_tests.jl
+++ b/test/misc/aliasing_tests.jl
@@ -1,4 +1,4 @@
-@testitem "PolyAlgorithm Aliasing" begin
+@testitem "PolyAlgorithm Aliasing" tags=[:misc] begin
     using NonlinearProblemLibrary
 
     # Use a problem that the initial solvers cannot solve and cause the initial value to
diff --git a/test/misc/banded_matrices_tests.jl b/test/misc/banded_matrices_tests.jl
index d4b8fbee5..a9a2bf4ca 100644
--- a/test/misc/banded_matrices_tests.jl
+++ b/test/misc/banded_matrices_tests.jl
@@ -1,4 +1,4 @@
-@testitem "Banded Matrix vcat" begin
+@testitem "Banded Matrix vcat" tags=[:misc] begin
     using BandedMatrices, LinearAlgebra, SparseArrays
 
     b = BandedMatrix(Ones(5, 5), (1, 1))
diff --git a/test/misc/bruss_tests.jl b/test/misc/bruss_tests.jl
index 43bcc3eb8..21b7d792b 100644
--- a/test/misc/bruss_tests.jl
+++ b/test/misc/bruss_tests.jl
@@ -1,4 +1,4 @@
-@testitem "Brusselator 2D" begin
+@testitem "Brusselator 2D" tags=[:misc] begin
     using LinearAlgebra, SparseArrays, Symbolics
 
     const N = 32
diff --git a/test/misc/linsolve_switching_tests.jl b/test/misc/linsolve_switching_tests.jl
index 521e417cc..34d10e294 100644
--- a/test/misc/linsolve_switching_tests.jl
+++ b/test/misc/linsolve_switching_tests.jl
@@ -1,4 +1,4 @@
-@testitem "Singular Systems -- Auto Linear Solve Switching" begin
+@testitem "Singular Systems -- Auto Linear Solve Switching" tags=[:misc] begin
     using LinearSolve, NonlinearSolve
 
     function f!(du, u, p)
diff --git a/test/misc/matrix_resizing_tests.jl b/test/misc/matrix_resizing_tests.jl
index 9210e3e6b..2aa8151ce 100644
--- a/test/misc/matrix_resizing_tests.jl
+++ b/test/misc/matrix_resizing_tests.jl
@@ -1,4 +1,4 @@
-@testitem "Out-of-place Matrix Resizing" begin
+@testitem "Out-of-place Matrix Resizing" tags=[:misc] begin
     using StableRNGs
 
     ff(u, p) = u .* u .- p
@@ -14,7 +14,7 @@
     end
 end
 
-@testitem "Inplace Matrix Resizing" begin
+@testitem "Inplace Matrix Resizing" tags=[:misc] begin
     using StableRNGs
 
     fiip(du, u, p) = (du .= u .* u .- p)
diff --git a/test/misc/noinit_caching_tests.jl b/test/misc/noinit_caching_tests.jl
index f9207d82f..4ab4138f8 100644
--- a/test/misc/noinit_caching_tests.jl
+++ b/test/misc/noinit_caching_tests.jl
@@ -1,4 +1,4 @@
-@testitem "NoInit Caching" begin
+@testitem "NoInit Caching" tags=[:misc] begin
     using LinearAlgebra
     import NLsolve, NLSolvers
 
diff --git a/test/misc/polyalg_tests.jl b/test/misc/polyalg_tests.jl
index c2864759c..88be7e8cc 100644
--- a/test/misc/polyalg_tests.jl
+++ b/test/misc/polyalg_tests.jl
@@ -1,4 +1,4 @@
-@testitem "Basic PolyAlgorithms" begin
+@testitem "Basic PolyAlgorithms" tags=[:misc] begin
     f(u, p) = u .* u .- 2
     u0 = [1.0, 1.0]
     probN = NonlinearProblem{false}(f, u0)
@@ -55,7 +55,7 @@
     end
 end
 
-@testitem "Testing #153 Singular Exception" begin
+@testitem "Testing #153 Singular Exception" tags=[:misc] begin
     # https://github.com/SciML/NonlinearSolve.jl/issues/153
     function f(du, u, p)
         s1, s1s2, s2 = u
@@ -71,7 +71,7 @@ end
     @test SciMLBase.successful_retcode(sol)
 end
 
-@testitem "PolyAlgorithms Autodiff" begin
+@testitem "PolyAlgorithms Autodiff" tags=[:misc] begin
     cache = zeros(2)
     function f(du, u, p)
         cache .= u .* u
@@ -104,7 +104,7 @@ end
     @test SciMLBase.successful_retcode(sol)
 end
 
-@testitem "Simple Scalar Problem #187" begin
+@testitem "Simple Scalar Problem #187" tags=[:misc] begin
     using NaNMath
 
     # https://github.com/SciML/NonlinearSolve.jl/issues/187
@@ -126,7 +126,7 @@ end
 # Shooting Problem: Taken from BoundaryValueDiffEq.jl
 # Testing for Complex Valued Root Finding. For Complex valued inputs we drop some of the
 # algorithms which dont support those.
-@testitem "Complex Valued Problems: Single-Shooting" begin
+@testitem "Complex Valued Problems: Single-Shooting" tags=[:misc] begin
     using OrdinaryDiffEq
 
     function ode_func!(du, u, p, t)
@@ -153,7 +153,7 @@ end
         prob, RobustMultiNewton(eltype(prob.u0)); abstol = 1e-19, maxiters = 10)
 end
 
-@testitem "No AD" begin
+@testitem "No AD" tags=[:misc] begin
     no_ad_fast = FastShortcutNonlinearPolyalg(autodiff = AutoFiniteDiff())
     no_ad_robust = RobustMultiNewton(autodiff = AutoFiniteDiff())
     no_ad_algs = Set([no_ad_fast, no_ad_robust, no_ad_fast.algs..., no_ad_robust.algs...])
@@ -226,7 +226,7 @@ end
 export f1_infeasible!, f1_infeasible
 end
 
-@testitem "[IIP] Infeasible" setup=[InfeasibleFunction] begin
+@testitem "[IIP] Infeasible" setup=[InfeasibleFunction] tags=[:misc] begin
     u0 = [0.0, 0.0, 0.0]
     prob = NonlinearProblem(f1_infeasible!, u0)
     sol = solve(prob)
@@ -236,7 +236,7 @@ end
     @inferred solve(prob)
 end
 
-@testitem "[OOP] Infeasible" setup=[InfeasibleFunction] begin
+@testitem "[OOP] Infeasible" setup=[InfeasibleFunction] tags=[:misc] begin
     using LinearAlgebra, StaticArrays
 
     u0 = [0.0, 0.0, 0.0]
diff --git a/test/misc/qa_tests.jl b/test/misc/qa_tests.jl
index 28cf1dbf4..79706a376 100644
--- a/test/misc/qa_tests.jl
+++ b/test/misc/qa_tests.jl
@@ -1,4 +1,4 @@
-@testitem "Aqua" begin
+@testitem "Aqua" tags=[:misc] begin
     using NonlinearSolve, SimpleNonlinearSolve, Aqua
 
     Aqua.find_persistent_tasks_deps(NonlinearSolve)
diff --git a/test/runtests.jl b/test/runtests.jl
index a3f786689..421a96f7b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,16 +1,10 @@
-using ReTestItems
+using ReTestItems, CUDA
 
-const GROUP = get(ENV, "GROUP", "All")
+const GROUP = get(ENV, "GROUP", CUDA.functional() ? "All" : "All")
 
-if GROUP == "All" || GROUP == "Core"
-    ReTestItems.runtests(joinpath(@__DIR__, "core/"), joinpath(@__DIR__, "misc/"),
-        joinpath(@__DIR__, "wrappers/"))
-end
-
-if GROUP == "Downstream"
-    include("downstream/downstream_tests.jl")
-end
-
-if GROUP == "GPU"
-    ReTestItems.runtests(joinpath(@__DIR__, "gpu/"))
+if GROUP == "All"
+    ReTestItems.runtests(@__DIR__)
+else
+    tags = [Symbol(lowercase(GROUP))]
+    ReTestItems.runtests(@__DIR__; tags)
 end
diff --git a/test/wrappers/fixedpoint_tests.jl b/test/wrappers/fixedpoint_tests.jl
index 3d5713844..a3ec497c7 100644
--- a/test/wrappers/fixedpoint_tests.jl
+++ b/test/wrappers/fixedpoint_tests.jl
@@ -5,7 +5,7 @@ import SIAMFANLEquations, FixedPointAcceleration, SpeedMapping, NLsolve
 end
 
 # Simple Scalar Problem
-@testitem "Simple Scalar Problem" setup=[WrapperFixedPointImports] begin
+@testitem "Simple Scalar Problem" setup=[WrapperFixedPointImports] tags=[:wrappers] begin
     f1(x, p) = cos(x) - x
     prob = NonlinearProblem(f1, 1.1)
 
@@ -22,7 +22,7 @@ end
 end
 
 # Simple Vector Problem
-@testitem "Simple Vector Problem" setup=[WrapperFixedPointImports] begin
+@testitem "Simple Vector Problem" setup=[WrapperFixedPointImports] tags=[:wrappers] begin
     f2(x, p) = cos.(x) .- x
     prob = NonlinearProblem(f2, [1.1, 1.1])
 
@@ -41,7 +41,7 @@ end
 
 # Fixed Point for Power Method
 # Taken from https://github.com/NicolasL-S/SpeedMapping.jl/blob/95951db8f8a4457093090e18802ad382db1c76da/test/runtests.jl
-@testitem "Power Method" setup=[WrapperFixedPointImports] begin
+@testitem "Power Method" setup=[WrapperFixedPointImports] tags=[:wrappers] begin
     C = [1 2 3; 4 5 6; 7 8 9]
     A = C + C'
     B = Hermitian(ones(10) * ones(10)' .* im + Diagonal(1:10))
diff --git a/test/wrappers/nlls_tests.jl b/test/wrappers/nlls_tests.jl
index a65ab5ec9..53cea758d 100644
--- a/test/wrappers/nlls_tests.jl
+++ b/test/wrappers/nlls_tests.jl
@@ -28,7 +28,7 @@ end
 export loss_function, θ_init, y_target, true_function, x, θ_true
 end
 
-@testitem "LeastSquaresOptim.jl" setup=[WrapperNLLSSetup] begin
+@testitem "LeastSquaresOptim.jl" setup=[WrapperNLLSSetup] tags=[:wrappers] begin
     prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x)
     prob_iip = NonlinearLeastSquaresProblem(
         NonlinearFunction(loss_function; resid_prototype = zero(y_target)), θ_init, x)
@@ -46,7 +46,7 @@ end
     end
 end
 
-@testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Provided" setup=[WrapperNLLSSetup] begin
+@testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Provided" setup=[WrapperNLLSSetup] tags=[:wrappers] begin
     function jac!(J, θ, p)
         resid = zeros(length(p))
         ForwardDiff.jacobian!(J, (resid, θ) -> loss_function(resid, θ, p), resid, θ)
@@ -70,14 +70,14 @@ end
             NonlinearFunction{false}(loss_function; jac), θ_init, x)]
 
     solvers = Any[FastLevenbergMarquardtJL(linsolve) for linsolve in (:cholesky, :qr)]
-    push!(solvers, CMINPACK())
+    Sys.isapple() || push!(solvers, CMINPACK())
     for solver in solvers, prob in probs
         sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
         @test maximum(abs, sol.resid) < 1e-6
     end
 end
 
-@testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Not Provided" setup=[WrapperNLLSSetup] begin
+@testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Not Provided" setup=[WrapperNLLSSetup] tags=[:wrappers] begin
     probs = [
         NonlinearLeastSquaresProblem(
             NonlinearFunction{true}(loss_function; resid_prototype = zero(y_target)),
@@ -90,7 +90,8 @@ end
     solvers = vec(Any[FastLevenbergMarquardtJL(linsolve; autodiff)
                       for linsolve in (:cholesky, :qr),
     autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff())])
-    append!(solvers, [CMINPACK(; method) for method in (:auto, :lm, :lmdif)])
+    Sys.isapple() ||
+        append!(solvers, [CMINPACK(; method) for method in (:auto, :lm, :lmdif)])
 
     for solver in solvers, prob in probs
         sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
@@ -98,7 +99,7 @@ end
     end
 end
 
-@testitem "FastLevenbergMarquardt.jl + StaticArrays" setup=[WrapperNLLSSetup] begin
+@testitem "FastLevenbergMarquardt.jl + StaticArrays" setup=[WrapperNLLSSetup] tags=[:wrappers] begin
     x_sa = SA[-1.0, -0.5, 0.0, 0.5, 1.0]
 
     const y_target_sa = true_function(x_sa, θ_true)
diff --git a/test/wrappers/rootfind_tests.jl b/test/wrappers/rootfind_tests.jl
index dcee9ceba..9568575e0 100644
--- a/test/wrappers/rootfind_tests.jl
+++ b/test/wrappers/rootfind_tests.jl
@@ -6,7 +6,7 @@ import NLSolvers, NLsolve, SIAMFANLEquations, MINPACK
 export NLSolvers
 end
 
-@testitem "Steady State Problems" setup=[WrapperRootfindImports] begin
+@testitem "Steady State Problems" setup=[WrapperRootfindImports] tags=[:wrappers] begin
     # IIP Tests
     function f_iip(du, u, p, t)
         du[1] = 2 - 2u[1]
@@ -18,6 +18,7 @@ end
     for alg in [
         NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())),
         NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+        alg isa CMINPACK && Sys.isapple() && continue
         sol = solve(prob_iip, alg)
         @test SciMLBase.successful_retcode(sol.retcode)
         @test maximum(abs, sol.resid) < 1e-6
@@ -31,13 +32,14 @@ end
     for alg in [
         NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())),
         NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+        alg isa CMINPACK && Sys.isapple() && continue
         sol = solve(prob_oop, alg)
         @test SciMLBase.successful_retcode(sol.retcode)
         @test maximum(abs, sol.resid) < 1e-6
     end
 end
 
-@testitem "Nonlinear Root Finding Problems" setup=[WrapperRootfindImports] begin
+@testitem "Nonlinear Root Finding Problems" setup=[WrapperRootfindImports] tags=[:wrappers] begin
     # IIP Tests
     function f_iip(du, u, p)
         du[1] = 2 - 2u[1]
@@ -49,6 +51,7 @@ end
     for alg in [
         NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())),
         NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+        alg isa CMINPACK && Sys.isapple() && continue
         local sol
         sol = solve(prob_iip, alg)
         @test SciMLBase.successful_retcode(sol.retcode)
@@ -62,6 +65,7 @@ end
     for alg in [
         NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())),
         NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+        alg isa CMINPACK && Sys.isapple() && continue
         local sol
         sol = solve(prob_oop, alg)
         @test SciMLBase.successful_retcode(sol.retcode)
@@ -78,6 +82,8 @@ end
             SIAMFANLEquationsJL(; method = :pseudotransient),
             SIAMFANLEquationsJL(; method = :secant)]
 
+        alg isa CMINPACK && Sys.isapple() && continue
+
         sol = solve(prob_tol, alg, abstol = tol)
         @test abs(sol.u[1] - sqrt(2)) < tol
     end