Merge remote-tracking branch 'origin/dev-0.3/gradient-limiter-refinem…

…ent' into dev-0.3-main

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)

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