FDS Source: Automatically search for HT3D mesh neighborhoods

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -2445,7 +2445,7 @@ \subsubsection{Limitations}
 \begin{enumerate}
 \item {\ct HT3D} cannot be applied to an exterior boundary; it must be applied to an {\ct OBST} for which at least one face in each coordinate direction is exposed.
 \item Avoid contact between 3-D and 1-D solids. If two sides of a 3-D solid touch 1-D solids, there will be no lateral heat conduction computed in that particular direction.
-\item If your 3-D obstruction extends beyond meshes that abut, add the parameter {\ct NEIGHBOR\_SEPARATION\_DISTANCE} to the {\ct MISC} line. Any mesh within this distance of another mesh will share geometry information for use in the 3-D heat conduction calculation.
+\item If your 3-D obstruction extends beyond meshes that abut, FDS uses a special algorithm to identify all the meshes where this obstruction lives, and also those obstructions connected to it. However, if this algorithm fails to detect all the meshes and you recieve an error stating that there is a problem, add the parameter {\ct NEIGHBOR\_SEPARATION\_DISTANCE} to the {\ct MISC} line. Any mesh within this distance of another mesh will share geometry information for use in the 3-D heat conduction calculation. If you set this parameter to a value larger than the width of the computational domain, then all meshes will establish communication channels for exchanging boundary information.
 \item By default, the interior nodes are clustered near the surface and stretched out deeper within the solid. If you want to maintain uniform spacing, set {\ct CELL\_SIZE} on the {\ct SURF} line to indicate the desired interior node spacing. The {\ct CELL\_SIZE} is typically the same as the gas phase cells.
 \end{enumerate}
 

Source/read.f90

@@ -49,6 +49,14 @@ MODULE READ_INPUT
 TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
 TYPE(REACTION_TYPE), POINTER :: RN=>NULL()
 LOGICAL :: RETURN_BEFORE_STOP_FILE=.FALSE., RETURN_BEFORE_SIM_MODE=.FALSE.
+INTEGER :: N_HT3D_SURF_LINES=0,N_HT3D_OBST=0
+CHARACTER(LABEL_LENGTH), DIMENSION(20) :: HT3D_SURF_LIST='null'
+TYPE HT3D_OBST_TYPE
+   INTEGER :: GROUP_INDEX=0
+   REAL(EB) :: XS,XF,YS,YF,ZS,ZF
+END TYPE
+TYPE(HT3D_OBST_TYPE), DIMENSION(:), ALLOCATABLE, TARGET :: HT3D_OBST
+
 
 CONTAINS
 
@@ -103,6 +111,7 @@ SUBROUTINE READ_DATA(DT)
 CALL READ_MOVE    ; CALL CHECK_STOP_STATUS ; IF (STOP_STATUS/=NO_STOP) RETURN
 CALL READ_MULT    ; CALL CHECK_STOP_STATUS ; IF (STOP_STATUS/=NO_STOP) RETURN
 CALL READ_MESH    ; CALL CHECK_STOP_STATUS ; IF (STOP_STATUS/=NO_STOP) RETURN
+CALL PROC_MESH    ; CALL CHECK_STOP_STATUS ; IF (STOP_STATUS/=NO_STOP) RETURN
 CALL READ_TRAN    ; CALL CHECK_STOP_STATUS ; IF (STOP_STATUS/=NO_STOP) RETURN
 CALL READ_TIME(DT); CALL CHECK_STOP_STATUS ; IF (STOP_STATUS/=NO_STOP) RETURN
 CALL READ_PRES    ; CALL CHECK_STOP_STATUS ; IF (STOP_STATUS/=NO_STOP) RETURN
@@ -539,15 +548,13 @@ END FUNCTION READ_STOP
 
 SUBROUTINE READ_MESH
 
