diff --git a/Manifest.toml b/Manifest.toml
index 12672e0e..b40d3fd7 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -244,9 +244,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[MPI]]
 deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "Pkg", "Random", "Requires", "Serialization", "Sockets"]
-git-tree-sha1 = "e4549a5ced642a73bd8ba2dd1d3b1d30b3530d94"
+git-tree-sha1 = "340d8dc89e1c85a846d3f38ee294bfdd1684055a"
 uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195"
-version = "0.19.0"
+version = "0.19.1"
 
 [[MPICH_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
diff --git a/Project.toml b/Project.toml
index 96b551a0..f4981849 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.2.0"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
diff --git a/test/PLaplacianTests.jl b/test/PLaplacianTests.jl
index b53d3094..911cc50c 100644
--- a/test/PLaplacianTests.jl
+++ b/test/PLaplacianTests.jl
@@ -46,7 +46,10 @@ function main(parts,strategy)
   fill!(x,1)
   @test (norm(A*x-_A*x)+1) ≈ 1
 
-  nls = NLSolver(show_trace=i_am_main(parts), method=:newton)
+  # This leads to a dead lock since the printing of the trace seems to lead to collective operations
+  # show_trace = i_am_main(parts)
+  show_trace = true # The output in MPI will be ugly, but fixing this would require to edit NLSolvers package.
+  nls = NLSolver(show_trace=show_trace, method=:newton)
   solver = FESolver(nls)
   uh = solve(solver,op)
 
diff --git a/test/mpi/runtests.jl b/test/mpi/runtests.jl
index bf399f1e..8e150ba0 100644
--- a/test/mpi/runtests.jl
+++ b/test/mpi/runtests.jl
@@ -1,3 +1,21 @@
 module MPITests
 
+using MPI
+using Test
+
+mpidir = @__DIR__
+testdir = joinpath(mpidir,"..")
+repodir = joinpath(testdir,"..")
+
+function run_driver(procs,file)
+  mpiexec() do cmd
+    run(`$cmd -n $procs $(Base.julia_cmd()) --project=$repodir $(joinpath(mpidir,file))`)
+    @test true
+  end
+end
+
+run_driver(4,"runtests_np4.jl")
+run_driver(1,"runtests_np4.jl") # Check that the degenerated case works
+
+
 end # module
diff --git a/test/mpi/runtests_np4.jl b/test/mpi/runtests_np4.jl
new file mode 100644
index 00000000..f14ad468
--- /dev/null
+++ b/test/mpi/runtests_np4.jl
@@ -0,0 +1,52 @@
+module NP4
+# All test running on 4 procs here
+
+using PartitionedArrays
+const PArrays = PartitionedArrays
+using MPI
+
+if ! MPI.Initialized()
+  MPI.Init()
+end
+
+include("../GeometryTests.jl")
+include("../CellDataTests.jl")
+include("../FESpacesTests.jl")
+include("../MultiFieldTests.jl")
+include("../PoissonTests.jl")
+include("../PLaplacianTests.jl")
+
+if MPI.Comm_size(MPI.COMM_WORLD) == 4
+  parts = get_part_ids(mpi,(2,2))
+elseif MPI.Comm_size(MPI.COMM_WORLD) == 1
+  parts = get_part_ids(mpi,(1,1))
+else
+  error()
+end
+
+display(parts)
+
+t = PArrays.PTimer(parts,verbose=true)
+PArrays.tic!(t)
+
+GeometryTests.main(parts)
+PArrays.toc!(t,"Geometry")
+
+CellDataTests.main(parts)
+PArrays.toc!(t,"CellData")
+
+FESpacesTests.main(parts)
+PArrays.toc!(t,"FESpaces")
+
+MultiFieldTests.main(parts)
+PArrays.toc!(t,"MultiField")
+
+PoissonTests.main(parts)
+PArrays.toc!(t,"Poisson")
+
+PLaplacianTests.main(parts)
+PArrays.toc!(t,"PLaplacian")
+
+display(t)
+
+end #module