Merge pull request #13942 from marcosvanella/master

FDS Source : Clean PRES_ON_WHOLE_DOMAIN, ccib.f90.
firemodels · Dec 24, 2024 · 287cc5e · 287cc5e
2 parents e11400e + 40e3576
commit 287cc5e
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 193 deletions.
diff --git a/Source/ccib.f90 b/Source/ccib.f90
@@ -2469,29 +2469,9 @@ SUBROUTINE CC_PROJECT_VELOCITY(NM,DT,STORE_FLG)
          END SELECT
       ENDDO WALL_CELL_LOOP_1
 
-      IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
-         DO K=1,KBAR
-            DO J=1,JBAR
-               DO I=0,IBAR
-                  IF(FCVAR(I,J,K,CC_FGSC,IAXIS)==CC_SOLID) US(I,J,K) = 0._EB
-               ENDDO
-            ENDDO
-         ENDDO
-         DO K=1,KBAR
-            DO J=0,JBAR
-               DO I=1,IBAR
-                  IF(FCVAR(I,J,K,CC_FGSC,JAXIS)==CC_SOLID) VS(I,J,K) = 0._EB
-               ENDDO
-            ENDDO
-         ENDDO
-         DO K=0,KBAR
-            DO J=1,JBAR
-               DO I=1,IBAR
-                  IF(FCVAR(I,J,K,CC_FGSC,KAXIS)==CC_SOLID) WS(I,J,K) = 0._EB
-               ENDDO
-            ENDDO
-         ENDDO
-      ENDIF
+      WHERE(FCVAR(0:IBAR,1:JBAR,1:KBAR,CC_FGSC,IAXIS)==CC_SOLID) US(0:IBAR,1:JBAR,1:KBAR) = 0._EB
+      WHERE(FCVAR(1:IBAR,0:JBAR,1:KBAR,CC_FGSC,JAXIS)==CC_SOLID) VS(1:IBAR,0:JBAR,1:KBAR) = 0._EB
+      WHERE(FCVAR(1:IBAR,1:JBAR,0:KBAR,CC_FGSC,KAXIS)==CC_SOLID) WS(1:IBAR,1:JBAR,0:KBAR) = 0._EB
 
    ELSE PRED_CORR_IF
 
@@ -2631,29 +2611,9 @@ SUBROUTINE CC_PROJECT_VELOCITY(NM,DT,STORE_FLG)
 
       DEALLOCATE(U_STORE,V_STORE,W_STORE)
 
-      IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
-         DO K=1,KBAR
-            DO J=1,JBAR
-               DO I=0,IBAR
-                  IF(FCVAR(I,J,K,CC_FGSC,IAXIS)==CC_SOLID) U(I,J,K) = 0._EB
-               ENDDO
-            ENDDO
-         ENDDO
-         DO K=1,KBAR
-            DO J=0,JBAR
-               DO I=1,IBAR
-                  IF(FCVAR(I,J,K,CC_FGSC,JAXIS)==CC_SOLID) V(I,J,K) = 0._EB
-               ENDDO
-            ENDDO
-         ENDDO
-         DO K=0,KBAR
-            DO J=1,JBAR
-               DO I=1,IBAR
-                  IF(FCVAR(I,J,K,CC_FGSC,KAXIS)==CC_SOLID) W(I,J,K) = 0._EB
-               ENDDO
-            ENDDO
-         ENDDO
-      ENDIF
+      WHERE(FCVAR(0:IBAR,1:JBAR,1:KBAR,CC_FGSC,IAXIS)==CC_SOLID) U(0:IBAR,1:JBAR,1:KBAR) = 0._EB
+      WHERE(FCVAR(1:IBAR,0:JBAR,1:KBAR,CC_FGSC,JAXIS)==CC_SOLID) V(1:IBAR,0:JBAR,1:KBAR) = 0._EB
+      WHERE(FCVAR(1:IBAR,1:JBAR,0:KBAR,CC_FGSC,KAXIS)==CC_SOLID) W(1:IBAR,1:JBAR,0:KBAR) = 0._EB
 
    ENDIF PRED_CORR_IF
 
@@ -6810,7 +6770,7 @@ SUBROUTINE CC_H_INTERP
 ! Local Variables:
 REAL(EB), POINTER, DIMENSION(:,:,:) :: UP,VP,WP,HP
 INTEGER :: NM, ICC, NCELL, I, J ,K
