diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 1119b5c..c5a3a2b 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -12,7 +12,7 @@ permissions:
   contents: read
 
 jobs:
-  test:
+  Ariadne:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
     continue-on-error: ${{ matrix.version == 'nightly' }}
@@ -42,3 +42,27 @@ jobs:
       - uses: codecov/codecov-action@v1
         with:
           file: lcov.info
+  Theseus:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    continue-on-error: ${{ matrix.version == 'nightly' }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - "1.10"
+          - "1.11"
+          # - "nightly"
+        os:
+          - ubuntu-latest
+          - macOS-latest
+          - windows-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v2
diff --git a/libs/Theseus/test/Project.toml b/libs/Theseus/test/Project.toml
new file mode 100644
index 0000000..e9e8f23
--- /dev/null
+++ b/libs/Theseus/test/Project.toml
@@ -0,0 +1,10 @@
+[deps]
+Ariadne = "0be81120-40bf-4f8b-adf0-26103efb66f1"
+ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Theseus = "7f538e44-2768-4ef2-af90-2110c0378286"
+TrixiTest = "0a316866-cbd0-4425-8bcb-08103b2c1f26"
+
+[sources]
+Ariadne = {path = "../../.."}
+Theseus = {path = ".."}
\ No newline at end of file
diff --git a/libs/Theseus/test/runtests.jl b/libs/Theseus/test/runtests.jl
new file mode 100644
index 0000000..17523e1
--- /dev/null
+++ b/libs/Theseus/test/runtests.jl
@@ -0,0 +1,4 @@
+using Theseus
+using ParallelTestRunner
+
+runtests(Theseus, ARGS)
diff --git a/test/Project.toml b/test/Project.toml
index 46dcc3a..0fef7d5 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,8 +1,9 @@
 [deps]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Ariadne = "0be81120-40bf-4f8b-adf0-26103efb66f1"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Ariadne = "0be81120-40bf-4f8b-adf0-26103efb66f1"
+ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[sources.Ariadne]
-path = ".."
\ No newline at end of file
+[sources]
+Ariadne = {path = ".."}
diff --git a/test/runtests.jl b/test/runtests.jl
index 122be50..0abdd17 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,70 +1,4 @@
-using Test
 using Ariadne
+using ParallelTestRunner
 
-function F!(res, x, _)
-    res[1] = x[1]^2 + x[2]^2 - 2
-    return res[2] = exp(x[1] - 1) + x[2]^2 - 2
-end
-
-function F(x, p)
-    res = similar(x)
-    F!(res, x, p)
-    return res
-end
-
-let x₀ = [2.0, 0.5]
-    x, stats = newton_krylov!(F!, x₀)
-    @test stats.solved
-end
-
-let x₀ = [3.0, 5.0]
-    x, stats = newton_krylov(F, x₀)
-    @test stats.solved
-end
-
-import Ariadne: JacobianOperator, BatchedJacobianOperator
-using Enzyme, LinearAlgebra
-
-@testset "Jacobian" begin
-    J_Enz = jacobian(Forward, x -> F(x, nothing), [3.0, 5.0]) |> only
-    J = JacobianOperator(F!, zeros(2), [3.0, 5.0], nothing)
-
-    @test size(J) == (2, 2)
-    @test length(J) == 4
-    @test eltype(J) == Float64
-
-    out = [NaN, NaN]
-    mul!(out, J, [1.0, 0.0])
-    @test out == [6.0, 7.38905609893065]
-
-    out = [NaN, NaN]
-    mul!(out, transpose(J), [1.0, 0.0])
-    @test out == [6.0, 10.0]
-
-    J_NK = collect(J)
-
-    @test J_NK == J_Enz
-
-    v = rand(2)
-    out = similar(v)
-    mul!(out, J, v)
-
-    @test out ≈ J_Enz * v
-
-    @test collect(transpose(J)) == transpose(collect(J))
-
-    # Batched
-    if VERSION >= v"1.11.0"
-        J = BatchedJacobianOperator{2}(F!, zeros(2), [3.0, 5.0], nothing)
-
-        V = [1.0 0.0; 0.0 1.0]
-        Out = similar(V)
-        mul!(Out, J, V)
-
-        @test Out == J_Enz
-
-        mul!(Out, transpose(J), V)
-        @test Out == J_Enz'
-        # @test Out == collect(transpose(J))
-    end
-end
+runtests(Ariadne, ARGS)
diff --git a/test/tests.jl b/test/tests.jl
new file mode 100644
index 0000000..122be50
--- /dev/null
+++ b/test/tests.jl
@@ -0,0 +1,70 @@
+using Test
+using Ariadne
+
+function F!(res, x, _)
+    res[1] = x[1]^2 + x[2]^2 - 2
+    return res[2] = exp(x[1] - 1) + x[2]^2 - 2
+end
+
+function F(x, p)
+    res = similar(x)
+    F!(res, x, p)
+    return res
+end
+
+let x₀ = [2.0, 0.5]
+    x, stats = newton_krylov!(F!, x₀)
+    @test stats.solved
+end
+
+let x₀ = [3.0, 5.0]
+    x, stats = newton_krylov(F, x₀)
+    @test stats.solved
+end
+
+import Ariadne: JacobianOperator, BatchedJacobianOperator
+using Enzyme, LinearAlgebra
+
+@testset "Jacobian" begin
+    J_Enz = jacobian(Forward, x -> F(x, nothing), [3.0, 5.0]) |> only
+    J = JacobianOperator(F!, zeros(2), [3.0, 5.0], nothing)
+
+    @test size(J) == (2, 2)
+    @test length(J) == 4
+    @test eltype(J) == Float64
+
+    out = [NaN, NaN]
+    mul!(out, J, [1.0, 0.0])
+    @test out == [6.0, 7.38905609893065]
+
+    out = [NaN, NaN]
+    mul!(out, transpose(J), [1.0, 0.0])
+    @test out == [6.0, 10.0]
+
+    J_NK = collect(J)
+
+    @test J_NK == J_Enz
+
+    v = rand(2)
+    out = similar(v)
+    mul!(out, J, v)
+
+    @test out ≈ J_Enz * v
+
+    @test collect(transpose(J)) == transpose(collect(J))
+
+    # Batched
+    if VERSION >= v"1.11.0"
+        J = BatchedJacobianOperator{2}(F!, zeros(2), [3.0, 5.0], nothing)
+
+        V = [1.0 0.0; 0.0 1.0]
+        Out = similar(V)
+        mul!(Out, J, V)
+
+        @test Out == J_Enz
+
+        mul!(Out, transpose(J), V)
+        @test Out == J_Enz'
+        # @test Out == collect(transpose(J))
+    end
+end