diff --git a/sopht_mpi/simulator/flow/flow_simulators_mpi_2d.py b/sopht_mpi/simulator/flow/flow_simulators_mpi_2d.py
index d937916..131c654 100644
--- a/sopht_mpi/simulator/flow/flow_simulators_mpi_2d.py
+++ b/sopht_mpi/simulator/flow/flow_simulators_mpi_2d.py
@@ -99,8 +99,12 @@ def init_mpi(self):
             real_t=self.real_t,
             rank_distribution=self.rank_distribution,
         )
+        # Unbounded (i.e. non periodic) 2D flow solver with mpi4py-fft (only slabs
+        # decomp) does not need full exchange in ghost communicator.
         self.mpi_ghost_exchange_communicator = MPIGhostCommunicator2D(
-            ghost_size=self.ghost_size, mpi_construct=self.mpi_construct
+            ghost_size=self.ghost_size,
+            mpi_construct=self.mpi_construct,
+            full_exchange=False,
         )
 
     def init_domain(self):
diff --git a/sopht_mpi/utils/mpi_utils_2d.py b/sopht_mpi/utils/mpi_utils_2d.py
index 239a745..72cff87 100644
--- a/sopht_mpi/utils/mpi_utils_2d.py
+++ b/sopht_mpi/utils/mpi_utils_2d.py
@@ -17,7 +17,7 @@ def __init__(
         self,
         grid_size_y,
         grid_size_x,
-        periodic_flag=False,
+        periodic_domain=False,
         real_t=np.float64,
         rank_distribution=None,
     ):
