diff --git a/hysop/backend/device/codegen/symbolic/expr.py b/hysop/backend/device/codegen/symbolic/expr.py
index ce7385f2f9d88ca68a08e5f975e22d998154ccf1..d6dca3eb0b79eaffca0f1467c92046e931809f85 100644
--- a/hysop/backend/device/codegen/symbolic/expr.py
+++ b/hysop/backend/device/codegen/symbolic/expr.py
@@ -9,7 +9,7 @@ from hysop.tools.types import check_instance, first_not_None, to_tuple, to_list
 from hysop.tools.numerics import is_fp, is_signed, is_unsigned, is_integer, is_complex
 
 from packaging import version
-if version.parse(sp.__version__) > version.parse("1.7"):
+if version.parse(sm.__version__) > version.parse("1.7"):
     from sympy.printing.c import C99CodePrinter
 else:
     from sympy.printing.ccode import C99CodePrinter
diff --git a/hysop/backend/device/opencl/opencl_printer.py b/hysop/backend/device/opencl/opencl_printer.py
index 7624591e5b350355509e85e053de77d1a987661c..96d95101af75dd11744c36ec38c02003d5bf78d6 100644
--- a/hysop/backend/device/opencl/opencl_printer.py
+++ b/hysop/backend/device/opencl/opencl_printer.py
@@ -2,7 +2,7 @@ from hysop.backend.device.opencl.opencl_types import OpenClTypeGen
 from hysop.tools.types import check_instance
 import sympy as sm
 from packaging import version
-if version.parse(sp.__version__) > version.parse("1.7"):
+if version.parse(sm.__version__) > version.parse("1.7"):
     from sympy.printing.c import C99CodePrinter
 else:
     from sympy.printing.ccode import C99CodePrinter
diff --git a/hysop/fields/tests/test_cartesian.py b/hysop/fields/tests/test_cartesian.py
index 6eb9e4fb033b5aa16550ef0e2df52a96cff0039f..f09eba05d4bd0cb62c8d2777e9b49582ac279d10 100644
--- a/hysop/fields/tests/test_cartesian.py
+++ b/hysop/fields/tests/test_cartesian.py
@@ -1,10 +1,13 @@
-import os, subprocess, sys, time
+import os
+import subprocess
+import sys
+import time
 import itertools as it
 import numpy as np
 
 from hysop import __ENABLE_LONG_TESTS__
 from hysop.constants import Backend, ExchangeMethod, GhostOperation, \
-                        GhostMask, DirectionLabels, BoundaryCondition
+    GhostMask, DirectionLabels, BoundaryCondition
 from hysop.tools.parameters import CartesianDiscretization
 from hysop.tools.numerics import is_integer, is_fp
 from hysop.tools.numpywrappers import npw
@@ -14,6 +17,7 @@ from hysop.topology.cartesian_topology import CartesianTopology, CartesianTopolo
 from hysop.testsenv import iter_clenv, test_context, domain_boundary_iterator
 from hysop.tools.numerics import is_fp, is_integer
 
+
 def __random_init(data, coords, component):
     shape = data.shape
     dtype = data.dtype
@@ -25,14 +29,17 @@ def __random_init(data, coords, component):
         msg = 'Unknown dtype {}.'.format(dtype)
         raise NotImplementedError(msg)
 
+
 def __zero_init(data, coords, component):
     data[...] = 0
 
+
 def __cst_init(cst):
-    def __init(data,coords,component):
+    def __init(data, coords, component):
         data[...] = cst
     return __init
 
+
 def test_serial_initialization_1d():
     print()
     print('test_serial_initialization_1d()')
@@ -42,20 +49,20 @@ def test_serial_initialization_1d():
 
     for (lbd, rbd) in domain_boundary_iterator(dim):
         domain = Box(dim=dim, lboundaries=lbd,
-                              rboundaries=rbd)
+                     rboundaries=rbd)
         F0 = Field(domain=domain, name='F0', nb_components=1)
         F1 = Field(domain=domain, name='F1', nb_components=2)