-USE COMP_FUNCTIONS, ONLY : SEARCH_INPUT_FILE
-INTEGER :: IJK(3),NM,NM2,CURRENT_MPI_PROCESS,MPI_PROCESS,RGB(3),N_MESH_NEW,N,II,JJ,KK,NMESHES_READ,NNN,JBAR_OLD_VALUE
-INTEGER, ALLOCATABLE, DIMENSION(:) :: NEIGHBOR_LIST
-LOGICAL :: OVERLAPPING_X,OVERLAPPING_Y,OVERLAPPING_Z,POSSIBLY_PERIODIC,CYLINDRICAL_OLD_VALUE,PERIODIC_FOUND_IN_FILE
+INTEGER :: IJK(3),NM,CURRENT_MPI_PROCESS,MPI_PROCESS,RGB(3),N_MESH_NEW,N,II,JJ,KK,NMESHES_READ,NNN,JBAR_OLD_VALUE
+LOGICAL :: CYLINDRICAL_OLD_VALUE
 REAL(EB) :: XB1,XB2,XB3,XB4,XB5,XB6
 CHARACTER(25) :: COLOR
 CHARACTER(LABEL_LENGTH) :: MULT_ID,TRNX_ID,TRNY_ID,TRNZ_ID
 NAMELIST /MESH/ CHECK_MESH_ALIGNMENT,COLOR,CYLINDRICAL,FYI,ID,IJK,MPI_PROCESS,MULT_ID,RGB,TRNX_ID,TRNY_ID,TRNZ_ID,XB
-TYPE (MESH_TYPE), POINTER :: M,M2
+TYPE (MESH_TYPE), POINTER :: M
 TYPE (MULTIPLIER_TYPE), POINTER :: MR
 
 NMESHES = 0
@@ -842,29 +849,117 @@ SUBROUTINE READ_MESH
    ENDIF
 ENDDO
 