@@ -27,6 +27,7 @@ def __init__(
         # Set the MPI dtype generator based on precision
         self.dtype_generator = MPI.FLOAT if real_t == np.float32 else MPI.DOUBLE
         # Setup MPI environment
+        self.periodic_domain = periodic_domain
         self.world = MPI.COMM_WORLD
         # Automatically create topologies
         if rank_distribution is None:
@@ -37,8 +38,10 @@ def __init__(
         else:
             self.rank_distribution = rank_distribution
         if 1 not in self.rank_distribution:
-            raise ValueError(
-                f"Rank distribution {self.rank_distribution} needs to be"
+            # Log warning here for more generic use.
+            # mpi4py-fft will take care of throwing errors if fft is invoked later.
+            logger.warning(
+                f"Rank distribution {self.rank_distribution} needs to be "
                 "aligned in at least one direction for fft"
             )
         self.grid_topology = MPI.Compute_dims(
@@ -60,7 +63,7 @@ def __init__(
 
         # Create Cartesian grid communicator
         self.grid = self.world.Create_cart(
-            self.grid_topology, periods=periodic_flag, reorder=False
+            self.grid_topology, periods=self.periodic_domain, reorder=False
         )
         # Determine neighbours in all directions
         self.previous_grid_along = np.zeros(self.grid_dim).astype(int)
@@ -84,13 +87,36 @@ def __init__(
 class MPIGhostCommunicator2D:
     """
     Class exclusive for ghost communication across ranks, initialises data types
-    that will be used for comm. in both blocking and non-blocking styles. Builds
-    dtypes based on ghost_size (determined from stencil width of the kernel)
-    This class wont be seen by the user, rather based on stencils we determine
-    the properties here.
+    for communication based on ghost_size (determined from stencil width of the kernel).
+
+    Here we refer the side and corners ghost cells following the terminologies for
+    polygons (https://en.wikipedia.org/wiki/Polygon) as below:
+    "The segments of a polygonal circuit are called its edges or sides. The points where
+    two edges meet are the polygon's vertices (singular: vertex) or corners."
+
+    v e e e e e v <- vertex (corner)
+    e x x x x x e
+    e x x x x x e <- edge (side)
+    e x x x x x e
+    v e e e e e v
+
+    (v) vertex (corner) ghost cell
+    (e) edge (side) ghost cell
+    (x) inner non-ghost cell
+
+    Note
+    ---
+    `full_exchange` allows for exchange both the edges (sides) and the vertex (corners)
+    ghost cells of the local domain (see illustration above). The full exchange mode is
+    required when Eulerian-to-Lagrangian (structured-to-unstructured) grid interpolation
+    is performed. In our 2D current unbounded flow solver, we won't need the
+    full exchange mode even when the interpolation is performed since mpi4py-fft only
+    for allows for slabs decomposition, and the corner cells corresponds to then the
+    corners of the unbounded domain. The full exchange mode may become useful if a
+    periodic domain flow solver is employed.
     """
 
-    def __init__(self, ghost_size, mpi_construct):
+    def __init__(self, ghost_size, mpi_construct, full_exchange=True):
         # extra width needed for kernel computation
         if ghost_size <= 0 and not isinstance(ghost_size, int):
             raise ValueError(
@@ -99,38 +125,83 @@ def __init__(self, ghost_size, mpi_construct):
             )
         self.ghost_size = ghost_size
         self.mpi_construct = mpi_construct
-        # define field_size variable for local field size (which includes ghost)
-        self.field_size = mpi_construct.local_grid_size + 2 * self.ghost_size
+        # exchange mode to include corners (forming full halo ring) or not
+        self.full_exchange = full_exchange
+        self.grid_coord = np.array(self.mpi_construct.grid.coords)
+
+        # Initialize data types
+        self.init_datatypes()
+
+        # Initialize requests list for non-blocking comm
+        self.comm_requests = []
 
-        # Set datatypes for ghost communication
-        # For row we use contiguous type
-        self.row_type = mpi_construct.dtype_generator.Create_contiguous(
-            count=self.ghost_size * self.field_size[1]
+        if self.full_exchange:
+            self.exchange_scalar_field_init = self.exchange_scalar_field_full_init
+        else:
+            self.exchange_scalar_field_init = self.exchange_scalar_field_edges_init
+
+    def init_datatypes(self):
+        """
+        Set datatypes for ghost communication
+        Note: this can be done more cleanly using Create_subarray, but using vector
+        approach here for illustration in a simpler 2d setting
+        """
+        # define field_size variable for local field size (without ghost)
+        self.field_size = self.mpi_construct.local_grid_size
+
+        # (1) For row data
+        self.row_type = self.mpi_construct.dtype_generator.Create_vector(
+            count=self.ghost_size,
+            blocklength=self.field_size[1],
+            stride=self.field_size[1] + 2 * self.ghost_size,
         )
         self.row_type.Commit()
-        # For column we use strided vector
-        self.column_type = mpi_construct.dtype_generator.Create_vector(
+
+        # (2) For column data
+        self.column_type = self.mpi_construct.dtype_generator.Create_vector(
             count=self.field_size[0],
             blocklength=self.ghost_size,
-            stride=self.field_size[1],
+            stride=self.field_size[1] + 2 * self.ghost_size,
         )
         self.column_type.Commit()
 
-        # Initialize requests list for non-blocking comm
-        self.comm_requests = []
+        # (3) For corner data (only if needed)
+        if self.full_exchange:
+            self.corner_type = self.mpi_construct.dtype_generator.Create_vector(
+                count=self.ghost_size,
+                blocklength=self.ghost_size,
+                stride=self.field_size[1] + 2 * self.ghost_size,
+            )
+            self.corner_type.Commit()
+
+    def _get_diagonally_shifted_coord_rank(self, coord_shift):
+        """Helper function to get diagonally shifted coords"""
+        shifted_coord = self.grid_coord + np.array(coord_shift)
+        if not self.mpi_construct.periodic_domain and (
+            np.any(shifted_coord >= self.mpi_construct.grid_topology)
+            or np.any(shifted_coord < 0)
+        ):
+            # The shifted coord is out of non-periodic domain
+            rank = MPI.PROC_NULL
+        else:
+            # Periodic domain is automatically taken care of in mpi cartersian grid
+            rank = self.mpi_construct.grid.Get_cart_rank(shifted_coord)
+        return rank
 
-    def exchange_scalar_field_init(self, local_field):
+    def exchange_scalar_field_edges_init(self, local_field):
         """
-        Exchange scalar field ghost data between neighbors.
+        Exchange scalar field ghost data on edges (sides) between neighbors.
         """
         # Lines below to make code more literal
         y_axis = 0
         x_axis = 1
+        ghost_rows_offset = self.ghost_size * local_field.shape[1]
+
         # Along Y: send to previous block, receive from next block
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
                 (
-                    local_field[self.ghost_size : 2 * self.ghost_size, :],
+                    local_field.ravel()[ghost_rows_offset + self.ghost_size :],
                     1,
                     self.row_type,
                 ),
@@ -140,7 +211,7 @@ def exchange_scalar_field_init(self, local_field):
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
                 (
-                    local_field[-self.ghost_size : local_field.shape[0], :],
+                    local_field.ravel()[-ghost_rows_offset + self.ghost_size :],
                     1,
                     self.row_type,
                 ),
@@ -152,7 +223,7 @@ def exchange_scalar_field_init(self, local_field):
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
                 (
-                    local_field[-2 * self.ghost_size : -self.ghost_size, :],
+                    local_field.ravel()[-2 * ghost_rows_offset + self.ghost_size :],
                     1,
                     self.row_type,
                 ),
@@ -162,7 +233,7 @@ def exchange_scalar_field_init(self, local_field):
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
                 (
-                    local_field[0 : self.ghost_size, :],
+                    local_field.ravel()[self.ghost_size :],
                     1,
                     self.row_type,
                 ),
@@ -174,7 +245,7 @@ def exchange_scalar_field_init(self, local_field):
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
                 (
-                    local_field.ravel()[self.ghost_size :],
+                    local_field.ravel()[ghost_rows_offset + self.ghost_size :],
                     1,
                     self.column_type,
                 ),
@@ -184,7 +255,9 @@ def exchange_scalar_field_init(self, local_field):
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
                 (
-                    local_field.ravel()[local_field.shape[1] - self.ghost_size :],
+                    local_field.ravel()[
+                        ghost_rows_offset + local_field.shape[1] - self.ghost_size :
+                    ],
                     1,
                     self.column_type,
                 ),
@@ -196,7 +269,9 @@ def exchange_scalar_field_init(self, local_field):
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
                 (
-                    local_field.ravel()[local_field.shape[1] - 2 * self.ghost_size :],
+                    local_field.ravel()[
+                        ghost_rows_offset + local_field.shape[1] - 2 * self.ghost_size :
+                    ],
                     1,
                     self.column_type,
                 ),
@@ -206,7 +281,7 @@ def exchange_scalar_field_init(self, local_field):
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
                 (
-                    local_field.ravel()[0:],
+                    local_field.ravel()[ghost_rows_offset:],
                     1,
                     self.column_type,
                 ),
@@ -214,6 +289,120 @@ def exchange_scalar_field_init(self, local_field):
             )
         )
 
+    def exchange_scalar_field_vertices_init(self, local_field):
+        """
+        Exchange scalar field ghost data on vertices (corners) between neighbors.
+        """
+        ghost_rows_offset = self.ghost_size * local_field.shape[1]
+        # Exchange corner ghost cells (four corners, so 4 exchanges here)
+        # (1) Update top right corner
+        # send bottom left corner to corresponding diagonal block (-1, -1),
+        # receive as top right corner from corresponding diagonal block (+1, +1)
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (
+                    local_field.ravel()[ghost_rows_offset + self.ghost_size :],
+                    1,
+                    self.corner_type,
+                ),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (
+                    local_field.ravel()[
+                        -(self.ghost_size - 1) * local_field.shape[1]
+                        - self.ghost_size :
+                    ],
+                    1,
+                    self.corner_type,
+                ),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 1]),
+            )
+        )
+        # (2) Update bottom left corner
+        # send top right corner to corresponding block (+1, +1),
+        # receive as bottom left corner from corresponding diagonal block (-1, -1)
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (
+                    local_field.ravel()[
+                        -(2 * self.ghost_size - 1) * local_field.shape[1]
+                        - 2 * self.ghost_size :
+                    ],
+                    1,
+                    self.corner_type,
+                ),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (
+                    local_field.ravel()[0:],
+                    1,
+                    self.corner_type,
+                ),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, -1]),
+            )
+        )
+        # (3) Update top left corner
+        # send bottom right corner to corresponding diagonal block (-1, +1),
+        # receive as top left corner from corresponding diagonal block (+1, -1)
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (
+                    local_field.ravel()[
+                        ghost_rows_offset + local_field.shape[1] - 2 * self.ghost_size :
+                    ],
+                    1,
+                    self.corner_type,
+                ),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (
+                    local_field.ravel()[-ghost_rows_offset:],
+                    1,
+                    self.corner_type,
+                ),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, -1]),
+            )
+        )
+        # (4) Update bottom right corner
+        # send top left corner to corresponding diagonal block (+1, -1),
+        # receive as bottom right corner from corresponding diagonal block (-1, +1)
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (
+                    local_field.ravel()[-2 * ghost_rows_offset + self.ghost_size :],
+                    1,
+                    self.corner_type,
+                ),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (
+                    local_field.ravel()[local_field.shape[1] - self.ghost_size :],
+                    1,
+                    self.corner_type,
+                ),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 1]),
+            )
+        )
+
+    def exchange_scalar_field_full_init(self, local_field):
+        """
+        Exchange scalar field ghost data (a full halo ring) between neighbors.
+        """
+        self.exchange_scalar_field_edges_init(local_field)
+        self.exchange_scalar_field_vertices_init(local_field)
+
     def exchange_vector_field_init(self, local_vector_field):
         self.exchange_scalar_field_init(
             local_field=local_vector_field[VectorField.x_axis_idx()]
diff --git a/sopht_mpi/utils/mpi_utils_3d.py b/sopht_mpi/utils/mpi_utils_3d.py
index 02a7637..5c75e47 100644
--- a/sopht_mpi/utils/mpi_utils_3d.py
+++ b/sopht_mpi/utils/mpi_utils_3d.py
@@ -16,7 +16,7 @@ def __init__(
         grid_size_z,
         grid_size_y,
         grid_size_x,
-        periodic_flag=False,
+        periodic_domain=False,
         real_t=np.float64,
         rank_distribution=None,
     ):
@@ -26,6 +26,7 @@ def __init__(
         # Set the MPI dtype generator based on precision
         self.dtype_generator = MPI.FLOAT if real_t == np.float32 else MPI.DOUBLE
         # Setup MPI environment
+        self.periodic_domain = periodic_domain
         self.world = MPI.COMM_WORLD
         # Automatically create topologies
         if rank_distribution is None:
@@ -36,8 +37,10 @@ def __init__(
         else:
             self.rank_distribution = rank_distribution
         if 1 not in self.rank_distribution:
-            raise ValueError(
-                f"Rank distribution {self.rank_distribution} needs to be"
+            # Log warning here for more generic use.
+            # mpi4py-fft will take care of throwing errors if fft is invoked later.
+            logger.warning(
+                f"Rank distribution {self.rank_distribution} needs to be "
                 "aligned in at least one direction for fft"
             )
         self.grid_topology = MPI.Compute_dims(
@@ -59,7 +62,7 @@ def __init__(
 
         # Create Cartesian grid communicator
         self.grid = self.world.Create_cart(
-            self.grid_topology, periods=periodic_flag, reorder=False
+            self.grid_topology, periods=self.periodic_domain, reorder=False
         )
         # Determine neighbours in all directions
         self.previous_grid_along = np.zeros(self.grid_dim).astype(int)
@@ -83,13 +86,36 @@ def __init__(
 class MPIGhostCommunicator3D:
     """
     Class exclusive for ghost communication across ranks, initialises data types
-    that will be used for comm. in both blocking and non-blocking styles. Builds
-    dtypes based on ghost_size (determined from stencil width of the kernel)
-    This class wont be seen by the user, rather based on stencils we determine
-    the properties here.
+    for communication based on ghost_size (determined from stencil width of the kernel).
+
+    Here we refer the ghost cells following the terminologies for cubes
+    (https://en.wikipedia.org/wiki/Cube): "The has 6 faces, 12 edges, and 8 vertices."
+    Consider a face of the cube in the ghost zone as illustrated below.
+
+    v e e e e e v <- vertex (corner)
+    e f f f f f e
+    e f f f f f e <- edge (side)
+    e f f f f f e
+    v e e e e e v
+
+    (v) vertex ghost cell
+    (e) edge ghost cell
+    (f) face ghost cell
+
+    Note
+    ---
+    `full_exchange` allows for exchanges ghost cells on the vertices, edges and faces of
+    the local domain (see illustration above). The full exchange mode is required only
+    when Eulerian-to-Lagrangian (structured-to-unstructured) grid interpolation is
+    performed. Similar to the 2D ghost communicator, we won't need the full exchange
+    mode even when the interpolation is performed if the domain decomposition mode
+    is slabs (unless a periodic flow simulator is employed). In the case where domain
+    decomposition mode is pencils (in 3D, mpi4py-fft allows for both slabs and pencils),
+    full exchange is needed if the structured-to-unstructured grid interpolation is
+    performed.
     """
 
-    def __init__(self, ghost_size, mpi_construct):
+    def __init__(self, ghost_size, mpi_construct, full_exchange=True):
         # extra width needed for kernel computation
         if ghost_size <= 0 and not isinstance(ghost_size, int):
             raise ValueError(
@@ -98,188 +124,1003 @@ def __init__(self, ghost_size, mpi_construct):
             )
         self.ghost_size = ghost_size
         self.mpi_construct = mpi_construct
+        self.full_exchange = full_exchange
+        self.grid_coord = np.array(self.mpi_construct.grid.coords)
+
+        # Initialize data types
+        self.init_datatypes()
+
+        # Initialize requests list for non-blocking comm
+        self.comm_requests = []
+
+        if self.full_exchange:
+            self.exchange_scalar_field_init = self.exchange_scalar_field_full_init
+        else:
+            self.exchange_scalar_field_init = self.exchange_scalar_field_faces_init
+
+    def init_datatypes(self):
+        """
+        Set datatypes for ghost communication on 6 faces, 12 edges and 8 vertices.
+        Each datatype will have a send and recv, giving a total of
+        (6 + 12 + 8) * 2 = 52 datatypes.
+
+        Note
+        ---
+        Here we are using the Create_subarray approach, each type for send/recv needs to
+        be initialized based on their starting index location. Alternatively, we can use
+        Create_vector approach which is more involved, but perhaps faster way.
+        Keeping subarray approach for now for its readibility and easy implementation.
+
+        For simplicity and readability, we use coordinate system (z, y, x) as used in
+        grid topology to define direction for send and recv. For example, when sending
+        along the (0, -1, +1) direction (i.e. along XY-plane in the positive X and
+        negative Y direction), we denote the data type name as `send_along_0_ny_px_type`
+        where 0, ny and px denotes 0 in Z, negative in Y, and positive in X coord shift,
+        respectively.
+        """
         # define field_size variable for local field size (which includes ghost)
-        self.field_size = mpi_construct.local_grid_size + 2 * self.ghost_size
-
-        # Set datatypes for ghost communication
-        # Note: these can be written in a more involved, but perhaps faster way.
-        # Keeping this for now for its readibility and easy implementation.
-        # Using the Create_subarray approach, each type for sending / receiving
-        # needs to be initialized based on their starting index location. In
-        # each dimension, we have 2 ghost layers to be sent (to next & prev) and
-        # 2 corresponding receiving layers (from next & prev). This amounts to
-        # (2 type for send/recv) * (2 dir along each dim) * (3 dim) = 12 type
-        # Along X (next)
-        self.send_next_along_x_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.field_size_with_ghost = (
+            self.mpi_construct.local_grid_size + 2 * self.ghost_size
+        )
+        self.field_size = self.mpi_construct.local_grid_size
+
+        # (1) For faces data
+        # Since there are 6 faces on a cube, each with a send and recv datatypes, we
+        # need 6 * 2 = 12 data types
+        # Comm. along +X direction: send to (0, 0, +1) with recv from (0, 0, -1) block
+        self.send_to_0_0_px_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.field_size[0], self.field_size[1], self.ghost_size],
-            starts=[0, 0, self.field_size[2] - 2 * self.ghost_size],
+            starts=[
+                self.ghost_size,
+                self.ghost_size,
+                self.field_size_with_ghost[2] - 2 * self.ghost_size,
+            ],
         )
-        self.send_next_along_x_type.Commit()
-        self.recv_next_along_x_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.send_to_0_0_px_type.Commit()
+        self.recv_from_0_0_nx_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.field_size[0], self.field_size[1], self.ghost_size],
-            starts=[0, 0, self.field_size[2] - self.ghost_size],
+            starts=[self.ghost_size, self.ghost_size, 0],
         )
-        self.recv_next_along_x_type.Commit()
-        # Along X (prev)
-        self.send_previous_along_x_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.recv_from_0_0_nx_type.Commit()
+        # Comm. along -X direction: send to (-1, 0, 0) with recv from (+1, 0, 0) block
+        self.send_to_0_0_nx_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.field_size[0], self.field_size[1], self.ghost_size],
-            starts=[0, 0, self.ghost_size],
+            starts=[self.ghost_size, self.ghost_size, self.ghost_size],
         )
-        self.send_previous_along_x_type.Commit()
-        self.recv_previous_along_x_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.send_to_0_0_nx_type.Commit()
+        self.recv_from_0_0_px_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.field_size[0], self.field_size[1], self.ghost_size],
-            starts=[0, 0, 0],
+            starts=[
+                self.ghost_size,
+                self.ghost_size,
+                self.field_size_with_ghost[2] - self.ghost_size,
+            ],
         )
