diff --git a/prototyping.jl b/prototyping.jl
index dfd40faa..9e46288b 100644
--- a/prototyping.jl
+++ b/prototyping.jl
@@ -29,6 +29,7 @@ phi = assign(phi,
 ) 
 
 gradPhi = Grad{Orthogonal}(phi)
+gradPhi = Grad{Midpoint}(phi)
 phif = FaceScalarField(mesh)
 
 solvers = (
@@ -54,4 +55,16 @@ scatter(cell_xyz.⋅Ref([1,0,0]), gradPhi.result.x.values)
 
 XCALibre.Solvers.limit_gradient!(gradPhi, phi, config)
 
-scatter(cell_xyz.⋅Ref([1,0,0]), gradPhi.result.x.values)
+scatter!(cell_xyz.⋅Ref([1,0,0]), gradPhi.result.x.values, alpha=0.25)
+
+fRange = mesh.boundaries[3].IDs_range
+cIDs = mesh.boundary_cellsID[fRange]
+b1_xyz = getproperty.(mesh.cells[cIDs], :centre)
+scatter!(b1_xyz.⋅Ref([1,0,0]), gradPhi.result.x.values[cIDs], alpha=0.50)
+
+vtkWriter = initialise_writer(mesh)
+
+write_vtk("uncorrected", mesh, vtkWriter, ("gradPhi", gradPhi.result.x)) #, Ux, Uy, Uz, p)
+write_vtk("uncorrected_orthogonal", mesh, vtkWriter, ("gradPhi", gradPhi.result.x)) #, Ux, Uy, Uz, p)
+write_vtk("corrected", mesh, vtkWriter, ("gradPhi", gradPhi.result.x)) #, Ux, Uy, Uz, p)
+write_vtk("corrected_orthogonal", mesh, vtkWriter, ("gradPhi", gradPhi.result.x)) #, Ux, Uy, Uz, p)
diff --git a/src/Solvers/Solvers_2_PISO.jl b/src/Solvers/Solvers_2_PISO.jl
index 23c12856..31441dd1 100644
--- a/src/Solvers/Solvers_2_PISO.jl
+++ b/src/Solvers/Solvers_2_PISO.jl
@@ -318,35 +318,43 @@ end
     end
 
     # g0 = ∇F[cID]
-    g0 = SVector{3}(x[cID] , y[cID] , z[cID])
+    grad0 = SVector{3}(x[cID] , y[cID] , z[cID])
 
     cc = cell.centre
     limiter = 1
     limiterf = 1
     for fi ∈ faces_range 
         fID = cell_faces[fi]
+        nID = cell_neighbours[fi]
         face = faces[fID]
-        nID = face.ownerCells[2]
-        phiN = F[nID]
+        cellN = cells[nID]
+        # nID = face.ownerCells[2]
+        # phiN = F[nID]
         normal = face.normal
         nsign = cell_nsign[fi]
-        # na = nsign*normal
-
-        fc = face.centre 
-        rf = fc - cc
-        δϕ = g0⋅rf
-        phif = phiP + δϕ
-        diff = phif - phiP
-        if diff > 0
-            limiterf = min(1, (phiN - phiP)/diff)
-        elseif diff < 0
-            # limiterf = min(1, (phiN - phiP)/diff)
-        elseif diff == 0
+        na = nsign*normal
+
+        # r = fc - cc
+        # fc = face.centre
+
+        nc = cellN.centre
+        r = nc - cc
+        δϕ = r⋅grad0
+
+        # rn = (nc - cc) ⋅ na
+        # gradn = grad0⋅na
+        # δϕ = rn* gradn
+        if δϕ > 0
+            limiterf = min(1, (phiMax - phiP)/δϕ)
+        elseif δϕ < 0
+            limiterf = min(1, (phiMin - phiP)/δϕ)
+        else
             limiterf = 1
         end
         limiter = min(limiterf, limiter)
     end
-    x.values[cID] = limiter*g0[1]
-    y.values[cID] = limiter*g0[2]
-    z.values[cID] = limiter*g0[3]
+    grad0 *= limiter
+    x.values[cID] = grad0[1]
+    y.values[cID] = grad0[2]
+    z.values[cID] = grad0[3]
 end
\ No newline at end of file