FDS Source: Define CFACE RDN from linked-cell bounding box normal pro…

Source/ccib.f90

@@ -54,6 +54,7 @@ MODULE CC_SCALARS
 INTEGER, ALLOCATABLE, DIMENSION(:)   :: JM_MAT_Z
 
 REAL(EB), ALLOCATABLE, DIMENSION(:)  :: F_Z, RZ_Z, RZ_ZS, P_0_CV, TMP_0_CV, RHO_0_CV, ZCEN_CV
+REAL(EB), ALLOCATABLE, DIMENSION(:,:):: DELTA_UNKZ
 REAL(EB), ALLOCATABLE, DIMENSION(:,:):: F_Z0, RZ_Z0
 
 REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP
@@ -587,9 +588,10 @@ SUBROUTINE CHECK_CFLVN_LINKED_CELLS(NM,DT,UVWMAX,R_DX2,MUTRM)
 
 ! Local variables:
 INTEGER :: I,J,K,ICC,JCC,IROW,IMAX,X1AXIS,ILH,IRC,IFC,IFACE,IFC2,IFACE2,ICFA,LOWHIGH
-REAL(EB):: MU_TMP,MURDN,AF,VELN,CFLMAX_TMP,VNMAX_TMP,TWOD_FCT
+REAL(EB):: MU_TMP,MURDN,AF,VELN,CFLMAX_TMP,VNMAX_TMP,TWOD_FCT,MU_AVG
 INTEGER,  ALLOCATABLE, DIMENSION(:,:) :: IJKT
 REAL(EB), ALLOCATABLE, DIMENSION(:)   :: UVWA, MUV, MURA, VOL, DIVG
+INTEGER, PARAMETER :: VN_DXN_METHOD = 2 ! 1: face-based sum, 2: face-based sum + dx/dy/dz in large CVs.
 
 IF(NUNKZ_LOC(NM)<1) RETURN
 
@@ -776,9 +778,25 @@ SUBROUTINE CHECK_CFLVN_LINKED_CELLS(NM,DT,UVWMAX,R_DX2,MUTRM)
 ENDIF
 
 IF (CHECK_VN) THEN
-   DO IROW=UNKZ_ILC(NM)+1,UNKZ_ILC(NM)+NUNKZ_LOC(NM)
-      MURA(IROW) = MURA(IROW)/VOL(IROW)
-   ENDDO
+   IF (VN_DXN_METHOD==2) THEN
+      DO IROW=UNKZ_ILC(NM)+1,UNKZ_ILC(NM)+NUNKZ_LOC(NM)
+         I = IJKT(IAXIS,IROW); J = IJKT(JAXIS,IROW); K = IJKT(KAXIS,IROW)
+         IF (VOL(IROW) >= DEFAULT_VOLFRAC_LINK*DX(I)*DY(J)*DZ(K)) THEN ! For large linked CVs relax VN criterion.
+            MU_AVG = MUV(IROW)/VOL(IROW)
+            IF (TWO_D) THEN
+               MURA(IROW) = 2._EB*MU_AVG*(1._EB/DX(I)**2 + 1._EB/DZ(K)**2)
+            ELSE
+               MURA(IROW) = 2._EB*MU_AVG*(1._EB/DX(I)**2 + 1._EB/DY(J)**2 + 1._EB/DZ(K)**2)
+            ENDIF
+         ELSE
+            MURA(IROW) = MURA(IROW)/VOL(IROW)
+         ENDIF
+      ENDDO
+   ELSE
+      DO IROW=UNKZ_ILC(NM)+1,UNKZ_ILC(NM)+NUNKZ_LOC(NM)
+         MURA(IROW) = MURA(IROW)/VOL(IROW)
+      ENDDO
+   ENDIF
    IMAX = UNKZ_ILC(NM)+MAXLOC(MURA,DIM=1)
    VNMAX_TMP = DT * MURA(IMAX)
    IF(VNMAX_TMP > VN) THEN
@@ -7133,57 +7151,77 @@ SUBROUTINE SET_CFACES_P1_RDN
 INTEGER :: ICF, IFACE
 INTEGER :: ICC, JCC, I, J, K
 INTEGER :: IFACE_CELL, ICF_CELL, IROW