-        self.recv_previous_along_x_type.Commit()
-        # Along Y (next)
-        self.send_next_along_y_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.recv_from_0_0_px_type.Commit()
+
+        # Comm. along +Y direction: send to (0, +1, 0) with recv from (0, -1, 0) block
+        self.send_to_0_py_0_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.field_size[0], self.ghost_size, self.field_size[2]],
-            starts=[0, self.field_size[1] - 2 * self.ghost_size, 0],
+            starts=[
+                self.ghost_size,
+                self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                self.ghost_size,
+            ],
         )
-        self.send_next_along_y_type.Commit()
-        self.recv_next_along_y_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.send_to_0_py_0_type.Commit()
+        self.recv_from_0_ny_0_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.field_size[0], self.ghost_size, self.field_size[2]],
-            starts=[0, self.field_size[1] - self.ghost_size, 0],
+            starts=[self.ghost_size, 0, self.ghost_size],
         )
-        self.recv_next_along_y_type.Commit()
-        # Along Y (prev)
-        self.send_previous_along_y_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.recv_from_0_ny_0_type.Commit()
+        # Comm. along -Y direction: send to (0, -1, 0) with recv from (0, +1, 0) block
+        self.send_to_0_ny_0_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.field_size[0], self.ghost_size, self.field_size[2]],
-            starts=[0, self.ghost_size, 0],
+            starts=[self.ghost_size, self.ghost_size, self.ghost_size],
         )
-        self.send_previous_along_y_type.Commit()
-        self.recv_previous_along_y_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.send_to_0_ny_0_type.Commit()
+        self.recv_from_0_py_0_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.field_size[0], self.ghost_size, self.field_size[2]],
-            starts=[0, 0, 0],
+            starts=[
+                self.ghost_size,
+                self.field_size_with_ghost[1] - self.ghost_size,
+                self.ghost_size,
+            ],
         )
-        self.recv_previous_along_y_type.Commit()
+        self.recv_from_0_py_0_type.Commit()
 
-        # Along Z (next)
-        self.send_next_along_z_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        # Comm. along +Z direction: send to (+1, 0, 0) with recv from (-1, 0, 0) block
+        self.send_to_pz_0_0_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.ghost_size, self.field_size[1], self.field_size[2]],
-            starts=[self.field_size[0] - 2 * self.ghost_size, 0, 0],
+            starts=[
+                self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                self.ghost_size,
+                self.ghost_size,
+            ],
         )
-        self.send_next_along_z_type.Commit()
-        self.recv_next_along_z_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.send_to_pz_0_0_type.Commit()
+        self.recv_from_nz_0_0_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.ghost_size, self.field_size[1], self.field_size[2]],
-            starts=[self.field_size[0] - self.ghost_size, 0, 0],
+            starts=[0, self.ghost_size, self.ghost_size],
         )
-        self.recv_next_along_z_type.Commit()
-        # Along Z (prev)
-        self.send_previous_along_z_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.recv_from_nz_0_0_type.Commit()
+        # Comm. along -Z direction: send to (-1, 0, 0) with recv from (+1, 0, 0) block
+        self.send_to_nz_0_0_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.ghost_size, self.field_size[1], self.field_size[2]],
-            starts=[self.ghost_size, 0, 0],
+            starts=[self.ghost_size, self.ghost_size, self.ghost_size],
         )
