diff --git a/benchmark/ConstitutiveModelsBenchmark/ViscousModelsBenchmarks.jl b/benchmark/ConstitutiveModelsBenchmark/ViscousModelsBenchmarks.jl
new file mode 100644
index 0000000..f897ed7
--- /dev/null
+++ b/benchmark/ConstitutiveModelsBenchmark/ViscousModelsBenchmarks.jl
@@ -0,0 +1,23 @@
+using Gridap.TensorValues
+using HyperFEM.PhysicalModels
+using BenchmarkTools
+
+
+function benchmark_viscous_model()
+  elasto = NeoHookean3D(λ=1e6, μ=1e3)
+  visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=1e3), 10.)
+  model = GeneralizedMaxwell(elasto, visco)
+  Ψ, ∂Ψu, ∂Ψuu = model(Δt = 1e-2)
+  F = TensorValue(1.:9...) * 1e-3 + I3
+  Fn = TensorValue(1.:9...) * 5e-4 + I3
+  Uvn = TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3
+  J = det(F)
+  Uvn *= J^(-1/3)
+  λvn = 1e-3
+  Avn = VectorValue(Uvn.data..., λvn)
+  SUITE["Constitutive models"]["Visco-elastic Ψ"] = @benchmarkable $Ψ($F, $Fn, $Avn)
+  SUITE["Constitutive models"]["Visco-elastic ∂Ψu"] = @benchmarkable $∂Ψu($F, $Fn, $Avn)
+  SUITE["Constitutive models"]["Visco-elastic ∂Ψuu"] = @benchmarkable $∂Ψuu($F, $Fn, $Avn)
+end
+
+benchmark_viscous_model()
diff --git a/benchmark/ConstitutiveModelsBenchmark/benchmarks.jl b/benchmark/ConstitutiveModelsBenchmark/benchmarks.jl
new file mode 100644
index 0000000..7f41120
--- /dev/null
+++ b/benchmark/ConstitutiveModelsBenchmark/benchmarks.jl
@@ -0,0 +1,4 @@
+
+SUITE["Constitutive models"] = BenchmarkGroup()
+
+include("ViscousModelsBenchmarks.jl")
diff --git a/test/TestTensorAlgebra/TensorAlgebraBenchmarks.jl b/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl
similarity index 73%
rename from test/TestTensorAlgebra/TensorAlgebraBenchmarks.jl
rename to benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl
index 749098c..66342fc 100644
--- a/test/TestTensorAlgebra/TensorAlgebraBenchmarks.jl
+++ b/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl
@@ -42,8 +42,5 @@ function _δδ_λ_2D(λ::Float64)
     λ)
 end
 
-@benchmark _δδ_μ_2D(1.0)
-@benchmark 1.0 * ((δᵢₖδⱼₗ2D + δᵢₗδⱼₖ2D))
-
-@benchmark _δδ_λ_2D(1.0)
-@benchmark 1.0 * δᵢⱼδₖₗ2D
+SUITE["Tensor algebra"]["δδ_μ_2d"] = @benchmarkable δᵢₖδⱼₗ2D + δᵢₗδⱼₖ2D
+SUITE["Tensor algebra"]["δδ_λ_2d"] = @benchmarkable 1.0 * δᵢⱼδₖₗ2D
diff --git a/benchmark/TensorAlgebraBenchmarks/benchmarks.jl b/benchmark/TensorAlgebraBenchmarks/benchmarks.jl
new file mode 100644
index 0000000..8efe105
--- /dev/null
+++ b/benchmark/TensorAlgebraBenchmarks/benchmarks.jl
@@ -0,0 +1,2 @@
+
+include("TensorAlgebraBenchmarks.jl")
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
new file mode 100644
index 0000000..679b1a1
--- /dev/null
+++ b/benchmark/benchmarks.jl
@@ -0,0 +1,8 @@
+using BenchmarkTools
+using HyperFEM
+
+const SUITE = BenchmarkGroup()
+
+include("TensorAlgebraBenchmarks/benchmarks.jl")
+
+include("ConstitutiveModelsBenchmark/benchmarks.jl")
diff --git a/test/TestConstitutiveModels/ViscousModelsBenchmarks.jl b/test/TestConstitutiveModels/ViscousModelsBenchmarks.jl
deleted file mode 100644
index 1acf758..0000000
--- a/test/TestConstitutiveModels/ViscousModelsBenchmarks.jl
+++ /dev/null
@@ -1,28 +0,0 @@
-using Gridap.TensorValues
-using HyperFEM.PhysicalModels
-using BenchmarkTools
-
-
-function benchmark_viscous_model(model)
-  Ψ, ∂Ψu, ∂Ψuu = model(Δt = 1e-2)
-  F = TensorValue(1.:9...) * 1e-3 + I3
-  Fn = TensorValue(1.:9...) * 5e-4 + I3
-  Uvn = TensorValue(1.,2.,3.,2.,4.,5.,3.,5.,6.) * 2e-4 + I3
-  J = det(F)
-  Uvn *= J^(-1/3)
-  λvn = 1e-3
-  Avn = VectorValue(Uvn.data..., λvn)
-  print("Ψ(F, Fn, Avn)    |")
-  @btime $Ψ($F, $Fn, $Avn)
-  print("∂Ψu(F, Fn, Avn)  |")
-  @btime $∂Ψu($F, $Fn, $Avn)
-  print("∂Ψuu(F, Fn, Avn) |")
-  @btime $∂Ψuu($F, $Fn, $Avn)
-end
-
-
-elasto = NeoHookean3D(λ=λ, μ=μ)
-visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ1)
-visco_elastic = GeneralizedMaxwell(elasto, visco)
-benchmark_viscous_model(visco);
-benchmark_viscous_model(visco_elastic);