diff --git a/prototyping.jl b/prototyping.jl
new file mode 100644
index 00000000..9e46288b
--- /dev/null
+++ b/prototyping.jl
@@ -0,0 +1,70 @@
+using XCALibre
+using LinearAlgebra
+using Plots
+
+grids_dir = pkgdir(XCALibre, "examples/0_GRIDS")
+# trig, trig40, trig100, quad, quad40, quad100
+grid = "trig40.unv"
+mesh_file = joinpath(grids_dir, grid)
+
+mesh = UNV2D_mesh(mesh_file, scale=0.001)
+
+phi = ScalarField(mesh)
+
+cell_xyz = getproperty.(mesh.cells, :centre)
+
+phi.values
+
+sin.(cell_xyz.⋅Ref([1,0,0])*π)
+
+phi.values .= sin.(cell_xyz.⋅Ref([1,0,0])*π)
+
+scatter(cell_xyz.⋅Ref([1,0,0]), phi.values)
+
+phi = assign(phi,
+    Dirichlet(:inlet, 0.0),
+    Dirichlet(:outlet, 0.0),
+    Neumann(:top, 0.0),
+    Neumann(:bottom, 0.0)
+) 
+
+gradPhi = Grad{Orthogonal}(phi)
+gradPhi = Grad{Midpoint}(phi)
+phif = FaceScalarField(mesh)
+
+solvers = (
+    phi = set_solver(
+        phi;
+        solver      = BicgstabSolver, # BicgstabSolver, GmresSolver
+        preconditioner = Jacobi(),
+        convergence = 1e-7,
+        relax = 1.0
+    ),
+)
+
+schemes = (phi = set_schemes(),)
+runtime = set_runtime(iterations=1, write_interval=-1, time_step=1)
+hardware = set_hardware(backend=CPU(), workgroup=1024)
+
+config = Configuration(
+    solvers=solvers, schemes=schemes, runtime=runtime, hardware=hardware)
+
+grad!(gradPhi, phif, phi, phi.BCs, 0, config)
+
+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, 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 a0304c11..31441dd1 100644
--- a/src/Solvers/Solvers_2_PISO.jl
+++ b/src/Solvers/Solvers_2_PISO.jl
@@ -256,47 +256,105 @@ function limit_gradient!(∇F::Grad{S,FF,R,I,M}, F, config) where {S,FF,R<:Tenso
     KernelAbstractions.synchronize(backend)
 end
 
+# @kernel function _limit_gradient!(x, y, z, F, cells, cell_neighbours, cell_faces, cell_nsign, faces, minPhi, maxPhi)
+#     cID = @index(Global)
+
+#     cell = cells[cID]
+#     faces_range = cell.faces_range
+#     phiP = F[cID]
+#     maxPhi = minPhi = phiP
+
+#     for fi ∈ faces_range
+#         nID = cell_neighbours[fi]
+#         phiN = F[nID]
+#         maxPhi = max(phiN, maxPhi)
+#         minPhi = min(phiN, minPhi)
+#     end
+
+#     # g0 = ∇F[cID]
+#     g0 = SVector{3}(x[cID] , y[cID] , z[cID])
+
+#     cc = cell.centre
+
+#     for fi ∈ faces_range 
+#         fID = cell_faces[fi]
+#         face = faces[fID]
+#         nID = face.ownerCells[2]
+#         # phiN = F[nID]
+#         normal = face.normal
+#         nsign = cell_nsign[fi]
+#         na = nsign*normal
+
+#         fc = face.centre 
+#         cc_fc = fc - cc
+#         n0 = cc_fc/norm(cc_fc)
+#         gn = g0⋅n0
+#         δϕ = g0⋅cc_fc
+#         gτ = g0 - gn*n0
+#         if (maxPhi > phiP) && (δϕ > maxPhi - phiP)
+#             g0 = gτ + na*(maxPhi - phiP)
+#         elseif (minPhi < phiP) && (δϕ < minPhi - phiP)
+#             g0 = gτ + na*(minPhi - phiP)
+#         end            
+#     end
+#     x.values[cID] = g0[1]
+#     y.values[cID] = g0[2]
+#     z.values[cID] = g0[3]
+# end
+
 @kernel function _limit_gradient!(x, y, z, F, cells, cell_neighbours, cell_faces, cell_nsign, faces, minPhi, maxPhi)
     cID = @index(Global)
 
     cell = cells[cID]
     faces_range = cell.faces_range
     phiP = F[cID]
-
+    phiMax = phiMin = phiP
+ 
     for fi ∈ faces_range
         nID = cell_neighbours[fi]
         phiN = F[nID]
-        maxPhi = max(phiN, maxPhi)
-        minPhi = min(phiN, minPhi)
+        phiMax = max(phiN, phiMax)
+        phiMin = min(phiN, phiMin)
     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]
+        cellN = cells[nID]
+        # nID = face.ownerCells[2]
         # phiN = F[nID]
         normal = face.normal
         nsign = cell_nsign[fi]
         na = nsign*normal
 
-        fc = face.centre 
-        cc_fc = fc - cc
-        n0 = cc_fc/norm(cc_fc)
-        gn = g0⋅n0
-        δϕ = g0⋅cc_fc
-        gτ = g0 - gn*n0
-        if (maxPhi > phiP) && (δϕ > maxPhi - phiP)
-            g0 = gτ + na*(maxPhi - phiP)
-        elseif (minPhi < phiP) && (δϕ < minPhi - phiP)
-            g0 = gτ + na*(minPhi - phiP)
-        end            
+        # 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] = g0[1]
-    y.values[cID] = g0[2]
-    z.values[cID] = 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