-        self.send_previous_along_z_type.Commit()
-        self.recv_previous_along_z_type = mpi_construct.dtype_generator.Create_subarray(
-            sizes=self.field_size,
+        self.send_to_nz_0_0_type.Commit()
+        self.recv_from_pz_0_0_type = self.mpi_construct.dtype_generator.Create_subarray(
+            sizes=self.field_size_with_ghost,
             subsizes=[self.ghost_size, self.field_size[1], self.field_size[2]],
-            starts=[0, 0, 0],
+            starts=[
+                self.field_size_with_ghost[0] - self.ghost_size,
+                self.ghost_size,
+                self.ghost_size,
+            ],
         )
-        self.recv_previous_along_z_type.Commit()
+        self.recv_from_pz_0_0_type.Commit()
 
-        # Initialize requests list for non-blocking comm
-        self.comm_requests = []
+        if self.full_exchange:
+            # (2) For edges
+            # Since there are 12 edges on a cube, each with a send and recv datatypes,
+            # we need 12 * 2 = 24 datatypes
+            # Comm. along +Y, +X: send to (0, +1, +1) with recv from (0, -1, -1) block
+            self.send_to_0_py_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.field_size[0], self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                        self.field_size_with_ghost[2] - 2 * self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_0_py_px_type.Commit()
+            self.recv_from_0_ny_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.field_size[0], self.ghost_size, self.ghost_size],
+                    starts=[self.ghost_size, 0, 0],
+                )
+            )
+            self.recv_from_0_ny_nx_type.Commit()
+            # Comm. along -Y, +X: send to (0, -1, +1) with recv from (0, +1, -1) block
+            self.send_to_0_ny_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.field_size[0], self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.ghost_size,
+                        self.field_size_with_ghost[2] - 2 * self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_0_ny_px_type.Commit()
+            self.recv_from_0_py_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.field_size[0], self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.field_size_with_ghost[1] - self.ghost_size,
+                        0,
+                    ],
+                )
+            )
+            self.recv_from_0_py_nx_type.Commit()
+            # Comm. along +Y, -X: send to (0, +1, -1) with recv from (0, -1, +1) block
+            self.send_to_0_py_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.field_size[0], self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_0_py_nx_type.Commit()
+            self.recv_from_0_ny_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.field_size[0], self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        0,
+                        self.field_size_with_ghost[2] - self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_0_ny_px_type.Commit()
+            # Comm. along -Y, -X: send to (0, -1, -1) with recv from (0, +1, +1) block
+            self.send_to_0_ny_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.field_size[0], self.ghost_size, self.ghost_size],
+                    starts=[self.ghost_size, self.ghost_size, self.ghost_size],
+                )
+            )
+            self.send_to_0_ny_nx_type.Commit()
+            self.recv_from_0_py_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.field_size[0], self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.field_size_with_ghost[1] - self.ghost_size,
+                        self.field_size_with_ghost[2] - self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_0_py_px_type.Commit()
+
+            # Comm. along +Z, +X: send to (+1, 0, +1) with recv from (-1, 0, -1) block
+            self.send_to_pz_0_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.field_size[1], self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                        self.ghost_size,
+                        self.field_size_with_ghost[2] - 2 * self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_pz_0_px_type.Commit()
+            self.recv_from_nz_0_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.field_size[1], self.ghost_size],
+                    starts=[0, self.ghost_size, 0],
+                )
+            )
+            self.recv_from_nz_0_nx_type.Commit()
+            # Comm. along -Z, +X: send to (-1, 0, +1) with recv from (+1, 0, -1) block
+            self.send_to_nz_0_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.field_size[1], self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.ghost_size,
+                        self.field_size_with_ghost[2] - 2 * self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_nz_0_px_type.Commit()
+            self.recv_from_pz_0_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.field_size[1], self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - self.ghost_size,
+                        self.ghost_size,
+                        0,
+                    ],
+                )
+            )
+            self.recv_from_pz_0_nx_type.Commit()
+            # Comm. along +Z, -X: send to (+1, 0, -1) with recv from (-1, 0, +1) block
+            self.send_to_pz_0_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.field_size[1], self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                        self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_pz_0_nx_type.Commit()
+            self.recv_from_nz_0_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.field_size[1], self.ghost_size],
+                    starts=[
+                        0,
+                        self.ghost_size,
+                        self.field_size_with_ghost[2] - self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_nz_0_px_type.Commit()
+            # Comm. along -Z, -X: send to (-1, 0, -1) with recv from (+1, 0, +1) block
+            self.send_to_nz_0_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.field_size[1], self.ghost_size],
+                    starts=[self.ghost_size, self.ghost_size, self.ghost_size],
+                )
+            )
+            self.send_to_nz_0_nx_type.Commit()
+            self.recv_from_pz_0_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.field_size[1], self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - self.ghost_size,
+                        self.ghost_size,
+                        self.field_size_with_ghost[2] - self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_pz_0_px_type.Commit()
+
+            # Comm. along +Z, +Y: send to (+1, +1, 0) with recv from (-1, -1, 0) block
+            self.send_to_pz_py_0_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.field_size[2]],
+                    starts=[
+                        self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                        self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_pz_py_0_type.Commit()
+            self.recv_from_nz_ny_0_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.field_size[2]],
+                    starts=[0, 0, self.ghost_size],
+                )
+            )
+            self.recv_from_nz_ny_0_type.Commit()
+            # Comm. along -Z, +Y: send to (-1, +1, 0) with recv from (+1, -1, 0) block
+            self.send_to_nz_py_0_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.field_size[2]],
+                    starts=[
+                        self.ghost_size,
+                        self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_nz_py_0_type.Commit()
+            self.recv_from_pz_ny_0_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.field_size[2]],
+                    starts=[
+                        self.field_size_with_ghost[0] - self.ghost_size,
+                        0,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_pz_ny_0_type.Commit()
+            # Comm. along +Z, -Y: send to (+1, -1, 0) with recv from (-1, +1, 0) block
+            self.send_to_pz_ny_0_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.field_size[2]],
+                    starts=[
+                        self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                        self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_pz_ny_0_type.Commit()
+            self.recv_from_nz_py_0_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.field_size[2]],
+                    starts=[
+                        0,
+                        self.field_size_with_ghost[1] - self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_nz_py_0_type.Commit()
+            # Comm. along -Z, -Y: send to (-1, -1, 0) with recv from (+1, +1, 0) block
+            self.send_to_nz_ny_0_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.field_size[2]],
+                    starts=[self.ghost_size, self.ghost_size, self.ghost_size],
+                )
+            )
+            self.send_to_nz_ny_0_type.Commit()
+            self.recv_from_pz_py_0_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.field_size[2]],
+                    starts=[
+                        self.field_size_with_ghost[0] - self.ghost_size,
+                        self.field_size_with_ghost[1] - self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_pz_py_0_type.Commit()
+
+            # (3) For vertices
+            # Since there are 8 vertices on a cube, each with a send and recv datatypes,
+            # we need 8 * 2 = 16 datatypes
+            # Comm. along +Z, +Y, +X: send to (+1, +1, +1) with recv from (-1, -1, -1) block
+            self.send_to_pz_py_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                        self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                        self.field_size_with_ghost[2] - 2 * self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_pz_py_px_type.Commit()
+            self.recv_from_nz_ny_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[0, 0, 0],
+                )
+            )
+            self.recv_from_nz_ny_nx_type.Commit()
+            # Comm. along -Z, +Y, +X: send to (-1, +1, +1) with recv from (+1, -1, -1) block
+            self.send_to_nz_py_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                        self.field_size_with_ghost[2] - 2 * self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_nz_py_px_type.Commit()
+            self.recv_from_pz_ny_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[self.field_size_with_ghost[0] - self.ghost_size, 0, 0],
+                )
+            )
+            self.recv_from_pz_ny_nx_type.Commit()
+            # Comm. along +Z, -Y, +X: send to (+1, -1, +1) with recv from (-1, +1, -1) block
+            self.send_to_pz_ny_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                        self.ghost_size,
+                        self.field_size_with_ghost[2] - 2 * self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_pz_ny_px_type.Commit()
+            self.recv_from_nz_py_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[0, self.field_size_with_ghost[1] - self.ghost_size, 0],
+                )
+            )
+            self.recv_from_nz_py_nx_type.Commit()
+            # Comm. along +Z, +Y, -X: send to (+1, +1, -1) with recv from (-1, -1, +1) block
+            self.send_to_pz_py_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                        self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_pz_py_nx_type.Commit()
+            self.recv_from_nz_ny_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[0, 0, self.field_size_with_ghost[2] - self.ghost_size],
+                )
+            )
+            self.recv_from_nz_ny_px_type.Commit()
+            # Comm. along -Z, -Y, +X: send to (-1, -1, +1) with recv from (+1, +1, -1) block
+            self.send_to_nz_ny_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.ghost_size,
+                        self.field_size_with_ghost[2] - 2 * self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_nz_ny_px_type.Commit()
+            self.recv_from_pz_py_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - self.ghost_size,
+                        self.field_size_with_ghost[1] - self.ghost_size,
+                        0,
+                    ],
+                )
+            )
+            self.recv_from_pz_py_nx_type.Commit()
+            # Comm. along -Z, +Y, -X: send to (-1, +1, -1) with recv from (+1, -1, +1) block
+            self.send_to_nz_py_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.ghost_size,
+                        self.field_size_with_ghost[1] - 2 * self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_nz_py_nx_type.Commit()
+            self.recv_from_pz_ny_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - self.ghost_size,
+                        0,
+                        self.field_size_with_ghost[2] - self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_pz_ny_px_type.Commit()
+            # Comm. along +Z, -Y, -X: send to (+1, -1, -1) with recv from (-1, +1, +1) block
+            self.send_to_pz_ny_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - 2 * self.ghost_size,
+                        self.ghost_size,
+                        self.ghost_size,
+                    ],
+                )
+            )
+            self.send_to_pz_ny_nx_type.Commit()
+            self.recv_from_nz_py_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        0,
+                        self.field_size_with_ghost[1] - self.ghost_size,
+                        self.field_size_with_ghost[2] - self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_nz_py_px_type.Commit()
+            # Comm. along -Z, -Y, -X: send to (-1, -1, -1) with recv from (+1, +1, +1) block
+            self.send_to_nz_ny_nx_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[self.ghost_size, self.ghost_size, self.ghost_size],
+                )
+            )
+            self.send_to_nz_ny_nx_type.Commit()
+            self.recv_from_pz_py_px_type = (
+                self.mpi_construct.dtype_generator.Create_subarray(
+                    sizes=self.field_size_with_ghost,
+                    subsizes=[self.ghost_size, self.ghost_size, self.ghost_size],
+                    starts=[
+                        self.field_size_with_ghost[0] - self.ghost_size,
+                        self.field_size_with_ghost[1] - self.ghost_size,
+                        self.field_size_with_ghost[2] - self.ghost_size,
+                    ],
+                )
+            )
+            self.recv_from_pz_py_px_type.Commit()
 