+INTEGER :: IROW_DELTA
 REAL(EB):: AREAI
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: DXN_UNKZ_LOC, AREA_UNKZ_LOC, VOL_UNKZ_LOC
+INTEGER, PARAMETER :: RDN_METHOD = 2 ! 1: VOL/SOLIDAREA, 2: bbox projection
+REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: XYZMIN_UNKZ, XYZMAX_UNKZ
+REAL(EB) :: DELTA_X, DELTA_Y, DELTA_Z, DXN, NVEC_ABS(MAX_DIM)
+TYPE(CFACE_TYPE), POINTER :: CFA
+TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
 
 ! ALLOCATE local arrays
 ALLOCATE(DXN_UNKZ_LOC(1:NUNKZ_LOCAL));  DXN_UNKZ_LOC(:)  = 0._EB
 ALLOCATE(AREA_UNKZ_LOC(1:NUNKZ_LOCAL)); AREA_UNKZ_LOC(:) = 0._EB
 ALLOCATE(VOL_UNKZ_LOC(1:NUNKZ_LOCAL));  VOL_UNKZ_LOC(:)  = 0._EB
+IF (RDN_METHOD==1) THEN
+   ! Main Loop:
+   MESH_LOOP_01 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
 
-! Main Loop:
-MESH_LOOP_01 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
-
-   CALL POINT_TO_MESH(NM)
+      CALL POINT_TO_MESH(NM)
 
-   ! Do a volume weighted average of distance to wall from linked cells, if one of them is a regular cell use 1/2 the
-   ! distance of corner to corner sqrt(DX^2+DY^2+DZ^2).
-   ! 1. Regular GASPHASE cells within the cc-region:
-   DO K=1,KBAR
-      DO J=1,JBAR
-         DO I=1,IBAR
-            IF (CCVAR(I,J,K,CC_UNKZ) <= 0 ) CYCLE ! Drop if regular gas cell has not been assigned unknown number.
-            IROW = CCVAR(I,J,K,CC_UNKZ) - UNKZ_IND(NM_START) ! All row indexes must refer to ind_loc.
-            VOL_UNKZ_LOC(IROW) = VOL_UNKZ_LOC(IROW) + (DX(I)*DY(J)*DZ(K))
+      ! Do a volume weighted average of distance to wall from linked cells, if one of them is a regular cell use 1/2 the
+      ! distance of corner to corner sqrt(DX^2+DY^2+DZ^2).
+      ! 1. Regular GASPHASE cells within the cc-region:
+      DO K=1,KBAR
+         DO J=1,JBAR
+            DO I=1,IBAR
+               IF (CCVAR(I,J,K,CC_UNKZ) <= 0 ) CYCLE ! Drop if regular gas cell has not been assigned unknown number.
+               IROW = CCVAR(I,J,K,CC_UNKZ) - UNKZ_IND(NM_START) ! All row indexes must refer to ind_loc.
+               VOL_UNKZ_LOC(IROW) = VOL_UNKZ_LOC(IROW) + (DX(I)*DY(J)*DZ(K))
+            ENDDO
          ENDDO
       ENDDO
-   ENDDO
-   ! 2. Cut-cells:
-   DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
-      CC => CUT_CELL(ICC); I = CC%IJK(IAXIS); J = CC%IJK(JAXIS); K = CC%IJK(KAXIS)
-      IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
-      DO JCC=1,CC%NCELL
-         IROW = CC%UNKZ(JCC) - UNKZ_IND(NM_START)
-         ! Mean INBOUNDARY cut-face distance to this cut-cell center, projected to cut-face normal:
-         AREAI = 0._EB
-         DO ICF_CELL=1,CC%CCELEM(1,JCC)
-            IFACE_CELL = CC%CCELEM(ICF_CELL+1,JCC)
-            IF (CC%FACE_LIST(1,IFACE_CELL) /= CC_FTYPE_CFINB) CYCLE
-            ! Indexes of INBOUNDARY cutface on CUT_FACE:
-            ICF   = CC%FACE_LIST(4,IFACE_CELL)
-            IFACE = CC%FACE_LIST(5,IFACE_CELL)
-            ! Area sum:
-            AREAI = AREAI + CUT_FACE(ICF)%AREA(IFACE)
+      ! 2. Cut-cells:
+      DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
+         CC => CUT_CELL(ICC); I = CC%IJK(IAXIS); J = CC%IJK(JAXIS); K = CC%IJK(KAXIS)
+         IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
+         DO JCC=1,CC%NCELL
+            IROW = CC%UNKZ(JCC) - UNKZ_IND(NM_START)
+            ! Mean INBOUNDARY cut-face distance to this cut-cell center, projected to cut-face normal:
+            AREAI = 0._EB
+            DO ICF_CELL=1,CC%CCELEM(1,JCC)
+               IFACE_CELL = CC%CCELEM(ICF_CELL+1,JCC)
+               IF (CC%FACE_LIST(1,IFACE_CELL) /= CC_FTYPE_CFINB) CYCLE
+               ! Indexes of INBOUNDARY cutface on CUT_FACE:
+               ICF   = CC%FACE_LIST(4,IFACE_CELL)
+               IFACE = CC%FACE_LIST(5,IFACE_CELL)
+               ! Area sum:
+               AREAI = AREAI + CUT_FACE(ICF)%AREA(IFACE)
+            ENDDO
+            AREA_UNKZ_LOC(IROW) = AREA_UNKZ_LOC(IROW) + AREAI
+            VOL_UNKZ_LOC(IROW) = VOL_UNKZ_LOC(IROW) + CC%VOLUME(JCC)
          ENDDO
