Merge pull request #14707 from drjfloyd/master

FDS Source: Fix terrain VENT error and SMV write error for GEOM

Author: drjfloyd

Source/dump.f90

@@ -2349,6 +2349,10 @@ SUBROUTINE WRITE_SMOKEVIEW_FILE
          ENDIF
       ENDIF
 
+      ! Render GEOM vents using SURF properties (that is, do not draw GEOM vents in Smokeview)
+
+      IF (VT%GEOM) VT%TYPE_INDICATOR = -2
+
    ENDDO VENT_LOOP
 
    ! Look for interpolated mesh boundaries and ensure that Smokeview leaves these blank (VENT_INDICES=-1).

Source/geom.f90

@@ -4121,11 +4121,13 @@ SUBROUTINE GET_EXT_INB_CUTFACES_TO_CFACE
 
 ! Local Variables:
 INTEGER :: ICF, CFACE_INDEX_LOCAL, SURF_INDEX
+INTEGER :: IVENT
 REAL(EB):: ADDMAT(IAXIS:KAXIS,LOW_IND:HIGH_IND)
 
 ! GET_CUTCELLS_VERBOSE variables:
 INTEGER, ALLOCATABLE, DIMENSION(:) :: NCFACE_BY_MESH
 
+TYPE(VENTS_TYPE), POINTER :: VT
 TYPE(CFACE_TYPE), POINTER :: CFA
 
 IF(GET_CUTCELLS_VERBOSE) CALL CPU_TIME(CPUTIME_START)
@@ -4308,10 +4310,50 @@ SUBROUTINE GET_EXT_INB_CUTFACES_TO_CFACE
    MESHES(NM)%N_INTERNAL_CFACE_CELLS   = CFACE_INDEX_LOCAL - MESHES(NM)%INTERNAL_CFACE_CELLS_LB
 ENDDO MESH_LOOP_1
 
+! Second loop, apply VENTS to change SURF_ID associated with CFACEs:
+MESH_LOOP_2 : DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
+   CALL POINT_TO_MESH(NM)
+
+   ! ! Currently : Modify CFACE SURF_INDEX with VENT information: This needs more development.
+
+   VENT_LOOP : DO IVENT=1,MESHES(NM)%N_VENT
+      VT => VENTS(IVENT)
+      IF(.NOT.VT%GEOM) CYCLE VENT_LOOP ! Do not apply vent to Geometries.
+
+      ! This test is a simplified test for VENTS changing the CFACE SURF_ID to VENT SURF_ID for all CFACEs whose
+      ! centroid locations lay within the frame of the IOR grid aligned VENT:
+      ADDMAT = 0._EB;
+      SELECT CASE(ABS(VT%IOR))
+      CASE(IAXIS)
+         ADDMAT(IAXIS,LOW_IND)  = -(XF_MAX-XS_MIN) ! -DX(VT%I1) Set normal size to 2 times domain size.
+         ADDMAT(IAXIS,HIGH_IND) =  (XF_MAX-XS_MIN) !  DX(VT%I2) XF_MAX, etc. defined in cons.f90.
+      CASE(JAXIS)
+         ADDMAT(JAXIS,LOW_IND)  = -(YF_MAX-YS_MIN) ! -DY(VT%J1)
+         ADDMAT(JAXIS,HIGH_IND) =  (YF_MAX-YS_MIN) !  DY(VT%J2)
+      CASE(KAXIS)
+         ADDMAT(KAXIS,LOW_IND)  = -(ZF_MAX-ZS_MIN) ! -DZ(VT%K1)
+         ADDMAT(KAXIS,HIGH_IND) =  (ZF_MAX-ZS_MIN) !  DZ(VT%K2)
+      END SELECT
+      ! CFACE Loop to modify SURF_INDEX in INTERNAL_CFACE_CELLS:
+      CFACE_LOOP_2 : DO CFACE_INDEX_LOCAL=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
+         CFA => CFACE(CFACE_INDEX_LOCAL)
+         BC => BOUNDARY_COORD(CFA%BC_INDEX)
+         IF (BC%X < X(VT%I1)+ADDMAT(IAXIS,LOW_IND )) CYCLE CFACE_LOOP_2
+         IF (BC%X > X(VT%I2)+ADDMAT(IAXIS,HIGH_IND)) CYCLE CFACE_LOOP_2
+         IF (BC%Y < Y(VT%J1)+ADDMAT(JAXIS,LOW_IND )) CYCLE CFACE_LOOP_2
+         IF (BC%Y > Y(VT%J2)+ADDMAT(JAXIS,HIGH_IND)) CYCLE CFACE_LOOP_2
+         IF (BC%Z < Z(VT%K1)+ADDMAT(KAXIS,LOW_IND )) CYCLE CFACE_LOOP_2
+         IF (BC%Z > Z(VT%K2)+ADDMAT(KAXIS,HIGH_IND)) CYCLE CFACE_LOOP_2
+         CFA%VENT_INDEX = IVENT
+         CFA%SURF_INDEX = VT%SURF_INDEX
+      ENDDO CFACE_LOOP_2
+   ENDDO VENT_LOOP
+ENDDO MESH_LOOP_2
 ! - At this pont all final values of SURF_INDEX have been given to CFACEs.
 
 ! Third loop, 1. Compute final FDS area integrals by SURF_ID and GEOM.
