Merge pull request #101 from mberto79/HM/normal-aligned-interpolation-weights

Face interpolation weights using face normal distances

Author: mberto79

CHANGELOG.md

@@ -19,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * +50x speed improvement for the method `construct_periodic` and also more robust algorithm used [#100](@ref)
 * Implementation to correct mass flux uses matrix coefficients directly for better stability when using periodic boundary conditions [#100](@ref)
 * New method to enforce matrix symmetry of scalar model equations when the only term is a laplacian [#100](@ref)
+* Change calculation of face interpolation weights to use face normal aligned weights, this is more physical than the current method using face-based distances (in preparation for formal support for non-orthogonality correction)[#101](@ref)
 
 ### Breaking
 * No breaking changes

src/Discretise/Discretise_1_schemes.jl

@@ -109,7 +109,8 @@ end
     )  where {F,P,I}
 
     w = face.weight
-    signbit(ns) ? w = one(w) - w : w
+    # signbit(ns) ? w = one(w) - w : w
+    w = 0.5 + ns*(w - 0.5)
 
     # Calculate link coefficients
     ap = term.sign*(term.flux[fID]*ns)

src/Discretise/boundary_conditions/periodic.jl

@@ -326,13 +326,23 @@ end
     pfID = bc.value.face_map[i] # id of periodic face 
     pface = faces[pfID]
     pcellID = pface.ownerCells[1]
+    C1 = cell.centre
+    C2 = cells[pcellID].centre - transform.distance
+    Cf = face.centre
+    n = face.normal
+
+    Pf = Cf - C1
+    PN = C2 - C1 
 
-    w = pface.delta/(face.delta + pface.delta)
+    wn = (Pf⋅n)/(PN⋅n)
+    w = one(wn) - wn
+    # w = pface.delta/(face.delta + pface.delta)
+    # wn = one(w) - w
 
     # Calculate link coefficients
     ap = term.sign*(term.flux[fID])
     ac = ap*w
-    an = ap*(one(w) - w)
+    an = ap*wn
 
     NN = spindex(rowptr, colval, pcellID, pcellID)
     NP = spindex(rowptr, colval, pcellID, cellID)
@@ -395,15 +405,21 @@ end
     pfID = bc.value.face_map[i] # id of periodic face 
     pface = faces[pfID]
     pcellID = pface.ownerCells[1]
+    C1 = cell.centre
+    C2 = cells[pcellID].centre - transform.distance
+    Cf = face.centre
+    n = face.normal
 
+    Pf = Cf - C1
+    PN = C2 - C1 
 
-    # Calculate interpoloation weight
-    w = pface.delta/(face.delta + pface.delta)
+    wn = (Pf⋅n)/(PN⋅n)
+    w = one(wn) - wn
 
     # Calculate link coefficients
     mdot = term.sign*(term.flux[fID])
     acLinear = mdot*w 
-    anLinear = mdot*(one(w) - w)
+    anLinear = mdot*wn
     acUpwind = max(mdot, 0.0) 
     anUpwind = -max(-mdot, 0.0)
     ac = 0.75*acLinear + 0.25*acUpwind

src/Discretise/boundary_conditions/periodic_interpolation.jl

@@ -14,19 +14,24 @@ end
         pfID = BC.value.face_map[i] # id of periodic face
         pface = faces[pfID]
         pcID = pface.ownerCells[1]
-        # pcell = cells[pcID]
+        pcell = cells[pcID]
         face = faces[fID]
-        cID = boundary_cellsID[fID]
+        # cID = boundary_cellsID[fID]
+        cID = face.ownerCells[1]
+        cell = cells[cID]
 
-        delta1 = face.delta #*norm(face.e ⋅ face.normal)
-        delta2 = pface.delta #*norm(pface.e ⋅ pface.normal)
-        delta = delta1 + delta2
-        weight = delta2/delta
+        # delta1 = face.delta #*norm(face.e ⋅ face.normal)
+        # delta2 = pface.delta #*norm(pface.e ⋅ pface.normal)
+        # delta = delta1 + delta2
+        # w = delta2/delta
 
-        one_minus_weight = one(weight) - weight
+        Pf = face.centre - cell.centre
+        PN = (pcell.centre - transform.distance) - cell.centre
+        normal = face.normal
+        wn = (Pf⋅normal)/(PN⋅normal)
+        w = one(wn) - wn
 
-        # phif_values[fID] = 0.5*(phi_values[cID] + phi_values[pcID]) # linear interpolation
-        phifi =  weight*phi[cID] + one_minus_weight*phi[pcID]
+        phifi =  w*phi[cID] + wn*phi[pcID]
         phif[fID] = phifi
         phif[pfID] = phifi
     end
@@ -38,21 +43,29 @@ end
     boundary_cellsID, time, fID)
     @inbounds begin 
         i = fID - BC.IDs_range.start + 1
-        (; transform) = BC.value
+        (; transform ) = BC.value
         (; faces, cells) = psif.mesh
         pfID = BC.value.face_map[i] # id of periodic face
         pface = faces[pfID]
         pcID = pface.ownerCells[1]
+        pcell = cells[pcID]
         face = faces[fID]
-        cID = boundary_cellsID[fID]
+        # cID = boundary_cellsID[fID]
+        cID = face.ownerCells[1]
+        cell = cells[cID]
+
+        # delta1 = face.delta #*norm(face.e ⋅ face.normal)
+        # delta2 = pface.delta #*norm(pface.e ⋅ pface.normal)
+        # delta = delta1 + delta2
+        # w = delta2/delta
 
-        delta1 = face.delta #*norm(face.e ⋅ face.normal)
-        delta2 = pface.delta #*norm(pface.e ⋅ pface.normal)
-        delta = delta1 + delta2
-        w = delta2/delta
-        one_w = one(w) - w
+        Pf = face.centre - cell.centre
+        PN = (pcell.centre - transform.distance) - cell.centre
+        normal = face.normal
+        wn = (Pf⋅normal)/(PN⋅normal)
+        w = one(wn) - wn
 
-        psifi = w*psi[cID] + one_w*psi[pcID] # linear interpolation 
+        psifi = w*psi[cID] + wn*psi[pcID] # linear interpolation 
         psif[fID] = psifi
         psif[pfID] = psifi
     end

src/Mesh/Mesh_1_functions.jl

@@ -16,7 +16,9 @@ _get_backend(mesh) = get_backend(mesh.cells)
 # C2F1 = distance vector from cell2 centre to face centre
 # C1C2 = distance vector from cell1 to cell2
 weight_delta_e(C1F1, C2F1, C1C2, normal) = begin
-    weight = norm(C2F1)/(norm(C1F1) + norm(C2F1))
+    # weight = norm(C2F1)/(norm(C1F1) + norm(C2F1)) # face-distance based
+    wi = (C1F1⋅normal)/(C1C2⋅normal)
+    weight = one(wi) - wi # normal aligned interpolation weight
     delta = norm(C1C2)
     e = C1C2/delta
     return weight, delta, e