-    def exchange_scalar_field_init(self, local_field):
+    def _get_diagonally_shifted_coord_rank(self, coord_shift):
+        """Helper function to get diagonally shifted coords"""
+        shifted_coord = self.grid_coord + np.array(coord_shift)
+        if not self.mpi_construct.periodic_domain and (
+            np.any(shifted_coord >= self.mpi_construct.grid_topology)
+            or np.any(shifted_coord < 0)
+        ):
+            # The shifted coord is out of non-periodic domain
+            rank = MPI.PROC_NULL
+        else:
+            # Periodic domain is automatically taken care of in mpi cartersian grid
+            rank = self.mpi_construct.grid.Get_cart_rank(shifted_coord)
+        return rank
+
+    def exchange_scalar_field_faces_init(self, local_field):
         """
-        Exchange ghost data between neighbors.
+        Exchange scalar field ghost data on faces between neighbors.
         """
         # Lines below to make code more literal
         z_axis = 0
         y_axis = 1
         x_axis = 2
-        # Along X: send to previous block, receive from next block
+
+        # Comm. along +X direction: send to (0, 0, +1) with recv from (0, 0, -1) block
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
-                (local_field, self.send_previous_along_x_type),
+                (local_field, self.send_to_0_0_px_type),
+                dest=self.mpi_construct.next_grid_along[x_axis],
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_0_0_nx_type),
+                source=self.mpi_construct.previous_grid_along[x_axis],
+            )
+        )
+        # Comm. along -X direction: send to (0, 0, -1) with recv from (0, 0, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_0_0_nx_type),
                 dest=self.mpi_construct.previous_grid_along[x_axis],
             )
         )
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
-                (local_field, self.recv_next_along_x_type),
+                (local_field, self.recv_from_0_0_px_type),
                 source=self.mpi_construct.next_grid_along[x_axis],
             )
         )
-        # Along X: send to next block, receive from previous block
+        # Comm. along +Y direction: send to (0, +1, 0) with recv from (0, -1, 0) block
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
-                (local_field, self.send_next_along_x_type),
-                dest=self.mpi_construct.next_grid_along[x_axis],
+                (local_field, self.send_to_0_py_0_type),
+                dest=self.mpi_construct.next_grid_along[y_axis],
             )
         )
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
-                (local_field, self.recv_previous_along_x_type),
-                source=self.mpi_construct.previous_grid_along[x_axis],
+                (local_field, self.recv_from_0_ny_0_type),
+                source=self.mpi_construct.previous_grid_along[y_axis],
             )
         )
-
-        # Along Y: send to previous block, receive from next block
+        # Comm. along -Y direction: send to (0, -1, 0) with recv from (0, +1, 0) block
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
-                (local_field, self.send_previous_along_y_type),
+                (local_field, self.send_to_0_ny_0_type),
                 dest=self.mpi_construct.previous_grid_along[y_axis],
             )
         )
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
-                (local_field, self.recv_next_along_y_type),
+                (local_field, self.recv_from_0_py_0_type),
                 source=self.mpi_construct.next_grid_along[y_axis],
             )
         )
-        # Along Y: send to next block, receive from previous block
+        # Comm. along +Z direction: send to (+1, 0, 0) with recv from (-1, 0, 0) block
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
-                (local_field, self.send_next_along_y_type),
-                dest=self.mpi_construct.next_grid_along[y_axis],
+                (local_field, self.send_to_pz_0_0_type),
+                dest=self.mpi_construct.next_grid_along[z_axis],
             )
         )
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
-                (local_field, self.recv_previous_along_y_type),
-                source=self.mpi_construct.previous_grid_along[y_axis],
+                (local_field, self.recv_from_nz_0_0_type),
+                source=self.mpi_construct.previous_grid_along[z_axis],
             )
         )
-
-        # Along Z: send to previous block, receive from next block
+        # Comm. along -Z direction: send to (-1, 0, 0) with recv from (+1, 0, 0) block
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
-                (local_field, self.send_previous_along_z_type),
+                (local_field, self.send_to_nz_0_0_type),
                 dest=self.mpi_construct.previous_grid_along[z_axis],
             )
         )
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
-                (local_field, self.recv_next_along_z_type),
+                (local_field, self.recv_from_pz_0_0_type),
                 source=self.mpi_construct.next_grid_along[z_axis],
             )
         )
-        # Along Z: send to next block, receive from previous block
+
+    def exchange_scalar_field_edges_init(self, local_field):
+        """
+        Exchange scalar field ghost data on edges between neighbors.
+        """
+        # Comm. along +Y, +X: send to (0, +1, +1) with recv from (0, -1, -1) block
         self.comm_requests.append(
             self.mpi_construct.grid.Isend(
-                (local_field, self.send_next_along_z_type),
-                dest=self.mpi_construct.next_grid_along[z_axis],
+                (local_field, self.send_to_0_py_px_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[0, 1, 1]),
             )
         )
         self.comm_requests.append(
             self.mpi_construct.grid.Irecv(
-                (local_field, self.recv_previous_along_z_type),
-                source=self.mpi_construct.previous_grid_along[z_axis],
+                (local_field, self.recv_from_0_ny_nx_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[0, -1, -1]),
+            )
+        )
+        # Comm. along -Y, +X: send to (0, -1, +1) with recv from (0, +1, -1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_0_ny_px_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[0, -1, 1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_0_py_nx_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[0, 1, -1]),
+            )
+        )
+        # Comm. along +Y, -X: send to (0, +1, -1) with recv from (0, -1, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_0_py_nx_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[0, 1, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_0_ny_px_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[0, -1, 1]),
+            )
+        )
+        # Comm. along -Y, -X: send to (0, -1, -1) with recv from (0, +1, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_0_ny_nx_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[0, -1, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_0_py_px_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[0, 1, 1]),
+            )
+        )
+        # Comm. along +Z, +X: send to (+1, 0, +1) with recv from (-1, 0, -1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_pz_0_px_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 0, 1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_nz_0_nx_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 0, -1]),
+            )
+        )
+        # Comm. along -Z, +X: send to (-1, 0, +1) with recv from (+1, 0, -1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_nz_0_px_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 0, 1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_pz_0_nx_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 0, -1]),
+            )
+        )
+        # Comm. along +Z, -X: send to (+1, 0, -1) with recv from (-1, 0, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_pz_0_nx_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 0, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_nz_0_px_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 0, 1]),
+            )
+        )
+        # Comm. along -Z, -X: send to (-1, 0, -1) with recv from (+1, 0, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_nz_0_nx_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 0, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_pz_0_px_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 0, 1]),
+            )
+        )
+        # Comm. along +Z, +Y: send to (+1, +1, 0) with recv from (-1, -1, 0) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_pz_py_0_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 1, 0]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_nz_ny_0_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, -1, 0]),
+            )
+        )
+        # Comm. along -Z, +Y: send to (-1, +1, 0) with recv from (+1, -1, 0) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_nz_py_0_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 1, 0]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_pz_ny_0_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, -1, 0]),
+            )
+        )
+        # Comm. along +Z, -Y: send to (+1, -1, 0) with recv from (-1, +1, 0) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_pz_ny_0_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, -1, 0]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_nz_py_0_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 1, 0]),
+            )
+        )
+        # Comm. along -Z, -Y: send to (-1, -1, 0) with recv from (+1, +1, 0) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_nz_ny_0_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, -1, 0]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_pz_py_0_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 1, 0]),
+            )
+        )
+
+    def exchange_scalar_field_vertices_init(self, local_field):
+        """
+        Exchange scalar field ghost data on vertices between neighbors.
+        """
+        # Comm. along +Z, +Y, +X: send to (+1, +1, +1) with recv from (-1, -1, -1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_pz_py_px_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 1, 1]),
             )
         )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_nz_ny_nx_type),