-! Determine mesh neighbors. MESH_SEPARATION_DISTANCE is a very small length
-! used to determine if there are periodic boundaries. NEIGHBOR_SEPARATION_DISANCE
-! is the distance beyond which no information or message passing is assumed
-! between the meshes. Its value is deliberately complicated to avoid having two
-! meshes separated by exactly that distance.
+END SUBROUTINE READ_MESH
+
+
+!> \brief Determine mesh neighborhoods for MPI communications
+
+SUBROUTINE PROC_MESH
+
+USE COMP_FUNCTIONS, ONLY : SEARCH_INPUT_FILE
+INTEGER :: NM,NM2,I,II,III,N_GROUPS=0
+INTEGER, ALLOCATABLE, DIMENSION(:) :: NEIGHBOR_LIST
+LOGICAL :: OVERLAPPING_X,OVERLAPPING_Y,OVERLAPPING_Z,POSSIBLY_PERIODIC,PERIODIC_FOUND_IN_FILE
+TYPE (MESH_TYPE), POINTER :: M,M2
+TYPE (HT3D_OBST_TYPE), POINTER :: HO,HO2,HO3
+TYPE MESH_COMM_TYPE
+   INTEGER :: N_MESHES=0
+   INTEGER, DIMENSION(100) :: LIST=0
+END TYPE MESH_COMM_TYPE
+TYPE(MESH_COMM_TYPE), ALLOCATABLE, DIMENSION(:), TARGET :: MESH_COMM
+TYPE(MESH_COMM_TYPE), POINTER :: MC
+
+! Read the SURF lines and make a list, HT3D_SURF_LIST, of those that are HT3D or VARIABLE_THICKNESS
+
+CALL READ_SURF(QUICK_READ=.TRUE.)
+
+! If there are HT3D solids, determine the indices of the meshes that contain connected HT3D OBSTs
+
+IF_HT3D: IF (N_HT3D_SURF_LINES>0) THEN
+
+   ! Read the OBST lines and make a list, HT3D_OBST, of those that have HT3D or VARIABLE_THICKNESS SURF lines
+
+   ALLOCATE(HT3D_OBST(100))  ! Initial allocation -- this array can be reallocated if needed
+   CALL READ_OBST(QUICK_READ=.TRUE.)
+
+   ! Assign each HT3D_OBST an integer, GROUP_INDEX, which indicates which group of obstructions it is connected to.
+   ! N_GROUPS is the number of these blocks of connected OBSTs.
+
+   N_GROUPS = 1
+
+   HT3D_OBST_LOOP: DO I=1,N_HT3D_OBST
+      HO => HT3D_OBST(I)
+      IF (I==1) THEN
+         HO%GROUP_INDEX = 1
+      ELSE
+         HO%GROUP_INDEX = 1000000
+         DO II=I-1,1,-1
+            HO2 => HT3D_OBST(II)
+            IF (HO%XS>HO2%XF.OR.HO%XF<HO2%XS.OR.HO%YS>HO2%YF.OR.HO%YF<HO2%YS.OR.HO%ZS>HO2%ZF.OR.HO%ZF<HO2%ZS) CYCLE
+            IF (HO2%GROUP_INDEX>HO%GROUP_INDEX) THEN  ! If current OBST is connected to one already assigned a GROUP_INDEX, 
+                                                      ! and that GROUP_INDEX is greater than that of the OBST, reduce the 
+                                                      ! N_GROUPS and GROUP_INDEXs of OBSTs with higher values
+               N_GROUPS = N_GROUPS - 1
+               DO III=1,I-1
+                  HO3 => HT3D_OBST(III)
+                  IF (HO3%GROUP_INDEX>HO%GROUP_INDEX) HO3%GROUP_INDEX = HO3%GROUP_INDEX - 1
+               ENDDO
+            ENDIF
+            HO%GROUP_INDEX = MIN(HO%GROUP_INDEX,HO2%GROUP_INDEX)
+            IF (HO%GROUP_INDEX==N_GROUPS) CYCLE HT3D_OBST_LOOP
+         ENDDO
+         IF (HO%GROUP_INDEX>100000) THEN  ! The current OBST is not connected to any previous OBST. Assign it a new GROUP_INDEX.
+            N_GROUPS = N_GROUPS + 1
+            HO%GROUP_INDEX = N_GROUPS
+         ENDIF
+      ENDIF
+   ENDDO HT3D_OBST_LOOP
+
+   ! Now determine all the MESH indices corresponding to each group of connected OBSTs
+
+   ALLOCATE(MESH_COMM(N_GROUPS))
+
+   DO I=1,N_HT3D_OBST
+      HO => HT3D_OBST(I)
+      MC => MESH_COMM(HO%GROUP_INDEX)
+      MESH_LOOP: DO NM=1,NMESHES
+         M => MESHES(NM)
+         IF (HO%XS>M%XF.OR.HO%XF<M%XS.OR.HO%YS>M%YF.OR.HO%YF<M%YS.OR.HO%ZS>M%ZF.OR.HO%ZF<M%ZS) CYCLE MESH_LOOP
+         DO II=1,MC%N_MESHES
+            IF (NM==MC%LIST(II)) CYCLE MESH_LOOP
+         ENDDO
+         MC%N_MESHES = MC%N_MESHES + 1
+         MC%LIST(MC%N_MESHES) = NM
+      ENDDO MESH_LOOP
+   ENDDO
+
+   DEALLOCATE(HT3D_OBST)
+
+ENDIF IF_HT3D
+
+! MESH_SEPARATION_DISTANCE is a very small length used to determine if there are periodic boundaries. NEIGHBOR_SEPARATION_DISANCE
+! is the distance beyond which no information or message passing is assumed between the meshes. Its value is deliberately 
+! complicated to avoid having two meshes separated by exactly that same distance.
 
 MESH_SEPARATION_DISTANCE = MIN(1.E-3_EB,0.05_EB*CHARACTERISTIC_CELL_SIZE)
 IF (NEIGHBOR_SEPARATION_DISTANCE<0._EB) NEIGHBOR_SEPARATION_DISTANCE = 4.56789_EB*CHARACTERISTIC_CELL_SIZE
 
-
 ! Search through the input file for any mention of the word PERIODIC. If not found, this simplifies neighbor selection.
 
 CALL SEARCH_INPUT_FILE(LU_INPUT,'PERIODIC',PERIODIC_FOUND_IN_FILE)
 
-! For MESHES controlled by a given MPI process, only allocate other MESHES that
-! are "close" as defined by the two parameters above.
+! For MESHES controlled by a given MPI process, only allocate other MESHES that are "close" as defined by the two parameters above.
 
 ALLOCATE(NEIGHBOR_LIST(10000))
 
 DO NM=1,NMESHES