-REAL(EB):: U_IBM, V_IBM, W_IBM, VCRT
+REAL(EB):: VCRT
 LOGICAL :: VOLFLG
 
 ! This routine interpolates H to cut cells/Cartesian cells at the end of step.
@@ -6866,52 +6826,12 @@ SUBROUTINE CC_H_INTERP
    ENDDO ICC_LOOP
 
    ! Finally set HP to zero inside immersed solids:
-   IF (.NOT.PRES_ON_WHOLE_DOMAIN) THEN
-   DO K=0,KBP1
-     DO J=0,JBP1
-        DO I=0,IBP1
-           IF (MESHES(NM)%CCVAR(I,J,K,CC_CGSC) /= CC_SOLID) CYCLE
-           HP(I,J,K) = 0._EB
-        ENDDO
-     ENDDO
-   ENDDO
-   ENDIF
-
-   ! In case of .NOT. PRES_ON_WHOLE_DOMAIN set velocities on solid faces to zero:
-   IF (.NOT.PRES_ON_WHOLE_DOMAIN) THEN
-   ! Force U velocities in CC_SOLID faces to zero
-   U_IBM = 0._EB ! Body doesn't move.
-   DO K=1,KBAR
-      DO J=1,JBAR
-         DO I=0,IBAR
-            IF (MESHES(NM)%FCVAR(I,J,K,CC_FGSC,IAXIS) /= CC_SOLID ) CYCLE
-            UP(I,J,K) = U_IBM
-         ENDDO
-      ENDDO
-   ENDDO
-
-   ! Force V velocities in CC_SOLID faces to zero
-   V_IBM = 0._EB ! Body doesn't move.
-   DO K=1,KBAR
-      DO J=0,JBAR
-         DO I=1,IBAR
-            IF (MESHES(NM)%FCVAR(I,J,K,CC_FGSC,JAXIS) /= CC_SOLID ) CYCLE
-            VP(I,J,K) = V_IBM
-         ENDDO
-      ENDDO
-   ENDDO
+   WHERE(MESHES(NM)%CCVAR(0:IBP1,0:JBP1,0:KBP1,CC_CGSC)==CC_SOLID) HP(0:IBP1,0:JBP1,0:KBP1) = 0._EB
 
-   ! Force W velocities in CC_SOLID faces to zero
-   W_IBM = 0._EB ! Body doesn't move.
-   DO K=0,KBAR
-      DO J=1,JBAR
-         DO I=1,IBAR
-            IF (MESHES(NM)%FCVAR(I,J,K,CC_FGSC,KAXIS) /= CC_SOLID ) CYCLE
-            WP(I,J,K) = W_IBM
-         ENDDO
-      ENDDO
-   ENDDO
-   ENDIF
+   ! Set velocities on solid faces to zero:
+   WHERE(MESHES(NM)%FCVAR(0:IBAR,1:JBAR,1:KBAR,CC_FGSC,IAXIS)==CC_SOLID) UP(0:IBAR,1:JBAR,1:KBAR) = 0._EB
+   WHERE(MESHES(NM)%FCVAR(1:IBAR,0:JBAR,1:KBAR,CC_FGSC,JAXIS)==CC_SOLID) VP(1:IBAR,0:JBAR,1:KBAR) = 0._EB
+   WHERE(MESHES(NM)%FCVAR(1:IBAR,1:JBAR,0:KBAR,CC_FGSC,KAXIS)==CC_SOLID) WP(1:IBAR,1:JBAR,0:KBAR) = 0._EB
 
    NULLIFY(UP,VP,WP,HP)
 
@@ -7356,11 +7276,9 @@ SUBROUTINE INIT_CUTCELL_DATA(T,DT,FIRST_CALL)
 LOGICAL, INTENT(IN)  :: FIRST_CALL
 
 ! Local Variables:
-INTEGER :: NM,I,J,K,N,ICC,JCC,X1AXIS,NFACE,ICF,IFACE
+INTEGER :: NM,I,J,K,N,ICC,JCC,X1AXIS,NFACE,ICF
 REAL(EB) TMP_CC,RHO_CC,AREAT,VEL_CF !,Z_CC,TMP_0_CC,P_0_CC
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_CC
-INTEGER :: INDADD, INDF, IFC, IFACE2, ICFC
-REAL(EB):: FSCU, AREATOT
 INTEGER :: IW,IROW_LOC
 
 REAL(EB) :: TNOW