+                source=self._get_diagonally_shifted_coord_rank(
+                    coord_shift=[-1, -1, -1]
+                ),
+            )
+        )
+        # Comm. along -Z, +Y, +X: send to (-1, +1, +1) with recv from (+1, -1, -1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_nz_py_px_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 1, 1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_pz_ny_nx_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, -1, -1]),
+            )
+        )
+        # Comm. along +Z, -Y, +X: send to (+1, -1, +1) with recv from (-1, +1, -1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_pz_ny_px_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, -1, 1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_nz_py_nx_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 1, -1]),
+            )
+        )
+        # Comm. along +Z, +Y, -X: send to (+1, +1, -1) with recv from (-1, -1, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_pz_py_nx_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 1, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_nz_ny_px_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, -1, 1]),
+            )
+        )
+        # Comm. along -Z, -Y, +X: send to (-1, -1, +1) with recv from (+1, +1, -1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_nz_ny_px_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, -1, 1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_pz_py_nx_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 1, -1]),
+            )
+        )
+        # Comm. along -Z, +Y, -X: send to (-1, +1, -1) with recv from (+1, -1, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_nz_py_nx_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 1, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_pz_ny_px_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, -1, 1]),
+            )
+        )
+        # Comm. along +Z, -Y, -X: send to (+1, -1, -1) with recv from (-1, +1, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_pz_ny_nx_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[1, -1, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_nz_py_px_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, 1, 1]),
+            )
+        )
+        # Comm. along -Z, -Y, -X: send to (-1, -1, -1) with recv from (+1, +1, +1) block
+        self.comm_requests.append(
+            self.mpi_construct.grid.Isend(
+                (local_field, self.send_to_nz_ny_nx_type),
+                dest=self._get_diagonally_shifted_coord_rank(coord_shift=[-1, -1, -1]),
+            )
+        )
+        self.comm_requests.append(
+            self.mpi_construct.grid.Irecv(
+                (local_field, self.recv_from_pz_py_px_type),
+                source=self._get_diagonally_shifted_coord_rank(coord_shift=[1, 1, 1]),
+            )
+        )
+
+    def exchange_scalar_field_full_init(self, local_field):
+        """
+        Exchange scalar field ghost data including all faces, edges and vertices between
+        neighbors.
+        """
+        self.exchange_scalar_field_faces_init(local_field)
+        self.exchange_scalar_field_edges_init(local_field)
+        self.exchange_scalar_field_vertices_init(local_field)
 
     def exchange_vector_field_init(self, local_vector_field):
         self.exchange_scalar_field_init(
diff --git a/tests/test_utils/test_mpi_utils_2d.py b/tests/test_utils/test_mpi_utils_2d.py
index 050b38c..288868c 100644
--- a/tests/test_utils/test_mpi_utils_2d.py
+++ b/tests/test_utils/test_mpi_utils_2d.py
@@ -12,7 +12,7 @@
 
 
 @pytest.mark.mpi(group="MPI_utils", min_size=4)
-@pytest.mark.parametrize("ghost_size", [1, 2, 3])
+@pytest.mark.parametrize("ghost_size", [1, 2])
 @pytest.mark.parametrize("precision", ["single", "double"])
 @pytest.mark.parametrize("rank_distribution", [(1, 0), (0, 1)])
 @pytest.mark.parametrize("aspect_ratio", [(1, 1), (1, 2), (2, 1)])
@@ -20,7 +20,7 @@
 def test_mpi_field_gather_scatter(
     ghost_size, precision, rank_distribution, aspect_ratio, master_rank
 ):
-    n_values = 32
+    n_values = 8
     real_t = get_real_t(precision)
     mpi_construct = MPIConstruct2D(
         grid_size_y=n_values * aspect_ratio[0],
@@ -79,25 +79,32 @@ def test_mpi_field_gather_scatter(
 
 
 @pytest.mark.mpi(group="MPI_utils", min_size=4)
-@pytest.mark.parametrize("ghost_size", [1, 2, 3])
+@pytest.mark.parametrize("ghost_size", [1, 2])
 @pytest.mark.parametrize("precision", ["single", "double"])
-@pytest.mark.parametrize("rank_distribution", [(1, 0), (0, 1)])
+@pytest.mark.parametrize("rank_distribution", [(1, 0), (0, 1), (2, 2)])
 @pytest.mark.parametrize("aspect_ratio", [(1, 1), (1, 2), (2, 1)])
+@pytest.mark.parametrize("full_exchange", [True, False])
 def test_mpi_ghost_communication(
-    ghost_size, precision, rank_distribution, aspect_ratio
+    ghost_size,
+    precision,
+    rank_distribution,
+    aspect_ratio,
+    full_exchange,
 ):
-    n_values = 32
+    n_values = 8
     real_t = get_real_t(precision)
     mpi_construct = MPIConstruct2D(
         grid_size_y=n_values * aspect_ratio[0],
         grid_size_x=n_values * aspect_ratio[1],
-        periodic_flag=True,
+        periodic_domain=True,
         real_t=real_t,
         rank_distribution=rank_distribution,
     )
     # extra width needed for kernel computation
     mpi_ghost_exchange_communicator = MPIGhostCommunicator2D(
-        ghost_size=ghost_size, mpi_construct=mpi_construct
+        ghost_size=ghost_size,
+        mpi_construct=mpi_construct,
+        full_exchange=full_exchange,
     )
     # Set internal field to manufactured values
     np.random.seed(0)
@@ -160,6 +167,62 @@ def test_mpi_ghost_communication(
         local_vector_field[:, ghost_size:-ghost_size, 0:ghost_size],
     )
 
+    if full_exchange:
+        # Check bottom left (ghost cell) == top right (inner cell)
+        np.testing.assert_allclose(
+            local_scalar_field[:ghost_size, :ghost_size],  # bottom left (ghost cell)
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size, -2 * ghost_size : -ghost_size
+            ],  # top right (inner cell)
+        )
+        np.testing.assert_allclose(
+            local_vector_field[:, :ghost_size, :ghost_size],  # bottom left (ghost cell)
+            local_vector_field[
+                :, -2 * ghost_size : -ghost_size, -2 * ghost_size : -ghost_size
+            ],  # top right (inner cell)
+        )
+        # Check bottom right (ghost cell) == top left (inner cell)
+        np.testing.assert_allclose(
+            local_scalar_field[:ghost_size, -ghost_size:],  # bottom right (ghost cell)
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size, ghost_size : 2 * ghost_size
+            ],  # top left (inner cell)
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :, :ghost_size, -ghost_size:
+            ],  # bottom right (ghost cell)
+            local_vector_field[
+                :, -2 * ghost_size : -ghost_size, ghost_size : 2 * ghost_size
+            ],  # top left (inner cell)
+        )
+        # Check top left (ghost cell) == bottom right (inner cell)
+        np.testing.assert_allclose(
+            local_scalar_field[-ghost_size:, :ghost_size],  # top left (ghost cell)
+            local_scalar_field[
+                ghost_size : 2 * ghost_size, -2 * ghost_size : -ghost_size
+            ],  # bottom right (inner cell)
+        )
+        np.testing.assert_allclose(
+            local_vector_field[:, -ghost_size:, :ghost_size],  # top left (ghost cell)
+            local_vector_field[
+                :, ghost_size : 2 * ghost_size, -2 * ghost_size : -ghost_size
+            ],  # bottom right (inner cell)
+        )
+        # Check top right (ghost cell) == bottom left (inner cell)
+        np.testing.assert_allclose(
+            local_scalar_field[-ghost_size:, -ghost_size:],  # top right (ghost cell)
+            local_scalar_field[
+                ghost_size : 2 * ghost_size, ghost_size : 2 * ghost_size
+            ],  # bottom left (inner cell)
+        )
+        np.testing.assert_allclose(
+            local_vector_field[:, -ghost_size:, -ghost_size:],  # top right (ghost cell)
+            local_vector_field[
+                :, ghost_size : 2 * ghost_size, ghost_size : 2 * ghost_size
+            ],  # bottom left (inner cell)
+        )
+
 
 @pytest.mark.mpi(group="MPI_utils", min_size=4)
 @pytest.mark.parametrize("precision", ["single", "double"])