-!             2. Compute input areas by SURF_ID and GEOM. First sum over GEOM FACES SURF_IDs.
+!             2. Compute input areas by SURF_ID and GEOM. First sum over GEOM FACES SURF_IDs,
+!                then VENTS input surfaces are assigned to corresponding GEOMs and SURF_IDs if present (VENTs take precedence).
 IF(N_GEOMETRY>0) THEN
    ALLOCATE(FDS_AREA_GEOM(0:N_SURF,N_GEOMETRY)); FDS_AREA_GEOM = 0._EB
 ENDIF
@@ -8835,16 +8877,15 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
    INS_INB_COND_1 : IF (IS_INB) THEN
       B1%VEL_ERR_NEW=CF%VEL(IFACE) - 0._EB ! Assumes zero velocity of solid.
 
-      IBOD           = CF%BODTRI(1,IFACE)
-      IWSEL          = CF%BODTRI(2,IFACE)
-
       ! Normal to cut-face:
       V2(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(2,IFACE))-CF%XYZCEN(IAXIS:KAXIS,IFACE)
       V3(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(3,IFACE))-CF%XYZCEN(IAXIS:KAXIS,IFACE)
       CALL CROSS_PRODUCT(BC%NVEC(IAXIS:KAXIS),V2,V3)
       IF(NORM2(BC%NVEC)>TWO_EPSILON_EB .AND. CF%CFACE_ORIGIN(IFACE)==BLOCKED_SPLIT_CELL) THEN
          BC%NVEC(IAXIS:KAXIS) = BC%NVEC(IAXIS:KAXIS)/NORM2(BC%NVEC)
       ELSE
+         IBOD =CF%BODTRI(1,IFACE)
+         IWSEL=CF%BODTRI(2,IFACE)
          BC%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)
       ENDIF
       X1AXIS = MAXLOC(ABS(BC%NVEC(IAXIS:KAXIS)),DIM=1)
@@ -8878,6 +8919,7 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
          ! External mesh boundary CFACES inherit the underlaying WALL type.
          CFA%BOUNDARY_TYPE = WC%BOUNDARY_TYPE
          CFA%NODE_INDEX    = SURFACE(WC%SURF_INDEX)%NODE_INDEX
+         CFA%VENT_INDEX    = WC%VENT_INDEX
 
          BC%II = WC_BC%II
          BC%JJ = WC_BC%JJ
@@ -23714,10 +23756,10 @@ SUBROUTINE READ_GEOM
 
 DEALLOCATE(GEOM_LINE)
 
-CC_IBM = .TRUE.
-
 IF( (T_END-T_BEGIN) < TWO_EPSILON_EB) RETURN
 