+
    M => MESHES(NM)
    M%N_NEIGHBORING_MESHES = 0
    NEIGHBOR_LIST = 0
+
+   ! Add adjacent meshes to the neighborhood of MESH NM
+
    DO NM2=1,NMESHES
       M2 => MESHES(NM2)
       OVERLAPPING_X = .TRUE.
@@ -888,16 +983,31 @@ SUBROUTINE READ_MESH
       M%N_NEIGHBORING_MESHES = M%N_NEIGHBORING_MESHES + 1
       NEIGHBOR_LIST(M%N_NEIGHBORING_MESHES) = NM2
    ENDDO
+
+   ! Add meshes containing the HT3D_OBST groups to the neighborhood of MESH NM
+
+   DO I=1,N_GROUPS
+      MC => MESH_COMM(I)
+      IF (ANY(MC%LIST==NM,DIM=1)) THEN
+         DO II=1,MC%N_MESHES
+            IF (ANY(NEIGHBOR_LIST==MC%LIST(II))) CYCLE
+            M%N_NEIGHBORING_MESHES = M%N_NEIGHBORING_MESHES + 1
+            NEIGHBOR_LIST(M%N_NEIGHBORING_MESHES) = MC%LIST(II)
+         ENDDO
+      ENDIF
+   ENDDO
+
+   ! Save the list of neighboring meshes into an array
+
    ALLOCATE(M%NEIGHBORING_MESH(M%N_NEIGHBORING_MESHES))
    DO I=1,M%N_NEIGHBORING_MESHES
       M%NEIGHBORING_MESH(I) = NEIGHBOR_LIST(I)
    ENDDO
 ENDDO
-DEALLOCATE(NEIGHBOR_LIST)
 
-REWIND(LU_INPUT) ; INPUT_FILE_LINE_NUMBER = 0
+DEALLOCATE(NEIGHBOR_LIST)
 
-END SUBROUTINE READ_MESH
+END SUBROUTINE PROC_MESH
 
 
 !> \brief Read the TRAN namelist lines and compute the polynomial transform function for the vertical coordinate
@@ -7051,10 +7161,11 @@ END SUBROUTINE PROC_MATL
 
 !> \brief Read the SURF namelist lines
 
-SUBROUTINE READ_SURF
+SUBROUTINE READ_SURF(QUICK_READ)
 
 USE MATH_FUNCTIONS, ONLY : GET_RAMP_INDEX
 USE DEVICE_VARIABLES, ONLY : PROPERTY_TYPE
+LOGICAL, INTENT(IN), OPTIONAL :: QUICK_READ
 CHARACTER(LABEL_LENGTH) :: PART_ID,RAMP_MF(MAX_SPECIES),RAMP_Q,RAMP_V,RAMP_T,RAMP_T_I,MATL_ID(MAX_LAYERS,MAX_MATERIALS),&
                  PROFILE,BACKING,GEOMETRY,RAMP_EF,RAMP_PART,NAME_LIST(MAX_MATERIALS*MAX_LAYERS),&
                  SPEC_ID(MAX_SPECIES),RAMP_TMP_BACK,RAMP_TMP_GAS_BACK,RAMP_TMP_GAS_FRONT,&
@@ -7130,18 +7241,28 @@ SUBROUTINE READ_SURF
    CALL CHECKREAD('SURF',LU_INPUT,IOS)  ; IF (STOP_STATUS==SETUP_STOP) RETURN
    IF (IOS==1) EXIT COUNT_SURF_LOOP
    ID = 'null'
+   HT3D = .FALSE.
+   VARIABLE_THICKNESS = .FALSE.
    READ(LU_INPUT,SURF,ERR=34,IOSTAT=IOS)
    IF (ID=='null') THEN
       WRITE(MESSAGE,'(A)') 'ERROR: SURF line must have an ID'
       CALL SHUTDOWN(MESSAGE) ; RETURN
    ENDIF