@@ -167,7 +230,7 @@ def test_mpi_ghost_communication(
 @pytest.mark.parametrize("aspect_ratio", [(1, 1), (1, 2), (2, 1)])
 def test_mpi_lagrangian_field_map(precision, rank_distribution, aspect_ratio):
     # Eulerian grid stuff
-    n_values = 32
+    n_values = 8
     real_t = get_real_t(precision)
     grid_size_y = n_values * aspect_ratio[0]
     grid_size_x = n_values * aspect_ratio[1]
@@ -247,7 +310,7 @@ def test_mpi_lagrangian_field_gather_scatter(
     precision, rank_distribution, aspect_ratio
 ):
     # Eulerian grid stuff
-    n_values = 32
+    n_values = 8
     real_t = get_real_t(precision)
     grid_size_y = n_values * aspect_ratio[0]
     grid_size_x = n_values * aspect_ratio[1]
diff --git a/tests/test_utils/test_mpi_utils_3d.py b/tests/test_utils/test_mpi_utils_3d.py
index caad257..3a1f77f 100644
--- a/tests/test_utils/test_mpi_utils_3d.py
+++ b/tests/test_utils/test_mpi_utils_3d.py
@@ -12,7 +12,7 @@
 
 
 @pytest.mark.mpi(group="MPI_utils", min_size=4)
-@pytest.mark.parametrize("ghost_size", [1, 2, 3])
+@pytest.mark.parametrize("ghost_size", [1, 2])
 @pytest.mark.parametrize("precision", ["single", "double"])
 @pytest.mark.parametrize(
     "rank_distribution",
@@ -26,7 +26,7 @@
 def test_mpi_field_gather_scatter(
     ghost_size, precision, rank_distribution, aspect_ratio, master_rank
 ):
-    n_values = 32
+    n_values = 8
     real_t = get_real_t(precision)
     mpi_construct = MPIConstruct3D(
         grid_size_z=n_values * aspect_ratio[0],
@@ -94,7 +94,7 @@ def test_mpi_field_gather_scatter(
 
 
 @pytest.mark.mpi(group="MPI_utils", min_size=4)
-@pytest.mark.parametrize("ghost_size", [1, 2, 3])
+@pytest.mark.parametrize("ghost_size", [1, 2])
 @pytest.mark.parametrize("precision", ["single", "double"])
 @pytest.mark.parametrize(
     "rank_distribution",
@@ -104,22 +104,23 @@ def test_mpi_field_gather_scatter(
     "aspect_ratio",
     [(1, 1, 1), (1, 1, 2), (1, 2, 1), (2, 1, 1), (1, 2, 2), (2, 1, 2), (2, 2, 1)],
 )
+@pytest.mark.parametrize("full_exchange", [True, False])
 def test_mpi_ghost_communication(
-    ghost_size, precision, rank_distribution, aspect_ratio
+    ghost_size, precision, rank_distribution, aspect_ratio, full_exchange
 ):
-    n_values = 32
+    n_values = 8
     real_t = get_real_t(precision)
     mpi_construct = MPIConstruct3D(
         grid_size_z=n_values * aspect_ratio[0],
         grid_size_y=n_values * aspect_ratio[1],
         grid_size_x=n_values * aspect_ratio[2],
-        periodic_flag=True,
+        periodic_domain=True,
         real_t=real_t,
         rank_distribution=rank_distribution,
     )
     # extra width needed for kernel computation
     mpi_ghost_exchange_communicator = MPIGhostCommunicator3D(
-        ghost_size=ghost_size, mpi_construct=mpi_construct
+        ghost_size=ghost_size, mpi_construct=mpi_construct, full_exchange=full_exchange
     )
     # Set internal field to manufactured values
     np.random.seed(0)
@@ -141,8 +142,8 @@ def test_mpi_ghost_communication(
     mpi_ghost_exchange_communicator.exchange_finalise()
 
     # check if comm. done rightly!
-    # Test scalar field
-    # Along X: comm with previous block
+    # (1) Test faces
+    # Comm. along (0, 0, -X)
     np.testing.assert_allclose(
         local_scalar_field[
             ghost_size:-ghost_size, ghost_size:-ghost_size, ghost_size : 2 * ghost_size
@@ -153,118 +154,100 @@ def test_mpi_ghost_communication(
             -ghost_size : local_scalar_field.shape[2],
         ],
     )
-    # Along X: comm with next block
     np.testing.assert_allclose(
-        local_scalar_field[
+        local_vector_field[
+            :,
             ghost_size:-ghost_size,
             ghost_size:-ghost_size,
-            -2 * ghost_size : -ghost_size,
-        ],
-        local_scalar_field[
-            ghost_size:-ghost_size, ghost_size:-ghost_size, 0:ghost_size
-        ],
-    )
-    # Along Y: comm with previous block
-    np.testing.assert_allclose(
-        local_scalar_field[
-            ghost_size:-ghost_size, ghost_size : 2 * ghost_size, ghost_size:-ghost_size
+            ghost_size : 2 * ghost_size,
         ],
-        local_scalar_field[
+        local_vector_field[
+            :,
             ghost_size:-ghost_size,
-            -ghost_size : local_scalar_field.shape[1],
             ghost_size:-ghost_size,
+            -ghost_size : local_vector_field.shape[3],
         ],
     )
-    # Along Y: comm with next block
+    # Comm. along (0, 0, +X)
     np.testing.assert_allclose(
         local_scalar_field[
             ghost_size:-ghost_size,
-            -2 * ghost_size : -ghost_size,
             ghost_size:-ghost_size,
+            -2 * ghost_size : -ghost_size,
         ],
         local_scalar_field[
-            ghost_size:-ghost_size, 0:ghost_size, ghost_size:-ghost_size
+            ghost_size:-ghost_size, ghost_size:-ghost_size, 0:ghost_size
         ],
     )
-    # Along Z: comm with previous block
     np.testing.assert_allclose(
-        local_scalar_field[
-            ghost_size : 2 * ghost_size, ghost_size:-ghost_size, ghost_size:-ghost_size
-        ],
-        local_scalar_field[
-            -ghost_size : local_scalar_field.shape[0],
+        local_vector_field[
+            :,
             ghost_size:-ghost_size,
             ghost_size:-ghost_size,
+            -2 * ghost_size : -ghost_size,
+        ],
+        local_vector_field[
+            :, ghost_size:-ghost_size, ghost_size:-ghost_size, 0:ghost_size
         ],
     )
-    # Along Z: comm with next block
+    # Comm. along (0, -Y, 0)
     np.testing.assert_allclose(
         local_scalar_field[
-            -2 * ghost_size : -ghost_size,
-            ghost_size:-ghost_size,
-            ghost_size:-ghost_size,
+            ghost_size:-ghost_size, ghost_size : 2 * ghost_size, ghost_size:-ghost_size
         ],
         local_scalar_field[
-            0:ghost_size, ghost_size:-ghost_size, ghost_size:-ghost_size
+            ghost_size:-ghost_size,
+            -ghost_size : local_scalar_field.shape[1],
+            ghost_size:-ghost_size,
         ],
     )
-
-    # Test vector field
-    # Along X: comm with previous block
     np.testing.assert_allclose(
         local_vector_field[
             :,
             ghost_size:-ghost_size,
-            ghost_size:-ghost_size,
             ghost_size : 2 * ghost_size,
+            ghost_size:-ghost_size,
         ],
         local_vector_field[
             :,
             ghost_size:-ghost_size,
+            -ghost_size : local_vector_field.shape[2],
             ghost_size:-ghost_size,
-            -ghost_size : local_vector_field.shape[3],
         ],
     )
-    # Along X: comm with next block
+    # Comm. along (0, +Y, 0)
     np.testing.assert_allclose(
-        local_vector_field[
-            :,
-            ghost_size:-ghost_size,
+        local_scalar_field[
             ghost_size:-ghost_size,
             -2 * ghost_size : -ghost_size,
+            ghost_size:-ghost_size,
         ],
-        local_vector_field[
-            :, ghost_size:-ghost_size, ghost_size:-ghost_size, 0:ghost_size
+        local_scalar_field[
+            ghost_size:-ghost_size, 0:ghost_size, ghost_size:-ghost_size
         ],
     )
-    # Along Y: comm with previous block
     np.testing.assert_allclose(
         local_vector_field[
             :,
             ghost_size:-ghost_size,
-            ghost_size : 2 * ghost_size,
+            -2 * ghost_size : -ghost_size,
             ghost_size:-ghost_size,
         ],
         local_vector_field[
-            :,
-            ghost_size:-ghost_size,
-            -ghost_size : local_vector_field.shape[2],
-            ghost_size:-ghost_size,
+            :, ghost_size:-ghost_size, 0:ghost_size, ghost_size:-ghost_size
         ],
     )
-    # Along Y: comm with next block
+    # Comm. along (-Z, 0, 0)
     np.testing.assert_allclose(
-        local_vector_field[
-            :,
+        local_scalar_field[
+            ghost_size : 2 * ghost_size, ghost_size:-ghost_size, ghost_size:-ghost_size
+        ],
+        local_scalar_field[
+            -ghost_size : local_scalar_field.shape[0],
             ghost_size:-ghost_size,
-            -2 * ghost_size : -ghost_size,
             ghost_size:-ghost_size,
         ],
-        local_vector_field[
-            :, ghost_size:-ghost_size, 0:ghost_size, ghost_size:-ghost_size
-        ],
     )
-    # Along Z: comm with previous block
     np.testing.assert_allclose(
         local_vector_field[
             :,
@@ -279,7 +262,17 @@ def test_mpi_ghost_communication(
             ghost_size:-ghost_size,
         ],
     )
-    # Along Z: comm with next block
+    # Comm. along (+Z, 0, 0)
+    np.testing.assert_allclose(
+        local_scalar_field[
+            -2 * ghost_size : -ghost_size,
+            ghost_size:-ghost_size,
+            ghost_size:-ghost_size,
+        ],
+        local_scalar_field[
+            0:ghost_size, ghost_size:-ghost_size, ghost_size:-ghost_size
+        ],
+    )
     np.testing.assert_allclose(
         local_vector_field[
             :,
@@ -292,6 +285,502 @@ def test_mpi_ghost_communication(
         ],
     )
 
+    if full_exchange:
+        # (2) Test edges
+        # Comm. along (0, +Y, +X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size:-ghost_size,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_scalar_field[ghost_size:-ghost_size, 0:ghost_size, 0:ghost_size],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size:-ghost_size,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_vector_field[:, ghost_size:-ghost_size, 0:ghost_size, 0:ghost_size],
+        )
+        # Comm. along (0, -Y, +X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size:-ghost_size,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_scalar_field[
+                ghost_size:-ghost_size,
+                -ghost_size : local_scalar_field.shape[1],
+                0:ghost_size,
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size:-ghost_size,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_vector_field[
+                :,
+                ghost_size:-ghost_size,
+                -ghost_size : local_vector_field.shape[2],
+                0:ghost_size,
+            ],
+        )
+        # Comm. along (0, +Y, -X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size:-ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_scalar_field[
+                ghost_size:-ghost_size,
+                0:ghost_size,
+                -ghost_size : local_scalar_field.shape[2],
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size:-ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_vector_field[
+                :,
+                ghost_size:-ghost_size,
+                0:ghost_size,
+                -ghost_size : local_vector_field.shape[3],
+            ],
+        )
+        # Comm. along (0, -Y, -X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size:-ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_scalar_field[
+                ghost_size:-ghost_size,
+                -ghost_size : local_scalar_field.shape[1],
+                -ghost_size : local_scalar_field.shape[2],
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size:-ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_vector_field[
+                :,
+                ghost_size:-ghost_size,
+                -ghost_size : local_vector_field.shape[2],
+                -ghost_size : local_vector_field.shape[3],
+            ],
+        )
+
+        # Comm. along (+Z, 0, +X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size,
+                ghost_size:-ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_scalar_field[0:ghost_size, ghost_size:-ghost_size, 0:ghost_size],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                -2 * ghost_size : -ghost_size,
+                ghost_size:-ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_vector_field[:, 0:ghost_size, ghost_size:-ghost_size, 0:ghost_size],
+        )
+        # Comm. along (-Z, 0, +X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size : 2 * ghost_size,
+                ghost_size:-ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_scalar_field[
+                -ghost_size : local_scalar_field.shape[0],
+                ghost_size:-ghost_size,
+                0:ghost_size,
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size : 2 * ghost_size,
+                ghost_size:-ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_vector_field[
+                :,
+                -ghost_size : local_vector_field.shape[1],
+                ghost_size:-ghost_size,
+                0:ghost_size,
+            ],
+        )
+        # Comm. along (+Z, 0, -X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size,
+                ghost_size:-ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_scalar_field[
+                0:ghost_size,
+                ghost_size:-ghost_size,
+                -ghost_size : local_scalar_field.shape[2],
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                -2 * ghost_size : -ghost_size,
+                ghost_size:-ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_vector_field[
+                :,
+                0:ghost_size,
+                ghost_size:-ghost_size,
+                -ghost_size : local_vector_field.shape[3],
+            ],
+        )
+        # Comm. along (-Z, 0, -X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size : 2 * ghost_size,
+                ghost_size:-ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_scalar_field[
+                -ghost_size : local_scalar_field.shape[0],
+                ghost_size:-ghost_size,
+                -ghost_size : local_scalar_field.shape[2],
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size : 2 * ghost_size,
+                ghost_size:-ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_vector_field[
+                :,
+                -ghost_size : local_vector_field.shape[1],
+                ghost_size:-ghost_size,
+                -ghost_size : local_vector_field.shape[3],
+            ],
+        )
+
+        # Comm. along (+Z, +Y, 0)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size:-ghost_size,
+            ],
+            local_scalar_field[0:ghost_size, 0:ghost_size, ghost_size:-ghost_size],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size:-ghost_size,
+            ],
+            local_vector_field[:, 0:ghost_size, 0:ghost_size, ghost_size:-ghost_size],
+        )
+        # Comm. along (-Z, +Y, 0)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size:-ghost_size,
+            ],
+            local_scalar_field[
+                -ghost_size : local_scalar_field.shape[0],
+                0:ghost_size,
+                ghost_size:-ghost_size,
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size:-ghost_size,
+            ],
+            local_vector_field[
+                :,
+                -ghost_size : local_vector_field.shape[1],
+                0:ghost_size,
+                ghost_size:-ghost_size,
+            ],
+        )
+        # Comm. along (+Z, -Y, 0)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size:-ghost_size,
+            ],
+            local_scalar_field[
+                0:ghost_size,
+                -ghost_size : local_scalar_field.shape[1],
+                ghost_size:-ghost_size,
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size:-ghost_size,
+            ],
+            local_vector_field[
+                :,
+                0:ghost_size,
+                -ghost_size : local_vector_field.shape[2],
+                ghost_size:-ghost_size,
+            ],
+        )
+        # Comm. along (-Z, -Y, 0)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size:-ghost_size,
+            ],
+            local_scalar_field[
+                -ghost_size : local_scalar_field.shape[0],
+                -ghost_size : local_scalar_field.shape[1],
+                ghost_size:-ghost_size,
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size:-ghost_size,
+            ],
+            local_vector_field[
+                :,
+                -ghost_size : local_vector_field.shape[1],
+                -ghost_size : local_vector_field.shape[2],
+                ghost_size:-ghost_size,
+            ],
+        )
+
+        # (3) Test vertices
+        # Comm. along (+Z, +Y, +X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_scalar_field[0:ghost_size, 0:ghost_size, 0:ghost_size],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_vector_field[:, 0:ghost_size, 0:ghost_size, 0:ghost_size],
+        )
+        # Comm. along (-Z, +Y, +X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_scalar_field[
+                -ghost_size : local_scalar_field.shape[0], 0:ghost_size, 0:ghost_size
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_vector_field[
+                :, -ghost_size : local_vector_field.shape[1], 0:ghost_size, 0:ghost_size
+            ],
+        )
+        # Comm. along (+Z, -Y, +X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_scalar_field[
+                0:ghost_size, -ghost_size : local_scalar_field.shape[1], 0:ghost_size
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_vector_field[
+                :, 0:ghost_size, -ghost_size : local_vector_field.shape[2], 0:ghost_size
+            ],
+        )
+        # Comm. along (+Z, +Y, -X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_scalar_field[
+                0:ghost_size, 0:ghost_size, -ghost_size : local_scalar_field.shape[2]
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                -2 * ghost_size : -ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_vector_field[
+                :, 0:ghost_size, 0:ghost_size, -ghost_size : local_vector_field.shape[3]
+            ],
+        )
+        # Comm. along (-Z, -Y, +X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_scalar_field[
+                -ghost_size : local_scalar_field.shape[0],
+                -ghost_size : local_scalar_field.shape[1],
+                0:ghost_size,
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+            ],
+            local_vector_field[
+                :,
+                -ghost_size : local_vector_field.shape[1],
+                -ghost_size : local_vector_field.shape[2],
+                0:ghost_size,
+            ],
+        )
+        # Comm. along (-Z, +Y, -X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_scalar_field[
+                -ghost_size : local_scalar_field.shape[0],
+                0:ghost_size,
+                -ghost_size : local_scalar_field.shape[2],
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size : 2 * ghost_size,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_vector_field[
+                :,
+                -ghost_size : local_vector_field.shape[1],
+                0:ghost_size,
+                -ghost_size : local_vector_field.shape[3],
+            ],
+        )
+        # Comm. along (+Z, -Y, -X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_scalar_field[
+                0:ghost_size,
+                -ghost_size : local_scalar_field.shape[1],
+                -ghost_size : local_scalar_field.shape[2],
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                -2 * ghost_size : -ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_vector_field[
+                :,
+                0:ghost_size,
+                -ghost_size : local_vector_field.shape[2],
+                -ghost_size : local_vector_field.shape[3],
+            ],
+        )
+        # Comm. along (-Z, -Y, -X)
+        np.testing.assert_allclose(
+            local_scalar_field[
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_scalar_field[
+                -ghost_size : local_scalar_field.shape[0],
+                -ghost_size : local_scalar_field.shape[1],
+                -ghost_size : local_scalar_field.shape[2],
+            ],
+        )
+        np.testing.assert_allclose(
+            local_vector_field[
+                :,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+                ghost_size : 2 * ghost_size,
+            ],
+            local_vector_field[
+                :,
+                -ghost_size : local_vector_field.shape[1],
+                -ghost_size : local_vector_field.shape[2],
+                -ghost_size : local_vector_field.shape[3],
+            ],
+        )
+
 
 @pytest.mark.mpi(group="MPI_utils", min_size=4)
 @pytest.mark.parametrize("precision", ["single", "double"])
@@ -305,7 +794,7 @@ def test_mpi_ghost_communication(
 )
 def test_mpi_lagrangian_field_map(precision, rank_distribution, aspect_ratio):
     # Eulerian grid stuff
-    n_values = 32
+    n_values = 8
     real_t = get_real_t(precision)
     grid_size_z = n_values * aspect_ratio[0]
     grid_size_y = n_values * aspect_ratio[1]
@@ -400,7 +889,7 @@ def test_mpi_lagrangian_field_gather_scatter(
     precision, rank_distribution, aspect_ratio
 ):
     # Eulerian grid stuff
-    n_values = 32
+    n_values = 8
     real_t = get_real_t(precision)
     grid_size_z = n_values * aspect_ratio[0]
     grid_size_y = n_values * aspect_ratio[1]