+CC_IBM = .TRUE.
+
 ! If unstructured projection defined set Pressure solver on unstructured grid.
 IF (PRES_FLAG/=UGLMAT_FLAG) THEN
    PRES_METHOD = 'ULMAT'

Source/hvac.f90

@@ -1038,7 +1038,7 @@ SUBROUTINE PROC_HVAC
 
 ENDDO NODE_LOOP
 
-GEOM_IF: IF (CC_IBM .AND. ANY(SURFACE%NODE_INDEX>0)) THEN
+GEOM_IF: IF (N_GEOMETRY > 0 .AND. ANY(SURFACE%NODE_INDEX>0)) THEN
    CF_AREA = 0._EB
    CF_X = 0._EB
    CF_Y = 0._EB

Source/init.f90

@@ -2854,6 +2854,7 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR)
 VENT_SEARCH_LOOP: DO N=1,M%N_VENT
 
    VT => M%VENTS(N)
+   IF (VT%GEOM) CYCLE VENT_SEARCH_LOOP
    IF (OBST_INDEX>0) THEN
       IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY)           CYCLE VENT_SEARCH_LOOP
       IF (.NOT.M%OBSTRUCTION(OBST_INDEX)%ALLOW_VENT) CYCLE VENT_SEARCH_LOOP

Source/read.f90

@@ -11844,14 +11844,14 @@ SUBROUTINE READ_VENT
 CHARACTER(LABEL_LENGTH) :: ID,DEVC_ID,CTRL_ID,SURF_ID,PRESSURE_RAMP,TMP_EXTERIOR_RAMP,MULT_ID,OBST_ID
 CHARACTER(25) :: COLOR
 TYPE(MULTIPLIER_TYPE), POINTER :: MR
-LOGICAL :: REJECT_VENT,OUTLINE,SOLID_FOUND,AREA_ADJUST,BLOCKED
+LOGICAL :: REJECT_VENT,OUTLINE,GEOM,SOLID_FOUND,AREA_ADJUST,BLOCKED
 TYPE IMPLICIT_VENT_TYPE
    REAL(EB) :: XB(6)
    INTEGER, DIMENSION(3) :: RGB=-1
    CHARACTER(LABEL_LENGTH) :: MB='null',SURF_ID='null',ID='null'
 END TYPE
 TYPE(IMPLICIT_VENT_TYPE), ALLOCATABLE, DIMENSION(:) :: IMPLICIT_VENT
-NAMELIST /VENT/ AREA_ADJUST,COLOR,CTRL_ID,DB,DEVC_ID,DYNAMIC_PRESSURE,FYI,ID,IOR,L_EDDY,L_EDDY_IJ, &
+NAMELIST /VENT/ AREA_ADJUST,COLOR,CTRL_ID,DB,DEVC_ID,DYNAMIC_PRESSURE,FYI,GEOM,ID,IOR,L_EDDY,L_EDDY_IJ, &
                 MB,MULT_ID,N_EDDY,OBST_ID,OUTLINE,PBX,PBY,PBZ,PRESSURE_RAMP,RADIUS,REYNOLDS_STRESS, &
                 RGB,SPREAD_RATE,SURF_ID,TEXTURE_ORIGIN,TMP_EXTERIOR,TMP_EXTERIOR_RAMP,TRANSPARENCY, &
                 UVW,VEL_RMS,XB,XYZ
@@ -12134,6 +12134,10 @@ SUBROUTINE READ_VENT
                   IF (BLOCKED) REJECT_VENT = .TRUE.
                ENDIF
 
+               ! If the VENT is on a GEOM do not reject (further setup in READ_GEOM)
+
+               IF (GEOM .AND. .NOT.(TERRAIN_CASE .AND. ALL(XB(1:6)>-1.01E6_EB))) REJECT_VENT = .FALSE.
+
                ! If the VENT is rejected, cycle
 
                IF (REJECT_VENT) THEN
@@ -12352,6 +12356,8 @@ SUBROUTINE READ_VENT
                   VT%UVW = VT%UVW/SQRT(VT%UVW(1)**2+VT%UVW(2)**2+VT%UVW(3)**2)
                ENDIF
 