-         AREA_UNKZ_LOC(IROW) = AREA_UNKZ_LOC(IROW) + AREAI
-         VOL_UNKZ_LOC(IROW) = VOL_UNKZ_LOC(IROW) + CC%VOLUME(JCC)
       ENDDO
-   ENDDO
 
-ENDDO MESH_LOOP_01
+   ENDDO MESH_LOOP_01
 
-! Compute average DXN of all linked cells:
-DXN_UNKZ_LOC = VOL_UNKZ_LOC / (AREA_UNKZ_LOC + TWENTY_EPSILON_EB)
+   ! Compute average DXN of all linked cells:
+   DXN_UNKZ_LOC = VOL_UNKZ_LOC / (AREA_UNKZ_LOC + TWO_EPSILON_EB)
+ELSE
+   ALLOCATE(XYZMIN_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL))
+   ALLOCATE(XYZMAX_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL))
+   CALL GET_LINKED_CELLS_BBOX(XYZMIN_UNKZ,XYZMAX_UNKZ)
+   IF (ALLOCATED(DELTA_UNKZ)) THEN
+      IF (SIZE(DELTA_UNKZ,DIM=2) /= NUNKZ_LOCAL) DEALLOCATE(DELTA_UNKZ)
+   ENDIF
+   IF (.NOT.ALLOCATED(DELTA_UNKZ)) ALLOCATE(DELTA_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL))
+   DO IROW_DELTA=1,NUNKZ_LOCAL
+      DELTA_UNKZ(IAXIS,IROW_DELTA) = MAX(0._EB,XYZMAX_UNKZ(IAXIS,IROW_DELTA)-XYZMIN_UNKZ(IAXIS,IROW_DELTA))
+      DELTA_UNKZ(JAXIS,IROW_DELTA) = MAX(0._EB,XYZMAX_UNKZ(JAXIS,IROW_DELTA)-XYZMIN_UNKZ(JAXIS,IROW_DELTA))
+      DELTA_UNKZ(KAXIS,IROW_DELTA) = MAX(0._EB,XYZMAX_UNKZ(KAXIS,IROW_DELTA)-XYZMIN_UNKZ(KAXIS,IROW_DELTA))
+   ENDDO
+ENDIF
 
 ! Finally Define B1%RDN:
 MESH_LOOP_02 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
@@ -7196,15 +7234,123 @@ SUBROUTINE SET_CFACES_P1_RDN
          ICC = CF%CELL_LIST(2,LOW_IND,IFACE)
          JCC = CF%CELL_LIST(3,LOW_IND,IFACE)
          IROW = CUT_CELL(ICC)%UNKZ(JCC) - UNKZ_IND(NM_START)
-         BOUNDARY_PROP1(CFACE(CF%CFACE_INDEX(IFACE))%B1_INDEX)%RDN = 1._EB/DXN_UNKZ_LOC(IROW)
+         SELECT CASE (RDN_METHOD)
+         CASE (1)
+            BOUNDARY_PROP1(CFACE(CF%CFACE_INDEX(IFACE))%B1_INDEX)%RDN = 1._EB/DXN_UNKZ_LOC(IROW)
+         CASE (2)
+            CFA => CFACE(CF%CFACE_INDEX(IFACE))
+            BC  => BOUNDARY_COORD(CFA%BC_INDEX)
+            DELTA_X = XYZMAX_UNKZ(IAXIS,IROW) - XYZMIN_UNKZ(IAXIS,IROW)
+            DELTA_Y = XYZMAX_UNKZ(JAXIS,IROW) - XYZMIN_UNKZ(JAXIS,IROW)
+            DELTA_Z = XYZMAX_UNKZ(KAXIS,IROW) - XYZMIN_UNKZ(KAXIS,IROW)
+            NVEC_ABS(IAXIS:KAXIS) = ABS(BC%NVEC(IAXIS:KAXIS))
+            DXN = DELTA_X*NVEC_ABS(IAXIS) + DELTA_Y*NVEC_ABS(JAXIS) + DELTA_Z*NVEC_ABS(KAXIS)
+            BOUNDARY_PROP1(CFA%B1_INDEX)%RDN = 1._EB/(DXN + TWO_EPSILON_EB)
+         END SELECT
       ENDDO
    ENDDO
 ENDDO MESH_LOOP_02