+   IF (PRESENT(QUICK_READ) .AND. (HT3D .OR. VARIABLE_THICKNESS)) THEN
+      N_HT3D_SURF_LINES = N_HT3D_SURF_LINES + 1
+      HT3D_SURF_LIST(N_HT3D_SURF_LINES) = ID
+   ENDIF
    N_SURF = N_SURF + 1
    34 IF (IOS>0) THEN
          WRITE(MESSAGE,'(A,I0,A,I0)') 'ERROR: Problem with SURF number ',N_SURF+1,', line number ',INPUT_FILE_LINE_NUMBER
          CALL SHUTDOWN(MESSAGE) ; RETURN
       ENDIF
 ENDDO COUNT_SURF_LOOP
 
+! Special case where SURF lines are scanned, looking for presence of HT3D
+
+IF (PRESENT(QUICK_READ)) RETURN
+
 ! Allocate the SURFACE derived type, leaving space for SURF entries not defined explicitly by the user
 
 N_SURF_RESERVED = 9
@@ -9532,16 +9653,17 @@ END SUBROUTINE READ_TABL
 
 !> \brief Read the OBSTruction namelist lines
 
-SUBROUTINE READ_OBST
+SUBROUTINE READ_OBST(QUICK_READ)
 
 USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_CELL
 USE GEOMETRY_FUNCTIONS, ONLY: BLOCK_CELL,CIRCLE_CELL_INTERSECTION_AREA,SEARCH_OTHER_MESHES
 USE COMPLEX_GEOMETRY, ONLY: INTERSECT_CONE_AABB,INTERSECT_CYLINDER_AABB,INTERSECT_SPHERE_AABB,INTERSECT_OBB_AABB,ROTATION_MATRIX
 USE MISC_FUNCTIONS, ONLY: PROCESS_MESH_NEIGHBORHOOD
+LOGICAL, INTENT(IN), OPTIONAL :: QUICK_READ
 TYPE(OBSTRUCTION_TYPE), POINTER :: OB2,OBT
 TYPE(MULTIPLIER_TYPE), POINTER :: MR
 TYPE(OBSTRUCTION_TYPE), DIMENSION(:), ALLOCATABLE, TARGET :: TEMP_OBSTRUCTION
-INTEGER :: NM,NOM,NOM2,N_OBST_O,IC,N,NN,NNN,NNNN,N_NEW_OBST,RGB(3),N_OBST_DIM,II,JJ,KK,MULT_INDEX,SHAPE_TYPE,IIO,JJO,KKO,IOR
+INTEGER :: I,NM,NOM,NOM2,N_OBST_O,IC,N,NN,NNN,NNNN,N_NEW_OBST,RGB(3),N_OBST_DIM,II,JJ,KK,MULT_INDEX,SHAPE_TYPE,IIO,JJO,KKO,IOR
 CHARACTER(LABEL_LENGTH) :: ID,DEVC_ID,SHAPE,SURF_ID,SURF_ID_INTERIOR,SURF_IDS(3),SURF_ID6(6),CTRL_ID,MULT_ID,MATL_ID(MAX_MATERIALS)
 CHARACTER(25) :: COLOR
 LOGICAL :: OVERLAY,IS_INTERSECT
@@ -9556,11 +9678,11 @@ SUBROUTINE READ_OBST
 
 MESH_LOOP: DO NM=1,NMESHES
 
-   M => MESHES(NM)
-
-   IF (.NOT.PROCESS_MESH_NEIGHBORHOOD(NM)) CYCLE MESH_LOOP
-
-   CALL POINT_TO_MESH(NM)
+   IF (.NOT.PRESENT(QUICK_READ)) THEN
+      M => MESHES(NM)
+      IF (.NOT.PROCESS_MESH_NEIGHBORHOOD(NM)) CYCLE MESH_LOOP
+      CALL POINT_TO_MESH(NM)
+   ENDIF
 
    ! Count OBST lines
 
@@ -9571,6 +9693,9 @@ SUBROUTINE READ_OBST
       CALL CHECKREAD('OBST',LU_INPUT,IOS)  ; IF (STOP_STATUS==SETUP_STOP) RETURN
       IF (IOS==1) EXIT COUNT_OBST_LOOP
       MULT_ID = 'null'
+      SURF_ID = 'null'
+      SURF_IDS = 'null'
+      SURF_ID6 = 'null'
       READ(LU_INPUT,NML=OBST,END=1,ERR=2,IOSTAT=IOS)
       CALL CHECK_XB(XB)
       MULT_INDEX = -1