+               VT%GEOM = GEOM
+
             ENDDO I_MULT_LOOP
          ENDDO J_MULT_LOOP
       ENDDO K_MULT_LOOP
@@ -12437,7 +12443,7 @@ SUBROUTINE READ_VENT
          ENDDO
       ENDIF
 
-      IF (VT%IOR==0) THEN
+      IF (VT%IOR==0 .AND. .NOT. VT%GEOM) THEN
          WRITE(MESSAGE,'(3A)')  'ERROR(818): VENT ',TRIM(VT%ID),' requires an orientation index, IOR.'
          CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.) ; RETURN
       ENDIF
@@ -12485,9 +12491,17 @@ SUBROUTINE READ_VENT
 
       ! Check UVW
 
-      IF (ABS(VT%UVW(ABS(VT%IOR))) < TWO_EPSILON_EB) THEN
-         WRITE(MESSAGE,'(3A)')  'ERROR(821): VENT ',TRIM(VT%ID),' cannot have normal component of UVW equal to 0.'
-         CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.) ; RETURN
+      IF (.NOT.VT%GEOM) THEN
+         IF (ABS(VT%UVW(ABS(VT%IOR))) < TWO_EPSILON_EB) THEN
+            WRITE(MESSAGE,'(3A)')  'ERROR(821): VENT ',TRIM(VT%ID),' cannot have normal component of UVW equal to 0.'
+            CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.) ; RETURN
+         ENDIF
+      ENDIF
+
+      ! Special treatment for coloring GEOM surface with HVAC_BOUNDARY
+
+      IF (VT%GEOM .AND. VT%BOUNDARY_TYPE==HVAC_BOUNDARY) THEN
+         SURFACE(VT%SURF_INDEX)%RGB = (/0,128,0/) ! green
       ENDIF
 
    ENDDO VENT_LOOP_2
@@ -12560,6 +12574,7 @@ SUBROUTINE SET_VENT_DEFAULTS
 DB                = 'null'
 DEVC_ID           = 'null'
 DYNAMIC_PRESSURE  = 0._EB
+GEOM              = .FALSE.
 ID                = 'null'
 IOR               = 0
 L_EDDY            = 0._EB

Source/type.f90

@@ -1339,6 +1339,7 @@ MODULE TYPES
    INTEGER :: BR_INDEX=0                 !< Derived type carrying angular-specific radiation intensities
    INTEGER :: SURF_INDEX=0
    INTEGER :: NODE_INDEX=0
+   INTEGER :: VENT_INDEX=0
    INTEGER :: BOUNDARY_TYPE=0
    INTEGER :: CUT_FACE_IND1=-11          !< First index pointing to CUT_FACE array for this CFACE.
    INTEGER :: CUT_FACE_IND2=-11          !< Second index pointing to CUT_FACE array for this CFACE.
@@ -1496,7 +1497,7 @@ MODULE TYPES
                X1_ORIG=0._EB,X2_ORIG=0._EB,Y1_ORIG=0._EB,Y2_ORIG=0._EB,Z1_ORIG=0._EB,Z2_ORIG=0._EB, &
                X0=-9.E6_EB,Y0=-9.E6_EB,Z0=-9.E6_EB,FIRE_SPREAD_RATE,UNDIVIDED_INPUT_AREA=0._EB,INPUT_AREA=0._EB,&
                TMP_EXTERIOR=-1000._EB,DYNAMIC_PRESSURE=0._EB,UVW(3)=-1.E12_EB,RADIUS=-1._EB
-   LOGICAL :: ACTIVATED=.TRUE.,AREA_ADJUST=.TRUE.,DRAW=.TRUE.
+   LOGICAL :: ACTIVATED=.TRUE.,GEOM=.FALSE.,AREA_ADJUST=.TRUE.,DRAW=.TRUE.
    CHARACTER(LABEL_LENGTH) :: DEVC_ID='null',CTRL_ID='null',ID='null'
    ! turbulent inflow (experimental)
    INTEGER :: N_EDDY=0