diff --git a/tests/mpi_routines.py b/tests/mpi_routines.py
index 74f6ab6..9ea5487 100644
--- a/tests/mpi_routines.py
+++ b/tests/mpi_routines.py
@@ -1,97 +1,139 @@
 # -*- coding: utf-8 -*-
 # @file
 # @section LICENSE
 #
 # Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
 # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Affero General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU Affero General Public License for more details.
 #
 # You should have received a copy of the GNU Affero General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 from __future__ import print_function
 
 import tamaas as tm
 import numpy as np
 
 from mpi4py import MPI
 
 
-def make_surface():
+def make_surface(N):
     spectrum = tm.Isopowerlaw2D()
     spectrum.q0 = 2
     spectrum.q1 = 2
     spectrum.q2 = 16
     spectrum.hurst = 0.8
 
     generator = tm.SurfaceGeneratorRandomPhase2D()
-    generator.setSizes([128, 128])
+    generator.setSizes(N)
     generator.setSpectrum(spectrum)
     generator.random_seed = 0
 
     return generator.buildSurface()
 
 
 def mpi_surface_generator():
+    N = [128, 128]
     tm.set_log_level(tm.LogLevel.debug)
 
     seq_surface = None
     comm = MPI.COMM_WORLD
     print('[{}] {}'.format(comm.rank, comm.size))
 
     with tm.sequential():
         if comm.rank == 0:
-            seq_surface = make_surface()
+            seq_surface = make_surface(N)
 
-    surface = make_surface()
+    surface = make_surface(N)
     print('[{}] {}'.format(comm.rank, surface.shape))
 
     recv = comm.gather(surface, root=0)
 
     if comm.rank == 0:
         gsurface = np.concatenate(recv)
 
         if False:
             import matplotlib.pyplot as plt
             plt.imshow(seq_surface)
             plt.colorbar()
             plt.figure()
             plt.imshow(gsurface)
             plt.colorbar()
             plt.show()
 
         assert np.all(seq_surface == gsurface)
 
 
 def mpi_model_creation():
     N = [20, 50, 50]
     S = [1., 1., 1.]
     comm = MPI.COMM_WORLD
 
     def get_discretizations(model):
         return model.getDiscretization(), model.getBoundaryDiscretization()
 
     model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S[1:], N[1:])
     n, bn = get_discretizations(model)
     n[0] = comm.allreduce(n[0])
     bn[0] = comm.allreduce(bn[0])
     assert n == N[1:] and bn == N[1:]
 
     model = tm.ModelFactory.createModel(tm.model_type.volume_2d, S, N)
     n, bn = get_discretizations(model)
     n[1] = comm.allreduce(n[1])
     bn[0] = comm.allreduce(bn[0])
     assert n == N and bn == N[1:]
 
 
+def mpi_polonski_keer():
+    N = [512, 512]
+    S = [1., 1.]
+    load = 0.1
+    comm = MPI.COMM_WORLD
+
+    seq_tractions = None
+    rms = 0
+    with tm.sequential():
+        if comm.rank == 0:
+            model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N)
+            surface = make_surface(N)
+            rms = tm.Statistics2D.computeSpectralRMSSlope(surface)
+            surface /= rms
+            solver = tm.PolonskyKeerRey(model, surface, 1e-14)
+            solver.solve(load)
+            seq_tractions = np.copy(model['traction'])
+
+    rms = comm.bcast(rms, root=0)
+
+    model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N)
+    surface = make_surface(N) / rms
+    solver = tm.PolonskyKeerRey(model, surface, 1e-14)
+    solver.solve(load)
+    tractions = model['traction']
+
+    recv = comm.gather(tractions, root=0)
+
+    if comm.rank == 0:
+        tractions = np.concatenate(recv)
+
+        if False:
+            import matplotlib.pyplot as plt
+            plt.imshow(seq_tractions - tractions)
+            plt.colorbar()
+            plt.show()
+
+        assert np.linalg.norm(seq_tractions - tractions) / seq_tractions.size \
+            < 1e-7
+
+
 if __name__ == '__main__':
     mpi_surface_generator()