-        F2 = Field(domain=domain, name='F2', shape=(2,2))
+        F2 = Field(domain=domain, name='F2', shape=(2, 2))
         print('[{}]'.format(F0.format_boundaries()))
 
         discretization = CartesianDiscretization(npts, nghosts,
-                                            lboundaries=F0.lboundaries,
-                                            rboundaries=F0.rboundaries)
+                                                 lboundaries=F0.lboundaries,
+                                                 rboundaries=F0.rboundaries)
 
         topo0 = CartesianTopology(domain=domain, discretization=discretization,
-                        backend=Backend.HOST)
+                                  backend=Backend.HOST)
         topos = (topo0,) + tuple(CartesianTopology(domain=domain, discretization=discretization,
-            backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
+                                                   backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
 
         assert all(t.mesh.global_lboundaries == t.mesh.local_lboundaries == F0.lboundaries for t in topos)
         assert all(t.mesh.global_rboundaries == t.mesh.local_rboundaries == F0.rboundaries for t in topos)
@@ -83,24 +90,24 @@ def test_serial_initialization_1d():
             Lx, = tuple(discretization.lboundaries)
             Rx, = tuple(discretization.rboundaries)
             try:
-                for i,d in enumerate(data):
+                for i, d in enumerate(data):
                     b = d.get().handle
-                    assert b.shape==(Nx+2*Gx,)
-                    if (Lx==BoundaryCondition.PERIODIC):
+                    assert b.shape == (Nx+2*Gx,)
+                    if (Lx == BoundaryCondition.PERIODIC):
                         assert (b[Gx:2*Gx] == b[Gx+Nx:]).all(), b
-                    elif (Lx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                    elif (Lx == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
                         assert (b[Gx] == 0).all(), b
                         assert (b[:Gx] == -b[Gx+1:2*Gx+1][::-1]).all(), b
-                    elif (Lx==BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                    elif (Lx == BoundaryCondition.HOMOGENEOUS_NEUMANN):
                         assert (b[:Gx] == +b[Gx+1:2*Gx+1][::-1]).all(), b
                     else:
                         raise NotImplementedError('Unknown boundary condition {}.'.format(Lx))
-                    if (Rx==BoundaryCondition.PERIODIC):
+                    if (Rx == BoundaryCondition.PERIODIC):
                         assert (b[:Gx] == b[Nx:Gx+Nx]).all(), b
-                    elif (Rx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                    elif (Rx == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
                         assert (b[Nx+Gx-1] == 0).all(), b
                         assert (b[Nx-1:Nx+Gx-1] == -b[Nx+Gx:][::-1]).all(), b
-                    elif (Rx==BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                    elif (Rx == BoundaryCondition.HOMOGENEOUS_NEUMANN):
                         assert (b[Nx-1:Nx+Gx-1] == +b[Nx+Gx:][::-1]).all(), b
                     else:
                         raise NotImplementedError('Unknown boundary condition {}.'.format(Rx))
@@ -114,23 +121,23 @@ def test_serial_initialization_2d():
     print()
     print('test_serial_initialization_2d()')
     dim = 2
-    npts = (4,8)
-    nghosts = (1,2)
+    npts = (4, 8)
+    nghosts = (1, 2)
     for (lbd, rbd) in domain_boundary_iterator(dim):
         domain = Box(dim=dim, lboundaries=lbd,
-                              rboundaries=rbd)
+                     rboundaries=rbd)
         F0 = Field(domain=domain, name='F0', nb_components=1)
         F1 = Field(domain=domain, name='F1', nb_components=2)
         print('[{}]'.format(F0.format_boundaries()))
 
         discretization = CartesianDiscretization(npts, nghosts,
-                                            lboundaries=F0.lboundaries,
-                                            rboundaries=F0.rboundaries)
+                                                 lboundaries=F0.lboundaries,
+                                                 rboundaries=F0.rboundaries)
 
         topo0 = CartesianTopology(domain=domain, discretization=discretization,
-                        backend=Backend.HOST)
+                                  backend=Backend.HOST)
         topos = (topo0,) + tuple(CartesianTopology(domain=domain, discretization=discretization,
-            backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
+                                                   backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
 
         assert all(np.all(t.mesh.local_lboundaries == F0.lboundaries) for t in topos)
         assert all(np.all(t.mesh.local_rboundaries == F0.rboundaries) for t in topos)
@@ -146,24 +153,24 @@ def test_serial_initialization_2d():
             assert all(np.all(d.local_lboundaries == F0.lboundaries) for d in dfields)
             assert all(np.all(d.local_rboundaries == F0.rboundaries) for d in dfields)
 
-            Ny,Nx = npts
-            Gy,Gx = nghosts
-            Ly,Lx = tuple(discretization.lboundaries)
-            Ry,Rx = tuple(discretization.rboundaries)
-            Xo = (slice(0,Gx),
-                   slice(Gx,Gx+Nx),
-                   slice(Gx+Nx,None))
-            Yo = (slice(0,Gy),
-                   slice(Gy,Gy+Ny),
-                   slice(Gy+Ny,None))
+            Ny, Nx = npts
+            Gy, Gx = nghosts
+            Ly, Lx = tuple(discretization.lboundaries)
+            Ry, Rx = tuple(discretization.rboundaries)
+            Xo = (slice(0, Gx),
+                  slice(Gx, Gx+Nx),
+                  slice(Gx+Nx, None))
+            Yo = (slice(0, Gy),
+                  slice(Gy, Gy+Ny),
+                  slice(Gy+Ny, None))
             Xi = (slice(Gx, 2*Gx),
-                   slice(2*Gx,Nx),
-                   slice(Nx, Gx+Nx))
+                  slice(2*Gx, Nx),
+                  slice(Nx, Gx+Nx))
             Yi = (slice(Gy, 2*Gy),
-                   slice(2*Gy,Ny),
-                   slice(Ny, Gy+Ny))
-            Ix = (slice(None,None,+1), slice(None,None,-1))
-            Iy = (slice(None,None,-1), slice(None,None,+1))
+                  slice(2*Gy, Ny),
+                  slice(Ny, Gy+Ny))
+            Ix = (slice(None, None, +1), slice(None, None, -1))
+            Iy = (slice(None, None, -1), slice(None, None, +1))
             data = dF0.data + dF1.data
 
             for ghost_mask in GhostMask.all:
@@ -174,8 +181,8 @@ def test_serial_initialization_2d():
                 sys.stdout.flush()
 
                 if ghost_mask is GhostMask.FULL:
-                    Fx = slice(None,None)
-                    Fy = slice(None,None)
+                    Fx = slice(None, None)
+                    Fy = slice(None, None)
                 elif ghost_mask is GhostMask.CROSS:
                     # we exclude exterior ghosts because pattern is CROSS
                     # we exclude interior ghost because of boundary clashes:
@@ -188,82 +195,82 @@ def test_serial_initialization_2d():
                     raise NotImplementedError(ghost_mask)
 
                 try:
-                    for (i,d) in enumerate(data):
+                    for (i, d) in enumerate(data):
                         b = d.get().handle
 
-                        assert b.shape==(Ny+2*Gy,Nx+2*Gx)
-                        assert b[Yo[1],Xo[1]].shape == npts
+                        assert b.shape == (Ny+2*Gy, Nx+2*Gx)
+                        assert b[Yo[1], Xo[1]].shape == npts
 
-                        assert b[Yo[0],Xo[0]].shape == nghosts
-                        assert b[Yo[0],Xo[2]].shape == nghosts
-                        assert b[Yo[2],Xo[0]].shape == nghosts
-                        assert b[Yo[2],Xo[2]].shape == nghosts
+                        assert b[Yo[0], Xo[0]].shape == nghosts
+                        assert b[Yo[0], Xo[2]].shape == nghosts
+                        assert b[Yo[2], Xo[0]].shape == nghosts
+                        assert b[Yo[2], Xo[2]].shape == nghosts
 
-                        assert b[Yi[0],Xi[0]].shape == nghosts
-                        assert b[Yi[0],Xi[2]].shape == nghosts
-                        assert b[Yi[2],Xi[0]].shape == nghosts
-                        assert b[Yi[2],Xi[2]].shape == nghosts
+                        assert b[Yi[0], Xi[0]].shape == nghosts
+                        assert b[Yi[0], Xi[2]].shape == nghosts
+                        assert b[Yi[2], Xi[0]].shape == nghosts
+                        assert b[Yi[2], Xi[2]].shape == nghosts
 
                         if ghost_mask is GhostMask.FULL:
-                            assert b[Fy,Fx].shape == b.shape
+                            assert b[Fy, Fx].shape == b.shape
                         elif ghost_mask is GhostMask.CROSS:
-                            assert b[Fy,Fx].shape == (Ny-2*Gy, Nx-2*Gx)
+                            assert b[Fy, Fx].shape == (Ny-2*Gy, Nx-2*Gx)
                         else:
                             raise NotImplementedError(ghost_mask)
 
-                        if (Lx==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Fy,Xo[0]] == b[Fy,Xi[2]]), '\n'+str(d)
-                        elif (Lx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Fy,Gx] == 0).all(), '\n'+str(d)
-                            assert (b[Fy,:Gx] == -b[Fy,Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
-                        elif (Lx==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[Fy,:Gx] == +b[Fy,Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
+                        if (Lx == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Fy, Xo[0]] == b[Fy, Xi[2]]), '\n'+str(d)
+                        elif (Lx == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Fy, Gx] == 0).all(), '\n'+str(d)
+                            assert (b[Fy, :Gx] == -b[Fy, Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
+                        elif (Lx == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[Fy, :Gx] == +b[Fy, Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Lx))
 
-                        if (Rx==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Fy,Xo[2]] == b[Fy,Xi[0]]), '\n'+str(d)
-                        elif (Rx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Fy,Nx+Gx-1] == 0).all(), '\n'+str(d)
-                            assert (b[Fy,Nx-1:Nx+Gx-1] == -b[Fy,Nx+Gx:][Ix]).all(), '\n'+str(d)
-                        elif (Rx==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[Fy,Nx-1:Nx+Gx-1] == +b[Fy,Nx+Gx:][Ix]).all(), '\n'+str(d)
+                        if (Rx == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Fy, Xo[2]] == b[Fy, Xi[0]]), '\n'+str(d)
+                        elif (Rx == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Fy, Nx+Gx-1] == 0).all(), '\n'+str(d)
+                            assert (b[Fy, Nx-1:Nx+Gx-1] == -b[Fy, Nx+Gx:][Ix]).all(), '\n'+str(d)
+                        elif (Rx == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[Fy, Nx-1:Nx+Gx-1] == +b[Fy, Nx+Gx:][Ix]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Rx))
 
-                        if (Ly==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Yo[0],Fx] == b[Yi[2],Fx]), '\n'+str(d)
-                        elif (Ly==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Gy,Fx] == 0).all(), '\n'+str(d)
-                            assert (b[:Gy,Fx] == -b[Gy+1:2*Gy+1,Fx][Iy]).all(), '\n'+str(d)
-                        elif (Ly==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[:Gy,Fx] == +b[Gy+1:2*Gy+1,Fx][Iy]).all(), '\n'+str(d)
+                        if (Ly == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Yo[0], Fx] == b[Yi[2], Fx]), '\n'+str(d)
+                        elif (Ly == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Gy, Fx] == 0).all(), '\n'+str(d)
+                            assert (b[:Gy, Fx] == -b[Gy+1:2*Gy+1, Fx][Iy]).all(), '\n'+str(d)
+                        elif (Ly == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[:Gy, Fx] == +b[Gy+1:2*Gy+1, Fx][Iy]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Ly))
 