@@ -7599,57 +7517,8 @@ SUBROUTINE INIT_CUTCELL_DATA(T,DT,FIRST_CALL)
          ENDDO
       ENDDO
 
-      ! Then INBOUNDARY cut-faces:
-      ICC_LOOP : DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
-         I      = CUT_CELL(ICC)%IJK(IAXIS)
-         J      = CUT_CELL(ICC)%IJK(JAXIS)
-         K      = CUT_CELL(ICC)%IJK(KAXIS)
-         FSCU = 0._EB
-
-         ! Loop on cells neighbors and test if they are of type CC_SOLID, if so
-         ! Add to velocity flux:
-         ! X faces
-         DO INDADD=-1,1,2
-            INDF = I - 1 + (INDADD+1)/2
-            IF( FCVAR(INDF,J,K,CC_FGSC,IAXIS) /= CC_SOLID) CYCLE
-            FSCU = FSCU + REAL(INDADD,EB)*U(INDF,J,K)*DY(J)*DZ(K)
-         ENDDO
-         ! Y faces
-         DO INDADD=-1,1,2
-            INDF = J - 1 + (INDADD+1)/2
-            IF( FCVAR(I,INDF,K,CC_FGSC,JAXIS) /= CC_SOLID ) CYCLE
-            FSCU = FSCU + REAL(INDADD,EB)*V(I,INDF,K)*DX(I)*DZ(K)
-         ENDDO
-         ! Z faces
-         DO INDADD=-1,1,2
-            INDF = K - 1 + (INDADD+1)/2
-            IF( FCVAR(I,J,INDF,CC_FGSC,KAXIS) /= CC_SOLID ) CYCLE
-            FSCU = FSCU + REAL(INDADD,EB)*W(I,J,INDF)*DX(I)*DY(J)
-         ENDDO
-
-         ! Now Define total area of INBOUNDARY cut-faces:
-         ICF=CCVAR(I,J,K,CC_IDCF);
-         ICF_COND : IF (ICF > 0) THEN
-            NFACE = CUT_FACE(ICF)%NFACE
-            AREATOT = SUM ( CUT_FACE(ICF)%AREA(1:NFACE) )
-            DO JCC =1,CUT_CELL(ICC)%NCELL
-               IFC_LOOP : DO IFC=1,CUT_CELL(ICC)%CCELEM(1,JCC)
-                  IFACE = CUT_CELL(ICC)%CCELEM(IFC+1,JCC)
-                  IF (CUT_CELL(ICC)%FACE_LIST(1,IFACE) == CC_FTYPE_CFINB) THEN
-                     IFACE2  = CUT_CELL(ICC)%FACE_LIST(5,IFACE)
-                     ICFC    = CUT_FACE(ICF)%CFACE_INDEX(IFACE2)
-                     IF(PRES_ON_WHOLE_DOMAIN) THEN
-                        CUT_FACE(ICF)%VELS(IFACE2)  = 1._EB/AREATOT*FSCU ! +ve into the solid Velocity error
-                        CUT_FACE(ICF)%VEL( IFACE2)  = 1._EB/AREATOT*FSCU
-                     ELSE
-                        CUT_FACE(ICF)%VELS(IFACE2)  = 0._EB
-                        CUT_FACE(ICF)%VEL( IFACE2)  = 0._EB
-                     ENDIF
-                  ENDIF
-               ENDDO IFC_LOOP
-            ENDDO
-         ENDIF ICF_COND
-      ENDDO ICC_LOOP
+      ! INBOUNDARY cut-faces are initialized with 0._EB velocity, that will be changed in
+      ! CFACE_PREDICT_NORMAL_VELOCITY.
 
    ENDIF PERIODIC_TEST_COND
 
@@ -11590,7 +11459,7 @@ SUBROUTINE CC_DENSITY_EXPLICIT(T,DT)
       F_Z(:) = 0._EB
       CALL GET_EXPLICIT_ADVDIFFVECTOR_SCALAR_3D(N)
 
