From de721f423ee357dd2d41e7c0bef6a02268158151 Mon Sep 17 00:00:00 2001 From: ericvmueller Date: Fri, 6 Oct 2023 09:18:29 -0400 Subject: [PATCH 1/3] FDS Source: update NVEC for cut faces coming from split cell blocking; add CELL_BLOCK_ORIENTATION input --- Source/geom.f90 | 137 ++++++++++++++++++++++++++++-------------------- Source/type.f90 | 2 +- 2 files changed, 80 insertions(+), 59 deletions(-) diff --git a/Source/geom.f90 b/Source/geom.f90 index 919f11e4e44..2c6c93b0ab4 100644 --- a/Source/geom.f90 +++ b/Source/geom.f90 @@ -785,7 +785,7 @@ SUBROUTINE SET_CUTCELLS_3D LOGICAL, SAVE :: FIRST_CALL_ARG=.TRUE., FIRST_CALL_ARG2=.TRUE. -REAL(EB):: VERT_AUX(IAXIS:KAXIS) +REAL(EB):: VERT_AUX(IAXIS:KAXIS),CELL_BLOCK_ORIENTATION(IAXIS:KAXIS) INTEGER :: ING,INOD,IWSEL,IEL,FACE_AUX(NOD1:NOD3),VOL_AUX(NOD1:NOD4),N_SPCELLCF_TOT,N_SPCELL_TOT CHARACTER(100) :: FILENAME @@ -1380,14 +1380,16 @@ SUBROUTINE SET_CUTCELLS_3D DO ICC1=1,MESHES(NM)%N_CUTCELL_MESH+MESHES(NM)%N_GCCUTCELL_MESH CC=>MESHES(NM)%CUT_CELL(ICC1); IF(CC%NCELL<2) CYCLE ! Find if any GEOMETRY related to CC_INBOUNDARY faces has CELL_BLOCK_IOR>0: - CELL_BLOCK_IOR=0 + CELL_BLOCK_IOR=0; CELL_BLOCK_ORIENTATION = 0._EB NCELL_LOOP_1 : DO J=1,CC%NCELL DO I=2,CC%CCELEM(1,J)+1 IF(CC%FACE_LIST(1,CC%CCELEM(I,J))==CC_FTYPE_CFINB) THEN ICF=CC%FACE_LIST(4,CC%CCELEM(I,J)); JCF=CC%FACE_LIST(5,CC%CCELEM(I,J)) IG=MESHES(NM)%CUT_FACE(ICF)%BODTRI(1,JCF) IF(IG>0) THEN - IF (ABS(GEOMETRY(IG)%CELL_BLOCK_IOR)>0) THEN + IF (NORM2(GEOMETRY(IG)%CELL_BLOCK_ORIENTATION(IAXIS:KAXIS))>TWO_EPSILON_EB) THEN + CELL_BLOCK_ORIENTATION = GEOMETRY(IG)%CELL_BLOCK_ORIENTATION + ELSEIF (ABS(GEOMETRY(IG)%CELL_BLOCK_IOR)>0) THEN CELL_BLOCK_IOR = GEOMETRY(IG)%CELL_BLOCK_IOR EXIT NCELL_LOOP_1 ENDIF @@ -1395,32 +1397,35 @@ SUBROUTINE SET_CUTCELLS_3D ENDIF ENDDO ENDDO NCELL_LOOP_1 - ALLOCATE(VOLUME(1:CC%NCELL)) - ! Make search for double precision min/max unambiguous. - SELECT CASE(CELL_BLOCK_IOR) - CASE(-IAXIS,IAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(IAXIS,1:CC%NCELL) - CASE(-JAXIS,JAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(JAXIS,1:CC%NCELL) - CASE(-KAXIS,KAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(KAXIS,1:CC%NCELL) - CASE DEFAULT; VOLUME(1:CC%NCELL) = CC%VOLUME(1:CC%NCELL) - END SELECT - ! IF(NORM2(CELL_BLOCK_NVEC)>TWO_EPSILON_EB) THEN ! Block in the NVEC direction: - ! DO J=1,CC%NCELL; VOLUME(J) = DOT_PRODUCT(CELL_BLOCK_NVEC,CC%XYZCEN(IAXIS:KAXIS,J)); ENDDO - ! ENDIF - + ALLOCATE(VOLUME(1:CC%NCELL)); VOLUME(1:CC%NCELL) = CC%VOLUME(1:CC%NCELL) DO J=1,CC%NCELL; VOLUME(J) = VOLUME(J)*(1._EB+REAL(J-1,EB)*GEOMEPS); ENDDO - SELECT CASE(CELL_BLOCK_IOR) - CASE(-IAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1) - CASE( IAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1) - CASE(-JAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1) - CASE( JAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1) - CASE(-KAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1) - CASE( KAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1) - CASE DEFAULT; I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1) - END SELECT + I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1) + IF(NORM2(CELL_BLOCK_ORIENTATION)>TWO_EPSILON_EB) THEN + ! Cell Block Orientation: + DO J=1,CC%NCELL; VOLUME(J) = DOT_PRODUCT(CELL_BLOCK_ORIENTATION,CC%XYZCEN(IAXIS:KAXIS,J)); ENDDO + DO J=1,CC%NCELL; VOLUME(J) = VOLUME(J)*(1._EB+REAL(J-1,EB)*GEOMEPS); ENDDO + I=MINLOC(VOLUME(1:CC%NCELL),DIM=1) + ELSEIF(ABS(CELL_BLOCK_IOR)>0) THEN + ! Make search for double precision min/max unambiguous. + SELECT CASE(CELL_BLOCK_IOR) + CASE(-IAXIS,IAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(IAXIS,1:CC%NCELL) + CASE(-JAXIS,JAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(JAXIS,1:CC%NCELL) + CASE(-KAXIS,KAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(KAXIS,1:CC%NCELL) + END SELECT + DO J=1,CC%NCELL; VOLUME(J) = VOLUME(J)*(1._EB+REAL(J-1,EB)*GEOMEPS); ENDDO + SELECT CASE(CELL_BLOCK_IOR) + CASE(-IAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1) + CASE( IAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1) + CASE(-JAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1) + CASE( JAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1) + CASE(-KAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1) + CASE( KAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1) + END SELECT + ENDIF DEALLOCATE(VOLUME) NCELL_LOOP_2 : DO J=1,CC%NCELL IF(J==I) CYCLE NCELL_LOOP_2 - CC%NOADVANCE(J) = BLOCKED_SPLIT_CELL + IF(CC%NOADVANCE(J)==NOT_BLOCKED) CC%NOADVANCE(J) = BLOCKED_SPLIT_CELL ENDDO NCELL_LOOP_2 ENDDO ENDIF @@ -1473,7 +1478,7 @@ SUBROUTINE SET_CUTCELLS_3D END SELECT ENDDO IF(SUM_FACE>1 .OR. SUM_CCELL>0) CYCLE - CC%NOADVANCE(J)=BLOCKED_CAVITY_CELL + IF(CC%NOADVANCE(J)==NOT_BLOCKED) CC%NOADVANCE(J)=BLOCKED_CAVITY_CELL K=K+1 ENDDO ENDDO @@ -2455,7 +2460,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT IFACE = CC%CCELEM(IFC,JCC) IF( (CC%FACE_LIST(1,IFACE)==CC_FTYPE_RCGAS .OR. CC%FACE_LIST(1,IFACE)==CC_FTYPE_CFGAS) .AND. & CC%FACE_LIST(2,IFACE)==LOHIF .AND. CC%FACE_LIST(3,IFACE)==X1AXIS ) THEN - CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER + IF(CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER CYCLE JCC_LOOP_1 ENDIF ENDDO @@ -2528,7 +2533,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT IFACE = CC%CCELEM(IFC,JCC) IF( (CC%FACE_LIST(1,IFACE)==CC_FTYPE_RCGAS .OR. CC%FACE_LIST(1,IFACE)==CC_FTYPE_CFGAS) .AND. & CC%FACE_LIST(2,IFACE)==LOHIF .AND. CC%FACE_LIST(3,IFACE)==X1AXIS ) THEN - CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER + IF(CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER CYCLE JCC_LOOP_3 ENDIF ENDDO @@ -2602,7 +2607,8 @@ SUBROUTINE TAG_BLOCK_CELL(NM,II1,JJ1,KK1,NOM,IIO1,JJO1,KKO1,FINE_CELL) ! Here Test if cut-cells in II,KK,KK are blocked or not in IIO,JJO,KKO: IF(ICC>0) THEN IF(M2%CCVAR(IIO1,JJO1,KKO1,CC_CGSC)==CC_SOLID) THEN - M%CUT_CELL(ICC)%NOADVANCE(1:M%CUT_CELL(ICC)%NCELL) = BLOCKED_REFI_INTER + WHERE(M%CUT_CELL(ICC)%NOADVANCE(1:M%CUT_CELL(ICC)%NCELL)==NOT_BLOCKED) & + M%CUT_CELL(ICC)%NOADVANCE(1:M%CUT_CELL(ICC)%NCELL) = BLOCKED_REFI_INTER ELSE; CALL TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2) ENDIF ENDIF @@ -2921,7 +2927,7 @@ SUBROUTINE TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2) ! ENDIF ENDDO INBFC_LOC_LOOP ! Here set no ADVANCE if BLOCK_CELL=T: - IF(BLOCK_CELL) CC%NOADVANCE(JCC) = BLOCKED_REFI_INTER + IF(BLOCK_CELL .AND. CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC) = BLOCKED_REFI_INTER ENDDO JCC_LOOP ! IF(NM==1 .AND. ICC<30) CLOSE(LU_CCELL) @@ -4148,7 +4154,7 @@ SUBROUTINE GET_EXT_INB_CUTFACES_TO_CFACE MESHES(NM)%INTERNAL_CFACE_CELLS_LB = MESHES(NM)%N_EXTERNAL_CFACE_CELLS + MESHES(NM)%N_INTWALL_CFACE_CELLS ! Define pointers among CC_INBOUNDARY CUT_FACE and CFACE (N_INTERNAL_CFACE_CELLS): DO ICF=1,MESHES(NM)%N_CUTFACE_MESH - CF => CUT_FACE(ICF); IF(CF%STATUS /= CC_INBOUNDARY) CYCLE + CF => MESHES(NM)%CUT_FACE(ICF); IF(CF%STATUS /= CC_INBOUNDARY) CYCLE I = CF%IJK(IAXIS); J = CF%IJK(JAXIS); K = CF%IJK(KAXIS) ! Don't count INB cut-faces inside an OBST: IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE @@ -7715,6 +7721,7 @@ SUBROUTINE DROP_CUTFACE(NM,FTYPE,I,J,K,ILHF,X1AXIS,IFC,JFC) CF%BODTRI( :,IND(DUM)) = CF%BODTRI( :,DUM) CF%SURF_INDEX( IND(DUM)) = CF%SURF_INDEX( DUM) CF%BLK_TAG( IND(DUM)) = CF%BLK_TAG( DUM) + CF%CFACE_ORIGIN( IND(DUM)) = CF%CFACE_ORIGIN( DUM) CF%AREA_ADJUST( IND(DUM)) = CF%AREA_ADJUST( DUM) ENDIF DO ILH=LOW_IND,CT @@ -7743,6 +7750,7 @@ SUBROUTINE DROP_CUTFACE(NM,FTYPE,I,J,K,ILHF,X1AXIS,IFC,JFC) IF(FTYPE==CC_FTYPE_CFINB) THEN CF%BODTRI( :,CF%NFACE) = CC_UNDEFINED CF%SURF_INDEX( CF%NFACE) = CC_UNDEFINED + CF%CFACE_ORIGIN( CF%NFACE) = CC_UNDEFINED ENDIF DEALLOCATE(IND) ENDIF NFACE_IF_1 @@ -8330,7 +8338,7 @@ SUBROUTINE GET_CELL_LINK_INFO(NM) DO ICC=1,M%N_CUTCELL_MESH DO JCC=1,M%CUT_CELL(ICC)%NCELL IF(M%CUT_CELL(ICC)%IJK_LINK(1,JCC)==CC_SOLID) THEN - M%CUT_CELL(ICC)%NOADVANCE(JCC) = BLOCKED_UNLINK_CELL + IF(M%CUT_CELL(ICC)%NOADVANCE(JCC)==NOT_BLOCKED) M%CUT_CELL(ICC)%NOADVANCE(JCC) = BLOCKED_UNLINK_CELL M%CUT_CELL(ICC)%LINK_LEV(JCC) = CC_UNDEFINED M%CUT_CELL(ICC)%IJK_LINK(2:5,JCC)= CC_UNDEFINED ELSEIF(M%CUT_CELL(ICC)%LINK_LEV(JCC) < M%FINEST_LINK_LEV) THEN @@ -8433,6 +8441,7 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, USE GEOMETRY_FUNCTIONS, ONLY : SEARCH_OTHER_MESHES USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_STORAGE +USE MATH_FUNCTIONS, ONLY : CROSS_PRODUCT ! Routine that initializes new CFACE with index CFACE_INDEX. ! Geometry information for CFACE is loaded from MESHES(NM)%CUT_FACE(ICF)%AREA(IFACE), etc. @@ -8453,9 +8462,12 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, TYPE (WALL_TYPE), POINTER :: WC=>NULL() TYPE (MESH_TYPE), POINTER :: M TYPE (CFACE_TYPE), POINTER :: CFA +TYPE (CC_CUTFACE_TYPE), POINTER :: CF + M => MESHES(NM) SF=> SURFACE(SURF_INDEX) +CF=> CUT_FACE(ICF) STAGE_FLG_BRANCH : SELECT CASE(STAGE_FLG) @@ -8469,17 +8481,17 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, CFA%SURF_INDEX = SURF_INDEX - BC%X = CUT_FACE(ICF)%XYZCEN(IAXIS,IFACE) - BC%Y = CUT_FACE(ICF)%XYZCEN(JAXIS,IFACE) - BC%Z = CUT_FACE(ICF)%XYZCEN(KAXIS,IFACE) - CFA%AREA = CUT_FACE(ICF)%AREA(IFACE) + BC%X = CF%XYZCEN(IAXIS,IFACE) + BC%Y = CF%XYZCEN(JAXIS,IFACE) + BC%Z = CF%XYZCEN(KAXIS,IFACE) + CFA%AREA = CF%AREA(IFACE) ! Now populate cut-face information: CFA%CUT_FACE_IND1 = ICF CFA%CUT_FACE_IND2 = IFACE INS_INB_COND_1 : IF (IS_INB) THEN - CFA%VEL_ERR_NEW=CUT_FACE(ICF)%VEL(IFACE) - 0._EB ! Assumes zero veloc of solid. + CFA%VEL_ERR_NEW=CF%VEL(IFACE) - 0._EB ! Assumes zero veloc of solid. ! Check if fire spreads radially over this surface type IF (SF%FIRE_SPREAD_RATE>0._EB) THEN @@ -8489,21 +8501,28 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, ENDIF ! Normal to cut-face: - IBOD =CUT_FACE(ICF)%BODTRI(1,IFACE) - IWSEL=CUT_FACE(ICF)%BODTRI(2,IFACE) - CFA%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL) + 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(CFA%NVEC(IAXIS:KAXIS),V2,V3) + IF(NORM2(CFA%NVEC)>TWO_EPSILON_EB .AND. CF%CFACE_ORIGIN(IFACE)==BLOCKED_SPLIT_CELL) THEN + CFA%NVEC(IAXIS:KAXIS) = CFA%NVEC(IAXIS:KAXIS)/NORM2(CFA%NVEC) + ELSE + IBOD =CF%BODTRI(1,IFACE) + IWSEL=CF%BODTRI(2,IFACE) + CFA%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL) + ENDIF ! Boundary CFACES processed are defined of type SOLID_BOUNDARY CFA%BOUNDARY_TYPE = SOLID_BOUNDARY ! Might need to rethink this, but for the time being... - BC%II = CUT_FACE(ICF)%IJK(IAXIS) - BC%JJ = CUT_FACE(ICF)%IJK(JAXIS) - BC%KK = CUT_FACE(ICF)%IJK(KAXIS) + BC%II = CF%IJK(IAXIS) + BC%JJ = CF%IJK(JAXIS) + BC%KK = CF%IJK(KAXIS) - BC%IIG = CUT_FACE(ICF)%IJK(IAXIS) - BC%JJG = CUT_FACE(ICF)%IJK(JAXIS) - BC%KKG = CUT_FACE(ICF)%IJK(KAXIS) + BC%IIG = CF%IJK(IAXIS) + BC%JJG = CF%IJK(JAXIS) + BC%KKG = CF%IJK(KAXIS) ELSE INS_INB_COND_1 ! External mesh boundary CFACE @@ -8532,14 +8551,14 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, ENDIF ENDIF INS_INB_COND_1 - B1%AREA = CUT_FACE(ICF)%AREA(IFACE) ! Init to CFACE AREA. + B1%AREA = CF%AREA(IFACE) ! Init to CFACE AREA. CASE(INTEGER_TWO) ! Assign AREA_ADJUST for CFACE, BCs information for CFACE. CFA => M%CFACE(CFACE_INDEX) B1 => M%BOUNDARY_PROP1(CFA%B1_INDEX) ! First: Assign AREA_ADJUST for CFACEs. - B1%AREA_ADJUST = CUT_FACE(ICF)%AREA_ADJUST(IFACE) + B1%AREA_ADJUST = CF%AREA_ADJUST(IFACE) CASE(INTEGER_THREE) @@ -8551,11 +8570,11 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, ! Associated cut-cell location in CUT_CELL array. ! This CFACE initialization assumes TMP,RHO,ZZ have been initialized in cut-cell ICC,JCC. - ICC = CUT_FACE(ICF)%CELL_LIST(2,LOW_IND,IFACE) - JCC = CUT_FACE(ICF)%CELL_LIST(3,LOW_IND,IFACE) + ICC = CF%CELL_LIST(2,LOW_IND,IFACE) + JCC = CF%CELL_LIST(3,LOW_IND,IFACE) ! Set TMP_F to Surface value and rest to ambient in underlying cartesian cell. - CFA%TMP_G = TMP_0(CUT_FACE(ICF)%IJK(KAXIS)) + CFA%TMP_G = TMP_0(CF%IJK(KAXIS)) IF (SF%TMP_FRONT > 0._EB) THEN B1%TMP_F = SF%TMP_FRONT ELSE @@ -8574,7 +8593,7 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, ! Assign normal velocity to CFACE from SURF input: B1%U_NORMAL_0 = SF%VEL ! Assign normal velocity from VOLUME_FLOW : - IBOD =CUT_FACE(ICF)%BODTRI(1,IFACE) + IBOD =CF%BODTRI(1,IFACE) IF(IBOD>0 .AND. ABS(SF%VOLUME_FLOW)>=TWO_EPSILON_EB) B1%U_NORMAL_0 = SF%VOLUME_FLOW / FDS_AREA_GEOM(SURF_INDEX,IBOD) ! Assign normal velocity from MASS_FLUX_TOTAL : IF(ABS(SF%MASS_FLUX_TOTAL)>=TWO_EPSILON_EB) B1%U_NORMAL_0 = SF%MASS_FLUX_TOTAL / RHOA * B1%AREA_ADJUST @@ -8583,8 +8602,8 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, ! Case of exposed Backing we need to find CFACE_INDEX of BACK CFACE. IF (SF%BACKING==EXPOSED .AND. SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN - IG = CUT_FACE(ICF)%BODTRI(1,IFACE) - TRI = CUT_FACE(ICF)%BODTRI(2,IFACE) + IG = CF%BODTRI(1,IFACE) + TRI = CF%BODTRI(2,IFACE) XP(IAXIS:KAXIS) = (/ BC%X, BC%Y, BC%Z /) ! CFACE centroid location. RDIR(IAXIS:KAXIS)= - GEOMETRY(IG)%FACES_NORMAL(IAXIS:KAXIS,TRI) ! Normal into the body. TRI_LOOP : DO IWSEL=1,GEOMETRY(IG)%N_FACES @@ -8720,8 +8739,8 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, B1%U_NORMAL_0 = WC_B1%U_NORMAL_0 ! Here downscale velocity: - IF (IFACE==CUT_FACE(ICF)%NFACE) WC_B1%U_NORMAL_0 = & - WC_B1%U_NORMAL_0 * SUM(CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE))/WC_B1%AREA + IF (IFACE==CF%NFACE) WC_B1%U_NORMAL_0 = & + WC_B1%U_NORMAL_0 * SUM(CF%AREA(1:CF%NFACE))/WC_B1%AREA ! Vegetation T_IGN setup: B1%T_IGN = WC_B1%T_IGN @@ -21646,7 +21665,7 @@ SUBROUTINE READ_GEOM REAL(EB) :: SPHERE_ORIGIN(3),SPHERE_RADIUS,TEXTURE_ORIGIN(3),TEXTURE_SCALE(2),XB(6),DX,BOX_XYZ(3),& ZMIN,VOLUME,TXMIN,TXMAX,TYMIN,TYMAX,TX,TY,DV1(MAX_DIM),DV2(MAX_DIM),& NVECI(MAX_DIM),DXCEN(MAX_DIM),DOTI,TRANSPARENCY,CYLINDER_ORIGIN(3),CYLINDER_AXIS(3),& - CYLINDER_RADIUS,CYLINDER_LENGTH,EXTRUDE + CYLINDER_RADIUS,CYLINDER_LENGTH,EXTRUDE,CELL_BLOCK_ORIENTATION(3) INTEGER :: MAX_IDS=0,MAX_SURF_IDS=0,MAX_ZVALS=0,MAX_VERTS=0,MAX_FACES=0,MAX_VOLUS=0,MAX_POLY_VERTS=0,& N_VERTS,N_FACES,N_FACES_TEMP,N_VOLUS,N_ZVALS,N_SURF_ID,N_SURF_ID2,N_POLY_VERTS,& @@ -21693,7 +21712,7 @@ SUBROUTINE READ_GEOM INTEGER, PARAMETER :: DELTA_GEOM_LINE=1000 INTEGER :: GEOM_LINE_SIZE -NAMELIST /GEOM/ BNDF_GEOM,BINARY_FILE,CELL_BLOCK_IOR,COLOR,CYLINDER_ORIGIN,CYLINDER_AXIS,& +NAMELIST /GEOM/ BNDF_GEOM,BINARY_FILE,CELL_BLOCK_IOR,CELL_BLOCK_ORIENTATION,COLOR,CYLINDER_ORIGIN,CYLINDER_AXIS,& CYLINDER_RADIUS,CYLINDER_LENGTH,CYLINDER_NSEG_THETA,CYLINDER_NSEG_AXIS,& EXTRUDE,EXTEND_TERRAIN,FACES,ID,IJK,IS_TERRAIN,MOVE_ID,N_LAT,N_LEVELS,N_LONG,POLY,& RGB,SPHERE_ORIGIN,SPHERE_RADIUS,SPHERE_TYPE,SURF_ID,SURF_IDS,SURF_ID6,& @@ -21770,6 +21789,7 @@ SUBROUTINE READ_GEOM CALL COLOR2RGB(RGB,COLOR) ENDIF G%CELL_BLOCK_IOR = CELL_BLOCK_IOR + G%CELL_BLOCK_ORIENTATION = CELL_BLOCK_ORIENTATION G%RGB = RGB G%TRANSPARENCY = TRANSPARENCY N_VERTS=0 @@ -23792,6 +23812,7 @@ SUBROUTINE SET_GEOM_DEFAULTS BINARY_FILE = 'null' RGB=-1 CELL_BLOCK_IOR=0 + CELL_BLOCK_ORIENTATION = 0._EB END SUBROUTINE SET_GEOM_DEFAULTS diff --git a/Source/type.f90 b/Source/type.f90 index ef2fad583b2..d51cf4bd7fd 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -1027,7 +1027,7 @@ MODULE TYPES AZIM,ELEV,SCALE(3),XYZ(3),AZIM_DOT,ELEV_DOT,SCALE_DOT(3),XYZ_DOT(3),GAXIS(3),GROTATE,GROTATE_DOT,GROTATE_BASE,& XB(6),SPHERE_ORIGIN(3),SPHERE_RADIUS,TEXTURE_ORIGIN(3),TEXTURE_SCALE(2),MIN_LEDGE,MAX_LEDGE,MEAN_LEDGE,& GEOM_BOX(LOW_IND:HIGH_IND,IAXIS:KAXIS),TRANSPARENCY,CYLINDER_ORIGIN(3),CYLINDER_AXIS(3),& - CYLINDER_RADIUS,CYLINDER_LENGTH + CYLINDER_RADIUS,CYLINDER_LENGTH,CELL_BLOCK_ORIENTATION(3)=0._EB INTEGER, ALLOCATABLE,DIMENSION(:) :: FACES,VOLUS,SUB_GEOMS,SURFS,MATLS INTEGER, ALLOCATABLE,DIMENSION(:,:) :: EDGES,FACE_EDGES,EDGE_FACES REAL(EB),ALLOCATABLE,DIMENSION(:) :: FACES_AREA,VERTS_BASE,VERTS,TFACES,DAZIM,DELEV,ZVALS From de3e43ae00a7ec3632928e1f859357a72eb1aa8b Mon Sep 17 00:00:00 2001 From: ericvmueller Date: Fri, 6 Oct 2023 09:40:00 -0400 Subject: [PATCH 2/3] FDS Verification: use CELL_BLOCK_ORIENTATION in geom_rad_2.fds --- Verification/Radiation/geom_rad_2.fds | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Verification/Radiation/geom_rad_2.fds b/Verification/Radiation/geom_rad_2.fds index 179af6d48eb..2ae05bc6b81 100644 --- a/Verification/Radiation/geom_rad_2.fds +++ b/Verification/Radiation/geom_rad_2.fds @@ -23,8 +23,9 @@ &SURF ID='COLD', TMP_FRONT=20, COLOR='BLUE', HEAT_TRANSFER_COEFFICIENT=0, TAU_T=0, EMISSIVITY=1 / &MOVE ID='move', X0=0., Y0=0.02, Z0=0., AXIS=0.1,0.3,0.5, ROTATION_ANGLE=42.4, DZ=-0.5, DY=-0.5 / -&GEOM XB=-0.02,0.02,0.0,1.0,-0.02,1.0, SURF_ID6='INERT','HOT','INERT','INERT','INERT','INERT', MOVE_ID='move' / -&GEOM XB=-0.02,1.0,0.0,1.0,-0.02,0.02, SURF_ID6='INERT','INERT','INERT','INERT','INERT','COLD', MOVE_ID='move' / +&GEOM XB=-0.02,0.02,0.0,1.0,-0.02,1.0, SURF_ID6='INERT','HOT','INERT','INERT','INERT','INERT', MOVE_ID='move', CELL_BLOCK_ORIENTATION= -0.746,-0.592,0.305/ +&GEOM XB=-0.02,1.0,0.0,1.0,-0.02,0.02, SURF_ID6='INERT','INERT','INERT','INERT','INERT','COLD', MOVE_ID='move', CELL_BLOCK_ORIENTATION= -0.379,0, -0.925/ + &BNDF QUANTITY='INCIDENT HEAT FLUX', CELL_CENTERED=T / &BNDF QUANTITY='WALL TEMPERATURE', CELL_CENTERED=T / From 81e40cf4edd31c72e4712030f3e43610f2d26fa6 Mon Sep 17 00:00:00 2001 From: ericvmueller Date: Fri, 6 Oct 2023 11:12:42 -0400 Subject: [PATCH 3/3] FDS User Guide: add description of CELL_BLOCK_ORIENTATION --- Manuals/FDS_User_Guide/FDS_User_Guide.tex | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Manuals/FDS_User_Guide/FDS_User_Guide.tex b/Manuals/FDS_User_Guide/FDS_User_Guide.tex index b5e24421e30..2b5a5819ffe 100644 --- a/Manuals/FDS_User_Guide/FDS_User_Guide.tex +++ b/Manuals/FDS_User_Guide/FDS_User_Guide.tex @@ -1535,7 +1535,7 @@ \subsection{Reading Geometry Node Locations And Connectivity Data From Binary} \subsection{Handling Split Cells and Thin Geometries} \label{info:thin_geom} -One current limitation of the unstructured geometry algorithm is that a Cartesian cell cannot be split by a thin geometry or corner, potentially creating two background pressure zones within one Cartesian cell. By default, if a cell is split by a thin geometry, FDS will fill in the smaller of the two disjoint volumes. However, it may be desired to leave one side perfectly smooth while allowing the other to absorb the geometry change. For example, you want the inside wall of a tunnel or ceiling to be smooth but you are not as concerned with the outer portion for your problem. In this scenario, you can control which side of the geometry gets filled by setting {\ct CELL\_BLOCK\_IOR} on the {\ct GEOM} line. An example of how this works is shown in Fig.~\ref{figure:thin_geometry}. Currently, this functionality is limited to geometries where the simple {\ct IOR} can define the direction. For example, this would not work for a hollow sphere. +One current limitation of the unstructured geometry algorithm is that a Cartesian cell cannot be split by a thin geometry or corner, potentially creating two background pressure zones within one Cartesian cell. By default, if a cell is split by a thin geometry, FDS will fill in the smaller of the two disjoint volumes. However, it may be desired to leave one side perfectly smooth while allowing the other to absorb the geometry change. For example, you want the inside wall of a tunnel or ceiling to be smooth but you are not as concerned with the outer portion for your problem. In this scenario, you can control which side of the geometry gets filled by setting either {\ct CELL\_BLOCK\_IOR} or {\ct CELL\_BLOCK\_ORIENTATION} on the {\ct GEOM} line. The former is an integer indicating a Cartesian axis while the latter allows the specification of a vector in any direction. An example of how this works is shown in Fig.~\ref{figure:thin_geometry}. Currently, this functionality is limited to geometries where the simple {\ct IOR} can define the direction. For example, this would not work for a hollow sphere. \begin{figure} @@ -11549,6 +11549,7 @@ \section{\texorpdfstring{{\tt GEOM}}{GEOM} (Unstructured Geometry Parameters)} {\ct BINARY\_FILE} & Character & Section~\ref{subsec:readbin} & & {\ct 'null'} \\ \hline {\ct BNDF\_GEOM} & Logical & Section~\ref{info:BNDF} & & {\ct T} \\ \hline {\ct CELL\_BLOCK\_IOR} & Integer & Section~\ref{info:thin_geom} & & 0 \\ \hline +{\ct CELL\_BLOCK\_ORIENTATION} & Real Triplet & Section~\ref{info:thin_geom} & & 0.0,0.0,0.0 \\ \hline {\ct COLOR} & Character & Section~\ref{info:colors} & & {\ct 'null'} \\ \hline {\ct CYLINDER\_LENGTH} & Real & Section~\ref{info:GEOM_cyls} & m & \\ \hline {\ct CYLINDER\_RADIUS} & Real & Section~\ref{info:GEOM_cyls} & m & \\ \hline