From 772a428c00ddb205f22f40b39e6d9cfb5c718392 Mon Sep 17 00:00:00 2001
From: marcosvanella <marcos_vanella@hotmail.com>
Date: Fri, 16 Aug 2024 11:01:17 -0400
Subject: [PATCH] FDS Source : Add CFACEs advection velocity for vorticity term
 in momentum RHS.

---
 Source/ccib.f90 | 38 +++++++++++++++++++++++++++++++++++++-
 1 file changed, 37 insertions(+), 1 deletion(-)

diff --git a/Source/ccib.f90 b/Source/ccib.f90
index 02067cc5704..a55eacb4427 100644
--- a/Source/ccib.f90
+++ b/Source/ccib.f90
@@ -394,7 +394,12 @@ SUBROUTINE CUTFACE_VELOCITIES(NM,UU,VV,WW,CUTFACES)
 REAL(EB), POINTER, DIMENSION(:,:,:), INTENT(INOUT) :: UU,VV,WW
 LOGICAL, INTENT(IN) :: CUTFACES
 TYPE(CC_CUTFACE_TYPE), POINTER :: CF=>NULL()
-INTEGER :: ICF,I,J,K,X1AXIS
+TYPE(CFACE_TYPE), POINTER :: CFA=>NULL()
+TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
+TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
+INTEGER :: ICF,ICFA,JCF,I,J,K,X1AXIS
+REAL(EB):: AREA,VELN(3),PREDFCT
+
 
 CUTFACES_IF : IF (CUTFACES) THEN ! USE CUT_FACE(ICF)%VEL_CF
    DO ICF=1,MESHES(NM)%N_CUTFACE_MESH+MESHES(NM)%N_GCCUTFACE_MESH
@@ -410,6 +415,29 @@ SUBROUTINE CUTFACE_VELOCITIES(NM,UU,VV,WW,CUTFACES)
       END SELECT
    ENDDO
 
+   PREDFCT = 0._EB; IF(PREDICTOR) PREDFCT = 1._EB
+   ! CFACEs, set velocity in underlaying solid cartesian faces to be used in VELOCITY_FLUX:
+   DO ICF=1,MESHES(NM)%N_CUTFACE_MESH
+      CF => CUT_FACE(ICF); IF(CF%STATUS/=CC_INBOUNDARY) CYCLE
+      I = CF%IJK(IAXIS); J = CF%IJK(JAXIS); K = CF%IJK(KAXIS)
+      ! Area Average velocity for boundary CFACEs:
+      AREA = 0._EB; VELN(IAXIS:KAXIS) = 0._EB
+      DO JCF=1,CF%NFACE
+         ICFA=CF%CFACE_INDEX(JCF); IF(ICFA<1) CYCLE
+         CFA => CFACE(ICFA)
+         BC => BOUNDARY_COORD(CFA%BC_INDEX)
+         B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
+         AREA = AREA+CFA%AREA
+         VELN(IAXIS:KAXIS) = VELN(IAXIS:KAXIS) - &
+                             (PREDFCT*B1%U_NORMAL+(1._EB-PREDFCT)*B1%U_NORMAL_S)*CFA%AREA*BC%NVEC(IAXIS:KAXIS)
+      ENDDO
+      VELN(IAXIS:KAXIS) = VELN(IAXIS:KAXIS)/(AREA+TWO_EPSILON_EB)
+      ! Distribute into Solid cartesian faces:
+      WHERE(FCVAR(I-1:I,J,K,CC_FGSC,IAXIS)==CC_SOLID) UU(I-1:I,J,K) = VELN(IAXIS)
+      WHERE(FCVAR(I,J-1:J,K,CC_FGSC,JAXIS)==CC_SOLID) VV(I,J-1:J,K) = VELN(JAXIS)
+      WHERE(FCVAR(I,J,K-1:K,CC_FGSC,KAXIS)==CC_SOLID) WW(I,J,K-1:K) = VELN(KAXIS)
+   ENDDO
+
 ELSE CUTFACES_IF ! USE CUT_FACE(ICF)%VEL_CRT
    DO ICF=1,MESHES(NM)%N_CUTFACE_MESH+MESHES(NM)%N_GCCUTFACE_MESH
       CF => CUT_FACE(ICF); IF(CF%STATUS/=CC_GASPHASE) CYCLE
@@ -424,6 +452,14 @@ SUBROUTINE CUTFACE_VELOCITIES(NM,UU,VV,WW,CUTFACES)
       END SELECT
    ENDDO
 
+   DO ICF=1,MESHES(NM)%N_CUTFACE_MESH
+      CF => CUT_FACE(ICF); IF(CF%STATUS/=CC_INBOUNDARY) CYCLE
+      I = CF%IJK(IAXIS); J = CF%IJK(JAXIS); K = CF%IJK(KAXIS)
+      WHERE(FCVAR(I-1:I,J,K,CC_FGSC,IAXIS)==CC_SOLID) UU(I-1:I,J,K) = 0._EB
+      WHERE(FCVAR(I,J-1:J,K,CC_FGSC,JAXIS)==CC_SOLID) VV(I,J-1:J,K) = 0._EB
+      WHERE(FCVAR(I,J,K-1:K,CC_FGSC,KAXIS)==CC_SOLID) WW(I,J,K-1:K) = 0._EB
+   ENDDO
+
 ENDIF CUTFACES_IF
 RETURN
 END SUBROUTINE CUTFACE_VELOCITIES