-      ! Add Advective fluxes due to PRES_ON_WHOLE_DOMAIN for F_Z:
+      ! Add Advective fluxes for F_Z:
       CALL GET_ADVDIFFVECTOR_SCALAR_3D(N)
 
       ! Here add the reaction source term M_DOT_PPP, treated explicitly:
@@ -12144,9 +12013,6 @@ SUBROUTINE GET_EXPLICIT_ADVDIFFVECTOR_SCALAR_3D(N)
 
    ENDDO
 
-   ! Case of PRES_ON_WHOLE_DOMAIN, we have non-zero velocities on INBOUNDARY cut-faces:
-   ! Done on CALL GET_ADVDIFFVECTOR_SCALAR_3D(N)
-
    ! Then add (Del rho D Del Z)*dv computed on CCDIVERGENCE_PART_1:
    ! Loop over regular cells on CC region:
    DO K=1,KBAR
@@ -15800,7 +15666,7 @@ SUBROUTINE CC_CHECK_DIVERGENCE(T,DT,PREDVEL)
       ENDIF
    ENDDO ICC_LOOP
 
-   IF(.NOT.PRES_ON_WHOLE_DOMAIN .AND. STORE_CARTESIAN_DIVERGENCE) THEN
+   IF(STORE_CARTESIAN_DIVERGENCE) THEN
       DO K=1,KBAR
          DO J=1,JBAR
             DO I=1,IBAR
@@ -19774,16 +19640,11 @@ SUBROUTINE GET_CC_UNKH(I,J,K,IUNKH)
 INTEGER :: ICC
 
 IUNKH    = CC_UNDEFINED ! This is < 0.
-IF(.NOT.PRES_ON_WHOLE_DOMAIN ) THEN
-
-   ! Regular gas cell, taken care of before.
-   ! Check cut-cell:
-   ICC = CCVAR(I,J,K,CC_IDCC)
-   ! If theres is a cut-cell ICC then CUT_CELL(ICC)%UNKH(1) has been populated.
-   IF (ICC > 0) IUNKH = CUT_CELL(ICC)%UNKH(1)
-
-ENDIF
-
+! Regular gas cell, taken care of before.
+! Check cut-cell:
+ICC = CCVAR(I,J,K,CC_IDCC)
+! If theres is a cut-cell ICC then CUT_CELL(ICC)%UNKH(1) has been populated.
+IF (ICC > 0) IUNKH = CUT_CELL(ICC)%UNKH(1)
 
 RETURN
 END SUBROUTINE GET_CC_UNKH
@@ -19798,12 +19659,10 @@ SUBROUTINE GET_CC_IROW(I,J,K,IPZ,IROW)
 
 ! Local variable:
 INTEGER :: ICC
-IROW    = CC_UNDEFINED ! This is < 0.
-IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
-   ICC = CCVAR(I,J,K,CC_IDCC)
-   ! If theres is a cut-cell ICC then CUT_CELL(ICC)%UNKH(1) has been populated.
-   IF (ICC > 0) IROW     = CUT_CELL(ICC)%UNKH(1) - ZONE_SOLVE(IPZ)%UNKH_IND(NM_START)
-ENDIF
+IROW = CC_UNDEFINED ! This is < 0.
+ICC  = CCVAR(I,J,K,CC_IDCC)
+! If theres is a cut-cell ICC then CUT_CELL(ICC)%UNKH(1) has been populated.
+IF (ICC > 0) IROW     = CUT_CELL(ICC)%UNKH(1) - ZONE_SOLVE(IPZ)%UNKH_IND(NM_START)
 
 RETURN
 END SUBROUTINE GET_CC_IROW
@@ -21836,14 +21695,14 @@ SUBROUTINE NUMBER_UNKH_CUTCELLS(FLAG12,NM,IPZ,NUNKH_LC)
    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