-                        if (Ry==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Yo[2],Fx] == b[Yi[0],Fx]), '\n'+str(d)
-                        elif (Ry==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Ny+Gy-1,Fx] == 0).all(), '\n'+str(d)
-                            assert (b[Ny-1:Ny+Gy-1,Fx] == -b[Ny+Gy:,Fx][Iy]).all(), '\n'+str(d)
-                        elif (Ry==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[Ny-1:Ny+Gy-1,Fx] == +b[Ny+Gy:,Fx][Iy]).all(), '\n'+str(d)
+                        if (Ry == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Yo[2], Fx] == b[Yi[0], Fx]), '\n'+str(d)
+                        elif (Ry == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Ny+Gy-1, Fx] == 0).all(), '\n'+str(d)
+                            assert (b[Ny-1:Ny+Gy-1, Fx] == -b[Ny+Gy:, Fx][Iy]).all(), '\n'+str(d)
+                        elif (Ry == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[Ny-1:Ny+Gy-1, Fx] == +b[Ny+Gy:, Fx][Iy]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Ry))
 
                         if (ghost_mask is GhostMask.FULL):
-                            if (Lx==Ly==Rx==Ry==BoundaryCondition.PERIODIC):
-                                assert np.all(b[Yo[0],Xo[0]]==b[Yi[2],Xi[2]]), '\n'+str(d)
-                                assert np.all(b[Yo[2],Xo[0]]==b[Yi[0],Xi[2]]), '\n'+str(d)
-                                assert np.all(b[Yo[2],Xo[2]]==b[Yi[0],Xi[0]]), '\n'+str(d)
-                                assert np.all(b[Yo[0],Xo[2]]==b[Yi[2],Xi[0]]), '\n'+str(d)
+                            if (Lx == Ly == Rx == Ry == BoundaryCondition.PERIODIC):
+                                assert np.all(b[Yo[0], Xo[0]] == b[Yi[2], Xi[2]]), '\n'+str(d)
+                                assert np.all(b[Yo[2], Xo[0]] == b[Yi[0], Xi[2]]), '\n'+str(d)
+                                assert np.all(b[Yo[2], Xo[2]] == b[Yi[0], Xi[0]]), '\n'+str(d)
+                                assert np.all(b[Yo[0], Xo[2]] == b[Yi[2], Xi[0]]), '\n'+str(d)
                         elif (ghost_mask is GhostMask.CROSS):
-                            assert np.all(np.isnan(b[Yo[0],Xo[0]])), '\n'+str(d)
-                            assert np.all(np.isnan(b[Yo[2],Xo[0]])), '\n'+str(d)
-                            assert np.all(np.isnan(b[Yo[2],Xo[2]])), '\n'+str(d)
-                            assert np.all(np.isnan(b[Yo[0],Xo[2]])), '\n'+str(d)
+                            assert np.all(np.isnan(b[Yo[0], Xo[0]])), '\n'+str(d)
+                            assert np.all(np.isnan(b[Yo[2], Xo[0]])), '\n'+str(d)
+                            assert np.all(np.isnan(b[Yo[2], Xo[2]])), '\n'+str(d)
+                            assert np.all(np.isnan(b[Yo[0], Xo[2]])), '\n'+str(d)
                         else:
-                            msg='Unknown ghost mask {}.'.format(ghost_mask)
+                            msg = 'Unknown ghost mask {}.'.format(ghost_mask)
                             raise NotImplementedError(msg)
                         sys.stdout.write('.')
                         sys.stdout.flush()
@@ -275,23 +282,23 @@ def test_serial_initialization_3d():
     print()
     print('test_serial_initialization_3d()')
     dim = 3
-    npts = (8,5,5)
-    nghosts = (3,1,2)
+    npts = (8, 5, 5)
+    nghosts = (3, 1, 2)
     for (lbd, rbd) in domain_boundary_iterator(dim):
         domain = Box(dim=dim, lboundaries=lbd,
-                              rboundaries=rbd)
+                     rboundaries=rbd)
         F0 = Field(domain=domain, name='F0', nb_components=1)
         F1 = Field(domain=domain, name='F1', nb_components=3)
         print('[{}]'.format(F0.format_boundaries()))
 
         discretization = CartesianDiscretization(npts, nghosts,
-                                            lboundaries=F0.lboundaries,
-                                            rboundaries=F0.rboundaries)
+                                                 lboundaries=F0.lboundaries,
+                                                 rboundaries=F0.rboundaries)
 
         topo0 = CartesianTopology(domain=domain, discretization=discretization,
-                        backend=Backend.HOST)
+                                  backend=Backend.HOST)
         topos = (topo0,) + tuple(CartesianTopology(domain=domain, discretization=discretization,
-            backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
+                                                   backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
 
         assert all(np.all(t.mesh.local_lboundaries == F0.lboundaries) for t in topos)
         assert all(np.all(t.mesh.local_rboundaries == F0.rboundaries) for t in topos)
@@ -307,31 +314,31 @@ def test_serial_initialization_3d():
             assert all(np.all(d.local_lboundaries == F0.lboundaries) for d in dfields)
             assert all(np.all(d.local_rboundaries == F0.rboundaries) for d in dfields)
 
-            Nz,Ny,Nx = npts
-            Gz,Gy,Gx = nghosts
-            Lz,Ly,Lx = tuple(discretization.lboundaries)
-            Rz,Ry,Rx = tuple(discretization.rboundaries)
-            Xo = (slice(0,Gx),
-                   slice(Gx,Gx+Nx),
-                   slice(Gx+Nx,None))
-            Yo = (slice(0,Gy),
-                   slice(Gy,Gy+Ny),
-                   slice(Gy+Ny,None))
-            Zo = (slice(0,Gz),
-                   slice(Gz,Gz+Nz),
-                   slice(Gz+Nz,None))
+            Nz, Ny, Nx = npts
+            Gz, Gy, Gx = nghosts
+            Lz, Ly, Lx = tuple(discretization.lboundaries)
+            Rz, Ry, Rx = tuple(discretization.rboundaries)
+            Xo = (slice(0, Gx),
+                  slice(Gx, Gx+Nx),
+                  slice(Gx+Nx, None))
+            Yo = (slice(0, Gy),
+                  slice(Gy, Gy+Ny),
+                  slice(Gy+Ny, None))
+            Zo = (slice(0, Gz),
+                  slice(Gz, Gz+Nz),
+                  slice(Gz+Nz, None))
             Xi = (slice(Gx, 2*Gx),
-                   slice(2*Gx,Nx),
-                   slice(Nx, Gx+Nx))
+                  slice(2*Gx, Nx),
+                  slice(Nx, Gx+Nx))
             Yi = (slice(Gy, 2*Gy),
-                   slice(2*Gy,Ny),
-                   slice(Ny, Gy+Ny))
+                  slice(2*Gy, Ny),
+                  slice(Ny, Gy+Ny))
             Zi = (slice(Gz, 2*Gz),
-                   slice(2*Gz,Nz),
-                   slice(Nz, Gz+Nz))
-            Ix = (slice(None,None,+1), slice(None,None,+1), slice(None,None,-1))
-            Iy = (slice(None,None,+1), slice(None,None,-1), slice(None,None,+1))
-            Iz = (slice(None,None,-1), slice(None,None,+1), slice(None,None,+1))
+                  slice(2*Gz, Nz),
+                  slice(Nz, Gz+Nz))
+            Ix = (slice(None, None, +1), slice(None, None, +1), slice(None, None, -1))
+            Iy = (slice(None, None, +1), slice(None, None, -1), slice(None, None, +1))
+            Iz = (slice(None, None, -1), slice(None, None, +1), slice(None, None, +1))
             data = dF0.data + dF1.data
 
             for ghost_mask in GhostMask.all:
@@ -342,9 +349,9 @@ def test_serial_initialization_3d():
                 sys.stdout.flush()
 
                 if ghost_mask is GhostMask.FULL:
-                    Fx = slice(None,None)
-                    Fy = slice(None,None)
-                    Fz = slice(None,None)
+                    Fx = slice(None, None)
+                    Fy = slice(None, None)
+                    Fz = slice(None, None)
                 elif ghost_mask is GhostMask.CROSS:
                     # we exclude exterior ghosts because pattern is CROSS
                     # we exclude interior ghost because of boundary clashes:
@@ -358,117 +365,117 @@ def test_serial_initialization_3d():
                     raise NotImplementedError(ghost_mask)
 
                 try:
-                    for (i,d) in enumerate(data):
+                    for (i, d) in enumerate(data):
                         b = d.get().handle
-                        assert b.shape==(Nz+2*Gz,Ny+2*Gy,Nx+2*Gx)
-                        assert b[Zo[1],Yo[1],Xo[1]].shape == npts
-
-                        assert b[Zo[0],Yo[0],Xo[0]].shape == nghosts
-                        assert b[Zo[0],Yo[0],Xo[2]].shape == nghosts
-                        assert b[Zo[0],Yo[2],Xo[0]].shape == nghosts
-                        assert b[Zo[0],Yo[2],Xo[2]].shape == nghosts
-                        assert b[Zo[2],Yo[0],Xo[0]].shape == nghosts
-                        assert b[Zo[2],Yo[0],Xo[2]].shape == nghosts
-                        assert b[Zo[2],Yo[2],Xo[0]].shape == nghosts
-                        assert b[Zo[2],Yo[2],Xo[2]].shape == nghosts
-
-                        assert b[Zi[0],Yi[0],Xi[0]].shape == nghosts
-                        assert b[Zi[0],Yi[0],Xi[2]].shape == nghosts
-                        assert b[Zi[0],Yi[2],Xi[0]].shape == nghosts
-                        assert b[Zi[0],Yi[2],Xi[2]].shape == nghosts
-                        assert b[Zi[2],Yi[0],Xi[0]].shape == nghosts
-                        assert b[Zi[2],Yi[0],Xi[2]].shape == nghosts
-                        assert b[Zi[2],Yi[2],Xi[0]].shape == nghosts
-                        assert b[Zi[2],Yi[2],Xi[2]].shape == nghosts
+                        assert b.shape == (Nz+2*Gz, Ny+2*Gy, Nx+2*Gx)
+                        assert b[Zo[1], Yo[1], Xo[1]].shape == npts
+
+                        assert b[Zo[0], Yo[0], Xo[0]].shape == nghosts
+                        assert b[Zo[0], Yo[0], Xo[2]].shape == nghosts
+                        assert b[Zo[0], Yo[2], Xo[0]].shape == nghosts
+                        assert b[Zo[0], Yo[2], Xo[2]].shape == nghosts
+                        assert b[Zo[2], Yo[0], Xo[0]].shape == nghosts
+                        assert b[Zo[2], Yo[0], Xo[2]].shape == nghosts
+                        assert b[Zo[2], Yo[2], Xo[0]].shape == nghosts
+                        assert b[Zo[2], Yo[2], Xo[2]].shape == nghosts
+
+                        assert b[Zi[0], Yi[0], Xi[0]].shape == nghosts
+                        assert b[Zi[0], Yi[0], Xi[2]].shape == nghosts
+                        assert b[Zi[0], Yi[2], Xi[0]].shape == nghosts
+                        assert b[Zi[0], Yi[2], Xi[2]].shape == nghosts
+                        assert b[Zi[2], Yi[0], Xi[0]].shape == nghosts
+                        assert b[Zi[2], Yi[0], Xi[2]].shape == nghosts
+                        assert b[Zi[2], Yi[2], Xi[0]].shape == nghosts
+                        assert b[Zi[2], Yi[2], Xi[2]].shape == nghosts
 
                         if ghost_mask is GhostMask.FULL:
-                            assert b[Fz,Fy,Fx].shape == b.shape
+                            assert b[Fz, Fy, Fx].shape == b.shape
                         elif ghost_mask is GhostMask.CROSS:
-                            assert b[Fz,Fy,Fx].shape == (Nz-2*Gz,Ny-2*Gy, Nx-2*Gx)
+                            assert b[Fz, Fy, Fx].shape == (Nz-2*Gz, Ny-2*Gy, Nx-2*Gx)
                         else:
                             raise NotImplementedError(ghost_mask)
 
-                        if (Lx==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Fz,Fy,Xo[0]] == b[Fz,Fy,Xi[2]]), '\n'+str(d)
-                        elif (Lx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Fz,Fy,Gx] == 0).all(), '\n'+str(d)
-                            assert (b[Fz,Fy,:Gx] == -b[Fz,Fy,Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
-                        elif (Lx==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[Fz,Fy,:Gx] == +b[Fz,Fy,Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
+                        if (Lx == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Fz, Fy, Xo[0]] == b[Fz, Fy, Xi[2]]), '\n'+str(d)
+                        elif (Lx == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Fz, Fy, Gx] == 0).all(), '\n'+str(d)
+                            assert (b[Fz, Fy, :Gx] == -b[Fz, Fy, Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
+                        elif (Lx == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[Fz, Fy, :Gx] == +b[Fz, Fy, Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Lx))
 
-                        if (Rx==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Fz,Fy,Xo[2]] == b[Fz,Fy,Xi[0]]), '\n'+str(d)
-                        elif (Rx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Fz,Fy,Nx+Gx-1] == 0).all(), '\n'+str(d)
-                            assert (b[Fz,Fy,Nx-1:Nx+Gx-1] == -b[Fz,Fy,Nx+Gx:][Ix]).all(), '\n'+str(d)
-                        elif (Rx==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[Fz,Fy,Nx-1:Nx+Gx-1] == +b[Fz,Fy,Nx+Gx:][Ix]).all(), '\n'+str(d)
+                        if (Rx == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Fz, Fy, Xo[2]] == b[Fz, Fy, Xi[0]]), '\n'+str(d)
+                        elif (Rx == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Fz, Fy, Nx+Gx-1] == 0).all(), '\n'+str(d)
+                            assert (b[Fz, Fy, Nx-1:Nx+Gx-1] == -b[Fz, Fy, Nx+Gx:][Ix]).all(), '\n'+str(d)
+                        elif (Rx == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[Fz, Fy, Nx-1:Nx+Gx-1] == +b[Fz, Fy, Nx+Gx:][Ix]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Rx))
 
-                        if (Ly==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Fz,Yo[0],Fx] == b[Fz,Yi[2],Fx]), '\n'+str(d)
-                        elif (Ly==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Fz,Gy,Fx] == 0).all(), '\n'+str(d)
-                            assert (b[Fz,:Gy,Fx] == -b[Fz,Gy+1:2*Gy+1,Fx][Iy]).all(), '\n'+str(d)
-                        elif (Ly==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[Fz,:Gy,Fx] == +b[Fz,Gy+1:2*Gy+1,Fx][Iy]).all(), '\n'+str(d)
+                        if (Ly == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Fz, Yo[0], Fx] == b[Fz, Yi[2], Fx]), '\n'+str(d)
+                        elif (Ly == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Fz, Gy, Fx] == 0).all(), '\n'+str(d)
+                            assert (b[Fz, :Gy, Fx] == -b[Fz, Gy+1:2*Gy+1, Fx][Iy]).all(), '\n'+str(d)
+                        elif (Ly == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[Fz, :Gy, Fx] == +b[Fz, Gy+1:2*Gy+1, Fx][Iy]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Ly))
 