@@ -9609,6 +9734,22 @@ SUBROUTINE READ_OBST
                   XB6 = XB(6) + MR%DZ0 + II*MR%DXB(6)
                ENDIF
                N_OBST_O = N_OBST_O + 1
+               IF (PRESENT(QUICK_READ)) THEN
+                  DO I=1,N_HT3D_SURF_LINES
+                     IF (HT3D_SURF_LIST(I)==SURF_ID .OR. ANY(HT3D_SURF_LIST(I)==SURF_IDS) .OR. &
+                                                         ANY(HT3D_SURF_LIST(I)==SURF_ID6)) THEN
+                        IF (SIZE(HT3D_OBST)<=N_HT3D_OBST) CALL REALLOCATE_HT3D_OBST(SIZE(HT3D_OBST),SIZE(HT3D_OBST)+100)
+                        N_HT3D_OBST = N_HT3D_OBST + 1
+                        HT3D_OBST(N_HT3D_OBST)%XS = XB1
+                        HT3D_OBST(N_HT3D_OBST)%XF = XB2
+                        HT3D_OBST(N_HT3D_OBST)%YS = XB3
+                        HT3D_OBST(N_HT3D_OBST)%YF = XB4
+                        HT3D_OBST(N_HT3D_OBST)%ZS = XB5
+                        HT3D_OBST(N_HT3D_OBST)%ZF = XB6
+                     ENDIF
+                  ENDDO
+                  CYCLE I_MULT_LOOP2
+               ENDIF
                IF (XB1>M%XF+M%DX(M%IBAR) .OR. XB2<M%XS-M%DX(1) .OR. &
                    XB3>M%YF+M%DY(M%JBAR) .OR. XB4<M%YS-M%DY(1) .OR. &
                    XB5>M%ZF+M%DZ(M%KBAR) .OR. XB6<M%ZS-M%DZ(1)) CYCLE I_MULT_LOOP2
@@ -9623,6 +9764,8 @@ SUBROUTINE READ_OBST
    ENDDO COUNT_OBST_LOOP
    1 REWIND(LU_INPUT) ; INPUT_FILE_LINE_NUMBER = 0
 
+   IF (PRESENT(QUICK_READ)) RETURN
+
    ! Allocate OBSTRUCTION array
 
    ALLOCATE(M%OBSTRUCTION(0:N_OBST_DIM),STAT=IZERO)
@@ -15534,4 +15677,15 @@ SUBROUTINE CALC_H2O_HV
 END SUBROUTINE CALC_H2O_HV
 
 
+SUBROUTINE REALLOCATE_HT3D_OBST(N1,N2)
+
+TYPE(HT3D_OBST_TYPE), ALLOCATABLE, DIMENSION(:) :: DUMMY
+INTEGER, INTENT(IN) :: N1,N2
+
+ALLOCATE(DUMMY(1:N2))
+DUMMY(1:N1) = HT3D_OBST(1:N1)
+CALL MOVE_ALLOC(DUMMY,HT3D_OBST)
+
+END SUBROUTINE REALLOCATE_HT3D_OBST
+
 END MODULE READ_INPUT

Verification/Heat_Transfer/ht3d_energy_conservation_8.fds

@@ -5,8 +5,6 @@
 
 &TIME T_END=900. /                                   
 
-&MISC NEIGHBOR_SEPARATION_DISTANCE=3 / 
-
 &SURF ID='HOT', TMP_FRONT=500, COLOR='RED', HEAT_TRANSFER_COEFFICIENT=0 /
 &OBST XB=-1.40, -1.20, -0.10, 0.10, 00.40, 01.00, COLOR='RED', SURF_ID6='HOT','INERT','INERT','INERT','INERT','INERT' /
 

Verification/Heat_Transfer/ht3d_sphere_96.fds

@@ -5,7 +5,7 @@
 
 &TIME T_END=180., DT=0.5 /
 
-&MISC NEIGHBOR_SEPARATION_DISTANCE=10, SOLID_PHASE_ONLY=T /
+&MISC SOLID_PHASE_ONLY=T /
 
 &RADI RADIATION=F /