From 2a83c710b9a5729cb55dc453a54ab3e2046eb16b Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Mon, 1 Jul 2019 18:26:31 +0100
Subject: [PATCH] Update sGamma handling

---
 .../immersedBoundaryFvPatch.C                   | 17 ++++++++++-------
 .../include/immersedBoundaryCourantNo.H         |  4 ++--
 2 files changed, 12 insertions(+), 9 deletions(-)

diff --git a/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C b/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C
index 831c1436d..4ce76b9da 100644
--- a/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C
+++ b/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C
@@ -124,7 +124,7 @@ void Foam::immersedBoundaryFvPatch::updatePhi
     scalar rDeltaT = 1.0/deltaT;
 
     // Multiply the raw mesh motion flux with the masking function
-    scalarField sGamma =
+    scalarField ibAreaRatio =
         mag(ibPolyPatch_.correctedFaceAreas())/mag(mesh.faceAreas());
 
     // Scaling of internal mesh flux field should be done only for the current
@@ -140,7 +140,7 @@ void Foam::immersedBoundaryFvPatch::updatePhi
         const label faceI = deadFaces[dfI];
         if (mesh.isInternalFace(faceI))
         {
-            phiIn[faceI] *= sGamma[faceI];
+            phiIn[faceI] = scalar(0);
         }
     }
  
@@ -150,7 +150,7 @@ void Foam::immersedBoundaryFvPatch::updatePhi
         const label faceI = cutFaces[cfI];
         if (mesh.isInternalFace(faceI))
         {
-            phiIn[faceI] *= sGamma[faceI];
+            phiIn[faceI] *= ibAreaRatio[faceI];
         }
     }
 
@@ -160,7 +160,7 @@ void Foam::immersedBoundaryFvPatch::updatePhi
         if (!isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
         {
             phi.boundaryField()[patchI] *=
-                mesh.boundary()[patchI].patchSlice(sGamma);
+                mesh.boundary()[patchI].patchSlice(ibAreaRatio);
         }
     }
 
@@ -231,12 +231,11 @@ void Foam::immersedBoundaryFvPatch::updatePhi
             // Attempt to correct via old volume
             scalar corrOldVol = newVols[cellI] - divPhi[cellI]*deltaT;
 
-            // Info<< "Flux maneouvre for cell " << cellI << ": "
+            // Pout<< "Flux maneouvre for cell " << cellI << ": "
             //     << " error: " << magDivPhi[cellI]
             //     << " V: " << newVols[cellI]
             //     << " V0: " << oldVols[cellI]
-            //     << " divPhi: " << divPhi[cellI]
-            //     << endl;
+            //     << " divPhi: " << divPhi[cellI];
 
             if (corrOldVol < SMALL)
             {
@@ -248,6 +247,10 @@ void Foam::immersedBoundaryFvPatch::updatePhi
             {
                 oldVols[cellI] = corrOldVol;
             }
+
+            // scalar corrDivMeshPhi =
+            //     mag((newVols[cellI] - oldVols[cellI]) - divPhi[cellI]*deltaT);
+            // Pout<< " Corrected: " << corrDivMeshPhi << endl;
         }
     }
 }
diff --git a/src/immersedBoundary/immersedBoundary/include/immersedBoundaryCourantNo.H b/src/immersedBoundary/immersedBoundary/include/immersedBoundaryCourantNo.H
index 41707b63a..05c1c266d 100644
--- a/src/immersedBoundary/immersedBoundary/include/immersedBoundaryCourantNo.H
+++ b/src/immersedBoundary/immersedBoundary/include/immersedBoundaryCourantNo.H
@@ -36,7 +36,7 @@ scalar velMag = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField magPhi = pos(sGamma - 0.5)*mag(phi);
+    surfaceScalarField magPhi = sGamma*mag(phi);
 
     surfaceScalarField SfUfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*magPhi;
@@ -50,7 +50,7 @@ if (mesh.nInternalFaces())
     velMag = max(magPhi/mesh.magSf()).value();
 }
 
-Info<< "Courant Number mean: " << meanCoNum
+Info<< "Immersed boundary Courant Number mean: " << meanCoNum
     << " max: " << CoNum
     << " velocity magnitude: " << velMag
     << endl;