-                        if (Ry==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Fz,Yo[2],Fx] == b[Fz,Yi[0],Fx]), '\n'+str(d)
-                        elif (Ry==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Fz,Ny+Gy-1,Fx] == 0).all(), '\n'+str(d)
-                            assert (b[Fz,Ny-1:Ny+Gy-1,Fx] == -b[Fz,Ny+Gy:,Fx][Iy]).all(), '\n'+str(d)
-                        elif (Ry==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[Fz,Ny-1:Ny+Gy-1,Fx] == +b[Fz,Ny+Gy:,Fx][Iy]).all(), '\n'+str(d)
+                        if (Ry == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Fz, Yo[2], Fx] == b[Fz, Yi[0], Fx]), '\n'+str(d)
+                        elif (Ry == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Fz, Ny+Gy-1, Fx] == 0).all(), '\n'+str(d)
+                            assert (b[Fz, Ny-1:Ny+Gy-1, Fx] == -b[Fz, Ny+Gy:, Fx][Iy]).all(), '\n'+str(d)
+                        elif (Ry == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[Fz, Ny-1:Ny+Gy-1, Fx] == +b[Fz, Ny+Gy:, Fx][Iy]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Ry))
 
-                        if (Lz==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Zo[0],Fy,Fx] == b[Zi[2],Fy,Fx]), '\n'+str(d)
-                        elif (Lz==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Gz,Fy,Fx] == 0).all(), '\n'+str(d)
-                            assert (b[:Gz,Fy,Fx] == -b[Gz+1:2*Gz+1,Fy,Fx][Iz]).all(), '\n'+str(d)
-                        elif (Lz==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[:Gz,Fy,Fx] == +b[Gz+1:2*Gz+1,Fy,Fx][Iz]).all(), '\n'+str(d)
+                        if (Lz == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Zo[0], Fy, Fx] == b[Zi[2], Fy, Fx]), '\n'+str(d)
+                        elif (Lz == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Gz, Fy, Fx] == 0).all(), '\n'+str(d)
+                            assert (b[:Gz, Fy, Fx] == -b[Gz+1:2*Gz+1, Fy, Fx][Iz]).all(), '\n'+str(d)
+                        elif (Lz == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[:Gz, Fy, Fx] == +b[Gz+1:2*Gz+1, Fy, Fx][Iz]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Lz))
 
-                        if (Rz==BoundaryCondition.PERIODIC):
-                            assert np.all(b[Zo[2],Fy,Fx] == b[Zi[0],Fy,Fx]), '\n'+str(d)
-                        elif (Rz==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
-                            assert (b[Nz+Gz-1,Fy,Fx] == 0).all(), '\n'+str(d)
-                            assert (b[Nz-1:Nz+Gz-1,Fy,Fx] == -b[Nz+Gz:,Fy,Fx][Iz]).all(), '\n'+str(d)
-                        elif (Rz==BoundaryCondition.HOMOGENEOUS_NEUMANN):
-                            assert (b[Nz-1:Nz+Gz-1,Fy,Fx] == +b[Nz+Gz:,Fy,Fx][Iz]).all(), '\n'+str(d)
+                        if (Rz == BoundaryCondition.PERIODIC):
+                            assert np.all(b[Zo[2], Fy, Fx] == b[Zi[0], Fy, Fx]), '\n'+str(d)
+                        elif (Rz == BoundaryCondition.HOMOGENEOUS_DIRICHLET):
+                            assert (b[Nz+Gz-1, Fy, Fx] == 0).all(), '\n'+str(d)
+                            assert (b[Nz-1:Nz+Gz-1, Fy, Fx] == -b[Nz+Gz:, Fy, Fx][Iz]).all(), '\n'+str(d)
+                        elif (Rz == BoundaryCondition.HOMOGENEOUS_NEUMANN):
+                            assert (b[Nz-1:Nz+Gz-1, Fy, Fx] == +b[Nz+Gz:, Fy, Fx][Iz]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Rz))
 
                         if (ghost_mask is GhostMask.FULL):
-                            if (Lx==Ly==Lz==Rx==Ry==Rz==BoundaryCondition.PERIODIC):
-                                assert np.all(b[Zo[0],Yo[0],Xo[0]]==b[Zi[2],Yi[2],Xi[2]])
-                                assert np.all(b[Zo[2],Yo[0],Xo[0]]==b[Zi[0],Yi[2],Xi[2]])
-                                assert np.all(b[Zo[2],Yo[2],Xo[0]]==b[Zi[0],Yi[0],Xi[2]])
-                                assert np.all(b[Zo[0],Yo[2],Xo[0]]==b[Zi[2],Yi[0],Xi[2]])
-                                assert np.all(b[Zo[0],Yo[0],Xo[2]]==b[Zi[2],Yi[2],Xi[0]])
-                                assert np.all(b[Zo[2],Yo[0],Xo[2]]==b[Zi[0],Yi[2],Xi[0]])
-                                assert np.all(b[Zo[2],Yo[2],Xo[2]]==b[Zi[0],Yi[0],Xi[0]])
-                                assert np.all(b[Zo[0],Yo[2],Xo[2]]==b[Zi[2],Yi[0],Xi[0]])
+                            if (Lx == Ly == Lz == Rx == Ry == Rz == BoundaryCondition.PERIODIC):
+                                assert np.all(b[Zo[0], Yo[0], Xo[0]] == b[Zi[2], Yi[2], Xi[2]])
+                                assert np.all(b[Zo[2], Yo[0], Xo[0]] == b[Zi[0], Yi[2], Xi[2]])
+                                assert np.all(b[Zo[2], Yo[2], Xo[0]] == b[Zi[0], Yi[0], Xi[2]])
+                                assert np.all(b[Zo[0], Yo[2], Xo[0]] == b[Zi[2], Yi[0], Xi[2]])
+                                assert np.all(b[Zo[0], Yo[0], Xo[2]] == b[Zi[2], Yi[2], Xi[0]])
+                                assert np.all(b[Zo[2], Yo[0], Xo[2]] == b[Zi[0], Yi[2], Xi[0]])
+                                assert np.all(b[Zo[2], Yo[2], Xo[2]] == b[Zi[0], Yi[0], Xi[0]])
+                                assert np.all(b[Zo[0], Yo[2], Xo[2]] == b[Zi[2], Yi[0], Xi[0]])
                         elif (ghost_mask is GhostMask.CROSS):
-                            assert np.all(np.isnan(b[Zo[0],Yo[0],Xo[0]]))
-                            assert np.all(np.isnan(b[Zo[2],Yo[0],Xo[0]]))
-                            assert np.all(np.isnan(b[Zo[2],Yo[2],Xo[0]]))
-                            assert np.all(np.isnan(b[Zo[0],Yo[2],Xo[0]]))
-                            assert np.all(np.isnan(b[Zo[0],Yo[0],Xo[2]]))
-                            assert np.all(np.isnan(b[Zo[2],Yo[0],Xo[2]]))
-                            assert np.all(np.isnan(b[Zo[2],Yo[2],Xo[2]]))
-                            assert np.all(np.isnan(b[Zo[0],Yo[2],Xo[2]]))
+                            assert np.all(np.isnan(b[Zo[0], Yo[0], Xo[0]]))
+                            assert np.all(np.isnan(b[Zo[2], Yo[0], Xo[0]]))
+                            assert np.all(np.isnan(b[Zo[2], Yo[2], Xo[0]]))
+                            assert np.all(np.isnan(b[Zo[0], Yo[2], Xo[0]]))
+                            assert np.all(np.isnan(b[Zo[0], Yo[0], Xo[2]]))
+                            assert np.all(np.isnan(b[Zo[2], Yo[0], Xo[2]]))
+                            assert np.all(np.isnan(b[Zo[2], Yo[2], Xo[2]]))
+                            assert np.all(np.isnan(b[Zo[0], Yo[2], Xo[2]]))
                         else:
-                            msg='Unknown ghost mask {}.'.format(ghost_mask)
+                            msg = 'Unknown ghost mask {}.'.format(ghost_mask)
                             raise NotImplementedError(msg)
                         sys.stdout.write('.')
                         sys.stdout.flush()
@@ -481,6 +488,7 @@ def iter_backends():
     for cl_env in iter_clenv():
         yield (Backend.OPENCL, cl_env)
 
+
 def test_mpi_ghost_exchange_periodic(comm=None):
     if comm is None:
         from mpi4py import MPI
@@ -491,7 +499,7 @@ def test_mpi_ghost_exchange_periodic(comm=None):
               np.int16,  np.int32,  np.int64,
               np.uint16, np.uint32, np.uint64)
     assert size-1 < len(dtypes)
-    if rank==0:
+    if rank == 0:
         print()
         msg = '*** COMM_WORLD_SIZE {} ***'.format(size)
         print()
@@ -499,98 +507,99 @@ def test_mpi_ghost_exchange_periodic(comm=None):
         print(msg)
         print('*'*len(msg))
         print('test_mpi_ghost_exchange_periodic()'.format(size))
-    for dim in range(1,3+__ENABLE_LONG_TESTS__):
-        if rank==0:
+    for dim in range(1, 3+__ENABLE_LONG_TESTS__):
+        if rank == 0:
             print(('  >DIM={}'.format(dim)))
 
-        npts = (53,47,59,23)[:dim]
-        nghosts = (2,1,0,3)[:dim]
+        npts = (53, 47, 59, 23)[:dim]
+        nghosts = (2, 1, 0, 3)[:dim]
         discretization = CartesianDiscretization(npts, nghosts,
-                                        default_boundaries=True)
+                                                 default_boundaries=True)
         domain = Box(dim=dim)
 
         for dtype in dtypes[size-1:size]:
-            if rank==0:
+            if rank == 0:
                 print('    >DTYPE={}'.format(dtype))
 
             F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype)
             F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype)
-            F2 = Field(domain=domain, name='F2', shape=(2,2), dtype=dtype)
+            F2 = Field(domain=domain, name='F2', shape=(2, 2), dtype=dtype)
 
             for (backend, cl_env) in iter_backends():
-                if rank==0:
+                if rank == 0:
                     print('      >BACKEND.{}{}'.format(backend,
-                            '' if (cl_env is None) else '::{}.{}'.format(
-                                cl_env.platform.name.strip(), cl_env.device.name.strip())))
+                                                       '' if (cl_env is None) else '::{}.{}'.format(
+                                                           cl_env.platform.name.strip(), cl_env.device.name.strip())))
                 for shape in it.product(range(0, size+1), repeat=dim):
-                    if np.prod(shape, dtype=np.uint32)!=size:
+                    if np.prod(shape, dtype=np.uint32) != size:
                         continue
-                    if rank==0:
+                    if rank == 0:
                         print('         *cart_shape: {}'.format(shape))
                     topo = CartesianTopology(domain=domain, discretization=discretization,
-                                               backend=backend, cart_shape=shape, cl_env=cl_env)
-                    assert (topo.proc_shape==shape).all()
+                                             backend=backend, cart_shape=shape, cl_env=cl_env)
+                    assert (topo.proc_shape == shape).all()
 
-                    def ghost_base(i,d,rank,local_dir):
+                    def ghost_base(i, d, rank, local_dir):
                         return (i+1)*17 + (d+1)*13 + (local_dir+1)*11 + (rank+1)*7
-                    def ghost_vals(shape,dtype,i,d,rank,local_dir):
-                        if (shape is None) or (len(shape)==0):
+
+                    def ghost_vals(shape, dtype, i, d, rank, local_dir):
+                        if (shape is None) or (len(shape) == 0):
                             raise ValueError('Shape is None or an empty tuple: {}'.format(shape))
                         base = np.full(shape=shape, dtype=dtype, fill_value=ghost_base(i, d, rank,
-                                                                                local_dir))
-                        I = np.ix_(*tuple(np.arange(shape[direction], dtype=dtype) \
-                                                for direction in range(len(shape))))
+                                                                                       local_dir))
+                        I = np.ix_(*tuple(np.arange(shape[direction], dtype=dtype)
+                                          for direction in range(len(shape))))
                         return base + I[d]
 
                     for F in (F0, F1, F2):
                         dF = F.discretize(topo)
-                        if rank==0:
+                        if rank == 0:
                             print('          |{} COMPONENT(S)'.format(F.nb_components))
                         for exchange_method in (ExchangeMethod.ISEND_IRECV,
                                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V,
                                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_W):
-                            if rank==0:
+                            if rank == 0:
                                 print('             ExchangeMethod.{} |'.format(exchange_method), end=' ')
                             sys.stdout.flush()
                             for d in range(dim):
                                 dF.initialize(__zero_init, exchange_ghosts=False)
-                                if rank==0:
+                                if rank == 0:
                                     print(DirectionLabels[dim-1-d], end=' ')
                                 sys.stdout.flush()
                                 lghosts, rghosts, shape = dF.inner_ghost_slices[d]
                                 _lghosts, _rghosts, shape = dF.outer_ghost_slices[d]
 
                                 if (shape is not None):
-                                    for (i,data) in enumerate(dF.data):
-                                            data[lghosts] = ghost_vals(shape, dtype, i,d,rank,0)
-                                            data[_lghosts] = -10
-                                            data[rghosts] = ghost_vals(shape, dtype, i,d,rank,1)
-                                            data[_rghosts] = +10
+                                    for (i, data) in enumerate(dF.data):
+                                        data[lghosts] = ghost_vals(shape, dtype, i, d, rank, 0)
+                                        data[_lghosts] = -10
+                                        data[rghosts] = ghost_vals(shape, dtype, i, d, rank, 1)
+                                        data[_rghosts] = +10
 
                                 dF.exchange_ghosts(directions=d, exchange_method=exchange_method)
 
                                 lghosts, rghosts, shape = dF.outer_ghost_slices[d]
-                                left_rank, right_rank = topo.proc_neighbour_ranks[:,d]
+                                left_rank, right_rank = topo.proc_neighbour_ranks[:, d]
                                 if (left_rank == -1):
-                                    assert right_rank==-1
+                                    assert right_rank == -1
                                     left_rank, right_rank = rank, rank
 
                                 if (shape is not None):
-                                    for (i,data) in enumerate(dF.data):
+                                    for (i, data) in enumerate(dF.data):
                                         ldata = data[lghosts]
                                         rdata = data[rghosts]
 
                                         ldata = np.atleast_1d(ldata.get())
                                         target_vals = ghost_vals(ldata.shape, dtype, i, d,
-                                                left_rank,1)
+                                                                 left_rank, 1)
                                         assert np.allclose(ldata, target_vals), (rank,
-                                                                                    target_vals)
+                                                                                 target_vals)
                                         rdata = np.atleast_1d(rdata.get())
                                         target_vals = ghost_vals(rdata.shape, dtype, i, d,
-                                                                                right_rank, 0)
+                                                                 right_rank, 0)
                                         assert np.allclose(rdata, target_vals), (rank,
-                                                                                        target_vals)
-                            if rank==0:
+                                                                                 target_vals)
+                            if rank == 0:
                                 print()
 
 