-      IF(.NOT.PRES_ON_WHOLE_DOMAIN .AND. ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE
+      IF(ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE
       NUNKH_LC(NM) = NUNKH_LC(NM) + 1
       CUT_CELL(ICC)%UNKH(1:CC%NCELL) = NUNKH_LC(NM)
    ENDDO
 ELSE
    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
-      IF(.NOT.PRES_ON_WHOLE_DOMAIN .AND. ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE
+      IF(ZONE_SOLVE(PRESSURE_ZONE(I,J,K))%CONNECTED_ZONE_PARENT/=IPZ ) CYCLE
       DO JCC=1,CC%NCELL
          CUT_CELL(ICC)%UNKH(JCC) = CUT_CELL(ICC)%UNKH(JCC) + ZONE_SOLVE(IPZ)%UNKH_IND(NM)
       ENDDO
@@ -21885,15 +21744,10 @@ SUBROUTINE COPY_CC_HS_TO_UNKH(NM)
 
    ICC=CCVAR(II,JJ,KK,CC_IDCC)
 
-   IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
-      IF (ICC > 0) THEN ! Cut-cells on this guard-cell Cartesian cell.
-         MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = INT(OM%HS(IIO,JJO,KKO))
-      ELSE
-         MESHES(NM)%CCVAR(II,JJ,KK,CC_UNKH) = INT(OM%HS(IIO,JJO,KKO))
-      ENDIF
+   IF (ICC > 0) THEN ! Cut-cells on this guard-cell Cartesian cell.
+      MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = INT(OM%HS(IIO,JJO,KKO))
    ELSE
       MESHES(NM)%CCVAR(II,JJ,KK,CC_UNKH) = INT(OM%HS(IIO,JJO,KKO))
-      IF (ICC > 0) MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = MESHES(NM)%CCVAR(II,JJ,KK,CC_UNKH)
    ENDIF
 
 ENDDO EXTERNAL_WALL_LOOP
@@ -21920,21 +21774,12 @@ SUBROUTINE COPY_CC_UNKH_TO_HS(NM)
 ! Local Variables:
 INTEGER :: I,J,K,ICC
 
-IF(.NOT.PRES_ON_WHOLE_DOMAIN) THEN
-   DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
-      I = MESHES(NM)%CUT_CELL(ICC)%IJK(IAXIS)
-      J = MESHES(NM)%CUT_CELL(ICC)%IJK(JAXIS)
-      K = MESHES(NM)%CUT_CELL(ICC)%IJK(KAXIS)
-      HS(I,J,K)= REAL(MESHES(NM)%CUT_CELL(ICC)%UNKH(1),EB)
-   ENDDO
-ELSE
-   DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
-      I = MESHES(NM)%CUT_CELL(ICC)%IJK(IAXIS)
-      J = MESHES(NM)%CUT_CELL(ICC)%IJK(JAXIS)
-      K = MESHES(NM)%CUT_CELL(ICC)%IJK(KAXIS)
-      MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = INT(HS(I,J,K))
-   ENDDO
-ENDIF
+DO ICC=1,MESHES(NM)%N_CUTCELL_MESH
+   I = MESHES(NM)%CUT_CELL(ICC)%IJK(IAXIS)
+   J = MESHES(NM)%CUT_CELL(ICC)%IJK(JAXIS)
+   K = MESHES(NM)%CUT_CELL(ICC)%IJK(KAXIS)
+   HS(I,J,K)= REAL(MESHES(NM)%CUT_CELL(ICC)%UNKH(1),EB)
+ENDDO
 
 RETURN
 END SUBROUTINE COPY_CC_UNKH_TO_HS

diff --git a/Source/geom.f90 b/Source/geom.f90
@@ -4745,8 +4745,8 @@ SUBROUTINE GET_GEOM_TRIBIN
    DO X1AXIS=IAXIS,KAXIS
 
       ! Here reduce the X1_LOW to X1_HIGH distance to the smallest of FDS Mesh and connected meshes BBOX or Geometry:
-      MIN_MESHGEOM = MAX(MINMAX_MESHES( LOW_IND,X1AXIS),G%GEOM_BOX( LOW_IND,X1AXIS))
-      MAX_MESHGEOM = MIN(MINMAX_MESHES(HIGH_IND,X1AXIS),G%GEOM_BOX(HIGH_IND,X1AXIS))
+      MIN_MESHGEOM = MAX(MINMAX_MESHES( LOW_IND,X1AXIS),G%GEOM_BOX( LOW_IND,X1AXIS)-G%MEAN_LEDGE)
+      MAX_MESHGEOM = MIN(MINMAX_MESHES(HIGH_IND,X1AXIS),G%GEOM_BOX(HIGH_IND,X1AXIS)+G%MEAN_LEDGE)
       LX1 = MAX_MESHGEOM - MIN_MESHGEOM
 
       ! Define number of bins in direction X1AXIS: