more robust implementation of gradient limiter

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)

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