+IF (ALLOCATED(XYZMIN_UNKZ)) DEALLOCATE(XYZMIN_UNKZ, XYZMAX_UNKZ)
 DEALLOCATE(DXN_UNKZ_LOC, VOL_UNKZ_LOC, AREA_UNKZ_LOC)
 
 RETURN
 END SUBROUTINE SET_CFACES_P1_RDN
 
+! ---------------------- GET_LINKED_CELLS_BBOX --------------------------------
+
+SUBROUTINE GET_LINKED_CELLS_BBOX(XYZMIN_UNKZ,XYZMAX_UNKZ)
+
+! Local Variables:
+INTEGER :: I, J, K, ICC, JCC, IROW
+INTEGER :: IFACE, IFC, FTYPE, X1AXIS, LOWHIGH, ILH, IFC2, IFACE2, NVFACE, IPT, IVERT
+REAL(EB) :: XLO, XHI, YLO, YHI, ZLO, ZHI
+REAL(EB), INTENT(OUT) :: XYZMIN_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL)
+REAL(EB), INTENT(OUT) :: XYZMAX_UNKZ(IAXIS:KAXIS,1:NUNKZ_LOCAL)
+
+XYZMIN_UNKZ =  HUGE(1._EB); XYZMAX_UNKZ = -HUGE(1._EB)
+MESH_LOOP_01 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
+   CALL POINT_TO_MESH(NM)
+   DO K=1,KBAR
+      DO J=1,JBAR
+         DO I=1,IBAR
+            IF (CCVAR(I,J,K,CC_UNKZ) <= 0 ) CYCLE
+            IROW = CCVAR(I,J,K,CC_UNKZ) - UNKZ_IND(NM_START)
+            XLO = X(I-1); XHI = X(I)
+            YLO = Y(J-1); YHI = Y(J)
+            ZLO = Z(K-1); ZHI = Z(K)
+            XYZMIN_UNKZ(IAXIS,IROW) = MIN(XYZMIN_UNKZ(IAXIS,IROW),XLO)
+            XYZMIN_UNKZ(JAXIS,IROW) = MIN(XYZMIN_UNKZ(JAXIS,IROW),YLO)
+            XYZMIN_UNKZ(KAXIS,IROW) = MIN(XYZMIN_UNKZ(KAXIS,IROW),ZLO)
+            XYZMAX_UNKZ(IAXIS,IROW) = MAX(XYZMAX_UNKZ(IAXIS,IROW),XHI)
+            XYZMAX_UNKZ(JAXIS,IROW) = MAX(XYZMAX_UNKZ(JAXIS,IROW),YHI)
+            XYZMAX_UNKZ(KAXIS,IROW) = MAX(XYZMAX_UNKZ(KAXIS,IROW),ZHI)
+         ENDDO
+      ENDDO
+   ENDDO
+
+   DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
+      CC => CUT_CELL(ICC); I = CC%IJK(IAXIS); J = CC%IJK(JAXIS); K = CC%IJK(KAXIS)
+      IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
+      DO JCC=1,CC%NCELL
+         IROW = CC%UNKZ(JCC) - UNKZ_IND(NM_START)
+         XLO =  HUGE(1._EB); XHI = -HUGE(1._EB)
+         YLO =  HUGE(1._EB); YHI = -HUGE(1._EB)
+         ZLO =  HUGE(1._EB); ZHI = -HUGE(1._EB)
+         DO IFC=1,CC%CCELEM(1,JCC)
+            IFACE = CC%CCELEM(IFC+1,JCC)
+            FTYPE = CC%FACE_LIST(1,IFACE)
+            SELECT CASE (FTYPE)
+            CASE (CC_FTYPE_CFGAS,CC_FTYPE_CFINB)
+               IFC2   = CC%FACE_LIST(4,IFACE)
+               IFACE2 = CC%FACE_LIST(5,IFACE)
+               NVFACE = CUT_FACE(IFC2)%CFELEM(1,IFACE2)
+               DO IPT=1,NVFACE
+                  IVERT = CUT_FACE(IFC2)%CFELEM(IPT+1,IFACE2)
+                  XLO = MIN(XLO,CUT_FACE(IFC2)%XYZVERT(IAXIS,IVERT))
+                  XHI = MAX(XHI,CUT_FACE(IFC2)%XYZVERT(IAXIS,IVERT))
+                  YLO = MIN(YLO,CUT_FACE(IFC2)%XYZVERT(JAXIS,IVERT))
+                  YHI = MAX(YHI,CUT_FACE(IFC2)%XYZVERT(JAXIS,IVERT))
+                  ZLO = MIN(ZLO,CUT_FACE(IFC2)%XYZVERT(KAXIS,IVERT))
+                  ZHI = MAX(ZHI,CUT_FACE(IFC2)%XYZVERT(KAXIS,IVERT))
+               ENDDO
+            CASE (CC_FTYPE_RCGAS)
+               LOWHIGH = CC%FACE_LIST(2,IFACE)
+               X1AXIS  = CC%FACE_LIST(3,IFACE)
+               ILH     = LOWHIGH - 1
+               SELECT CASE (X1AXIS)
+               CASE (IAXIS)
+                  XLO = MIN(XLO,X(I-1+ILH)); XHI = MAX(XHI,X(I-1+ILH))
+                  YLO = MIN(YLO,Y(J-1));     YHI = MAX(YHI,Y(J))
+                  ZLO = MIN(ZLO,Z(K-1));     ZHI = MAX(ZHI,Z(K))
+               CASE (JAXIS)
+                  XLO = MIN(XLO,X(I-1));     XHI = MAX(XHI,X(I))
+                  YLO = MIN(YLO,Y(J-1+ILH)); YHI = MAX(YHI,Y(J-1+ILH))
+                  ZLO = MIN(ZLO,Z(K-1));     ZHI = MAX(ZHI,Z(K))
+               CASE (KAXIS)
+                  XLO = MIN(XLO,X(I-1));     XHI = MAX(XHI,X(I))
+                  YLO = MIN(YLO,Y(J-1));     YHI = MAX(YHI,Y(J))
+                  ZLO = MIN(ZLO,Z(K-1+ILH)); ZHI = MAX(ZHI,Z(K-1+ILH))
+               END SELECT
+            END SELECT
+         ENDDO
+         IF (XLO > XHI) THEN
+            XLO = X(I-1); XHI = X(I)
+            YLO = Y(J-1); YHI = Y(J)
+            ZLO = Z(K-1); ZHI = Z(K)
+         ENDIF
+         XYZMIN_UNKZ(IAXIS,IROW) = MIN(XYZMIN_UNKZ(IAXIS,IROW),XLO)
+         XYZMIN_UNKZ(JAXIS,IROW) = MIN(XYZMIN_UNKZ(JAXIS,IROW),YLO)
+         XYZMIN_UNKZ(KAXIS,IROW) = MIN(XYZMIN_UNKZ(KAXIS,IROW),ZLO)
+         XYZMAX_UNKZ(IAXIS,IROW) = MAX(XYZMAX_UNKZ(IAXIS,IROW),XHI)
+         XYZMAX_UNKZ(JAXIS,IROW) = MAX(XYZMAX_UNKZ(JAXIS,IROW),YHI)
+         XYZMAX_UNKZ(KAXIS,IROW) = MAX(XYZMAX_UNKZ(KAXIS,IROW),ZHI)
+      ENDDO
+   ENDDO
+ENDDO MESH_LOOP_01
+
+RETURN
+END SUBROUTINE GET_LINKED_CELLS_BBOX
+
 END SUBROUTINE CC_SET_DATA
 
 ! ------------------------------- CC_END_STEP --------------------------------

Source/cons.f90

@@ -697,6 +697,7 @@ MODULE GLOBAL_CONSTANTS
 
 ! Threshold factor for volume of cut-cells respect to volume of Cartesian cells:
 ! Currently used in the thermo div definition of cut-cells.
+REAL(EB), PARAMETER :: DEFAULT_VOLFRAC_LINK = 0.5_EB
 
 REAL(EB) :: CCVOL_LINK=0.95_EB
 LOGICAL  :: GET_CUTCELLS_VERBOSE=.FALSE.

Source/geom.f90

@@ -23781,7 +23781,7 @@ SUBROUTINE READ_GEOM
    PRES_FLAG   = ULMAT_FLAG
 ENDIF
 PRES_ON_WHOLE_DOMAIN = .FALSE.
-IF (ABS(CCVOL_LINK-0.95_EB)<TWENTY_EPSILON_EB) CCVOL_LINK = 0.5_EB
+IF (ABS(CCVOL_LINK-0.95_EB)<TWENTY_EPSILON_EB) CCVOL_LINK = DEFAULT_VOLFRAC_LINK
 IF (CHECK_POISSON) GLMAT_VERBOSE=.TRUE.
 
 CONTAINS