diff --git a/examples/CFD/cavity3d.py b/examples/CFD/cavity3d.py
index f7e35c8..0afa648 100644
--- a/examples/CFD/cavity3d.py
+++ b/examples/CFD/cavity3d.py
@@ -17,7 +17,7 @@
 # Use 8 CPU devices
 # os.environ["XLA_FLAGS"] = '--xla_force_host_platform_device_count=8'
 from src.models import BGKSim, KBCSim
-from src.lattice import LatticeD3Q27
+from src.lattice import LatticeD3Q19
 from src.utils import *
 from jax.config import config
 from src.boundary_conditions import *
@@ -31,6 +31,13 @@ def __init__(self, **kwargs):
         super().__init__(**kwargs)
 
     def set_boundary_conditions(self):
+
+        # apply inlet equilibrium boundary condition to the top wall
+        moving_wall = self.boundingBoxIndices['top']
+        vel_wall = np.zeros(moving_wall.shape, dtype=self.precisionPolicy.compute_dtype)
+        vel_wall[:, 0] = prescribed_vel
+        self.BCs.append(BounceBackHalfway(tuple(moving_wall.T), self.gridInfo, self.precisionPolicy, vel_wall))
+
         # concatenate the indices of the left, right, and bottom walls
         walls = np.concatenate(
             (self.boundingBoxIndices['left'], self.boundingBoxIndices['right'],
@@ -38,14 +45,7 @@ def set_boundary_conditions(self):
              self.boundingBoxIndices['bottom']))
         # apply bounce back boundary condition to the walls
         self.BCs.append(BounceBackHalfway(tuple(walls.T), self.gridInfo, self.precisionPolicy))
-
-        # apply inlet equilibrium boundary condition to the top wall
-        moving_wall = self.boundingBoxIndices['top']
-
-        rho_wall = np.ones((moving_wall.shape[0], 1), dtype=self.precisionPolicy.compute_dtype)
-        vel_wall = np.zeros(moving_wall.shape, dtype=self.precisionPolicy.compute_dtype)
-        vel_wall[:, 0] = prescribed_vel
-        self.BCs.append(BounceBackHalfway(tuple(moving_wall.T), self.gridInfo, self.precisionPolicy, vel_wall))
+        return
 
     def output_data(self, **kwargs):
         # 1: -1 to remove boundary voxels (not needed for visualization when using full-way bounce-back)
@@ -64,8 +64,10 @@ def output_data(self, **kwargs):
 
         # output profiles of velocity at mid-plane for benchmarking
         output_filename = "./profiles_" + f"{timestep:07d}.json"
-        ldc_ref_result = {'ux(x=y=0)': list(u[nx//2, ny//2, :, 0]/prescribed_vel),
-                          'uz(z=y=0)': list(u[:, ny//2, nz//2, 2]/prescribed_vel)}
+        ux_mid = 0.5*(u[nx//2, ny//2, :, 0] + u[nx//2+1, ny//2+1, :, 0])
+        uz_mid = 0.5*(u[:, ny//2, nz//2, 2] + u[:, ny//2+1, nz//2+1, 2])
+        ldc_ref_result = {'ux(x=y=0)': list(ux_mid/prescribed_vel),
+                          'uz(z=y=0)': list(uz_mid/prescribed_vel)}
         json.dump(ldc_ref_result, codecs.open(output_filename, 'w', encoding='utf-8'),
                 separators=(',', ':'),
                 sort_keys=True,
@@ -76,7 +78,7 @@ def output_data(self, **kwargs):
         # live_volume_randering(timestep, u_mag)
 
 if __name__ == '__main__':
-    lattice = LatticeD3Q27(precision)
+    lattice = LatticeD3Q19(precision)
 
     nx = 256
     ny = 256