@@ -603,7 +612,7 @@ def test_mpi_ghost_exchange_runtime(comm=None):
     size = comm.Get_size()
     dtype = np.float32
 
-    if rank==0:
+    if rank == 0:
         print()
         msg = '*** COMM_WORLD_SIZE {} ***'.format(size)
         print()
@@ -613,30 +622,30 @@ def test_mpi_ghost_exchange_runtime(comm=None):
         print('test_mpi_ghost_exchange_runtime()'.format(size))
 
     for dim in range(1, 3+__ENABLE_LONG_TESTS__):
-        if rank==0:
+        if rank == 0:
             sys.stdout.write('>DIM={}\n'.format(dim))
 
-        npts = (17,16,19)[:dim]
-        nghosts = (2,1,3)[:dim]
+        npts = (17, 16, 19)[:dim]
+        nghosts = (2, 1, 3)[:dim]
 
         for shape in it.product(range(0, size+1), repeat=dim):
-            if np.prod(shape, dtype=np.uint32)!=size:
+            if np.prod(shape, dtype=np.uint32) != size:
                 continue
-            if rank==0:
+            if rank == 0:
                 sys.stdout.write('  >CART SHAPE: {}\n'.format(shape))
 
             for (backend, cl_env) in iter_backends():
-                if rank==0:
+                if rank == 0:
                     sys.stdout.write('    >BACKEND.{:<7} '.format(str(backend)+':'))
 
                 def breakline(i):
-                    if (rank==0) and ((i+1)%63==0):
+                    if (rank == 0) and ((i+1) % 63 == 0):
                         sys.stdout.write('\n' + ' '*21)
                         sys.stdout.flush()
                         return True
                     return False
 
-                i=0
+                i = 0
                 brk = False
                 try:
                     for (lbd, rbd) in domain_boundary_iterator(dim):
@@ -645,32 +654,33 @@ def test_mpi_ghost_exchange_runtime(comm=None):
                         F = Field(domain=domain, name='F', nb_components=1, dtype=dtype)
 
                         discretization = CartesianDiscretization(npts, nghosts,
-                                                lboundaries=F.lboundaries,
-                                                rboundaries=F.rboundaries)
+                                                                 lboundaries=F.lboundaries,
+                                                                 rboundaries=F.rboundaries)
 
                         topo = CartesianTopology(domain=domain, discretization=discretization,
-                                                   backend=backend, cart_shape=shape, cl_env=cl_env)
-                        assert (topo.proc_shape==shape).all()
+                                                 backend=backend, cart_shape=shape, cl_env=cl_env)
+                        assert (topo.proc_shape == shape).all()
 
                         dF = F.discretize(topo)
+                        dF.initialize(__random_init, exchange_ghosts=False)
                         for exchange_method in (ExchangeMethod.ISEND_IRECV,
                                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V,
                                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_W):
-                                for d in range(dim):
-                                    dF.exchange_ghosts(directions=d, exchange_method=exchange_method)
-                                    if rank==0:
-                                        sys.stdout.write('.')
-                                        brk = breakline(i)
-                                        i+=1
-
-                                    dF.accumulate_ghosts(directions=d, exchange_method=exchange_method)
-                                    if rank==0:
-                                        sys.stdout.write('.')
-                                        brk = breakline(i)
-                                        i+=1
-                                        sys.stdout.flush()
+                            for d in range(dim):
+                                dF.exchange_ghosts(directions=d, exchange_method=exchange_method)
+                                if rank == 0:
+                                    sys.stdout.write('.')
+                                    brk = breakline(i)
+                                    i += 1
+
+                                dF.accumulate_ghosts(directions=d, exchange_method=exchange_method)
+                                if rank == 0:
+                                    sys.stdout.write('.')
+                                    brk = breakline(i)
+                                    i += 1
+                                    sys.stdout.flush()
                 finally:
-                    if (rank==0):
+                    if (rank == 0):
                         sys.stdout.write('\n')
                         sys.stdout.flush()
 
@@ -681,7 +691,7 @@ def test_mpi_ghost_accumulate_periodic(comm=None):
         comm = MPI.COMM_WORLD
     rank = comm.Get_rank()
     size = comm.Get_size()
-    if rank==0:
+    if rank == 0:
         msg = '*** COMM_WORLD_SIZE {} ***'.format(size)
         print()
         print('*'*len(msg))
@@ -694,106 +704,108 @@ def test_mpi_ghost_accumulate_periodic(comm=None):
               np.int16,  np.int32,  np.int64,
               np.uint16, np.uint32, np.uint64)
     assert size-1 < len(dtypes)
-    for dim in range(1,3+__ENABLE_LONG_TESTS__):
-        if rank==0:
+    for dim in range(1, 3+__ENABLE_LONG_TESTS__):
+        if rank == 0:
             print(('  >DIM={}'.format(dim)))
 
-        npts    = (53,57,51,49)[:dim]
-        nghosts = (1,3,0,2)[:dim]
+        npts = (53, 57, 51, 49)[:dim]
+        nghosts = (1, 3, 0, 2)[:dim]
         discretization = CartesianDiscretization(npts, nghosts,
-                                        default_boundaries=True)
+                                                 default_boundaries=True)
         domain = Box(dim=dim)
 
         for dtype in dtypes[size-1:size]:
-            if rank==0:
+            if rank == 0:
                 print('    >DTYPE={}'.format(dtype))
 
             F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype)
             F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype)
-            F2 = Field(domain=domain, name='F2', shape=(2,2), dtype=dtype)
+            F2 = Field(domain=domain, name='F2', shape=(2, 2), dtype=dtype)
 
             for (backend, cl_env) in iter_backends():
-                if rank==0:
+                if rank == 0:
                     print('      >BACKEND.{}{}'.format(backend,
-                            '' if (cl_env is None) else '::{}.{}'.format(
-                                cl_env.platform.name.strip(), cl_env.device.name.strip())))
-                for shape in it.product(range(0,size+1), repeat=dim):
-                    if np.prod(shape, dtype=np.uint32)!=size:
+                                                       '' if (cl_env is None) else '::{}.{}'.format(
+                                                           cl_env.platform.name.strip(), cl_env.device.name.strip())))
+                for shape in it.product(range(0, size+1), repeat=dim):
+                    if np.prod(shape, dtype=np.uint32) != size:
                         continue
-                    if rank==0:
+                    if rank == 0:
                         print('        *cart_shape: {}'.format(shape))
                     topo = CartesianTopology(domain=domain, discretization=discretization,
-                                               backend=backend, cart_shape=shape, cl_env=cl_env)
-                    assert (topo.proc_shape==shape).all()
+                                             backend=backend, cart_shape=shape, cl_env=cl_env)
+                    assert (topo.proc_shape == shape).all()
 
                     def ghost_base(rank, directions, displacements, i):
-                        disweight     = np.asarray([19,23,29,31], dtype=np.int32)
-                        dirweight     = np.asarray([37,41,43,51], dtype=np.int32)
-                        directions    = np.asarray(directions, dtype=np.int32) + 1
+                        disweight = np.asarray([19, 23, 29, 31], dtype=np.int32)
+                        dirweight = np.asarray([37, 41, 43, 51], dtype=np.int32)
+                        directions = np.asarray(directions, dtype=np.int32) + 1
                         displacements = np.asarray(displacements, dtype=np.int32) + 2
-                        tag =  (rank+1)*17 + (i+1)*13
+                        tag = (rank+1)*17 + (i+1)*13
                         tag += dirweight[:directions.size].dot(directions)
                         tag += disweight[:displacements.size].dot(displacements)
                         return tag
+
                     def ghost_vals(shape, dtype, rank, directions, displacements, i):
-                        if (shape is None) or (len(shape)==0):
+                        if (shape is None) or (len(shape) == 0):
                             raise ValueError('Shape is None or an empty tuple: {}'.format(shape))
-                        base_value = ghost_base(rank,directions,displacements,i)
+                        base_value = ghost_base(rank, directions, displacements, i)
                         vals = np.full(shape=shape, fill_value=base_value, dtype=dtype)
                         return vals
 
-                    for F in (F0,F1,F2):
+                    for F in (F0, F1, F2):
                         dF = F.discretize(topo)
-                        proc_ranks  = topo.proc_ranks
-                        proc_shape  = topo.proc_shape
+                        dF.initialize(__random_init, exchange_ghosts=False)
+                        proc_ranks = topo.proc_ranks
+                        proc_shape = topo.proc_shape
                         proc_coords = tuple(topo.proc_coords.tolist())
                         assert proc_ranks[proc_coords] == rank
                         all_inner_ghost_slices = dF.all_inner_ghost_slices
                         all_outer_ghost_slices = dF.all_outer_ghost_slices
-                        if rank==0:
+                        if rank == 0:
                             print('          |{} COMPONENT(S)'.format(F.nb_components))
                         for exchange_method in (ExchangeMethod.ISEND_IRECV,
                                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V,
                                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_W,):
-                            if rank==0:
+                            if rank == 0:
                                 print('             ExchangeMethod.{:<25} |'.format(str(exchange_method)), end=' ')
                             sys.stdout.flush()
 
                             # test one direction at a time
                             max_displacements = 1
-                            for ndirections in range(1,dim+1):
-                                all_displacements = tuple(it.product((-1,0,+1), repeat=ndirections))
-                                all_directions    = tuple(it.combinations(range(dim), ndirections))
-                                masks             = tuple(it.product((0,1), repeat=ndirections))
+                            for ndirections in range(1, dim+1):
+                                all_displacements = tuple(it.product((-1, 0, +1), repeat=ndirections))
+                                all_directions = tuple(it.combinations(range(dim), ndirections))
+                                masks = tuple(it.product((0, 1), repeat=ndirections))
                                 for directions in all_directions:
-                                    if rank==0:
+                                    if rank == 0:
                                         if directions:
                                             print(''.join(DirectionLabels[dim-1-d]
-                                                            for d in directions), end=' ')
+                                                          for d in directions), end=' ')
                                         else:
                                             print('--', end=' ')
                                         sys.stdout.flush()
 
-                                    for (i,data) in enumerate(dF.data):
+                                    for (i, data) in enumerate(dF.data):
                                         data[...] = (rank+1)*(i+1)
                                         for displacements in all_displacements:
-                                            if sum(d!=0 for d in displacements) == 0:
+                                            if sum(d != 0 for d in displacements) == 0:
                                                 continue
-                                            (iview,ishape) = all_inner_ghost_slices[ndirections][directions][displacements]
-                                            (oview,oshape) = all_outer_ghost_slices[ndirections][directions][displacements]
+                                            (iview, ishape) = all_inner_ghost_slices[ndirections][directions][displacements]
+                                            (oview, oshape) = all_outer_ghost_slices[ndirections][directions][displacements]
                                             if (oshape is not None):
                                                 assert (ishape is not None)
                                                 data[oview] = ghost_vals(oshape, dtype, rank, directions,
-                                                                                            displacements, i)
+                                                                         displacements, i)
 
                                     dF.accumulate_ghosts(directions=directions,
                                                          exchange_method=exchange_method)
 
-                                    for (i,data) in enumerate(dF.data):
+                                    for (i, data) in enumerate(dF.data):
                                         for displacements in all_displacements:
-                                            ndisplacements = sum(d!=0 for d in displacements)
-                                            (iview,ishape) = all_inner_ghost_slices[ndirections][directions][displacements]
-                                            (oview,oshape) = all_outer_ghost_slices[ndirections][directions][displacements]
+                                            ndisplacements = sum(d != 0 for d in displacements)
+                                            (iview, ishape) = all_inner_ghost_slices[ndirections][directions][displacements]
+                                            (oview, oshape) = all_outer_ghost_slices[ndirections][directions][displacements]
 
                                             if (ishape is None):
                                                 assert (oshape is None)
@@ -803,17 +815,17 @@ def test_mpi_ghost_accumulate_periodic(comm=None):
                                             assert np.array_equal(data[oview].shape, oshape)
 
                                             overlaping_neighbours = set(tuple((np.asarray(mask, dtype=np.int32)*displacements).tolist()) for mask in masks)
-                                            overlaping_neighbours = filter(lambda x: 0<sum(_!=0 for _ in x)<=max_displacements, overlaping_neighbours) #diagonals
+                                            overlaping_neighbours = filter(lambda x: 0 < sum(_ != 0 for _ in x) <= max_displacements, overlaping_neighbours)  # diagonals
                                             overlaping_neighbours = tuple(overlaping_neighbours)
 
                                             expected_values = (rank+1)*(i+1)
                                             for disp in overlaping_neighbours:
                                                 ncoords = ()
-                                                j=0
+                                                j = 0
                                                 for _ in range(dim):
                                                     if _ in directions:
                                                         ci = (proc_shape[_]+proc_coords[_]+disp[j]) % proc_shape[_]
-                                                        j+=1
+                                                        j += 1
                                                     else:
                                                         ci = proc_coords[_]
                                                     ncoords += (ci,)
@@ -823,12 +835,13 @@ def test_mpi_ghost_accumulate_periodic(comm=None):
                                             idata = data[iview].get()
                                             odata = data[oview].get()
                                             assert np.allclose(idata, expected_values)
-                                            if (ndisplacements>0):
+                                            if (ndisplacements > 0):
                                                 assert np.allclose(odata, ghost_vals(oshape, dtype, rank, directions, displacements, i))
                                 continue
-                            if rank==0:
+                            if rank == 0:
                                 print()
 
+
 if __name__ == '__main__':
     from mpi4py import MPI
     comm = MPI.COMM_WORLD
@@ -837,7 +850,7 @@ if __name__ == '__main__':
     from hysop.tools.warning import disable_hysop_warnings
 
     with test_context():
-        if (size==1):
+        if (size == 1):
             test_serial_initialization_1d()
             test_serial_initialization_2d()
             if __ENABLE_LONG_TESTS__:
diff --git a/hysop/tools/field_utils.py b/hysop/tools/field_utils.py
index bc4c7c52a517e8aa647c1d8c418812eb1f23a7c9..042a67150aa0047bd90d7a113942802295f48cc6 100644
--- a/hysop/tools/field_utils.py
+++ b/hysop/tools/field_utils.py
@@ -4,10 +4,11 @@ from hysop.tools.types import first_not_None, to_tuple
 from hysop.tools.sympy_utils import nabla, partial, subscript, subscripts, \
     exponent, exponents, xsymbol, get_derivative_variables
 
+import sympy as sm
 from sympy.printing.str import StrPrinter, StrReprPrinter
 from sympy.printing.latex import LatexPrinter
 from packaging import version
-if version.parse(sp.__version__) > version.parse("1.7"):
+if version.parse(sm.__version__) > version.parse("1.7"):
     from sympy.printing.c import C99CodePrinter
 else:
     from sympy.printing.ccode import C99CodePrinter